text:
module FockPackage

using LinearAlgebra: dot
using Printf: @printf,@sprintf
using ..Spatials: PID,AbstractBond,Bond,decompose
using ..DegreesOfFreedom: IID,Index,Internal,FilteredAttributes,IDFConfig,Table,OID,Operator,Operators,LaTeX,latexformat,coordpresent
using ..Terms: wildcard,constant,Subscript,Subscripts,subscriptexpr,Coupling,Couplings,@subscript,couplingcenters,Term,TermCouplings,TermAmplitude,TermModulate
using ...Interfaces: id,rank,kind
using ...Prerequisites: Float,delta,decimaltostr
using ...Mathematics.AlgebraOverFields: SimpleID,ID,Element
using ...Mathematics.VectorSpaces: VectorSpace,IsMultiIndexable,MultiIndexOrderStyle

import ..DegreesOfFreedom: script,otype,isHermitian
import ..Terms: couplingcenter,statistics,abbr,termfactor
import ...Interfaces: dims,inds,⊗,⋅,expand,expand!,permute

export ANNIHILATION,CREATION,MAJORANA,FID,FIndex,Fock
export usualfockindextotuple,nambufockindextotuple
export FockOperator,FOperator,BOperator,isnormalordered
export foptdefaultlatex,boptdefaultlatex
export FCID,FockCoupling
export σ⁰,σˣ,σʸ,σᶻ,σ⁺,σ⁻
export @fc_str,@σ⁰_str,@σˣ_str,@σʸ_str,@σᶻ_str,@σ⁺_str,@σ⁻_str
export Onsite,Hopping,Pairing
export hubbard,Hubbard
export interorbitalinterspin,InterOrbitalInterSpin
export interorbitalintraspin,InterOrbitalIntraSpin
export spinflip,SpinFlip
export pairhopping,PairHopping
export Coulomb
export FFockTerm,BFockTerm

"""
    ANNIHILATION

Indicate that the nambu index is ANNIHILATION.
"""
const ANNIHILATION=1
"""
    CREATION

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

"""
    MAJORANA

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

"""
    FID <: IID

The Fock id.
"""
struct FID <: IID
    orbital::Int
    spin::Int
    nambu::Int
    FID(orbital::Int,spin::Int,nambu::Int)=(@assert nambu ∈ (0,1,2) "FID error: wrong input nambu($nambu).";new(orbital,spin,nambu))
end
Base.fieldnames(::Type{FID})=(:orbital,:spin,:nambu)

"""
    FID(;orbital::Int=1,spin::Int=1,nambu::Int=ANNIHILATION)

Create a Fock id.
"""
FID(;orbital::Int=1,spin::Int=1,nambu::Int=ANNIHILATION)=FID(orbital,spin,nambu)

"""
    adjoint(fid::FID) -> FID

Get the adjoint of a Fock id.
"""
Base.adjoint(fid::FID)=FID(fid.orbital,fid.spin,fid.nambu==0 ? 0 : 3-fid.nambu)

"""
    Fock <: Internal{FID}

The Fock interanl degrees of freedom.
"""
struct Fock <: Internal{FID}
    atom::Int
    norbital::Int
    nspin::Int
    nnambu::Int
    Fock(atom::Int,norbital::Int,nspin::Int,nnambu::Int)=(@assert nnambu ∈ (1,2) "Fock error: wrong input nnambu($nnambu).";new(atom,norbital,nspin,nnambu))
end
IsMultiIndexable(::Type{Fock})=IsMultiIndexable(true)
MultiIndexOrderStyle(::Type{Fock})=MultiIndexOrderStyle('C')
dims(fock::Fock)=(fock.norbital,fock.nspin,fock.nnambu)
inds(fid::FID,::Fock)=(fid.orbital,fid.spin,fid.nambu)
FID(inds::NTuple{3,Int},::Fock)=FID(inds[1],inds[2],inds[3])

"""
    Fock(;atom::Int=1,norbital::Int=1,nspin::Int=2,nnambu::Int=2)

Construct a Fock degrees of freedom.
"""
Fock(;atom::Int=1,norbital::Int=1,nspin::Int=2,nnambu::Int=2)=Fock(atom,norbital,nspin,nnambu)

"""
    FIndex{S} <: Index{PID{S},FID}

The Fock index.
"""
struct FIndex{S} <: Index{PID{S},FID}
    scope::S
    site::Int
    orbital::Int
    spin::Int
    nambu::Int
    function FIndex(scope::S,site::Int,orbital::Int,spin::Int,nambu::Int) where S
        @assert nambu ∈ (0,1,2) "FIndex error: wrong input nambu($nambu)."
        new{S}(scope,site,orbital,spin,nambu)
    end
end
Base.fieldnames(::Type{<:FIndex})=(:scope,:site,:orbital,:spin,:nambu)

"""
    union(::Type{P},::Type{FID}) where P<:PID

Get the union type of `PID` and `FID`.
"""
Base.union(::Type{P},::Type{FID}) where P<:PID=FIndex{fieldtype(P,:scope)}

"""
    script(oid::OID{<:FIndex},::Val{:site}) -> Int
    script(oid::OID{<:FIndex},::Val{:orbital}) -> Int
    script(oid::OID{<:FIndex},::Val{:spinint}) -> Int
    script(oid::OID{<:FIndex},::Val{:spinsym}) -> String
    script(oid::OID{<:FIndex},::Val{:nambu}) -> String

Get the required script of an Fock oid.
"""
script(oid::OID{<:FIndex},::Val{:site})=oid.index.site
script(oid::OID{<:FIndex},::Val{:orbital})=oid.index.orbital
script(oid::OID{<:FIndex},::Val{:spinint})=oid.index.spin
script(oid::OID{<:FIndex},::Val{:spinsym})=oid.index.spin==1 ? "↓" : oid.index.spin==2 ? "↑" : error("script error: spin($oid.index.spin) not in (1,2).")
script(oid::OID{<:FIndex},::Val{:nambu})=oid.index.nambu==CREATION ? "\\dagger" : ""

"""
    angle(id::OID{<:FIndex},vectors::AbstractVector{<:AbstractVector{Float}},values::AbstractVector{Float}) -> Complex{Float}

Get the twist phase corresponding to a Fock oid.
"""
function Base.angle(id::OID{<:FIndex},vectors::AbstractVector{<:AbstractVector{Float}},values::AbstractVector{Float})
    phase=  length(vectors)==1 ? 2pi*dot(decompose(id.icoord,vectors[1]),values) :
            length(vectors)==2 ? 2pi*dot(decompose(id.icoord,vectors[1],vectors[2]),values) :
            length(vectors)==3 ? 2pi*dot(decompose(id.icoord,vectors[1],vectors[2],vectors[3]),values) :
            error("angle error: not supported number of input basis vectors.")
    id.index.nambu==ANNIHILATION ? phase : id.index.nambu==CREATION ? -phase : error("angle error: not supported Fock index.")
end

"""
    usualfockindextotuple

Indicate that the filtered attributes are `(:scope,:site,:orbital,:spin)` when converting a Fock index to tuple.
"""
const usualfockindextotuple=FilteredAttributes(:scope,:site,:orbital,:spin)
"""
    nambufockindextotuple

Indicate that the filtered attributes are `(:scope,:nambu,:site,:orbital,:spin)` when converting a Fock index to tuple.
"""
const nambufockindextotuple=FilteredAttributes(:scope,:nambu,:site,:orbital,:spin)

"""
    FockOperator{V<:Number,I<:ID{OID}} <: Operator{V,I}

Abstract type for all Fock operators.
"""
abstract type FockOperator{V<:Number,I<:ID{OID}} <: Operator{V,I} end

"""
    isnormalordered(opt::FockOperator) -> Bool

Judge whether a FockOperator is normal ordered.
"""
function isnormalordered(opt::FockOperator)
    flag=true
    for i=1:rank(opt)
        flag && opt.id[i].index.nambu==ANNIHILATION && (flag=false)
        flag || opt.id[i].index.nambu==CREATION && return false
    end
    return true
end

"""
    FOperator(value::Number,id::ID{OID}=ID())

Fermionic Fock operator.
"""
struct FOperator{V<:Number,I<:ID{OID}} <: FockOperator{V,I}
    value::V
    id::I
    FOperator(value::Number)=new{typeof(value),Tuple{}}(value,())
    FOperator(value::Number,id::ID{OID})=new{typeof(value),typeof(id)}(value,id)
end

"""
    foptdefaultlatex

The default LaTeX format of the oids of a fermionic operator.
"""
const foptdefaultlatex=LaTeX{(:nambu,),(:site,:orbital,:spinsym)}('c')
latexformat(FOperator,foptdefaultlatex)

"""
    statistics(opt::FOperator) -> Char
    statistics(::Type{<:FOperator}) -> Char

Get the statistics of FOperator.
"""
statistics(opt::FOperator)=opt|>typeof|>statistics
statistics(::Type{<:FOperator})='F'

"""
    permute(::Type{<:FOperator},id1::OID{<:FIndex},id2::OID{<:FIndex},::Any=nothing) -> Tuple{Vararg{FOperator}}

Permute two fermionic oid and get the result.
"""
function permute(::Type{<:FOperator},id1::OID{<:FIndex},id2::OID{<:FIndex},::Any=nothing)
    @assert id1.index≠id2.index || id1.rcoord≠id2.rcoord || id1.icoord≠id2.icoord "permute error: permuted ids should not be equal to each other."
    if id1.index'==id2.index && id1.rcoord==id2.rcoord && id1.icoord==id2.icoord
        return (FOperator(1),FOperator(-1,ID(id2,id1)))
    else
        return (FOperator(-1,ID(id2,id1)),)
    end
end

"""
    *(f1::FOperator,f2::FOperator) -> Union{Nothing,FOperator}

Get the multiplication of two fermionic Fock operators.
"""
function Base.:*(f1::FOperator,f2::FOperator)
    rank(f1)>0 && rank(f2)>0 && f1.id[end]==f2.id[1] && return nothing
    return invoke(*,Tuple{Element,Element},f1,f2)
end

"""
    BOperator(value::Number,id::ID{OID}=ID())

Bosonic Fock operator.
"""
struct BOperator{V<:Number,I<:ID{OID}} <: FockOperator{V,I}
    value::V
    id::I
    BOperator(value::Number)=new{typeof(value),Tuple{}}(value,())
    BOperator(value::Number,id::ID{OID})=new{typeof(value),typeof(id)}(value,id)
end

"""
    boptdefaultlatex

The default LaTeX format of the oids of a bosonic operator.
"""
const boptdefaultlatex=LaTeX{(:nambu,),(:site,:orbital,:spinsym)}('b')
latexformat(BOperator,boptdefaultlatex)

"""
    statistics(opt::BOperator)
    statistics(::Type{<:BOperator})

Get the statistics of BOperator.
"""
statistics(opt::BOperator)=opt|>typeof|>statistics
statistics(::Type{<:BOperator})='B'

"""
    permute(::Type{<:BOperator},id1::OID{<:FIndex},id2::OID{<:FIndex},::Any=nothing) -> Tuple{Vararg{BOperator}}

Permute two bosonic oid and get the result.
"""
function permute(::Type{<:BOperator},id1::OID{<:FIndex},id2::OID{<:FIndex},::Any=nothing)
    @assert id1.index≠id2.index || id1.rcoord≠id2.rcoord || id1.icoord≠id2.icoord "permute error: permuted ids should not be equal to each other."
    if id1.index'==id2.index && id1.rcoord==id2.rcoord && id1.icoord==id2.icoord
        if id1.index.nambu==CREATION
            return (BOperator(1),BOperator(1,ID(id2,id1)))
        else
            return (BOperator(-1),BOperator(1,ID(id2,id1)))
        end
    else
        return (BOperator(1,ID(id2,id1)),)
    end
end

"""
    FCID(;center=wildcard,atom=wildcard,orbital=wildcard,spin=wildcard,nambu=wildcard,obsub=wildcard,spsub=wildcard)

The id of a Fock coupling.
"""
struct FCID{C,A,O,S,N,OS,SS} <: SimpleID
    center::C
    atom::A
    orbital::O
    spin::S
    nambu::N
    obsub::OS
    spsub::SS
end
Base.fieldnames(::Type{<:FCID})=(:center,:atom,:orbital,:spin,:nambu,:obsub,:spsub)
FCID(;center=wildcard,atom=wildcard,orbital=wildcard,spin=wildcard,nambu=wildcard,obsub=wildcard,spsub=wildcard)=FCID(center,atom,orbital,spin,nambu,obsub,spsub)

"""
    FockCoupling(value::Number,id::ID{FCID},obsubscripts::Subscripts,spsubscripts::Subscripts)
    FockCoupling{N}(value::Number=1;
                    centers::Union{NTuple{N,Int},Nothing}=nothing,
                    atoms::Union{NTuple{N,Int},Nothing}=nothing,
                    orbitals::Union{NTuple{N,Int},Subscript,Nothing}=nothing,
                    spins::Union{NTuple{N,Int},Subscript,Nothing}=nothing,
                    nambus::Union{NTuple{N,Int},Nothing}=nothing) where N

Fock coupling.
"""
struct FockCoupling{V,I<:ID{FCID},OS<:Subscripts,SS<:Subscripts} <: Coupling{V,I}
    value::V
    id::I
    obsubscripts::OS
    spsubscripts::SS
    function FockCoupling(value::Number,id::ID{FCID},obsubscripts::Subscripts,spsubscripts::Subscripts)
        new{typeof(value),typeof(id),typeof(obsubscripts),typeof(spsubscripts)}(value,id,obsubscripts,spsubscripts)
    end
end
function FockCoupling{N}(   value::Number=1;
                            centers::Union{NTuple{N,Int},Nothing}=nothing,
                            atoms::Union{NTuple{N,Int},Nothing}=nothing,
                            orbitals::Union{NTuple{N,Int},Subscript,Nothing}=nothing,
                            spins::Union{NTuple{N,Int},Subscript,Nothing}=nothing,
                            nambus::Union{NTuple{N,Int},Nothing}=nothing) where N
    centers===nothing && (centers=ntuple(i->wildcard,N))
    atoms===nothing && (atoms=ntuple(i->wildcard,N))
    obsubscript=isa(orbitals,NTuple{N,Int}) ? Subscript(orbitals) : isa(orbitals,Nothing) ? Subscript{N}() : orbitals
    spsubscript=isa(spins,NTuple{N,Int}) ? Subscript(spins) : isa(spins,Nothing) ? Subscript{N}() : spins
    orbitals=obsubscript.opattern
    spins=spsubscript.opattern
    nambus==nothing && (nambus=ntuple(i->wildcard,N))
    obsubs=ntuple(i->obsubscript.identifier,N)
    spsubs=ntuple(i->spsubscript.identifier,N)
    return FockCoupling(value,ID(FCID,centers,atoms,orbitals,spins,nambus,obsubs,spsubs),Subscripts(obsubscript),Subscripts(spsubscript))
end

"""
    show(io::IO,fc::FockCoupling)

Show a Fock coupling.
"""
function Base.show(io::IO,fc::FockCoupling)
    @printf io "FockCoupling{%s}(value=%s" rank(fc) decimaltostr(fc.value)
    cache=[]
    for attrname in (:centers,:atoms,:orbitals,:spins,:nambus)
        any((attrvalue=getproperty(fc.id,attrname)).≠wildcard) && @printf io ",%s=(%s)" attrname join(attrvalue,",")
    end
    for attrname in (:obsubs,:spsubs)
        all((attrvalue=getproperty(fc.id,attrname)).==wildcard) || all((attrvalue.==constant)) || @printf io ",%s=(%s)" attrname join(attrvalue,",")
    end
    @printf io ")"
end

"""
    repr(fc::FockCoupling) -> String

Get the repr representation of a Fock coupling.
"""
function Base.repr(fc::FockCoupling)
    cache=[]
    flag=true
    for (attrname,abbr) in zip((:atoms,:orbitals,:spins,:nambus),("sl","ob","sp","ph"))
        any((attrvalue=getproperty(fc.id,attrname)).≠wildcard) && push!(cache,@sprintf "%s(%s)" abbr join(attrvalue,','))
    end
    result=decimaltostr(fc.value)
    length(cache)>0 && (flag=false;result=@sprintf "%s %s" result join(cache,"⊗"))
    any((centers=fc.id.centers).≠wildcard) && (flag=false;result=@sprintf "%s%s@(%s)" result (length(cache)>0 ? "" : " ") join(centers,','))
    obsubs,spsubs=fc.id.obsubs,fc.id.spsubs
    ((all(obsubs.==wildcard) || all(obsubs.==constant)) && (all(spsubs.==wildcard) || all(spsubs.==constant))) || (
            flag=false;
            result=@sprintf "%s with %s" result join(((@sprintf "(%s,%s)" obsub spsub) for (obsub,spsub) in zip(obsubs,spsubs))," && ")
            )
    flag && (result=@sprintf "%s {%s}" result rank(fc))
    return result
end

"""
    *(fc1::FockCoupling,fc2::FockCoupling) -> FockCoupling

Get the multiplication between two Fock couplings.
"""
function Base.:*(fc1::FockCoupling,fc2::FockCoupling)
    FockCoupling(fc1.value*fc2.value,fc1.id*fc2.id,fc1.obsubscripts*fc2.obsubscripts,fc1.spsubscripts*fc2.spsubscripts)
end

"""
    ⊗(fc1::FockCoupling,fc2::FockCoupling) -> FockCoupling

Get the direct product between two Fock couplings.
"""
function ⊗(fc1::FockCoupling,fc2::FockCoupling)
    @assert fc1|>rank==fc2|>rank "⊗ error: dismatched rank."
    value,attrvalues,subscripts,N=fc1.value*fc2.value,[],[],rank(fc1)
    wildcards=NTuple{N,Char}(wildcard for i=1:N)
    v1,v2=fc1.id.centers,fc2.id.centers; t1,t2=v1==wildcards,v2==wildcards
    push!(attrvalues,t1 && t2 ? wildcards : t1 ? v2 : t2 ? v1 : v1==v2 ? v1 : error("⊗ error: dismatched centers ($v1 and $v2)."))
    for (i,attrname) in enumerate((:atoms,:orbitals,:spins,:nambus,:obsubs,:spsubs))
        v1,v2=getproperty(fc1.id,attrname),getproperty(fc2.id,attrname)
        t1,t2=v1==wildcards,v2==wildcards
        if attrname ∉ (:obsubs,:spsubs)
            @assert (t1 || t2) "⊗ error: dismatched $attrname ($v1 and $v2)."
        else
            @assert t1 || t2 "⊗ error: dismatched $attrname ($v1 and $v2)."
            attrname==:obsubs && push!(subscripts,t1 ? fc2.obsubscripts : fc1.obsubscripts)
            attrname==:spsubs && push!(subscripts,t1 ? fc2.spsubscripts : fc1.spsubscripts)
        end
        push!(attrvalues,t1 ? v2 : v1)
    end
    return FockCoupling(value,ID(FCID,attrvalues...),subscripts...)
end

"""
    ⋅(fc1::FockCoupling,fc2::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 sublattice, orbital, spin and nambu spaces.
The dot product here is defined as the multiplication between such matrices.
"""
function ⋅(fc1::FockCoupling,fc2::FockCoupling)
    @assert rank(fc1)==rank(fc2)==2 "⋅ error: input fockcouplings must be rank-2."
    attrpairs=[]
    value=fc1.value*fc2.value
    wildcards=(wildcard,wildcard)
    v1,v2=fc1.id.centers,fc2.id.centers
    t1,t2=v1==wildcards,v2==wildcards
    push!(attrpairs,:centers=>(t1 && t2 ? nothing : t1 ? v2 : t2 ? v1 : v1==v2 ? v1 : error("⋅ error: dismatched centers ($v1 and $v2).")))
    for (i,attrname) in enumerate((:atoms,:orbitals,:spins,:nambus))
        v1,v2=getproperty(fc1.id,attrname),getproperty(fc2.id,attrname)
        if isa(v1,NTuple{2,Int}) && isa(v2,NTuple{2,Int})
            value=value*delta(v1[2],v2[1])
            push!(attrpairs,attrname=>(v1[1],v2[2]))
        else
            @assert v1==wildcards && v2==wildcards "⋅ error: dismatched $attrname ($v1 and $v2)."
        end
    end
    return FockCoupling{2}(value;attrpairs...)
end

"""
    expand(fc::FockCoupling,pid::PID,fock::Fock,kind::Union{Val{K},Nothing}=nothing) where K -> Union{FCExpand,Tuple{}}
    expand(fc::FockCoupling,pids::NTuple{R,PID},focks::NTuple{R,Fock},kind::Union{Val{K},Nothing}=nothing) where {R,K} -> Union{FCExpand,Tuple{}}

Expand a Fock coupling with the given set of point ids and Fock degrees of freedom.
"""
expand(fc::FockCoupling,pid::PID,fock::Fock,kind::Union{Val{K},Nothing}=nothing) where K=expand(fc,(pid,),(fock,),kind)
function expand(fc::FockCoupling,pids::NTuple{R,PID},focks::NTuple{R,Fock},kind::Union{Val{K},Nothing}=nothing) where {R,K}
    centers=couplingcenters(typeof(fc),fc.id.centers,Val(R))
    rpids=NTuple{rank(fc),eltype(pids)}(pids[centers[i]] for i=1:rank(fc))
    rfocks=NTuple{rank(fc),eltype(focks)}(focks[centers[i]] for i=1:rank(fc))
    nambus=fockcouplingnambus(kind,fc.id.nambus,NTuple{rank(fc),Int}(rfocks[i].nnambu for i=1:rank(fc)))
    for (i,atom) in enumerate(fc.id.atoms)
       isa(atom,Int) && (atom≠rfocks[i].atom) && return ()
    end
    obsbexpands=expand(fc.obsubscripts,NTuple{rank(fc),Int}(rfocks[i].norbital for i=1:rank(fc)))|>collect
    spsbexpands=expand(fc.spsubscripts,NTuple{rank(fc),Int}(rfocks[i].nspin for i=1:rank(fc)))|>collect
    return FCExpand(fc.value,rpids,obsbexpands,spsbexpands,nambus)
end
couplingcenter(::Type{<:FockCoupling},i::Int,n::Int,::Val{2})=i<=n/2 ? 1 : 2
@generated function fockcouplingnambus(::Union{Val{K},Nothing},nambus::NTuple{N,Any},ranges::NTuple{N,Int}) where {N,K}
    exprs=[:(isa(nambus[$i],Int) ?
            ((ranges[$i]==1 && nambus[$i]==0) || (ranges[$i]==2 && 0<nambus[$i]<=2) ? nambus[$i] : error("fockcouplingnambus error: nambu out of range.")) :
            (ranges[$i]==1 ? 0 : ($i)%2==1 ? CREATION : ANNIHILATION)
            ) for i=1:N]
    return Expr(:tuple,exprs...)
end
struct FCExpand{V,N,S} <: VectorSpace{Tuple{V,NTuple{N,FIndex{S}}}}
    value::V
    pids::NTuple{N,PID{S}}
    obsbexpands::Vector{NTuple{N,Int}}
    spsbexpands::Vector{NTuple{N,Int}}
    nambus::NTuple{N,Int}
end
IsMultiIndexable(::Type{<:FCExpand})=IsMultiIndexable(true)
MultiIndexOrderStyle(::Type{<:FCExpand})=MultiIndexOrderStyle('C')
dims(fce::FCExpand)=(length(fce.obsbexpands),length(fce.spsbexpands))
@generated function Tuple(inds::NTuple{2,Int},fce::FCExpand{V,N}) where {V,N}
    exprs=[:(FIndex(fce.pids[$i],FID(fce.obsbexpands[inds[1]][$i],fce.spsbexpands[inds[2]][$i],fce.nambus[$i]))) for i=1:N]
    return Expr(:tuple,:(fce.value),Expr(:tuple,exprs...))
end

"""
    σ⁰(mode::String;centers::Union{NTuple{2,Int},Nothing}=nothing) -> Couplings{ID{FCID,2},FockCoupling{Int,<:ID{FCID,2}}}

The Pauli matrix σ⁰, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
function σ⁰(mode::String;centers::Union{NTuple{2,Int},Nothing}=nothing)
    @assert (mode=="sp" || mode=="ob" || mode=="sl" || mode=="ph") "σ⁰ error: not supported mode($mode)."
    attrname=mode=="sp" ? :spins : mode=="ob" ? :orbitals : mode=="sl" ? :atoms : :nambus
    (attrval1,attrval2)=mode=="ph" ? ((ANNIHILATION,CREATION),(CREATION,ANNIHILATION)) : ((1,1),(2,2))
    FockCoupling{2}(1;attrname=>attrval1,:centers=>centers)+FockCoupling{2}(1;attrname=>attrval2,:centers=>centers)
end

"""
    σˣ(mode::String;centers::Union{NTuple{2,Int},Nothing}=nothing) -> Couplings{ID{FCID,2},FockCoupling{Int,<:ID{FCID,2}}}

The Pauli matrix σˣ, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
function σˣ(mode::String;centers::Union{NTuple{2,Int},Nothing}=nothing)
    @assert (mode=="sp" || mode=="ob" || mode=="sl" || mode=="ph") "σˣ error: not supported mode($mode)."
    attrname=mode=="sp" ? :spins : mode=="ob" ? :orbitals : mode=="sl" ? :atoms : :nambus
    (attrval1,attrval2)=mode=="ph" ? ((ANNIHILATION,ANNIHILATION),(CREATION,CREATION)) : ((1,2),(2,1))
    FockCoupling{2}(1;attrname=>attrval1,:centers=>centers)+FockCoupling{2}(1;attrname=>attrval2,:centers=>centers)
end

"""
    σʸ(mode::String;centers::Union{NTuple{2,Int},Nothing}=nothing) -> Couplings{ID{FCID,2},FockCoupling{Int,<:ID{FCID,2}}}

The Pauli matrix σʸ, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
function σʸ(mode::String;centers::Union{NTuple{2,Int},Nothing}=nothing)
    @assert (mode=="sp" || mode=="ob" || mode=="sl" || mode=="ph") "σʸ error: not supported mode($mode)."
    attrname=mode=="sp" ? :spins : mode=="ob" ? :orbitals : mode=="sl" ? :atoms : :nambus
    (attrval1,attrval2)=mode=="ph" ? ((ANNIHILATION,ANNIHILATION),(CREATION,CREATION)) : ((1,2),(2,1))
    FockCoupling{2}(1im;attrname=>attrval1,:centers=>centers)+FockCoupling{2}(-1im;attrname=>attrval2,:centers=>centers)
end

"""
    σᶻ(mode::String;centers::Union{NTuple{2,Int},Nothing}=nothing) -> Couplings{ID{FCID,2},FockCoupling{Int,<:ID{FCID,2}}}

The Pauli matrix σᶻ, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
function σᶻ(mode::String;centers::Union{NTuple{2,Int},Nothing}=nothing)
    @assert (mode=="sp" || mode=="ob" || mode=="sl" || mode=="ph") "σᶻ error: not supported mode($mode)."
    attrname=mode=="sp" ? :spins : mode=="ob" ? :orbitals : mode=="sl" ? :atoms : :nambus
    (attrval1,attrval2)=mode=="ph" ? ((ANNIHILATION,CREATION),(CREATION,ANNIHILATION)) : ((1,1),(2,2))
    FockCoupling{2}(-1;attrname=>attrval1,:centers=>centers)+FockCoupling{2}(1;attrname=>attrval2,:centers=>centers)
end

"""
    σ⁺(mode::String;centers::Union{NTuple{2,Int},Nothing}=nothing) -> Couplings{ID{FCID,2},FockCoupling{Int,<:ID{FCID,2}}}

The Pauli matrix σ⁺, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
function σ⁺(mode::String;centers::Union{NTuple{2,Int},Nothing}=nothing)
    @assert (mode=="sp" || mode=="ob" || mode=="sl" || mode=="ph") "σ⁺ error: not supported mode($mode)."
    attrname=mode=="sp" ? :spins : mode=="ob" ? :orbitals : mode=="sl" ? :atoms : :nambus
    attrval=mode=="ph" ? (CREATION,CREATION) : (2,1)
    Couplings(FockCoupling{2}(1;attrname=>attrval,:centers=>centers))
end

"""
    σ⁻(mode::String;centers::Union{NTuple{2,Int},Nothing}=nothing) -> Couplings{ID{FCID,2},FockCoupling{Int,<:ID{FCID,2}}}

The Pauli matrix σ⁻, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
function σ⁻(mode::String;centers::Union{NTuple{2,Int},Nothing}=nothing)
    @assert (mode=="sp" || mode=="ob" || mode=="sl" || mode=="ph") "σ⁻ error: not supported mode($mode)."
    attrname=mode=="sp" ? :spins : mode=="ob" ? :orbitals : mode=="sl" ? :atoms : :nambus
    attrval=mode=="ph" ? (ANNIHILATION,ANNIHILATION) : (1,2)
    Couplings(FockCoupling{2}(1;attrname=>attrval,:centers=>centers))
end

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

Construct a FockCoupling from a literal string.
"""
macro fc_str(str)
    ps=split(str," with ")
    conditions=length(ps)==2 ? fcconditions(ps[2]) : nothing
    ps=split(ps[1],' ')
    coeff=eval(Meta.parse(ps[1]))
    if ps[2][1]=='{' && ps[2][3]=='}'
        N=parse(Int,ps[2][2])
        return FockCoupling{N}(coeff)
    else
        attrpairs=[]
        components=split(ps[2],'@')
        centers=length(components)==2 ? couplingcenters(components[2]) : nothing
        push!(attrpairs,:centers=>centers)
        N=centers===nothing ? nothing : length(centers)
        count=0
        if length(components[1])>0
            for component in split(components[1],'⊗')
                attrname,attrvalue=fccomponent(component)
                N===nothing && (N=length(attrvalue))
                if isa(attrvalue,Expr)
                    @assert attrname==:orbitals || attrname==:spins "@fc_str error: wrong input pattern."
                    @assert N==length(attrvalue.args) "@fc_str error: dismatched ranks."
                    count=count+1
                    condition=conditions===nothing ? :nothing : conditions[count]
                    push!(attrpairs,attrname=>eval(subscriptexpr(attrvalue,condition)))
                else
                    @assert N==length(attrvalue) "@fc_str error: dismatched ranks."
                    push!(attrpairs,attrname=>attrvalue)
                end
            end
        end
        return FockCoupling{N}(coeff;attrpairs...)
    end
end
function fcconditions(str::AbstractString)
    conditions=Meta.parse(str)
    conditions=conditions.head==:tuple ? conditions.args : [conditions]
    return [condition=="*" ? :nothing : condition for condition in conditions]
end
function fccomponent(str::AbstractString)
    @assert str[3]=='(' && str[end]==')' "fccomponent error: wrong input pattern."
    attrname=str[1:2]=="sl" ? :atoms : str[1:2]=="ob" ? :orbitals : str[1:2]=="sp" ? :spins : str[1:2]=="ph" ? :nambus : error("fccomponent error: wrong input pattern.")
    expr=Meta.parse(str[3:end])
    attrvalue=isa(expr,Expr) ? (all(isa(arg,Int) for arg in expr.args) ? Tuple(expr.args) : expr) : (expr,)
    return attrname=>attrvalue
end

σᵅsplit(str::AbstractString)=(ps=split(str,'@'); length(ps)==1 ? (ps[1],nothing) : length(ps)==2 ? (ps[1],couplingcenters(ps[2])) : "σᵅsplit error: wrong input pattern.")
σᵅname(mode::AbstractString)=mode=="sp" ? :spins : mode=="ob" ? :orbitals : mode=="sl" ? :atoms : mode=="ph" ? :nambus : error("σᵅname error: wrong input mode.")

"""
    σ⁰"sp"/σ⁰"sp@(c₁,c₂)" -> Couplings
    σ⁰"ob"/σ⁰"ob@(c₁,c₂)" -> Couplings
    σ⁰"sl"/σ⁰"sl@(c₁,c₂)" -> Couplings
    σ⁰"ph"/σ⁰"ph@(c₁,c₂)" -> Couplings

The Pauli matrix σ⁰, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
macro σ⁰_str(str::String)
    mode,centers=σᵅsplit(str)
    attrname=σᵅname(mode)
    (attrval1,attrval2)=mode=="ph" ? ((ANNIHILATION,CREATION),(CREATION,ANNIHILATION)) : ((1,1),(2,2))
    FockCoupling{2}(1;attrname=>attrval1,:centers=>centers)+FockCoupling{2}(1;attrname=>attrval2,:centers=>centers)
end

"""
    σˣ"sp"/σˣ"sp@(c₁,c₂)" -> Couplings
    σˣ"ob"/σˣ"ob@(c₁,c₂)" -> Couplings
    σˣ"sl"/σˣ"sl@(c₁,c₂)" -> Couplings
    σˣ"ph"/σˣ"ph@(c₁,c₂)" -> Couplings

The Pauli matrix σˣ, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
macro σˣ_str(str::String)
    mode,centers=σᵅsplit(str)
    attrname=σᵅname(mode)
    (attrval1,attrval2)=mode=="ph" ? ((ANNIHILATION,ANNIHILATION),(CREATION,CREATION)) : ((1,2),(2,1))
    FockCoupling{2}(1;attrname=>attrval1,:centers=>centers)+FockCoupling{2}(1;attrname=>attrval2,:centers=>centers)
end

"""
    σʸ"sp"/σʸ"sp@(c₁,c₂)" -> Couplings
    σʸ"ob"/σʸ"ob@(c₁,c₂)" -> Couplings
    σʸ"sl"/σʸ"sl@(c₁,c₂)" -> Couplings
    σʸ"ph"/σʸ"ph@(c₁,c₂)" -> Couplings

The Pauli matrix σʸ, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
macro σʸ_str(str::String)
    mode,centers=σᵅsplit(str)
    attrname=σᵅname(mode)
    (attrval1,attrval2)=mode=="ph" ? ((ANNIHILATION,ANNIHILATION),(CREATION,CREATION)) : ((1,2),(2,1))
    FockCoupling{2}(1im;attrname=>attrval1,:centers=>centers)+FockCoupling{2}(-1im;attrname=>attrval2,:centers=>centers)
end

"""
    σᶻ"sp"/σᶻ"sp@(c₁,c₂)" -> Couplings
    σᶻ"ob"/σᶻ"ob@(c₁,c₂)" -> Couplings
    σᶻ"sl"/σᶻ"sl@(c₁,c₂)" -> Couplings
    σᶻ"ph"/σᶻ"ph@(c₁,c₂)" -> Couplings

The Pauli matrix σᶻ, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
macro σᶻ_str(str::String)
    mode,centers=σᵅsplit(str)
    attrname=σᵅname(mode)
    (attrval1,attrval2)=mode=="ph" ? ((ANNIHILATION,CREATION),(CREATION,ANNIHILATION)) : ((1,1),(2,2))
    FockCoupling{2}(-1;attrname=>attrval1,:centers=>centers)+FockCoupling{2}(1;attrname=>attrval2,:centers=>centers)
end

"""
    σ⁺"sp"/σ⁺"sp@(c₁,c₂)" -> Couplings
    σ⁺"ob"/σ⁺"ob@(c₁,c₂)" -> Couplings
    σ⁺"sl"/σ⁺"sl@(c₁,c₂)" -> Couplings
    σ⁺"ph"/σ⁺"ph@(c₁,c₂)" -> Couplings

The Pauli matrix σ⁺, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
macro σ⁺_str(str::String)
    mode,centers=σᵅsplit(str)
    attrname=σᵅname(mode)
    attrval=mode=="ph" ? (CREATION,CREATION) : (2,1)
    Couplings(FockCoupling{2}(1;attrname=>attrval,:centers=>centers))
end

"""
    σ⁻"sp"/σ⁻"sp@(c₁,c₂)" -> Couplings
    σ⁻"ob"/σ⁻"ob@(c₁,c₂)" -> Couplings
    σ⁻"sl"/σ⁻"sl@(c₁,c₂)" -> Couplings
    σ⁻"ph"/σ⁻"ph@(c₁,c₂)" -> Couplings

The Pauli matrix σ⁻, which can act on the space of spins("sp"), orbitals("ob"), sublattices("sl") or particle-holes("ph").
"""
macro σ⁻_str(str::String)
    mode,centers=σᵅsplit(str)
    attrname=σᵅname(mode)
    attrval=mode=="ph" ? (ANNIHILATION,ANNIHILATION) : (1,2)
    Couplings(FockCoupling{2}(1;attrname=>attrval,:centers=>centers))
end

"""
    Onsite{ST}( id::Symbol,value::Any;
                couplings::Union{Function,Coupling,Couplings,Nothing}=nothing,
                amplitude::Union{Function,Nothing}=nothing,
                modulate::Union{Function,Bool}=false,
                ) where {ST}

Onsite term.

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

"""
    Hopping{ST}(id::Symbol,value::Any,bondkind::Int=1;
                couplings::Union{Function,Coupling,Couplings,Nothing}=nothing,
                amplitude::Union{Function,Nothing}=nothing,
                modulate::Union{Function,Bool}=false,
                ) where {ST}

Hopping term.

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

"""
    Pairing{ST}(id::Symbol,value::Any,bondkind::Int=0;
                couplings::Union{Function,Coupling,Couplings,Nothing}=nothing,
                amplitude::Union{Function,Nothing}=nothing,
                modulate::Union{Function,Bool}=false,
                ) where {ST}

Pairing term.

Type alias for `Term{statistics,:Pairing,2,id,V,Int,<:TermCouplings,<:TermAmplitude,<:Union{TermModulate,Nothing}}`.
"""
const Pairing{statistics,id,V,C<:TermCouplings,A<:TermAmplitude,M<:Union{TermModulate,Nothing}}=Term{statistics,:Pairing,2,id,V,Int,C,A,M}
function Pairing{ST}(id::Symbol,value::Any,bondkind::Int=0;
                    couplings::Union{Function,Coupling,Couplings,Nothing}=nothing,
                    amplitude::Union{Function,Nothing}=nothing,
                    modulate::Union{Function,Bool}=false,
                    ) where {ST}
    couplings=TermCouplings(couplings===nothing ? FockCoupling{2}() : couplings)
    Term{ST,:Pairing,2}(id,value,bondkind,couplings=couplings,amplitude=amplitude,modulate=modulate)
end
abbr(::Type{<:Pairing})=:pr
isHermitian(::Type{<:Pairing})=false
function fockcouplingnambus(::Val{:Pairing},nambus::NTuple{2,Any},ranges::NTuple{2,Int})
    @assert ranges==(2,2) "fockcouplingnambus error: ranges for Pairing terms must be (2,2)."
    (   isa(nambus[1],Int) ? (0<nambus[1]<=2 ? nambus[1] : error("fockcouplingnambus error: nambu out of range.")) : ANNIHILATION,
        isa(nambus[2],Int) ? (0<nambus[2]<=2 ? nambus[2] : error("fockcouplingnambus error: nambu out of range.")) : ANNIHILATION
        )
end
function expand!(   operators::Operators,term::Pairing,bond::AbstractBond,config::IDFConfig,
                    table::Union{Table,Nothing}=nothing,
                    half::Bool=false,
                    coord::Union{Val{true},Val{false}}=coordpresent
                    )
    argtypes=Tuple{Operators,Term,AbstractBond,IDFConfig,Union{Table,Nothing},Bool,Union{Val{true},Val{false}}}
    invoke(expand!,argtypes,operators,term,bond,config,table,half,coord)
    isa(bond,Bond) && invoke(expand!,argtypes,operators,term,bond|>reverse,config,table,half,coord)
    return operators
end

const hubbard=FockCoupling{4}(spins=(2,2,1,1),nambus=(CREATION,ANNIHILATION,CREATION,ANNIHILATION))
"""
    Hubbard{ST}(id::Symbol,value::Any;
                amplitude::Union{Function,Nothing}=nothing,
                modulate::Union{Function,Bool}=false,
                ) where {ST}

Hubbard term.

Type alias for `Term{statistics,:Hubbard,4,id,V,Int,<:TermCouplings,<:TermAmplitude,<:Union{TermModulate,Nothing}}`.
"""
const Hubbard{statistics,id,V,C<:TermCouplings,A<:TermAmplitude,M<:Union{TermModulate,Nothing}}=Term{statistics,:Hubbard,4,id,V,Int,C,A,M}
function Hubbard{ST}(id::Symbol,value::Any;
                    amplitude::Union{Function,Nothing}=nothing,
                    modulate::Union{Function,Bool}=false,
                    ) where {ST}
    @assert value==value' "Hubbard error: only real values are allowed."
    Term{ST,:Hubbard,4}(id,value,0,couplings=hubbard,amplitude=amplitude,modulate=modulate)
end
abbr(::Type{<:Hubbard})=:hb
isHermitian(::Type{<:Hubbard})=true

const interorbitalinterspin=FockCoupling{4}(orbitals=(@subscript (α,α,β,β) with α<β),spins=(@subscript (σ₁,σ₁,σ₂,σ₂) with σ₁≠σ₂),nambus=(CREATION,ANNIHILATION,CREATION,ANNIHILATION))
"""
    InterOrbitalInterSpin{ST}(  id::Symbol,value::Any;
                                amplitude::Union{Function,Nothing}=nothing,
                                modulate::Union{Function,Bool}=false,
                                ) where {ST}

Interorbital-interspin term.

Type alias for `Term{statistics,:InterOrbitalInterSpin,4,id,V,Int,<:TermCouplings,<:TermAmplitude,<:Union{TermModulate,Nothing}}`.
"""
const InterOrbitalInterSpin{statistics,id,V,C<:TermCouplings,A<:TermAmplitude,M<:Union{TermModulate,Nothing}}=Term{statistics,:InterOrbitalInterSpin,4,id,V,Int,C,A,M}
function InterOrbitalInterSpin{ST}( id::Symbol,value::Any;
                                    amplitude::Union{Function,Nothing}=nothing,
                                    modulate::Union{Function,Bool}=false,
                                    ) where {ST}
    @assert value==value' "InterOrbitalInterSpin error: only real values are allowed."
    Term{ST,:InterOrbitalInterSpin,4}(id,value,0,couplings=interorbitalinterspin,amplitude=amplitude,modulate=modulate)
end
abbr(::Type{<:InterOrbitalInterSpin})=:nons
isHermitian(::Type{<:InterOrbitalInterSpin})=true

const interorbitalintraspin=FockCoupling{4}(orbitals=(@subscript (α,α,β,β) with α<β),nambus=(CREATION,ANNIHILATION,CREATION,ANNIHILATION))
"""
    InterOrbitalIntraSpin{ST}(  id::Symbol,value::Any;
                                amplitude::Union{Function,Nothing}=nothing,
                                modulate::Union{Function,Bool}=false,
                                ) where {ST}

Interorbital-intraspin term.

Type alias for `Term{statistics,:InterOrbitalIntraSpin,4,id,V,Int,<:TermCouplings,<:TermAmplitude,<:Union{TermModulate,Nothing}}`.
"""
const InterOrbitalIntraSpin{statistics,id,V,C<:TermCouplings,A<:TermAmplitude,M<:Union{TermModulate,Nothing}}=Term{statistics,:InterOrbitalIntraSpin,4,id,V,Int,C,A,M}
function InterOrbitalIntraSpin{ST}( id::Symbol,value::Any;
                                    amplitude::Union{Function,Nothing}=nothing,
                                    modulate::Union{Function,Bool}=false,
                                    ) where {ST}
    @assert value==value' "InterOrbitalIntraSpin error: only real values are allowed."
    Term{ST,:InterOrbitalIntraSpin,4}(id,value,0,couplings=interorbitalintraspin,amplitude=amplitude,modulate=modulate)
end
abbr(::Type{<:InterOrbitalIntraSpin})=:noes
isHermitian(::Type{<:InterOrbitalIntraSpin})=true

const spinflip=FockCoupling{4}(orbitals=(@subscript (α,β,α,β) with α<β),spins=(2,1,1,2),nambus=(CREATION,CREATION,ANNIHILATION,ANNIHILATION))
"""
    SpinFlip{ST}(   id::Symbol,value::Any;
                    amplitude::Union{Function,Nothing}=nothing,
                    modulate::Union{Function,Bool}=false,
                    ) where {ST}

Spin-flip term.

Type alias for `Term{statistics,:SpinFlip,4,id,V,Int,<:TermCouplings,<:TermAmplitude,<:Union{TermModulate,Nothing}}`.
"""
const SpinFlip{statistics,id,V,C<:TermCouplings,A<:TermAmplitude,M<:Union{TermModulate,Nothing}}=Term{statistics,:SpinFlip,4,id,V,Int,C,A,M}
function SpinFlip{ST}(  id::Symbol,value::Any;
                        amplitude::Union{Function,Nothing}=nothing,
                        modulate::Union{Function,Bool}=false,
                        ) where {ST}
    Term{ST,:SpinFlip,4}(id,value,0,couplings=spinflip,amplitude=amplitude,modulate=modulate)
end
abbr(::Type{<:SpinFlip})=:sf
isHermitian(::Type{<:SpinFlip})=false

const pairhopping=FockCoupling{4}(orbitals=(@subscript (α,α,β,β) with α<β),spins=(2,1,1,2),nambus=(CREATION,CREATION,ANNIHILATION,ANNIHILATION))
"""
    PairHopping{ST}(    id::Symbol,value::Any;
                        amplitude::Union{Function,Nothing}=nothing,
                        modulate::Union{Function,Bool}=false,
                        ) where {ST}

Pair-hopping term.

Type alias for `Term{statistics,:PairHopping,4,id,V,Int,<:TermCouplings,<:TermAmplitude,<:Union{TermModulate,Nothing}}`.
"""
const PairHopping{statistics,id,V,C<:TermCouplings,A<:TermAmplitude,M<:Union{TermModulate,Nothing}}=Term{statistics,:PairHopping,4,id,V,Int,C,A,M}
function PairHopping{ST}(   id::Symbol,value::Any;
                            amplitude::Union{Function,Nothing}=nothing,
                            modulate::Union{Function,Bool}=false,
                            ) where {ST}
    Term{ST,:PairHopping,4}(id,value,0,couplings=pairhopping,amplitude=amplitude,modulate=modulate)
end
abbr(::Type{<:PairHopping})=:ph
isHermitian(::Type{<:PairHopping})=false

"""
    Coulomb{ST}(    id::Symbol,value::Any,bondkind::Int=1;
                    couplings::Union{Function,Coupling,Couplings,Nothing}=nothing,
                    amplitude::Union{Function,Nothing}=nothing,
                    modulate::Union{Function,Bool}=false,
                    ) where {ST}

Coulomb term.

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

"""
    FFockTerm

Fermionic Fock term types.
"""
const FFockTerm=Union{Onsite{'F'},Hopping{'F'},Pairing{'F'},Hubbard{'F'},InterOrbitalInterSpin{'F'},InterOrbitalIntraSpin{'F'},SpinFlip{'F'},PairHopping{'F'},Coulomb{'F'}}

"""
    BFockTerm

Bosonic Fock term types.
"""
const BFockTerm=Union{Onsite{'B'},Hopping{'B'},Pairing{'B'},Hubbard{'B'},InterOrbitalInterSpin{'B'},InterOrbitalIntraSpin{'B'},SpinFlip{'B'},PairHopping{'B'},Coulomb{'B'}}

"""
    otype(T::Type{<:FFockTerm},I::Type{<:OID})
    otype(T::Type{<:BFockTerm},I::Type{<:OID})

Get the operator type of a Fock term.
"""
otype(T::Type{<:FFockTerm},I::Type{<:OID})=FOperator{T|>valtype,ID{I,T|>rank}}
otype(T::Type{<:BFockTerm},I::Type{<:OID})=BOperator{T|>valtype,ID{I,T|>rank}}

end # module


text:
export inner,
       _abs_sum,
       to_matrix

# computes ‖a^T G‖₁
@inline function _abs_sum(a::AbstractVector, G::AbstractMatrix)
   n, p = size(G)
   N = promote_type(eltype(a), eltype(G))
   abs_sum = zero(N)
   @inbounds for j in 1:p
       aux = zero(N)
       @simd for i in 1:n
           aux += a[i] * G[i, j]
       end
       abs_sum += abs(aux)
   end
   return abs_sum
end

# computes ‖a^T G‖₁ for `a` being a sparse vector
@inline function _abs_sum(a::AbstractSparseVector, G::AbstractMatrix)
   return sum(abs, transpose(a) * G)
end

# computes ‖a^T G‖₁ for `a` having only one nonzero element
@inline function _abs_sum(a::SingleEntryVector, G::AbstractMatrix)
   p = size(G, 2)
   i = a.i
   v = abs(a.v)
   N = promote_type(eltype(a), eltype(G))
   abs_sum = zero(N)
   @inbounds for j in 1:p
       abs_sum += abs(G[i, j])
   end
   abs_sum *= v
   return abs_sum
end

"""
    inner(x::AbstractVector{N}, A::AbstractMatrix{N}, y::AbstractVector{N}
         ) where {N}

Compute the inner product ``xᵀ A y``.

### Input

- `x` -- vector on the left
- `A` -- matrix
- `y` -- vector on the right

### Output

The (scalar) result of the multiplication.
"""
function inner(x::AbstractVector{N}, A::AbstractMatrix{N}, y::AbstractVector{N}
              ) where {N}
    dot(x, A * y)
end

"""
    _vector_type(T)

Return a corresponding vector type with respect to type `T`.

### Input

- `T` -- vector or matrix type

### Output

A vector type that corresponds in some sense (see Notes below) to `T`.

### Notes

- If `T` is a sparse vector or a sparse matrix, the corresponding type is also
  a sparse vector.
- If `T` is a regular vector (i.e. `Vector`) or a regular matrix (i.e. `Matrix`),
  the corresponding type is also a regular vector.
- Otherwise, the corresponding type is a regular vector.
"""
function _vector_type end

"""
    _matrix_type(T)

Return a corresponding matrix type with respect to type `T`.

### Input

- `T` -- vector or matrix type

### Output

A matrix type that corresponds in some sense (see Notes below) to `T`.

### Notes

- If `T` is a sparse vector or a sparse matrix, the corresponding type is also
  a sparse matrix.
- If `T` is a regular vector (i.e. `Vector`) or a regular matrix (i.e. `Matrix`),
  the corresponding type is also a regular matrix.
- Otherwise, the corresponding type is a regular matrix.
"""
function _matrix_type end

_vector_type(::Type{<:AbstractSparseArray{T}}) where T = SparseVector{T, Int}
_vector_type(VT::Type{<:AbstractVector{T}}) where T = VT
_vector_type(::Type{<:AbstractMatrix{T}}) where T = Vector{T}

_matrix_type(::Type{<:AbstractVector{T}}) where T = Matrix{T}
_matrix_type(MT::Type{<:AbstractMatrix{T}}) where T = MT
_matrix_type(::Type{<:AbstractSparseVector{T}}) where T = SparseMatrixCSC{T, Int}

# matrix constructors
_matrix(m, n, MT::Type{<:AbstractMatrix{T}}) where T = Matrix{T}(undef, m, n)
_matrix(m, n, MT::Type{<:SparseMatrixCSC{T}}) where T = spzeros(T, m, n)

"""
    to_matrix(vectors::AbstractVector{VN},
              [m]=length(vectors[1]),
              [MT]=_matrix_type(VN)) where {VN}

### Input

- `vectors` -- list of vectors
- `m`       -- (optional; default: `length(vectors[1])`) number of rows
- `MT`      -- (optional; default: `_matrix_type(VN)`) type of target matrix

### Output

A matrix with the column vectors from `vectors` in the same order.
"""
function to_matrix(vectors::AbstractVector{VN},
                   m=length(vectors[1]),
                   mat_type=_matrix_type(VN)) where {VN}
    n = length(vectors)
    M = _matrix(m, n, mat_type)
    @inbounds for (j, vj) in enumerate(vectors)
        M[:, j] = vj
    end
    return M
end

# no-op
_similar_type(x::AbstractVector) = typeof(x)

function load_copy_finalize_static()

return quote
    _similar_type(x::StaticArrays.StaticArray) = StaticArrays.similar_type(x)
end # quote

end # end load_copy_finalize_static

ext: jl


text:
using Showoff
using Compat
using Compat.Test
using Compat.Dates

@testset "Internals" begin
    @test Showoff.@grisu_ccall(1, 2, 3) === nothing
    @test Showoff.grisu(1.0, Base.Grisu.SHORTEST, 2) == (1, 1, false, Base.Grisu.DIGITS)

    let x = [1.0, Inf, 2.0, NaN]
        @test Showoff.concrete_minimum(x) == 1.0
        @test Showoff.concrete_maximum(x) == 2.0
    end

    @test_throws ArgumentError Showoff.concrete_minimum([])
    @test_throws ArgumentError Showoff.concrete_maximum([])

    let x = [1.12345, 4.5678]
        @test Showoff.plain_precision_heuristic(x) == 5
        @test Showoff.scientific_precision_heuristic(x) == 6
    end
end

@testset "Formatting" begin
    @test Showoff.format_fixed(-10.0, 0) == "-10"
    @test Showoff.format_fixed(0.012345, 3) == "0.012"
    @test Showoff.format_fixed(Inf, 1) == "∞"
    @test Showoff.format_fixed(-Inf, 1) == "-∞"
    @test Showoff.format_fixed(NaN, 1) == "NaN"
    @test Showoff.format_fixed_scientific(0.0, 1, false) == "0"
    @test Showoff.format_fixed_scientific(Inf, 1, false) == "∞"
    @test Showoff.format_fixed_scientific(-Inf, 1, false) == "-∞"
    @test Showoff.format_fixed_scientific(NaN, 1, false) == "NaN"
    @test Showoff.format_fixed_scientific(0.012345678, 4, true) == "12.34568×10⁻³"
    @test Showoff.format_fixed_scientific(0.012345678, 4, false) == "1.234568×10⁻²"
    @test Showoff.format_fixed_scientific(-10.0, 4, false) == "-1.000×10¹"
end

@testset "Showoff" begin
    x = [1.12345, 4.5678]
    @test showoff(x) == ["1.12345", "4.56780"]
    @test showoff([0.0, 50000.0]) == ["0", "5×10⁴"]
    @test showoff(x, :plain) == ["1.12345", "4.56780"]
    @test showoff(x, :scientific) == ["1.12345×10⁰", "4.56780×10⁰"]
    @test showoff(x, :engineering) == ["1.12345×10⁰", "4.56780×10⁰"]
    @test showoff([DateTime("2017-04-11", "yyyy-mm-dd")]) == ["Apr 11, 2017"]
    @test showoff(["a", "b"]) == ["\"a\"", "\"b\""]
    @test_throws ArgumentError showoff(x, :nevergonnagiveyouup)
    @test_throws ArgumentError showoff([Inf, Inf, NaN])
end


text:
"""
    hadi1992_handle_singularity(S)

Perform the sub-algorithm of handling singularity defined in Hadi (1992).

# Arguments
- `S::Array{Float64, 2}`: A covariance matrix.

# Reference
Hadi, Ali S. "Identifying multiple outliers in multivariate data."
Journal of the Royal Statistical Society: Series B (Methodological) 54.3 (1992): 761-771.
 """
function hadi1992_handle_singularity(S::Array{Float64,2})::Array{Float64,2}
    p, _ = size(S)
    eigen_structure = eigen(S)
    values = eigen_structure.values
    vectors = eigen_structure.vectors
    lambda_s = find_minimum_nonzero(values)
    W = zeros(Float64, p, p)
    for i in 1:p
        @inbounds W[i, i] = 1 / max(values[i], lambda_s)
    end
    newS = vectors * W * vectors
    return newS
end

"""
    hadi1992(multivariateData)

Perform Hadi (1992) algorithm for a given multivariate data.

# Arguments
-`multivariate::Array{Float64, 2}`: Multivariate data.

# Examples
```julia-repl
julia> multidata = hcat(hbk.x1, hbk.x2, hbk.x3);

julia> hadi1992(multidata)
Dict{Any,Any} with 3 entries:
  "outliers"              => [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
  "criticial.chi.squared" => 7.81473
  "rth.robust.distance"   => 5.04541
```

# Reference
Hadi, Ali S. "Identifying multiple outliers in multivariate data."
Journal of the Royal Statistical Society: Series B (Methodological) 54.3 (1992): 761-771.
"""
function hadi1992(multivariateData::Array{Float64,2}; alpha=0.05)
    n, p = size(multivariateData)
    h = Int(floor((n + p + 1.0) / 2.0))
    chi_50_quantile = quantile(Chisq(p), 0.50)
    critical_quantile = quantile(Chisq(p), 1 - alpha)
    allindices = collect(1:n)

    # Step 0
    meds = coordinatwisemedians(multivariateData)
    Sm = (1.0 / (n - 1.0)) * (multivariateData .- meds')' * (multivariateData .- meds')
    mah0 = diag(mahalanobisSquaredMatrix(multivariateData, meanvector=meds, covmatrix=Sm))
    ordering_indices_mah0 = sortperm(mah0)
    best_indices_mah0 = ordering_indices_mah0[1:h]
    starting_data = multivariateData[best_indices_mah0, :]

    Cv = coordinatwisemedians(starting_data)
    Sv = (1.0 / (h - 1.0)) * (starting_data .- Cv')' * (starting_data .- Cv')
    mah1 = diag(mahalanobisSquaredMatrix(multivariateData, meanvector=Cv, covmatrix=Sv))
    ordering_indices_mah1 = sortperm(mah1)

    r = p + 1
    basic_subset_indices = []
    basic_subset = []
    sorted_mah1 = []

    while r < n
        cnpr = 1 + (r / (n - p))^2.0
        basic_subset_indices = ordering_indices_mah1[1:r]
        basic_subset = multivariateData[basic_subset_indices, :]
        Cb = applyColumns(mean, basic_subset)
        Sb = cov(basic_subset)

        r += 1
        cfactor = cnpr * sqrt(sort(mah1)[h]) / chi_50_quantile
        if det(cfactor * Sb) == 0
            @info "singular Sb case"
            newSb = hadi1992_handle_singularity(cfactor * Sb)
            mah1 = diag(mahalanobisSquaredMatrix(multivariateData, meanvector=Cb, covmatrix=newSb))
            ordering_indices_mah1 = sortperm(mah1)
            basic_subset_indices = ordering_indices_mah1[1:r]
        else
            mah1 = diag(mahalanobisSquaredMatrix(multivariateData, meanvector=Cb, covmatrix=(cfactor * Sb)))
            ordering_indices_mah1 = sortperm(mah1)
            basic_subset_indices = ordering_indices_mah1[1:r]
        end

        sorted_mah1 = sort(mah1)
        if sorted_mah1[r] >= critical_quantile
            break
        end
    end

    outlierset = setdiff(allindices, basic_subset_indices)

    result = Dict()
    result["outliers"] = sort(outlierset)
    result["criticial.chi.squared"] = critical_quantile
    result["rth.robust.distance"] = sorted_mah1[r - 1]
    return result
end

max_stars_repo_name: akadal/LinRegOutliers

max_issues_repo_name: akadal/LinRegOutliers

max_forks_repo_name: akadal/LinRegOutliers


text:
package ln

import "math"

struct Box
    Min::Float64
    Max::Float64
end

type Box struct {
	Min, Max Vector
}

func BoxForShapes(shapes []Shape) Box {
	if len(shapes) == 0 {
		return Box{}
	}
	box := shapes[0].BoundingBox()
	for _, shape := range shapes {
		box = box.Extend(shape.BoundingBox())
	}
	return box
}

func BoxForTriangles(shapes []*Triangle) Box {
	if len(shapes) == 0 {
		return Box{}
	}
	box := shapes[0].BoundingBox()
	for _, shape := range shapes {
		box = box.Extend(shape.BoundingBox())
	}
	return box
}

func BoxForVectors(vectors []Vector) Box {
	if len(vectors) == 0 {
		return Box{}
	}
	min := vectors[0]
	max := vectors[0]
	for _, v := range vectors {
		min = min.Min(v)
		max = max.Max(v)
	}
	return Box{min, max}
}

func (a Box) Anchor(anchor Vector) Vector {
	return a.Min.Add(a.Size().Mul(anchor))
}

func (a Box) Center() Vector {
	return a.Anchor(Vector{0.5, 0.5, 0.5})
}

func (a Box) Size() Vector {
	return a.Max.Sub(a.Min)
}

func (a Box) Contains(b Vector) bool {
	return a.Min.X <= b.X && a.Max.X >= b.X &&
		a.Min.Y <= b.Y && a.Max.Y >= b.Y &&
		a.Min.Z <= b.Z && a.Max.Z >= b.Z
}

func (a Box) Extend(b Box) Box {
	return Box{a.Min.Min(b.Min), a.Max.Max(b.Max)}
}

func (b *Box) Intersect(r Ray) (float64, float64) {
	x1 := (b.Min.X - r.Origin.X) / r.Direction.X
	y1 := (b.Min.Y - r.Origin.Y) / r.Direction.Y
	z1 := (b.Min.Z - r.Origin.Z) / r.Direction.Z
	x2 := (b.Max.X - r.Origin.X) / r.Direction.X
	y2 := (b.Max.Y - r.Origin.Y) / r.Direction.Y
	z2 := (b.Max.Z - r.Origin.Z) / r.Direction.Z
	if x1 > x2 {
		x1, x2 = x2, x1
	}
	if y1 > y2 {
		y1, y2 = y2, y1
	}
	if z1 > z2 {
		z1, z2 = z2, z1
	}
	t1 := math.Max(math.Max(x1, y1), z1)
	t2 := math.Min(math.Min(x2, y2), z2)
	return t1, t2
}

func (b *Box) Partition(axis Axis, point float64) (left, right bool) {
	switch axis {
	case AxisX:
		left = b.Min.X <= point
		right = b.Max.X >= point
	case AxisY:
		left = b.Min.Y <= point
		right = b.Max.Y >= point
	case AxisZ:
		left = b.Min.Z <= point
		right = b.Max.Z >= point
	}
	return
}


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Constants for the dipole model for the Earth geomagnetic field.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] http://wdc.kugi.kyoto-u.ac.jp/poles/polesexp.html
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# Location of the South pole (which lies in the North hemisphere) along the
# years [°]. This data was obtained from [1].
const _dipole_pole = [1900    +78.7   -68.8
                      1905    +78.7   -68.7
                      1910    +78.7   -68.7
                      1915    +78.6   -68.6
                      1920    +78.6   -68.4
                      1925    +78.6   -68.3
                      1930    +78.6   -68.3
                      1935    +78.6   -68.4
                      1940    +78.5   -68.5
                      1945    +78.5   -68.5
                      1950    +78.5   -68.8
                      1955    +78.5   -69.2
                      1960    +78.6   -69.5
                      1965    +78.6   -69.9
                      1970    +78.7   -70.2
                      1975    +78.8   -70.5
                      1980    +78.9   -70.8
                      1985    +79.0   -70.9
                      1990    +79.2   -71.1
                      1995    +79.4   -71.4
                      2000    +79.6   -71.6
                      2005    +79.8   -71.8
                      2010    +80.1   -72.2
                      2011    +80.1   -72.3
                      2012    +80.2   -72.4
                      2013    +80.3   -72.5
                      2014    +80.3   -72.5
                      2015    +80.4   -72.6
                      2016    +80.4   -72.7
                      2017    +80.5   -72.8
                      2018    +80.5   -73.0
                      2019    +80.6   -73.1
                      2020    +80.6   -73.2]

# Dipole momentum along the years [10²² A.m²]. This data was obtained from [1].
const _dipole_mag = [1900	8.32
                     1905	8.30
                     1910	8.27
                     1915	8.24
                     1920	8.20
                     1925	8.16
                     1930	8.13
                     1935	8.11
                     1940	8.09
                     1945	8.08
                     1950	8.06
                     1955	8.05
                     1960	8.03
                     1965	8.00
                     1970	7.97
                     1975	7.94
                     1980	7.91
                     1985	7.87
                     1990	7.84
                     1995	7.81
                     2000	7.79
                     2005	7.77
                     2010	7.75
                     2011	7.74
                     2012	7.74
                     2013	7.73
                     2014	7.73
                     2015	7.72
                     2016	7.72
                     2017	7.72
                     2018	7.71
                     2019	7.71
                     2020	7.70]

# Interpolations.
const _itp_dipole_lat = extrapolate(
    interpolate( (_dipole_pole[:,1],), _dipole_pole[:,2], Gridded(Linear())),
    Flat())
const _itp_dipole_lon = extrapolate(
    interpolate( (_dipole_pole[:,1],), _dipole_pole[:,3], Gridded(Linear())),
    Flat())
const _itp_dipole_mag = extrapolate(
    interpolate( (_dipole_mag[:,1],), _dipole_mag[:,2], Gridded(Linear())),
    Flat())


text:
# Atmos RHS

#=
Ekman column:
    ∂_t ρ =  ∇ (μ ∇ ρ - w ρ)
    ∂_t ρθ =  ∇ (μ ∇ ρθ - w ρθ)
    ∂_t u =  ∇ (μ ∇ u - w u) + (v - v_g)
    ∂_t v =  ∇ (μ ∇ v - w v) - (u - u_g)
    ∂_t w =  ∇ (μ ∇ w - w w) - g - c_p θ ∂_z Π

where
    Π = (p/p_0)^{R/c_p}

top BCs are insulating and impenetrable:
    ∂_t T = 0
    u = u_g
    v = v_g
    w = 0.0
    ∂_t ρ = 0

and bottom BCs use bulk formulae for surface fluxes of heat and momentum:
    ∂_t ρθ = F₃ = -Ch ρ ||u|| (T_sfc - ρθ / ρ)
    ∂_t u  = F₁ = -Cd u ||u||
    ∂_t v  = F₂ = -Cd v ||u||
    w = 0.0
    ∂_t ρ = 0

We also use this model to accumulate fluxes it calculates
    ∂_t F_accum = -(F₁, F₂, F₃)
=#
function ∑tendencies_atm!(dY, Y, (parameters, T_sfc), t)

    UnPack.@unpack Cd, f, ν, ug, vg, C_p, MSLP, R_d, R_m, C_v, grav = parameters

    (Yc, Yf, _) = Y.x
    (dYc, dYf, dF_sfc) = dY.x
    UnPack.@unpack ρ, u, v, ρθ = Yc
    UnPack.@unpack w = Yf
    dρ = dYc.ρ
    du = dYc.u
    dv = dYc.v
    dρθ = dYc.ρθ
    dw = dYf.w

    # auxiliary calculations of surface variables (can be replaced by functions which call exatrapolation or MO)
    If = Operators.InterpolateC2F(bottom = Operators.Extrapolate(), top = Operators.Extrapolate()) # XXX is this correct? (testing: all elems aside from boundaries change)
    u_1 = parent(If(u))[1]
    v_1 = parent(If(v))[1]
    ρ_1 = parent(If(ρ))[1]
    θ_1 = parent(If(ρθ))[1] / parent(If(ρ))[1]
    u_wind = sqrt(u_1^2 + v_1^2)

    # surface flux calculations
    F_sfc = -calculate_sfc_fluxes_energy(DryBulkFormulaWithRadiation(), parameters, T_sfc[1], θ_1, u_1, v_1, ρ_1, t) # W / m2
    dY.x[3] .= -F_sfc[1]

    # density (centers)
    gradc2f = Operators.GradientC2F()
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))

    If = Operators.InterpolateC2F(bottom = Operators.Extrapolate(), top = Operators.Extrapolate())
    @. dρ = gradf2c(-w * If(ρ)) # Eq. 4.11

    # potential temperature (centers)
    gradc2f = Operators.GradientC2F()
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(F_sfc), top = Operators.SetValue(0.0))

    @. dρθ = gradf2c(-w * If(ρθ) + ν * gradc2f(ρθ / ρ))
    # u velocity (centers)
    gradc2f = Operators.GradientC2F(top = Operators.SetValue(ug))
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(Cd * u_wind * u_1))

    A = Operators.AdvectionC2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))
    @. du = gradf2c(ν * gradc2f(u)) + f * (v - vg) - A(w, u)

    # v velocity (centers)
    gradc2f = Operators.GradientC2F(top = Operators.SetValue(vg))
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(Cd * u_wind * v_1))

    A = Operators.AdvectionC2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))
    @. dv = gradf2c(ν * gradc2f(v)) - f * (u - ug) - A(w, v)

    # w velocity (faces)
    gradc2f = Operators.GradientC2F()
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))

    B = Operators.SetBoundaryOperator(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))
    Π(ρθ) = C_p .* (R_d .* ρθ ./ MSLP) .^ (R_m ./ C_v) ## This should be R_d/ C_p (but needs stabilising...reported to original developer)
    @. dw = B(-(If(ρθ / ρ) * gradc2f(Π(ρθ))) - grav + gradc2f(ν * gradf2c(w)) - w * If(gradf2c(w)))
    return dY

end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
export get_eigenmodes, assumed_modes_solve

"""
    ```
    function get_eigenmodes(
        aircraft::Aircraft{Fg},
        n::Int64 = 6,
        MinvK::Union{AbstractMatrix, Nothing} = nothing;
        fixed_points::Vector{Fi} = Int64[],
        kwargs... # keyword arguments for Arpack
    ) where {Fg, Fi <: Int}
    ```

Get the `n` first eigenmodes for an aircraft,
in order of increasing frequency

Returns a matrix in which each column is an eigenmode,
and a vector of corresponding frequencies
"""
function get_eigenmodes(
    aircraft::Aircraft{Fg},
    n::Int64 = 10,
    MinvK::Union{AbstractMatrix, Nothing} = nothing;
    fixed_points::Vector{Fi} = Int64[],
    kwargs... # keyword arguments for Arpack
) where {Fg, Fi <: Int}

    MinvK = (
        isnothing(MinvK) ?
        get_elastic_matrix(aircraft; fixed_points = fixed_points) :
        MinvK
    )

    N = ndofs(aircraft) ÷ 2
    ϕ, ω = let A = MinvK[(N + 1):end, 1:N]
        #=
        for f in fixed_points
            for idof = (6 * (f - 1) + 1):(6 * f)
                A[idof, :] .= 0.0

                A[idof, idof] = - Fmax
            end
        end
        =#

        dropzeros!(A)

        ω, ϕ = eigs(A; nev = n, which = :SM, kwargs...)

        ϕ = [
            real.(ϕ);
            zeros(Fg, N, n)
        ]

        ϕ, sqrt.(abs.(ω))
    end

    ϕ, ω

end

"""
```
	function assumed_modes_solve(
		acft::Aircraft{Fg},
		x::AbstractVector,
		U∞::Fg;
		n_eig::Int64 = 6,
        transient_aerodynamics::Bool = true,
		kwargs...
	) where Fg
```

Solve aircraft aerostructural problem using Krylov subspace,
Proper Orthogonal Decomposition method

* `acft`: aircraft object
* `q`: variables of interest
* `U∞`: freestream velocity
* `n_eig`: number of eigenvalues to consider
* `transient_aerodynamics`: whether to use transient
aerodynamics for calculations
* `kwargs`: keyword arguments for `state_space`
"""
function assumed_modes_solve(
    acft::Aircraft{Fg},
    x::AbstractVector,
    U∞::Fg;
    n_eig::Int64 = 6,
    transient_aerodynamics::Bool = true,
    kwargs...
) where {Fg}

    fpoints = Int64[]
    for k in kwargs
        if k[1] == :fixed_points
            fpoints = k[2]
        end
    end

    ϕ, _ = get_eigenmodes(acft, n_eig; fixed_points = fpoints)

    r, fcon = state_space(
        acft,
        x,
        U∞;
        kwargs...
    )

    f = (q; Γt = nothing) -> begin
        qd, fcon = state_space(
            acft,
            q,
            U∞;
            ∂Γ!∂t = Γt,
            kwargs...
        )

        qd, fcon.Γ
    end
    df = qd -> begin
        qduals = [
            ForwardDiff.Dual(qq, qqd) for (qq, qqd) in zip(x, qd)
        ]

        qd_d, Γ_d = f(qduals)

        [
            v.partials[1] for v in qd_d
        ], [
            v.partials[1] for v in Γ_d
        ]
    end

    A = Matrix{Fg}(undef, ndofs(acft), n_eig)
    JV = similar(A)

    dΓ!dx = Matrix{Fg}(undef, length(acft.wing_strips), n_eig)

    for i = 1:n_eig
        A[:, i] .= ϕ[:, i]

        #=
        for j = 1:(i-1)
            A[:, i] .-= (A[:, i] ⋅ A[:, j]) .* A[:, j]
        end
        =#

        A[:, i] ./= norm(A[:, i])

        qd, Γ = df(A[:, i])

        JV[:, i] .= qd
        dΓ!dx[:, i] .= Γ
    end

    for i = 1:n_eig
        v = [
            zeros(Fg, ndofs(acft) ÷ 2);
            A[1:(ndofs(acft) ÷ 2), i]
        ]

        #=
        for j = 1:size(A, 2)
            v .-= (v ⋅ A[:, j]) .* A[:, j]
        end
        =#

        v ./= norm(v)

        qd, Γ = df(v)

        A = [A v]
        JV = [JV qd]
        dΓ!dx = [dΓ!dx Γ]
    end

    D = (
        transient_aerodynamics ?
        hcat(
            [
                f(x; Γt = dΓ!dx[:, i])[1] .- r for i = 1:size(dΓ!dx, 2)
            ]...
        ) :
        nothing
    )

    #=
    A, JV, D = map(
        mat -> mat[:, 1:n_eig] .+ mat[:, (n_eig + 1):end] .* 1im,
        [
            A, JV, D
        ]
    )
    =#

    Ainv = _linv_reg(A)

    eig = let (M1, M2) = (
        (
            transient_aerodynamics ?
            Matrix{Fg}(I, size(Ainv, 1), size(Ainv, 1)) .- Ainv * D :
            Matrix{Fg}(I, size(Ainv, 1), size(Ainv, 1))
        ),
        Ainv * JV
    )
        eigen(M2, M1)
    end

    A, (
        ϕ = eig.vectors,
        λ = eig.values,
        JA = JV,
        B = D
    )

end

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
using REPL
using REPL: LineEdit

function setup_keybinds()
    schedule(Task(setup_keybinds_background))
    return
end

function setup_keybinds_background()
    try
        setup_keybinds_impl()
    catch err
        @error(
            "Unexpected error from `setup_keybinds_impl`",
            exception = (err, catch_backtrace()),
        )
    end
end

function setup_keybinds_impl()
    # This is why we need https://github.com/JuliaLang/julia/pull/29896...
    for _ in 1:20
        try
            Base.active_repl.interface.modes[1].keymap_dict
            @goto ok
        catch
        end
        sleep(0.05)
    end
    @warn "Failed to wait for REPL"
    return
    @label ok

    trigger_key = CONFIG.trigger_key
    if trigger_key === nothing
        @debug "`trigger_key` is `nothing`; not defining a shortcut key"
        return
    end
    trigger_key::Char

    repl = Base.active_repl
    repl isa REPL.LineEditREPL || return
    insert_search = function(s, _...)
        if isempty(s) || position(LineEdit.buffer(s)) == 0
            LineEdit.edit_insert(s, "@search")
        else
            LineEdit.edit_insert(s, trigger_key)
        end
    end
    new_keymap = Dict{Any,Any}(trigger_key => insert_search)

    main_mode = repl.interface.modes[1]
    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict,
                                                  new_keymap)
    return
end

precompile(Tuple{typeof(setup_keybinds)})
precompile(Tuple{typeof(setup_keybinds_background)})


text:
import Statistics

"""
    simple_linear_regression(x::AbstractVector, y::AbstractVector)

Simple linear regression - given a set of two-dimensional points (x, y), use
the ordinary least squares method to find the best fit line of the form
y = a + b*x (where a and b are real numbers) and return the tuple (a, b).
"""
function simple_linear_regression(
        x::AbstractVector,
        y::AbstractVector,
        )::Tuple
    if length(x) != length(y)
        error("length(x) != length(y)")
    end
    if length(x) == 0
        error("length(x) == 0")
    end

    x_bar = Statistics.mean(x)
    y_bar = Statistics.mean(y)
    var_x = Statistics.var(x)
    cov_x_y = Statistics.cov(x,y)

    @assert(isfinite(x_bar))
    @assert(isfinite(y_bar))
    @assert(isfinite(var_x))
    @assert(isfinite(cov_x_y))

    coefficient = cov_x_y/var_x

    if isfinite(coefficient)
        intercept = y_bar - coefficient*x_bar
        @debug(
            string("Found best fit line: "),
            intercept,
            coefficient,
            )
    else
        @warn(
            string(
                "The best fit line does not have a finite slope. ",
                "I will ignore this result and will instead return ",
                "intercept = 0 and coefficient = 0",
                )
            )
        intercept = 0
        coefficient = 0
    end

    return intercept, coefficient
end


text:
function upsidedown{T}(tensor6::Array{T,6})
	upsidedowntensor = permutedims(tensor6, [1,2,3,4,6,5])
	return upsidedowntensor
end

function renormalizeTrotter!{T}(simulator::Quantum2dFractalSimulator{T}, dimM::Int)
	lattice = simulator.lattice
	tensorT, tensorTtilde = getTensorT(lattice)
 	tensorPx = getTensorP(lattice,1)
 	tensorPy = getTensorP(lattice,2)
 	tensorQ = getTensorQ(lattice)

 	if isodd(getTrotterCount(simulator))
 		whichisT = 1
		tensorT1 = tensorT
		tensorT2 = tensorTtilde
	else
 		whichisT = 2
		tensorT2 = tensorT
		tensorT1 = tensorTtilde
 	end

	nextTensorT, nextTensorTtilde, projectorUx, projectorUy = renormalizeTensorT_trotter(tensorT1, tensorT2, dimM , whichIsTensorT = whichisT)
	nextTensorTtilde = symmetrizeTensorTrotter(nextTensorTtilde)
 	nextTensorPy = renormalizeTensorPy_trotter(tensorPy, projectorUx, projectorUy, dimM)
 	nextTensorPx = renormalizeTensorPx_trotter(tensorPx, projectorUx, projectorUy, dimM)
 	nextTensorQ = renormalizeTensorQ_trotter(tensorQ, projectorUx, projectorUy, dimM)

	setTensorT!(lattice, nextTensorT)
	setTensorT!(lattice, nextTensorTtilde; tilde = true)
 	setTensorP!(lattice, 1, nextTensorPx)
 	setTensorP!(lattice, 2, nextTensorPy)
 	setTensorQ!(lattice, nextTensorQ)

end

function renormalizeTensorQ_trotter{T}(tensorQ::Array{T,4}, projectorUx::Array{T,3}, projectorUy::Array{T,3}, dimM::Int)
	@tensor begin
		dummy1[x1,z,i,y2,y] :=
		tensorQ[x1,y1,z,i] * projectorUy[y1,y2,y]
		dummy2[z,i,y2,y,x2,x] :=
		dummy1[x1,z,i,y2,y] * projectorUx[x1,x2,x]
		nextTensorQ[x,y,z,zp] :=
		dummy2[z,i,y2,y,x2,x] * tensorQ[x2,y2,i,zp]
	end
	return nextTensorQ
end

function renormalizeTensorPx_trotter{T}(tensorPx::Array{T,5}, projectorUx::Array{T,3}, projectorUy::Array{T,3}, dimM::Int)
	@tensor begin
		dummy1[y1p,s1,z,i,y2,y] :=
		tensorPx[y1,y1p,s1,z,i] * projectorUy[y1,y2,y]
		dummy2[y1p,s1,z,y,y2p,s2,zp] :=
		dummy1[y1p,s1,z,i,y2,y] * tensorPx[y2,y2p,s2,i,zp]
		dummy3[s1,z,y,s2,zp,yp] :=
		dummy2[y1p,s1,z,y,y2p,s2,zp] * projectorUy[y1p,y2p,yp]
		nextTensorPx[y,yp,s,z,zp] :=
		dummy3[s1,z,y,s2,zp,yp] * projectorUx[s1,s2,s]
	end
	return nextTensorPx
end

function renormalizeTensorPy_trotter{T}(tensorPy::Array{T,5}, projectorUx::Array{T,3}, projectorUy::Array{T,3}, dimM::Int)
	@tensor begin
		dummy1[x1,x1p,z,i,s2,s] :=
		tensorPy[x1,x1p,s1,z,i] * projectorUy[s1,s2,s]
 		dummy2[x1,x1p,z,s,x2,x2p,zp] :=
 		dummy1[x1,x1p,z,i,s2,s] * tensorPy[x2,x2p,s2,i,zp]
 		dummy3[x1p,z,s,x2p,zp,x] :=
 		dummy2[x1,x1p,z,s,x2,x2p,zp] * projectorUx[x1,x2,x]
 		nextTensorPy[x,xp,s,z,zp] :=
 		dummy3[x1p,z,s,x2p,zp,x] * projectorUx[x1p,x2p,xp]
 	end
 	return nextTensorPy
end

function renormalizeTensorT_trotter{T}(tensorT1::Array{T,6}, tensorT2::Array{T,6}, dimM::Int ; whichIsTensorT::Int = 1)
	 if whichIsTensorT == 1
		tensorT = tensorT1
		tensorTtilde = tensorT2
	 else
		tensorT = tensorT2
		tensorTtilde = tensorT1
	 end
	tensorU = getTensorU(tensorT, dimM)
	tensorV = getTensorV(tensorT, dimM)
	newTensorT = getNewTensorT_2dQ(tensorT, tensorT, tensorU, tensorV)
	newTensorTtilde = getNewTensorT_2dQ(tensorT1, tensorT2, tensorU, tensorV)
	return newTensorT, newTensorTtilde, tensorU, tensorV
end

#---------------------------------------------------------------------------------------
#---------------------------------------------------------------------------------------
function renormalizeSpace!{T}(simulator::Quantum2dFractalSimulator{T}, dimM::Int)
	lattice = simulator.lattice
	tensorT, tensorTtilde = getTensorT(lattice)
	# println(traceTensorTPeriodic(tensorTtilde), " over ", traceTensorTPeriodic(tensorT)) ###debug
	tensorPx = getTensorP(lattice,1)
	tensorQ = getTensorQ(lattice)

	tensorCore, tensorCoretilde, tensorUy =
	calculateCoreTensors(tensorT, tensorTtilde, dimM)
	tensorCore, normC = normalizeTensor(tensorCore)
	tensorCoretilde = tensorCoretilde ./ normC
	setNorm!(simulator, normC, getIndexOf(simulator, "c"))

	tensorLeg = calculateLegTensor(tensorPx, tensorUy, dimM)
	tensorLeg, normLy = normalizeAndSetNorm!(simulator, tensorLeg, "ly")
  		###### optional(later) extend leg
	tensorCaret = calculateCaretTensor(tensorQ, tensorUy, dimM)
	tensorCaret, normEy = normalizeAndSetNorm!(simulator, tensorCaret, "ey")

	tensorLegCore, tensorLegCoretilde = calculateTensorLegCore(tensorLeg, tensorCore, tensorCoretilde, dimM)
	tensorLegCore, normLc = normalizeAndSetNorm!(simulator, tensorLegCore, "lc")
	tensorCaretCore = calculateTensorCaretCore(tensorCaret, tensorCore, dimM)
	tensorCaretCore, normEc = normalizeAndSetNorm!(simulator, tensorCaretCore, "ec")

	tensorTau = rotateTensorS_xy(tensorLegCore)
	tensorTautilde = rotateTensorS_xy(tensorLegCoretilde)

 	# do the similar things
	# second stage
	# P = getTensorP(lattice,2)
	# Tau -> Core, tensorU
	# P, Q, tensorU -> leg, caret

	tensorP = getTensorP(lattice,2)
	tensorCore, tensorCoretilde, tensorUy =
	calculateCoreTensors(tensorTau, tensorTautilde, dimM, mirror = true)
	tensorCore, normLccl = normalizeAndSetNorm!(simulator, tensorCore, "lccl")
	tensorCoretilde = tensorCoretilde ./ normLccl
	tensorLeg = calculateLegTensor(tensorP, tensorUy, dimM)
	tensorLeg, normLx = normalizeAndSetNorm!(simulator, tensorLeg, "lx")
  		###### optional(later) extend leg
	tensorCaret = calculateCaretTensor(permutedims(tensorQ, [2,1,3,4]), tensorUy, dimM)
	tensorCaret, normEx = normalizeAndSetNorm!(simulator, tensorCaret, "ex")

	# calculateNextPQ:
	# input:	CaretCore, LegCore, leg, caret, tensorUy, dimM
	# output:	nextTensorPx, nextTensorQ
	## CareCore, LegCore from the prev stage
	## CaretCore, LegCore, leg, caret -> Py, Q
	##### 137,153 #noleg#
		# 1. tensorPi0: CaretCore + CoreLeg = eccl
		tensorPi0 = calculateWallCore(tensorCaretCore, mirrorTensorS(tensorLegCore), tensorUy, dimM)
		tensorPi0, normEccl = normalizeAndSetNorm!(simulator, tensorPi0, "eccl")
		# 2. tensorPi1: leg + tensorPi0 = eccll
			# tensorPi1 = calculatePi1(tensorPi0, leg, dimM)
			tensorPi1 = calculateLegWall(tensorLeg, tensorPi0, dimM)
		tensorPi1, normEccll = normalizeAndSetNorm!(simulator, tensorPi1, "eccll")
		# 3. nextTensorPy: tensorPi1 + leg = ecclll
			nextTensorPy = calculateWallLeg(tensorPi1, tensorLeg, dimM)
			nextTensorPy, normPy = normalizeAndSetNorm!(simulator, nextTensorPy, "py")
		# 4. nextTensorQ:  tensorPi1 + caret = ecclle
			nextTensorQ = calculateWallCaret(tensorPi1, tensorCaret, dimM)
			nextTensorQ, normQ = normalizeAndSetNorm!(simulator, nextTensorQ, "q")
	# Core, leg -> LegCore ( attach leg to y index)
	# Core = lccl, leg -> LegCore = lccll
	tensorLegCore, tensorLegCoretilde = calculateTensorLegCore(tensorLeg, tensorCore, tensorCoretilde, dimM)
	tensorLegCore, normLccll = normalizeAndSetNorm!(simulator, tensorLegCore, "lccll")
	tensorLegCoretilde = tensorLegCoretilde ./ normLccll
		## LegCore, leg (attach leg to the yp index) -> rotate -> T
		## LegCore = lccll + leg = lcclll = next tensorT
		tensorT, tensorTtilde = calculateTensorCoreLeg(tensorLegCore, tensorLegCoretilde, tensorLeg, dimM)
		tensorT, normT = normalizeAndSetNorm!(simulator, tensorT, "t")
		tensorTtilde = tensorTtilde ./ normT
		nextTensorT = rotateTensorS_xy(tensorT)
		nextTensorTtilde = rotateTensorS_xy(tensorTtilde)
		# symmetrize tensorTtilde
		nextTensorTtilde = symmetrizeTensorSpatial(nextTensorTtilde)
		## LegCore, caret -> Px (lccll + caret = nextTensorPx = lcclle)
		nextTensorPx = calculateTensorCoreCaret(tensorLegCore, tensorCaret, dimM)
		nextTensorPx, normPx = normalizeAndSetNorm!(simulator, nextTensorPx, "px")

		setTensorT!(lattice, nextTensorT)
		setTensorT!(lattice, nextTensorTtilde; tilde = true)
		setTensorP!(lattice, 1, nextTensorPx)
	setTensorP!(lattice, 2, nextTensorPy)
	setTensorQ!(lattice, nextTensorQ)

# 	# with mirror
# 	tensorPx = getTensorP(lattice,2)
# 	tensorUy = getTensorUy(dimM, "x", rotatedTensorS1, mirror = true)
# 	tensorUz = getTensorUz(dimM, "x", rotatedTensorS1, mirror = true)
# 	tensorCore2 = getTensorCore(rotatedTensorS1, mirrorTensorS(rotatedTensorS1), tensorUy, tensorUz)
# 	tensorCore2tilde = getTensorCore(rotatedTensorS1tilde, mirrorTensorS(rotatedTensorS1tilde), tensorUy, tensorUz)
#
#  	legProjectorZ = getLegProjectorZ(tensorPx, dimM)
#  	tensorLegY = constructLegY(tensorPx, tensorUy, legProjectorZ)
#  		###### optional(later) extend leg
#  	caretProjectorZ = getCaretProjectorZ(tensorQ, dimM)
#  	tensorCaretY = constructCaretY(tensorQ, tensorUy, caretProjectorZ)
#
#  	tensorSProjectorZ = getSProjectorZ(tensorLegY, tensorCore2, dimM)
#  	tensorS2 = constructTensorS(tensorLegY, tensorSProjectorZ, tensorCore2)
#  	tensorS2tilde = constructTensorS(tensorLegY, tensorSProjectorZ, tensorCore2tilde)
#
# 	# 1. T: S2 + leg
# 	# 2. Ttilde: S2tilde + leg
# 		#
# 	# 3. Py: S2 + caret
#
	# from CaretCore, LegCore, newLeg, newCaret
	# -> Py, Q
	# do it without rotating anything
	# 1. 	CaretCore + CoreLeg
# 	# 1. ScS1: Sc + mirror(S1)
# 	# 2. ScS1Leg: ScS1 + leg
# 	# 3. Px: ScS1Leg + leg
# 	# 4. Q: ScS1Leg + caret
#
# 	# rotate
# 	# store them
# 	# T, Ttilde
# 	# Px, Py
# 	# Q
#
# 	# from S and Sc
# 	###### from here
# 	# attach one more leg or caret
# 	# to make next tensor T

end

#---
# rotation on x,y plane

function rotateTensorS_xy{T}(tensorS::Array{T,6})
	return permutedims(tensorS,[3,4,2,1,5,6])
end

function mirrorTensorS{T}(tensorS::Array{T,6})
	return permutedims(tensorS,[2,1,3,4,5,6])
end

# symmetrizeTensorSpatial
function symmetrizeTensorSpatial{T}(tensor6::Array{T,6})
	rotated90Tensor = rotateTensorS_xy(tensor6)
	rotated180Tensor = rotateTensorS_xy(rotated90Tensor)
	rotated270Tensor = rotateTensorS_xy(rotated180Tensor)
	symmetrizedTensor = (tensor6 + rotated90Tensor + rotated180Tensor + rotated270Tensor) ./ 4.0
	return symmetrizedTensor
end

function symmetrizeTensorTrotter{T}(tensor6::Array{T,6})
	upsidedowntensor = upsidedown(tensor6)
	symmetrizedTensor = (tensor6 + upsidedowntensor) ./ 2.0
	return symmetrizedTensor
end

#----------------------------------------------------------------------
#----------------------------------------------------------------------
# functions mainly for the second part
#---
# functions for
# WallCaret
function calculateWallCaret{T}(tensorWall::Array{T,5}, tensorCaret::Array{T,3}, dimM::Int)
	projectorWallCaret = getWallCaretProjectorZ(tensorWall, tensorCaret, dimM)
	tensorWallCaret = constructTensorWallCaret(tensorWall, projectorWallCaret, tensorCaret)
	return tensorWallCaret
end

function getWallCaretProjectorZ{T}(tensorWall::Array{T,5}, tensorCaret::Array{T,3}, dimM::Int)
	tenmatMMd = getTenmatMMd_WallCaret(tensorWall, tensorCaret)
	projectorWallCaret, lambdaVector = getProjectorFromMMd(tenmatMMd, dimM)
	return projectorWallCaret
end

function getTenmatMMd_WallCaret{T}(tensorWall::Array{T,5}, tensorCaret::Array{T,3})
	@tensor begin
		dummy1[i,z1,j,a1] := tensorWall[x,i,s,z1,z1p] * tensorWall[x,j,s,a1,z1p]
		dummy2[i,z2,j,a2] := tensorCaret[i,z2,z2p] * tensorCaret[j,a2,z2p]
		tensorMMd[z1,z2,a1,a2] := dummy1[i,z1,j,a1] * dummy2[i,z2,j,a2]
	end
	tenmatMMd =tensor2tenmat(tensorMMd,[1,2],[3,4])
	return tenmatMMd
end

function constructTensorWallCaret{T}(tensorWall::Array{T,5}, projectorWallCaret::Array{T,3}, tensorCaret::Array{T,3})
	@tensor begin
		dummy1[i,z2p,z1,z] :=
		tensorCaret[i,z2,z2p] * projectorWallCaret[z1,z2,z]
		dummy2[i,z1,z,z1p,zp] :=
		dummy1[i,z2p,z1,z] * projectorWallCaret[z1p,z2p,zp]
		tensorWallCaret[x,y,z,zp] :=
		dummy2[i,z1,z,z1p,zp] * tensorWall[x,i,y,z1,z1p]
	end
	return tensorWallCaret
end

#---
# functions for
# WallLeg
function calculateWallLeg{T}(tensorWall::Array{T,5}, tensorLeg::Array{T,4}, dimM::Int)
	tensorWallLeg = calculateLegWall(permutedims(tensorLeg,[2,1,3,4]), permutedims(tensorWall,[2,1,3,4,5]), dimM)
	return tensorWallLeg
end

#---
# functions for
# LegWall

function calculateLegWall{T}(tensorLeg::Array{T,4}, tensorWall::Array{T,5}, dimM::Int)
	projectorLegWall = getLegWallProjectorZ(tensorLeg, tensorWall, dimM)
	tensorLegWall = constructTensorLegWall(tensorLeg, projectorLegWall, tensorWall)
	return tensorLegWall
end

function getLegWallProjectorZ{T}(tensorLeg::Array{T,4}, tensorWall::Array{T,5}, dimM::Int)
	tenmatMMd = getTenmatMMd_LegWall(tensorLeg, tensorWall)
	projectorLegWall, lambdaVector = getProjectorFromMMd(tenmatMMd, dimM)
	return projectorLegWall
end

function getTenmatMMd_LegWall{T}(tensorLeg::Array{T,4}, tensorWall::Array{T,5})
	@tensor begin
		dummy1[i,j,z1,a1] := tensorLeg[x,i,z1,z1p] * tensorLeg[x,j,a1,z1p]
		dummy2[i,z2,j,a2] := tensorWall[i,xp,s,z2,z2p] * tensorWall[j,xp,s,a2,z2p]
		tensorMMd[z1,z2,a1,a2] := dummy1[i,j,z1,a1] * dummy2[i,z2,j,a2]
	end
	tenmatMMd =tensor2tenmat(tensorMMd,[1,2],[3,4])
	return tenmatMMd
end

function constructTensorLegWall{T}(tensorLeg::Array{T,4}, projectorLegWall::Array{T,3}, tensorWall::Array{T,5})
	@tensor begin
		dummy1[x,i,z1p,z2,z] :=
		tensorLeg[x,i,z1,z1p] * projectorLegWall[z1,z2,z]
		dummy2[x,i,z2,z,z2p,zp] :=
		dummy1[x,i,z1p,z2,z] * projectorLegWall[z1p,z2p,zp]
		tensorLegWall[x,xp,s,z,zp] :=
		dummy2[x,i,z2,z,z2p,zp] * tensorWall[i,xp,s,z2,z2p]
	end
	return tensorLegWall
end

#---
# functions for
# WallCore

function calculateWallCore{T}(tensorWall::Array{T,5}, tensorCore::Array{T,6}, tensorUx::Array{T,3}, dimM::Int)
		projectorWallCore = getWallCoreProjectorZ(tensorWall, tensorCore, dimM)
		tensorWallCore = constructTensorWallCore(tensorWall, tensorUx, projectorWallCore, tensorCore)
		return tensorWallCore
end

function getWallCoreProjectorZ{T}(tensorWall::Array{T,5}, tensorCore::Array{T,6}, dimM::Int)
	tenmatMMd = getTenmatMMd_WallCore(tensorWall, tensorCore)
	projectorWallCore, lambdaVector = getProjectorFromMMd(tenmatMMd, dimM)
	return projectorWallCore
end

function getTenmatMMd_WallCore{T}(tensorWall::Array{T,5}, tensorCore::Array{T,6})
	@tensor begin
		dummy1[i,z1,j,a1] :=
		tensorWall[x1,x1p,i,z1,z1p] * tensorWall[x1,x1p,j,a1,z1p]
		dummy2[i,z2,j,a2] :=
		tensorCore[x2,x2p,i,y,z2,z2p] * tensorCore[x2,x2p,j,y,a2,z2p]
		tensorMMd[z1,z2,a1,a2] :=
		dummy1[i,z1,j,a1] * dummy2[i,z2,j,a2]
	end
	tenmatMMd =tensor2tenmat(tensorMMd,[1,2],[3,4])
	return tenmatMMd
end

function constructTensorWallCore{T}(tensorWall::Array{T,5}, tensorUx::Array{T,3}, projectorWallCore::Array{T,3}, tensorCore::Array{T,6})
	@tensor begin
		dummy1[x1,x1p,i,z1p,z2,z] :=
		tensorWall[x1,x1p,i,z1,z1p] * projectorWallCore[z1,z2,z]
		dummy2[x1,x1p,i,z2,z,z2p,zp] :=
		dummy1[x1,x1p,i,z1p,z2,z] * projectorWallCore[z1p,z2p,zp]
		dummy3[x1,x1p,z,zp,x2,x2p,s] :=
		dummy2[x1,x1p,i,z2,z,z2p,zp] * tensorCore[x2,x2p,i,s,z2,z2p]
		dummy4[x1p,z,zp,x2p,s,x] :=
		dummy3[x1,x1p,z,zp,x2,x2p,s] * tensorUx[x1,x2,x]
		tensorWallCore[x,xp,s,z,zp] :=
		dummy4[x1p,z,zp,x2p,s,x] * tensorUx[x1p,x2p,xp]
	end
	return tensorWallCore
end

#----------------------------------------------------------------------
# functions mainly for the first part

#---
# calculate
# CoreCaret (attach Caret to yp index)
function calculateTensorCoreCaret{T}(tensorCore::Array{T,6}, tensorCaret::Array{T,3}, dimM::Int)
	tensorCoreCaret = calculateTensorCaretCore(tensorCaret, permutedims(tensorCore,[1,2,4,3,5,6]), dimM)
	return tensorCoreCaret
end

#---
# get functions for
# CaretCore (attach Caret to y index)

function calculateTensorCaretCore{T}(tensorCaret::Array{T,3}, tensorCore::Array{T,6}, dimM::Int)
  	tensorScProjectorZ = getScProjectorZ(tensorCaret, tensorCore, dimM)
  	tensorSc = constructTensorSc(tensorCaret, tensorScProjectorZ, tensorCore)
	return tensorSc
end

function getScProjectorZ{T}(tensorCaretY::Array{T,3}, tensorCore::Array{T,6}, dimM)
	tenmatMMd = getTenmatMMd_caretCore(tensorCaretY, tensorCore)
	tensorScProjectorZ, lambdaVector = getProjectorFromMMd(tenmatMMd, dimM)
	return tensorScProjectorZ
end

function getTenmatMMd_caretCore{T}(tensorCaretY::Array{T,3}, tensorCore::Array{T,6})
	@tensor begin
		dummy1[i,z2,j,a2] :=
		tensorCore[x,xp,i,yp,z2,z2p] *
		tensorCore[x,xp,j,yp,a2,z2p]
		dummy2[z1,z1p,z2,j,a2] :=
		dummy1[i,z2,j,a2] *
		tensorCaretY[i,z1,z1p]
		tensorMMd[z1,z2,a1,a2] :=
		dummy2[z1,z1p,z2,j,a2] *
		tensorCaretY[j,a1,z1p]
	end
	tenmatMMd =tensor2tenmat(tensorMMd,[1,2],[3,4])
	return tenmatMMd
end

function constructTensorSc{T}(tensorCaretY::Array{T,3}, tensorScProjectorZ::Array{T,3}, tensorCore::Array{T,6})
	@tensor begin
		dummy1[x,xp,i,s,z2p,z1,z] :=
		tensorCore[x,xp,i,s,z2,z2p] *
		tensorScProjectorZ[z1,z2,z]
		dummy2[x,xp,s,z2p,z,z1p] :=
		dummy1[x,xp,i,s,z2p,z1,z] *
		tensorCaretY[i,z1,z1p]
		tensorSc[x,xp,s,z,zp] :=
		dummy2[x,xp,s,z2p,z,z1p] *
		tensorScProjectorZ[z1p,z2p,zp]
	end
	return tensorSc
end

#---
# calculate
# CoreLeg ( attach leg to the yp index )

function calculateTensorCoreLeg{T}(tensorCore::Array{T,6}, tensorCoretilde::Array{T,6}, tensorLeg::Array{T,4}, dimM::Int)
	tensorS, tensorStilde = calculateTensorLegCore(permutedims(tensorLeg,[2,1,3,4]), permutedims(tensorCore,[1,2,4,3,5,6]), permutedims(tensorCoretilde,[1,2,4,3,5,6]), dimM)
	return permutedims(tensorS,[1,2,4,3,5,6]), permutedims(tensorStilde,[1,2,4,3,5,6])
end

#---
# get functions for
# LegCore

function calculateTensorLegCore{T}(tensorLeg::Array{T,4}, tensorCore::Array{T,6}, tensorCoretilde::Array{T,6}, dimM::Int)
	# attach leg to the y index
  	tensorSProjectorZ = getSProjectorZ(tensorLeg, tensorCore, dimM)
  	tensorS = constructTensorS(tensorLeg, tensorSProjectorZ, tensorCore)
  	tensorStilde = constructTensorS(tensorLeg, tensorSProjectorZ, tensorCoretilde)
	return tensorS, tensorStilde
end

function getSProjectorZ{T}(tensorLegY::Array{T,4}, tensorCore::Array{T,6}, dimM::Int)
	tenmatMMd = getTenmatMMd_legCore(tensorLegY, tensorCore)
	tensorSProjectorZ,lambdaVector = getProjectorFromMMd(tenmatMMd, dimM)
	return tensorSProjectorZ
end

function getTenmatMMd_legCore{T}(tensorLegY::Array{T,4}, tensorCore::Array{T,6})
	@tensor begin
		dummy1[j,z2,k,a2] :=
		tensorCore[x,xp,j,yp,z2,z2p] *
		tensorCore[x,xp,k,yp,a2,z2p]
		dummy2[j,z1,k,a1] :=
		tensorLegY[y,j,z1,z1p] *
		tensorLegY[y,k,a1,z1p]
		tensorMMd[z1,z2,a1,a2] :=
		dummy1[j,z2,k,a2] *
		dummy2[j,z1,k,a1]
	end
	tenmatMMd =tensor2tenmat(tensorMMd,[1,2],[3,4])
	return tenmatMMd
end

function constructTensorS{T}(tensorLegY::Array{T,4}, tensorSProjectZ::Array{T,3}, tensorCore::Array{T,6})
	@tensor begin
		dummy1[x,xp,j,yp,z2p,z1,z] :=
		tensorCore[x,xp,j,yp,z2,z2p] *
		tensorSProjectZ[z1,z2,z]
		dummy2[x,xp,yp,z2p,z,y,z1p] :=
		dummy1[x,xp,j,yp,z2p,z1,z] *
		tensorLegY[y,j,z1,z1p]
		tensorS[x,xp,y,yp,z,zp] :=
		dummy2[x,xp,yp,z2p,z,y,z1p] *
		tensorSProjectZ[z1p,z2p,zp]
	end
	return tensorS
end

#---
# get functions for
# caret

function calculateCaretTensor{T}(tensorQ::Array{T,4}, tensorUy::Array{T,3}, dimM::Int)
  	caretProjectorZ = getCaretProjectorZ(tensorQ, dimM)
	tensorCaretY = constructCaretY(tensorQ, tensorUy, caretProjectorZ)
	return tensorCaretY
end

function getCaretProjectorZ{T}(tensorQ::Array{T,4}, dimM::Int)
	tenmatMMd = getTenmatMMd_caret(tensorQ)
	caretProjectorZ, lambdaVector = getProjectorFromMMd(tenmatMMd, dimM)
	return caretProjectorZ
end

function getTenmatMMd_caret{T}(tensorQ::Array{T,4})
	@tensor begin
		dummy1[z1,a1,s,t] :=
		tensorQ[s,y1,z1,z1p] *
		tensorQ[t,y1,a1,z1p]
		dummy2[z2,a2,s,t] :=
		tensorQ[s,y2,z2,z2p] *
		tensorQ[t,y2,a2,z2p]
		tensorMMd[z1,z2,a1,a2] :=
		dummy1[z1,a1,s,t] *
		dummy2[z2,a2,s,t]
	end
	tenmatMMd =tensor2tenmat(tensorMMd,[1,2],[3,4])
	return tenmatMMd
end

function constructCaretY{T}(tensorQ::Array{T,4}, tensorUy::Array{T,3}, caretProjectorZ::Array{T,3})
	@tensor begin
		dummy1[s,y1,z1p,z2,z] :=
		 tensorQ[s,y1,z1,z1p] *
		 caretProjectorZ[z1,z2,z]
		dummy2[s,y2,z2,z1p,zp] :=
		 tensorQ[s,y2,z2,z2p] *
		 caretProjectorZ[z1p,z2p,zp]
		dummy3[y1,y2,z,zp] :=
		 dummy1[s,y1,z1p,z2,z] *
		 dummy2[s,y2,z2,z1p,zp]
		tensorCaretY[y,z,zp] :=
		 dummy3[y1,y2,z,zp] *
		 tensorUy[y1,y2,y]
	end
	return tensorCaretY
end

#---
# get functions for
# leg

######

function calculateLegTensor{T}(tensorPx::Array{T,5}, tensorUy::Array{T,3}, dimM::Int)
  	legProjectorZ = getLegProjectorZ(tensorPx, dimM)
  	tensorLegY = constructLegY(tensorPx, tensorUy, legProjectorZ)
	return tensorLegY
end

function getLegProjectorZ{T}(tensorPx::Array{T,5}, dimM::Int)
	tenmatMMd = getTenmatMMd_leg(tensorPx)
	legProjectorZ,lambdaVector = getProjectorFromMMd(tenmatMMd, dimM)
	return legProjectorZ
end

function getTenmatMMd_leg{T}(tensorPx::Array{T,5})
	@tensor begin
		dummy1[z1,a1,s,t] :=
		tensorPx[y1,y1p,s,z1,z1p] *
		tensorPx[y1,y1p,t,a1,z1p]
		dummy2[z2,a2,s,t] :=
		tensorPx[y2p,y2,s,z2,z2p] *
		tensorPx[y2p,y2,t,a2,z2p]
		tensorMMd[z1,z2,a1,a2] :=
		dummy1[z1,a1,s,t] *
		dummy2[z2,a2,s,t]
	end
	tenmatMMd =tensor2tenmat(tensorMMd,[1,2],[3,4])
	return tenmatMMd
end

function constructLegY{T}(tensorPx::Array{T,5}, tensorUy::Array{T,3}, legProjectorZ::Array{T,3})
	@tensor begin
		dummy1[y1p,s,z1,z1p,y2,y] :=
		tensorPx[y1,y1p,s,z1,z1p] *
		tensorUy[y1,y2,y]
		dummy2[y2,s,z2,z2p,y1p,yp] :=
		tensorPx[y2p,y2,s,z2,z2p] *
		tensorUy[y1p,y2p,yp]
		dummy3[z1,z1p,y,z2,z2p,yp] :=
		dummy1[y1p,s,z1,z1p,y2,y] *
		dummy2[y2,s,z2,z2p,y1p,yp]
		dummy4[z1p,y,z2p,yp,z] :=
		dummy3[z1,z1p,y,z2,z2p,yp] *
		legProjectorZ[z1,z2,z]
		tensorLegY[y,yp,z,zp] :=
		dummy4[z1p,y,z2p,yp,z] *
		legProjectorZ[z1p,z2p,zp]
	end
	return tensorLegY
end

#---
# functions for Core Tensor

function calculateCoreTensors{T}(tensorT::Array{T,6}, tensorTtilde::Array{T,6}, dimM::Int; mirror::Bool = false)
 	tensorUy = getTensorUy(dimM, "x", tensorT, mirror = mirror)
 	tensorUz = getTensorUz(dimM, "x", tensorT, mirror = mirror)
	if mirror
		tensorT2 = mirrorTensorS(tensorT)
		tensorT2tilde = mirrorTensorS(tensorTtilde)
	else
		tensorT2 = tensorT
		tensorT2tilde = tensorTtilde
	end
  	tensorCore = getTensorCore(tensorT, tensorT2, tensorUy, tensorUz)
  	tensorCoretilde = getTensorCore(tensorT, tensorT2tilde, tensorUy, tensorUz)
	return tensorCore, tensorCoretilde, tensorUy
end

function getTensorUy{T}(dimM::Int, axis::AbstractString, tensorT::Array{T,6}; mirror::Bool = false)
	if axis == "x"
		if mirror
			tenmatMMd = getTenmatMMd_CoreCore(tensorT, mirrorTensorS(tensorT))
		else
			tenmatMMd = getTenmatMMd_CoreCore(tensorT)
		end
		tensorUy, lambdaVector = getProjectorFromMMd(tenmatMMd, dimM)
		return tensorUy
	end
end

function getTensorUz{T}(dimM::Int, axis::AbstractString, tensorT::Array{T,6}; mirror::Bool = false)
	tensorUz = getTensorUy(dimM, axis, permutedims(tensorT,[1,2,5,6,4,3]), mirror = mirror)
	return tensorUz
end

function getTenmatMMd_CoreCore{T}(tensorT::Array{T,6}, tensorT2::Array{T,6} = tensorT)
	@tensor begin
		dummy1[i,y1,j,a1] :=
		tensorT[x,i,y1,y1p,z1,z1p] *
		tensorT[x,j,a1,y1p,z1,z1p]
		dummy2[i,y2,j,a2] :=
		tensorT2[i,xp,y2,y2p,z2,z2p] *
		tensorT2[j,xp,a2,y2p,z2,z2p]
		tensorMMd[y1,y2,a1,a2] :=
		dummy1[i,y1,j,a1] *
		dummy2[i,y2,j,a2]
	end
	tenmatMMd = tensor2tenmat(tensorMMd,[1,2],[3,4])
	return tenmatMMd
end

function getTensorCore{T}(tensorT1::Array{T,6}, tensorT2::Array{T,6}, tensorUy::Array{T,3}, tensorUz::Array{T,3})
	newTensorT = getNewTensorT_2dQ(permutedims(tensorT1,[3,4,5,6,1,2]),permutedims(tensorT2,[3,4,5,6,1,2]), tensorUy, tensorUz)
	# from:	 simulator_quantum_2d_square_renormalize.jl
	return permutedims(newTensorT, [5,6,1,2,3,4])
end

#---
# functions to creat projector
function getProjectorFromMMd{T}(matrix::Array{T,2},dimM::Int)
	tensorU, lambdaVector, tensorUd = svd(matrix)
	truncatedU = truncMatrixU(tensorU,dimM)
	tensorU = matU2tenU(truncatedU)
	return tensorU, lambdaVector
end

function getProjectorFromMMd{T}(tenmatMMd::Tenmat{T}, dimM::Int)
	projectorMatrix = getMatrix(tenmatMMd)
	indicesTuple = tensorSize(tenmatMMd)[1:2]
	tensorU, lambdaVector, tensorUd = svd(projectorMatrix)
	truncatedU = truncMatrixU(tensorU,dimM)
	tensorU = matU2tenU_asym(truncatedU, indicesTuple)
	return tensorU, lambdaVector
end

function matU2tenU_asym(matrixU::Matrix, indicesTuple::Tuple{Int,Int})
	sizeTuple = (indicesTuple[1], indicesTuple[2], size(matrixU,2))
	tenmatU = Tenmat(matrixU,[1,2],[3],sizeTuple)
	tensorU = tenmat2tensor(tenmatU)
	return tensorU
end


text:
"""
arrow_read(file::String)::DataFrame

Read contents from a arrow file into a DataFrame

# Arguments
- `file::String`: path to the input file (arrow file)

# Example

    df = arrow_read("data.arrow")

"""
function arrow_read(file::String)::DataFrame
    DataFrame(Arrow.Table(file))
end

"""
    arrow_write(df::DataFrame,file::String)::Nothing

Write DataFrame to a arrow file.

# Arguments
- `df::DataFrame`: DataFrame
- `file::String`: filepath to which object shall be written

# Example

    arrow_write(df, "data.arrow")

"""
function arrow_write(df::DataFrame,file::String)::String
    if isfile(file)
        rm(file)
    end
    return Arrow.write(file,df)
end

max_stars_repo_name: mrchaos/fxrobot

max_issues_repo_name: mrchaos/fxrobot

max_forks_repo_name: mrchaos/fxrobot


text:
retrieve(c::Client, ::Type{VoiceRegion}, g::AbstractGuild) = get_guild_voice_regions(c, g.id)
retrieve(c::Client, ::Type{VoiceRegion}) = list_voice_regions(c)