text:
"""
`module JAC.SphericalTensor`
    ... a submodel of JAC that help deal with and simplify composed (single- and many-electron) spherical tensor operators,
        matrix elements and spherical amplitudes into a reduced matrix-element representation.
"""
module  SphericalTensor

    using   SymEngine,  ..AngularMomentum, ..Basics,  ..Defaults,  ..RacahAlgebra

    export  SphericalState, SphericalMatrixElement

    """
    `abstract type SphericalTensor.AbstractSphericalTensor`
        ... defines an abstract type to comprise various spherical (electron) operators and spherical fields.
    """
    abstract type  AbstractSphericalTensor                          end

    """
    `abstract type SphericalTensor.AbstractSphericalOperator`
        ... defines an abstract type to comprise various spherical (electron) operators.

        + CkOperator                ... to represent a C^(K) tensor operator.
        + CoulombOperator           ... to represent a V^(0, Coulomb) tensor operator.
        + DipoleOperator            ... to represent a D^(1) dipole operator.
        + MultipoleOperator         ... to represent a O^(L, multipole) operator.
        + TkOperator                ... to represent a general T^(K) tensor operator.
    """
    abstract type  AbstractSphericalOperator  <:  AbstractSphericalTensor  end

    """
    `struct  SphericalTensor.CkOperator  <: AbstractSphericalOperator`   ... to represent a C^(K) spherical tensor operator.

        + star          ::Bool      ... defines the complex-conjugate form, if true.
        + k             ::Basic     ... Tensor rank.
    """
    struct   CkOperator  <: AbstractSphericalOperator
        star            ::Bool
        k               ::Basic
    end

    """
    `SphericalTensor.CkOperator(k::Basic)`  ... constructor for just setting the rank (star = false).
    """
    function CkOperator(k::Basic)
    	CkOperator(false, k)
    end

    function Base.string(op::CkOperator)
        sa = "C^(" * string(op.k);        if op.star   sa = sa * "*"   end;      sa = sa * ")"
        return( sa )
    end

    function Base.show(io::IO, op::CkOperator)
        sa = string(op) * " spherical tensor operator.";       print(io, sa, "\n")
    end

    """
    `struct  SphericalTensor.CoulombOperator  <: AbstractSphericalOperator`   ... to represent a V^(0, Coulomb) spherical tensor operator.

        + star          ::Bool      ... defines the complex-conjugate form, if true.
    """
    struct   CoulombOperator  <: AbstractSphericalOperator
        star            ::Bool
    end

    """
    `SphericalTensor.CoulombOperator()`  ... constructor for just setting the dipole operator (star = false).
    """
    function CoulombOperator()
    	CoulombOperator(false)
    end

    function Base.string(op::CoulombOperator)
        sa = "V^(0, Coulomb";        if op.star   sa = sa * "*"   end;      sa = sa * ")"
        return( sa )
    end

    function Base.show(io::IO, op::CoulombOperator)
        sa = string(op) * " rank-0 Coulomb operator.";       print(io, sa, "\n")
    end

    """
    `struct  SphericalTensor.DipoleOperator  <: AbstractSphericalOperator`   ... to represent a D^(1) spherical tensor operator.

        + star          ::Bool      ... defines the complex-conjugate form, if true.
    """
    struct   DipoleOperator  <: AbstractSphericalOperator
        star            ::Bool
    end

    """
    `SphericalTensor.DipoleOperator()`  ... constructor for just setting the dipole operator (star = false).
    """
    function DipoleOperator()
    	DipoleOperator(false)
    end

    function Base.string(op::DipoleOperator)
        sa = "D^(1";        if op.star   sa = sa * "*"   end;      sa = sa * ")"
        return( sa )
    end

    function Base.show(io::IO, op::DipoleOperator)
        sa = string(op) * " spherical dipole operator.";       print(io, sa, "\n")
    end

    """
    `struct  SphericalTensor.TkOperator  <: AbstractSphericalOperator`   ... to represent a general rank-K T^(K) spherical tensor operator.

        + star          ::Bool      ... defines the complex-conjugate form, if true.
        + k             ::Basic     ... Tensor rank.
    """
    struct   TkOperator  <: AbstractSphericalOperator
        star            ::Bool
        k               ::Basic
    end

    """
    `SphericalTensor.TkOperator(k::Basic)`  ... constructor for just setting the rank (star = false).
    """
    function TkOperator(k::Basic)
    	TkOperator(false, k)
    end

    function Base.string(op::TkOperator)
        sa = "T^(" * string(op.k);        if op.star   sa = sa * "*"   end;      sa = sa * ")"
        return( sa )
    end

    function Base.show(io::IO, op::TkOperator)
        sa = string(op) * " spherical tensor operator.";       print(io, sa, "\n")
    end

    """
    `abstract type SphericalTensor.AbstractSphericalConstant  <:  AbstractSphericalTensor`
        ... defines an abstract type to comprise various (parametric-dependent) scalar constants, such as field amplitudes, etc.

        + ScalarConstant            ... to represent a scalar constant to comprise all all factors that do not
                                        affect the tensorial structure of matrix elements and amplitues.
    """
    abstract type  AbstractSphericalConstant   <:  AbstractSphericalTensor  end

    """
    `struct  SphericalTensor.ScalarConstant  <: AbstractSphericalConstant`
        ... to represent a scalar constant which does not affect the tensorial evaluation of matrix elements.

        + star          ::Bool      ... defines the complex-conjugate form, if true.
        + symbol        ::Basic     ... symbol of the constant.
    """
    struct   ScalarConstant  <: AbstractSphericalConstant
        star            ::Bool
        symbol          ::Basic
    end

    function Base.string(op::ScalarConstant)
        sa = string(op.symbol);        if op.star   sa = sa * "^(*)"   end
        return( sa )
    end

    function Base.show(io::IO, op::ScalarConstant)
        sa = string(op) * " scalar constant.";       print(io, sa, "\n")
    end

    """
    `abstract type SphericalTensor.AbstractSphericalField  <:  AbstractSphericalTensor`
        ... defines an abstract type to comprise various (external and electron-independent) field operators.

        + UVector                   ... to represent a u^(1) polarization (unit) vector.
        + TwistedU                  ... to represent a u^(1, twisted) unit vector.
    """
    abstract type  AbstractSphericalField   <:  AbstractSphericalTensor  end

    """
    `struct  SphericalTensor.UVector  <: AbstractSphericalField`   ... to represent a u^(1) polarization (unit) vector.

        + star          ::Bool      ... defines the complex-conjugate form, if true.
    """
    struct   UVector  <: AbstractSphericalField
        star            ::Bool
    end

    """
    `SphericalTensor.UVector()`  ... constructor for just setting the vector (star = false).
    """
    function UVector()
    	UVector(false)
    end

    function Base.string(op::UVector)
        sa = "u^(1";        if op.star   sa = sa * "*"   end;      sa = sa * ")"
        return( sa )
    end

    function Base.show(io::IO, op::UVector)
        sa = string(op) * " spherical polarization (unit) vector.";       print(io, sa, "\n")
    end

    """
    `struct  SphericalTensor.TensorProduct  <:  AbstractSphericalTensor`
        ... defines the coupled tensor product of two spherical tensors [T1^(k1) * T2^(k2)]^(k).

        + star          ::Bool      ... defines the complex-conjugate form, if true.
        + k             ::Basic     ... (total) tensor rank k.
        + tensor1       ::AbstractSphericalTensor  ... spherical tensor T1.
        + tensor2       ::AbstractSphericalTensor  ... spherical tensor T2.
    """
    struct   TensorProduct  <:  AbstractSphericalTensor
        star            ::Bool
        k               ::Basic
        tensor1         ::AbstractSphericalTensor
        tensor2         ::AbstractSphericalTensor
    end

    """
    `SphericalTensor.TensorProduct(k::Basic, tensor1::AbstractSphericalTensor, tensor2::AbstractSphericalTensor)`
        ... constructor for just setting the tensor product (star = false).
    """
    function TensorProduct(k::Basic, tensor1::AbstractSphericalTensor, tensor2::AbstractSphericalTensor)
    	TensorProduct(false, k, tensor1, tensor2)
    end

    function Base.string(op::TensorProduct)
        sa = "[" * string(op.tensor1) * " * " * string(op.tensor2) * "]^(" * string(op.k);        if op.star   sa = sa * "*"   end;      sa = sa * ")"
        return( sa )
    end

    function Base.show(io::IO, op::TensorProduct)
        sa = string(op) * " spherical tensor product.";       print(io, sa, "\n")
    end

    """
    `abstract type SphericalTensor.AbstractSphericalTensorComp`
        ... defines an abstract type to comprise various spherical (electron) operators and spherical field tensor components.
    """
    abstract type  AbstractSphericalTensorComp     end

    """
    `struct  SphericalTensor.SphericalOperatorComp  <: AbstractSphericalTensorComp`   ... to represent a spherical tensor component.

        + tensor        ::AbstractSphericalOperator   ... tensor
        + q             ::Basic                       ... q-component
    """
    struct   SphericalOperatorComp  <: AbstractSphericalTensorComp
        tensor          ::AbstractSphericalTensor
        q               ::Basic
    end

    function Base.string(comp::SphericalOperatorComp)
        sa = string(comp.tensor) * "_" * string(comp.q)
        return( sa )
    end

    function Base.show(io::IO, comp::SphericalOperatorComp)
        sa = string(comp) * " spherical operator component.";       print(io, sa, "\n")
    end

    """
    `struct  SphericalTensor.SphericalFieldComp  <: AbstractSphericalTensorComp`   ... to represent a spherical field component.

        + tensor        ::AbstractSphericalField    ... exernal tensor field
        + q             ::Basic                     ... q-component
    """
    struct   SphericalFieldComp  <: AbstractSphericalTensorComp
        tensor          ::AbstractSphericalField
        q               ::Basic
    end

    function Base.string(comp::SphericalFieldComp)
        sa = string(comp.tensor) * "_" * string(comp.q)
        return( sa )
    end

    function Base.show(io::IO, comp::SphericalFieldComp)
        sa = string(comp) * " spherical field component.";       print(io, sa, "\n")
    end

    """
    `struct  SphericalTensor.SphericalState`   ... to represent a  |jm > state (spherical representation).

        + star          ::Bool      ... defines the complex-conjugate form, if true.
        + j             ::Basic     ... total angular momentum
        + m             ::Basic     ... projection of total angular momentum
    """
    struct   SphericalState
        star            ::Bool
        j               ::Basic
        m               ::Basic
    end

    """
    `SphericalTensor.SphericalState(j::Basic, m::Basic)`  ... constructor for just setting the ket-vector (star = false).
    """
    function SphericalState(j::Basic, m::Basic)
    	SphericalState(false, j, m)
    end

    function Base.string(state::SphericalState; plain=false)
        if  plain   sa =       string(state.j) * " " * string(state.m);        if state.star   sa = sa * " *"   end
        else        sa = "|" * string(state.j) * " " * string(state.m);        if state.star   sa = sa * " *"   end;      sa = sa * ">"
        end
        return( sa )
    end

    function Base.show(io::IO, state::SphericalState)
        sa = string(state) * " spherical state.";       print(io, sa, "\n")
    end

    """
    `struct  SphericalTensor.SphericalMatrixElement`
        ... to represent a (full)  <ja ma|| T^(k)_q |jb mb>  matrix element of the spherical operator T^(k).

        + star             ::Bool                         ... defines the complex-conjugate form, if true.
        + leftState        ::SphericalState               ... (reduced) bra-state.
        + constant         ::ScalarConstant               ... a scalar constant that is treated independently.
        + tensor           ::AbstractSphericalTensor      ... spherical tensor.
        + q                ::Basic                        ... q-quantum number of tensor
        + rightState       ::SphericalState               ... (reduced) ket-state.
    """
    struct   SphericalMatrixElement
        star               ::Bool
        leftState          ::SphericalState
        constant           ::ScalarConstant
        tensor             ::AbstractSphericalTensor
        q                  ::Basic
        rightState         ::SphericalState
    end

    """
    `SphericalTensor.SphericalMatrixElement(leftState::SphericalState, tensor::AbstractSphericalTensor, q::Basic, rightState::SphericalState)`
        ... constructor for just setting the matrix element without constant and with star = false.
    """
    function SphericalMatrixElement(leftState::SphericalState, tensor::AbstractSphericalTensor, q::Basic, rightState::SphericalState)
    	SphericalMatrixElement(false, leftState, ScalarConstant(false, Basic(1)), tensor, q, rightState)
    end

    function Base.string(me::SphericalMatrixElement)
        sa = "<" * string(me.leftState, plain=true) * " | "
        if  me.constant.symbol != Basic(1)    sa = sa * string(me.constant)    end
        sa = sa * string(me.tensor) * "_" * string(me.q) * " | " * string(me.rightState, plain=true) * ">"
        if  me.star                    sa = sa * "^(*)"                 end
        return( sa )
    end

    function Base.show(io::IO, me::SphericalMatrixElement)
        sa = string(me) * " reduced matrix element.";       print(io, sa, "\n")
    end

    """
    `struct  SphericalTensor.ProductMatrixElement`
        ... to represent a (reduced) <ja|| T^(k) |jb> matrix element of the spherical operator T^(k).

        + leftState        ::SphericalState                    ... (spherical) bra-state.
        + components       ::Array{SphericalOperatorComp,1}    ... (product of) spherical tensors.
        + rightState       ::SphericalState                    ... (spherical) ket-state.
    """
    struct   ProductMatrixElement
        leftState          ::SphericalState
        components         ::Array{SphericalOperatorComp,1}
        rightState         ::SphericalState
    end

    function Base.string(me::ProductMatrixElement)
        sa = "<" * string(me.leftState, plain=true) * " | "
        for  comp in me.components    sa = sa * "  " * string(comp)      end
        sa = sa * " | " * string(me.rightState, plain=true) * "> "
        return( sa )
    end

    function Base.show(io::IO, me::ProductMatrixElement)
        sa = string(me) * " product matrix element.";       print(io, sa, "\n")
    end

    """
    `struct  SphericalTensor.SphericalAmplitude`
        ... to represent a (non-reduced, many-electron) spherical amplitude that may include one or several spherical matrix elements, spherical
            field operator components as well as  various symbols from Racah's algebra. In general, these spherical amplitudes need further to be
            simplified to (fully) reduced amplitudes prior to the calculations.

        + rex              ::RacahAlgebra.RacahExpression       ... contains all factors from Racah's algebra including summations and delta-factors
        + constants        ::Array{AbstractSphericalConstant,1} ... product (list) of scalar constants.
        + fields           ::Array{SphericalFieldComp,1}        ... product (list) of external field tensors.
        + productMes       ::Array{ProductMatrixElement,1}      ... product (list) of product matrix elements.
    """
    struct   SphericalAmplitude
        rex                ::RacahAlgebra.RacahExpression
        constants          ::Array{AbstractSphericalConstant,1}
        fields             ::Array{SphericalFieldComp,1}
        productMes         ::Array{ProductMatrixElement,1}
    end

    function Base.string(amp::SphericalAmplitude)
        sa =      "\nRacahExpression:    " * string(amp.rex )
        sa = sa * "\nConstants:        "
        for cst in amp.constants       sa = sa * "  " * string(cst)      end
        sa = sa * "\nField components: "
        for field in amp.fields        sa = sa * "  " * string(field)    end
        sa = sa * "\nproductMes:       "
        for rme in amp.productMes      sa = sa * "  " * string(rme)      end
        return( sa )
    end

    function Base.show(io::IO, amp::SphericalAmplitude)
        sa = "Spherical amplitude with: " * string(amp);       print(io, sa, "\n")
    end

    """
    `struct  SphericalTensor.ReducedState`   ... to represent a (reduced) |j> state (spherical representation).

        + j             ::Basic     ... total angular momentum
    """
    struct   ReducedState
        j               ::Basic
    end

    function Base.string(state::ReducedState; plain=false)
        if  plain      sa = string(state.j)    else    sa = "|" * string(state.j) * ">"    end
        return( sa )
    end

    function Base.show(io::IO, state::ReducedState)
        sa = string(state) * " reduced state.";       print(io, sa, "\n")
    end

    """
    `struct  SphericalTensor.ReducedMatrixElement`
        ... to represent a (reduced) <ja|| T^(k) |jb> matrix element of the spherical operator T^(k).

        + leftState        ::ReducedState                      ... (reduced) bra-state.
        + tensors          ::Array{AbstractSphericalTensor,1}  ... (product of) spherical tensors.
        + rightState       ::ReducedState                      ... (reduced) ket-state.
    """
    struct   ReducedMatrixElement
        leftState          ::ReducedState
        tensors            ::Array{AbstractSphericalTensor,1}
        rightState         ::ReducedState
    end

    function Base.string(me::ReducedMatrixElement)
        sa = "<" * string(me.leftState, plain=true) * " || "
        for  tn in me.tensors   sa = sa * " " * string(tn)      end
        sa = sa * " || " * string(me.rightState, plain=true) * ">"
        return( sa )
    end

    function Base.show(io::IO, me::ReducedMatrixElement)
        sa = string(me) * " reduced matrix element.";       print(io, sa, "\n")
    end

    """
    `struct  SphericalTensor.ReducedAmplitude`
        ... to represent a (reduced, many-electron) amplitude that may include one or several reduced matrix elements, spherical field operators
            and various symbols from Racah's algebra. In general, these reduced amplitudes need to be simplified prior to the calculations.

        + rex              ::RacahAlgebra.RacahExpression       ... contains all factors from Racah's algebra including summations and delta-factors
        + fields           ::Array{AbstractSphericalField,1}    ... product (list) of external field tensors.
        + reducedMes       ::Array{ReducedMatrixElement,1}      ... product (list) of reduced matrix elements.
    """
    struct   ReducedAmplitude
        rex                ::RacahAlgebra.RacahExpression
        fields             ::Array{AbstractSphericalField,1}
        reducedMes         ::Array{ReducedMatrixElement,1}
    end

    function Base.string(amp::ReducedAmplitude)
        sa =      "\nRacahExpression:   " * string(amp.rex )
        sa = sa * "\nFields:           "
        for field in amp.fields    sa = sa * " " * string(field)    end
        sa = sa * "\nreducedMes:       "
        for rme in amp.reducedMes  sa = sa * " " * string(rme)      end
        return( sa )
    end

    function Base.show(io::IO, amp::ReducedAmplitude)
        sa = "Reduced amplitude with: " * string(amp);       print(io, sa, "\n")
    end

    """
    `SphericalTensor.conjugate(me::SphericalTensor.SphericalMatrixElement)`
        ... conjugates the spherical matrix element by conjugating all parts of the tensor operator and by interchanging the
            left- and right-hand states. A conjme::SphericalTensor.SphericalMatrixElement is returned.
    """
    function conjugate(me::SphericalTensor.SphericalMatrixElement)
        starx       = !me.star
        leftStatex  = me.rightState
        constantx   = ScalarConstant( !me.constant.star,   me.constant.symbol)
        tensorx     = SphericalTensor.conjugate( me.tensor )
        qx          = me.q
        rightStatex = me.leftState
        conjme      = SphericalTensor.SphericalMatrixElement(starx, leftStatex, constantx, tensorx, qx, rightStatex)

        return( conjme )
    end

    """
    `SphericalTensor.conjugate(tn::SphericalTensor.TensorProduct)`
        ... conjugates the (abstract) spherical matrix element by conjugating all parts of the tensor operator and by interchanging
            the individual tensor components of each tensor product.
    """
    function conjugate(tn::SphericalTensor.TensorProduct)

        starx       = tn.star
        kx          = tn.k
        tensor1x    = SphericalTensor.conjugate( tn.tensor2 )
        tensor2x    = SphericalTensor.conjugate( tn.tensor1 )
        conjtn      = SphericalTensor.TensorProduct(starx, kx, tensor1x, tensor2x)

        return( conjtn )
    end

    conjugate(tn::SphericalTensor.CkOperator)      = SphericalTensor.CkOperator( !tn.star, tn.k)
    conjugate(tn::SphericalTensor.CoulombOperator) = SphericalTensor.CoulombOperator( !tn.star)
    conjugate(tn::SphericalTensor.DipoleOperator)  = SphericalTensor.DipoleOperator( !tn.star)
    conjugate(tn::SphericalTensor.TkOperator)      = SphericalTensor.TkOperator( !tn.star, tn.k)

    conjugate(tn::SphericalTensor.ScalarConstant)  = SphericalTensor.ScalarConstant( !tn.star, tn.symbol)
    conjugate(tn::SphericalTensor.UVector)         = SphericalTensor.UVector( !tn.star)

    """
    `SphericalTensor.expandSphericalMatrixElements(mes::Array{SphericalTensor.SphericalMatrixElement,1})`
        ... expand the (product of) spherical matrix elements into a reduced spherical amplitude; it separates the scalar constants from the
            field operators and the reduced (electron) matrix elements. The Wigner-Eckardt theorem is applied to reduce the spherical
            matrix elements and a Clebsch-Gordan expansion for each product of tensor operators. All additional magnetic quantum numbers,
            that are needed for this expansion, are derived from the corresponding tensor ranks. The Clebsch-Gordan/Wigner symbols are
            'collected' in a Racah expression for further simplification. A sphericalAmp::SphericalTensor.SphericalAmplitude is returned.
    """
    function expandSphericalMatrixElements(mes::Array{SphericalTensor.SphericalMatrixElement,1})
        rex         = RacahAlgebra.RacahExpression()
        constants   = AbstractSphericalConstant[]
        fieldComps  = SphericalFieldComp[]
        productMes  = ProductMatrixElement[]

        for  mex in mes
           # Always work with the non-conjugated form of spherical ME
           if  mex.star  me = SphericalTensor.conjugate(mex)    else   me = mex     end
           # 'Add' the constants
           if me.constant.symbol != Basic(1)    push!(constants, me.constant)   end
           # Decompose the spherical operator recusively
           wa, wb, wc, wd = SphericalTensor.expandSphericalTensorComponent(me.tensor, me.q)
           rex = rex * wa;    append!(constants, wb);    append!(fieldComps, wc)
           ##x # Apply Wigner-Eckardt theorem and append reduced matrix element
           ##x k     = SphericalTensor.getRank(me.tensor)
           ##x rex   = rex * SphericalTensor.getWignerEckardtFactor(me.leftState, k, me.q, me.rightState)
           if  me.leftState.star  ||  me.rightState.star   error("SphericalState() assumes non-conjugated spherical states.")     end
           @show wd
           prodMe = SphericalTensor.ProductMatrixElement( me.leftState, wd, me.rightState)
           push!(productMes, prodMe )
        end
        sphericalAmp = SphericalTensor.SphericalAmplitude(rex, constants, fieldComps, productMes)

        return( sphericalAmp )
    end

    """
    `SphericalTensor.expandSphericalTensorComponent(tn::SphericalTensor.AbstractSphericalTensorComp, q::Basic)`
        ... expand a coupled product of spherical tensors into its (uncoupled) tensor components by means of a Clebsch-Gordan expansion.
            All necessary magnetic quantum numbers are 'derived' from the corresponding tensor ranks: K -> qK, 2 -> q2,  0 -> 0.
            No expansion is made for elementary spherical tensors. The tensor components of this expansion is returned by (the tuple of)
            three variables: (rex::RacahExpression, fields::Array{SphericalTensor.AbstractSphericalField,1},
                              operators::Array{SphericalTensor.AbstractSphericalOperator,1}).
    """
    function expandSphericalTensorComponent(tn::SphericalTensor.AbstractSphericalTensor, q::Basic)
        rex           = RacahAlgebra.RacahExpression()
        constants     = SphericalTensor.AbstractSphericalConstant[]
        fieldComps    = SphericalTensor.SphericalFieldComp[]
        operatorComps = SphericalTensor.SphericalOperatorComp[]

        if  typeof(tn) in [SphericalTensor.CkOperator, SphericalTensor.CoulombOperator, SphericalTensor.DipoleOperator,
                           SphericalTensor.TkOperator]
            push!(operatorComps, SphericalOperatorComp(tn,q) )
        elseif   typeof(tn) in [SphericalTensor.ScalarConstant]
            push!(constants, tn )
        elseif   typeof(tn) in [SphericalTensor.UVector]
            push!(fieldComps, SphericalFieldComp(tn,q) )
        elseif   typeof(tn) in [SphericalTensor.TensorProduct]
            B0  = Basic(0)
            k1  = getRank(tn.tensor1);   if  k1 == B0 q1 = B0   else   q1 = Basic("q" * string(k1) * "_" *string(floor(10000*rand()))[1:2]) end
            k2  = getRank(tn.tensor2);   if  k2 == B0 q2 = B0   else   q2 = Basic("q" * string(k2) * "_" *string(floor(10000*rand()))[1:2]) end
            kk  = getRank(tn);           if  kk == B0 qq = B0   else   qq = Basic("q" * string(kk) * "_" *string(floor(10000*rand()))[1:2]) end

            rx1, cx1, fx1, ox1 = SphericalTensor.expandSphericalTensorComponent(tn.tensor1, q1)
            rx2, cx2, fx2, ox2 = SphericalTensor.expandSphericalTensorComponent(tn.tensor2, q2)
            rex = rex * rx1 * rx2 * RacahAlgebra.ClebschGordanExpansion(k1, q1, k2, q2, kk, qq)
            append!(constants,  cx1);      append!(constants,  cx2)
            append!(fieldComps, fx1);      append!(fieldComps, fx2)
            append!(operatorComps, ox1);   append!(operatorComps, ox2)
        else     error("stop a")
        end

        return( rex, constants, fieldComps, operatorComps )
    end

    getRank(tn::SphericalTensor.CkOperator)      = tn.k
    getRank(tn::SphericalTensor.CoulombOperator) = Basic(0)
    getRank(tn::SphericalTensor.DipoleOperator)  = Basic(1)
    getRank(tn::SphericalTensor.TkOperator)      = tn.k

    getRank(tn::SphericalTensor.ScalarConstant)  = Basic(0)
    getRank(tn::SphericalTensor.UVector)         = Basic(1)

    getRank(tn::SphericalTensor.TensorProduct)   = tn.k

    """
    `SphericalTensor.getWignerEckardtFactor(leftState::SphericalState, k::Basic, q::Basic, rightState::SphericalState)`
        ... Returns the pre-factor from the Wigner-Eckardt theorem. A rex::RacahAlgebra.RacahExpression is returned.
    """
    function getWignerEckardtFactor(leftState::SphericalState, k::Basic, q::Basic, rightState::SphericalState)
        rex = RacahAlgebra.ClebschGordan(rightState.j, rightState.m, k, q, leftState.j, leftState.j)
        return( rex )
    end

end # module


text:
###
# Lmul
####

# Diagonal multiplication never changes structure
similar(M::Lmul{<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.B, T, axes)
# equivalent to rescaling
function materialize!(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout}})
    M.B .= getindex_value(M.A.diag) .* M.B
    M.B
end

copy(M::Lmul{<:DiagonalLayout,<:DiagonalLayout}) = Diagonal(M.A.diag .* M.B.diag)
copy(M::Lmul{<:DiagonalLayout}) = M.A.diag .* M.B
copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout}}) = getindex_value(M.A.diag) * M.B
copy(M::Lmul{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = Diagonal(getindex_value(M.A.diag) * M.B.diag)

# Diagonal multiplication never changes structure
similar(M::Rmul{<:Any,<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.A, T, axes)
# equivalent to rescaling
function materialize!(M::Rmul{<:Any,<:DiagonalLayout{<:AbstractFillLayout}})
    M.A .= M.A .* getindex_value(M.B.diag)
    M.A
end

copy(M::Rmul{<:Any,<:DiagonalLayout}) = M.A .* permutedims(M.B.diag)
copy(M::Rmul{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) =  M.A .* getindex_value(M.B.diag)

copy(M::Ldiv{<:DiagonalLayout,<:DiagonalLayout}) = Diagonal(inv.(M.A.diag) .* M.B.diag)
copy(M::Ldiv{<:DiagonalLayout}) = inv.(M.A.diag) .* M.B
copy(M::Ldiv{<:DiagonalLayout{<:AbstractFillLayout}}) = inv(getindex_value(M.A.diag)) .* M.B
copy(M::Ldiv{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout}) = Diagonal(inv(getindex_value(M.A.diag)) .* M.B.diag)

copy(M::Rdiv{<:DiagonalLayout,<:DiagonalLayout}) = Diagonal(M.A.diag .* inv.(M.B.diag))
copy(M::Rdiv{<:Any,<:DiagonalLayout}) = M.A .* inv.(permutedims(M.B.diag))
copy(M::Rdiv{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) = M.A .* inv(getindex_value(M.B.diag))
copy(M::Rdiv{<:DiagonalLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = Diagonal(M.A.diag .* inv(getindex_value(M.B.diag)))

ext: jl


text:
type MonteCarloProblem{T,T2,T3} <: AbstractMonteCarloProblem
  prob::T
  prob_func::T2
  output_func::T3
end

MonteCarloProblem(prob::DEProblem;
                  output_func = (sol,i)-> sol,
                  prob_func= (prob,i)->prob) =
                  MonteCarloProblem(prob,prob_func,output_func)


text:
"""
    TextHashFunction(hash_function::Function, cardinality::Int)

The basic structure for performing text hashing: uses the `hash_function`
to generate feature vectors of length `cardinality`.

# Details
The hash trick is the use a hash function instead of a lexicon to determine the
columns of a DocumentTermMatrix-like encoding of the data. To produce a DTM for
a Corpus for which we do not have an existing lexicon, we need someway to map
the terms from each document into column indices. We use the now standard "Hash Trick"
in which we hash strings and then reduce the resulting integers modulo N, which
defines the numbers of columns we want our DTM to have. This amounts to
doing a non-linear dimensionality reduction with low probability that similar
terms hash to the same dimension.

To make things easier, we wrap Julia's hash functions in a new type,
TextHashFunction, which maintains information about the desired cardinality
of the hashes.

# References:
  * [The "Hash Trick" wiki page](https://en.wikipedia.org/wiki/Feature_hashing)
  * [Moody, John 1989](http://papers.nips.cc/paper/175-fast-learning-in-multi-resolution-hierarchies.pdf)

# Examples
```
julia> doc = StringDocument("this is a text")
       thf = TextHashFunction(hash, 13)
       hash_dtv(doc, thf, Float16)
13-element Array{Float16,1}:
 1.0
 1.0
 0.0
 0.0
 0.0
 0.0
 0.0
 2.0
 0.0
 0.0
 0.0
 0.0
 0.0
```
"""
struct TextHashFunction
    hash_function::Function
    cardinality::Int
end

TextHashFunction(cardinality::Int) = TextHashFunction(DEFAULT_HASH_FUNCTION, cardinality)

TextHashFunction() = TextHashFunction(DEFAULT_HASH_FUNCTION, DEFAULT_CARDINALITY)

cardinality(h::TextHashFunction) = h.cardinality

hash_function(h::TextHashFunction) = h.hash_function

function index_hash(s::AbstractString, h::TextHashFunction)
    return Int(rem(h.hash_function(s), h.cardinality)) + 1
end


text:
# Feynman and Hibbs, page 319.
athermal_values = [
    [3.00, 5.00, 7.00, 9.00, 11.00],
    [3.44, 4.02, 5.81, 9.85, 15.5],
    [2.55, 2.13, 1.6, 1.28, 1.15]
]

# Hellwarth & Biaggio (1999) data. Table 3. [α, β, [(v, w)...]]
thermal_values = [
    [1.0, 2.0, 3.0, 4.0, 5.0],
    [1.00, 1.25, 1.50, 1.75, 2.00, 2.50, 2.75, 3.00],
    [
        [(7.20, 6.50), (5.94, 5.31), (5.11, 4.54), (4.54, 4.01), (4.12, 3.63), (3.57, 3.13), (3.38, 2.97), (3.24, 2.84)],
        [(7.69, 6.20), (6.38, 5.03), (5.52, 4.29), (4.92, 3.77), (4.48, 3.40), (3.89, 2.93), (3.69, 2.77), (3.53, 2.65)],
        [(8.26, 5.87), (6.91, 4.73), (6.02, 4.01), (5.40, 3.51), (4.94, 3.16), (4.32, 2.70), (4.10, 2.55), (3.93, 2.44)],
        [(8.95, 5.51), (7.57, 4.41), (6.66, 3.70), (6.02, 3.22), (5.54, 2.88), (4.90, 2.45), (4.67, 2.31), (4.48, 2.21)],
        [(9.78, 5.13), (8.39, 4.05), (7.48, 3.37), (6.84, 2.92), (6.36, 2.60), (5.71, 2.20), (5.48, 2.07), (5.29, 1.97)]
    ]
]

@testset "variation" begin

    @testset "athermal" begin
        for (α, v, w) in zip(athermal_values...)
            @test all(all.(map((i,j) -> isapprox.(i,j, atol = 0.1), variation(α), (v, w))))
        end
    end

    @testset "thermal" begin
        for (α, vw) in zip(thermal_values[1], thermal_values[3])
            for (β, p) in zip(thermal_values[2], vw)
                @test all(all.(map((i,j) -> isapprox.(i,j, atol = 0.01), variation(α, β), p)))
            end
        end
    end
end


text:
"""
# Onda Corporation
# Hydrophone Calibration Data
#
Calibration_DATE	29-Sep-2020
#
# Hydrophone description
HYD_MFG	Onda
HYD_MODEL	HGL
HYD_APERTURE_NOM_UM	0200
HYD_SN	2322
HYD_POLARITY	Positive
#
# External preamplifier description
AMP_MFG	NONE
AMP_MODEL	NONE
AMP_SN	NONE
AMP_GAIN_NOM_DB	NONE
AMP_BW_MHZ	NONE
AMP_OUTPUT_IMP_NOM_OHMS	NONE
#
# Cable/adaptor between hydrophone and preamplifier
CA_BWT_HYD_AMP	NONE
#
# Environment
ELECTRICAL_LOAD	OpenCircuit
WATER_TEMP_DEGC	23.5
WATER_RESISTIVITY_MOHMS-CM	7.0
#
# Calibration data
# Freq(MHz) Sens(dB re. 1V/uPa) Sens(V/Pa) Sens(V^2cm^2/W) Cap(pF)
#
DATA_FIELDS	5
DATA_FIELD	FREQ_MHz
DATA_FIELD	SENS_DB
DATA_FIELD	SENS_VPERPA
DATA_FIELD	SENS_V2CM2PERW
DATA_FIELD	CAP_PF
# End of file header
HEADER_END	0
"""

@calibration :Onda_HGL0200_2322 [
     1.000	 -267.72	  4.110E-008	2.534E-005	      9.963
     1.050	 -268.05	  3.958E-008	2.350E-005	      9.953
     1.100	 -268.12	  3.925E-008	2.311E-005	      9.943
     1.150	 -268.47	  3.773E-008	2.136E-005	      9.933
     1.200	 -268.60	  3.715E-008	2.070E-005	      9.923
     1.250	 -268.89	  3.595E-008	1.939E-005	      9.914
     1.300	 -268.95	  3.568E-008	1.909E-005	      9.904
     1.350	 -269.05	  3.529E-008	1.868E-005	      9.895
     1.400	 -269.27	  3.441E-008	1.776E-005	      9.886
     1.450	 -269.80	  3.235E-008	1.569E-005	      9.876
     1.500	 -269.87	  3.209E-008	1.545E-005	      9.867
     1.550	 -269.70	  3.272E-008	1.606E-005	      9.858
     1.600	 -269.94	  3.185E-008	1.522E-005	      9.849
     1.650	 -270.10	  3.125E-008	1.465E-005	      9.840
     1.700	 -270.09	  3.131E-008	1.470E-005	      9.832
     1.750	 -270.23	  3.079E-008	1.422E-005	      9.823
     1.800	 -270.46	  3.001E-008	1.351E-005	      9.814
     1.850	 -270.70	  2.919E-008	1.278E-005	      9.806
     1.900	 -270.81	  2.881E-008	1.245E-005	      9.797
     1.950	 -270.83	  2.874E-008	1.239E-005	      9.789
     2.000	 -270.82	  2.876E-008	1.241E-005	      9.781
     2.050	 -270.77	  2.894E-008	1.256E-005	      9.772
     2.100	 -270.64	  2.937E-008	1.294E-005	      9.764
     2.150	 -270.65	  2.934E-008	1.291E-005	      9.756
     2.200	 -270.66	  2.932E-008	1.289E-005	      9.748
     2.250	 -270.68	  2.924E-008	1.282E-005	      9.740
     2.300	 -270.47	  2.995E-008	1.345E-005	      9.732
     2.350	 -270.48	  2.993E-008	1.344E-005	      9.725
     2.400	 -270.46	  3.000E-008	1.350E-005	      9.717
     2.450	 -270.47	  2.994E-008	1.345E-005	      9.709
     2.500	 -270.58	  2.959E-008	1.314E-005	      9.702
     2.550	 -270.51	  2.983E-008	1.335E-005	      9.695
     2.600	 -270.04	  3.149E-008	1.487E-005	      9.687
     2.650	 -269.35	  3.409E-008	1.743E-005	      9.680
     2.700	 -269.00	  3.548E-008	1.889E-005	      9.673
     2.750	 -268.89	  3.592E-008	1.935E-005	      9.666
     2.800	 -268.62	  3.708E-008	2.063E-005	      9.658
     2.850	 -268.68	  3.682E-008	2.033E-005	      9.651
     2.900	 -268.60	  3.715E-008	2.070E-005	      9.645
     2.950	 -268.48	  3.766E-008	2.127E-005	      9.638
     3.000	 -268.83	  3.618E-008	1.963E-005	      9.631
     3.050	 -269.18	  3.477E-008	1.813E-005	      9.624
     3.100	 -269.05	  3.528E-008	1.867E-005	      9.618
     3.150	 -268.81	  3.629E-008	1.975E-005	      9.611
     3.200	 -268.68	  3.681E-008	2.032E-005	      9.604
     3.250	 -268.67	  3.685E-008	2.037E-005	      9.598
     3.300	 -268.58	  3.725E-008	2.081E-005	      9.592
     3.350	 -268.51	  3.752E-008	2.112E-005	      9.585
     3.400	 -268.75	  3.650E-008	1.999E-005	      9.579
     3.450	 -268.99	  3.551E-008	1.891E-005	      9.573
     3.500	 -268.86	  3.607E-008	1.951E-005	      9.567
     3.550	 -268.56	  3.733E-008	2.090E-005	      9.560
     3.600	 -268.33	  3.832E-008	2.203E-005	      9.554
     3.650	 -268.35	  3.823E-008	2.192E-005	      9.548
     3.700	 -268.61	  3.711E-008	2.065E-005	      9.543
     3.750	 -268.92	  3.580E-008	1.922E-005	      9.537
     3.800	 -269.12	  3.500E-008	1.838E-005	      9.531
     3.850	 -269.18	  3.477E-008	1.813E-005	      9.525
     3.900	 -269.15	  3.487E-008	1.823E-005	      9.519
     3.950	 -269.21	  3.465E-008	1.801E-005	      9.514
     4.000	 -269.21	  3.464E-008	1.800E-005	      9.508
     4.050	 -269.11	  3.505E-008	1.843E-005	      9.503
     4.100	 -269.03	  3.536E-008	1.875E-005	      9.497
     4.150	 -269.01	  3.542E-008	1.882E-005	      9.492
     4.200	 -269.12	  3.499E-008	1.837E-005	      9.487
     4.250	 -269.30	  3.427E-008	1.762E-005	      9.481
     4.300	 -269.42	  3.382E-008	1.716E-005	      9.476
     4.350	 -269.45	  3.369E-008	1.702E-005	      9.471
     4.400	 -269.38	  3.394E-008	1.728E-005	      9.466
     4.450	 -269.22	  3.461E-008	1.797E-005	      9.461
     4.500	 -269.09	  3.511E-008	1.849E-005	      9.455
     4.550	 -269.04	  3.533E-008	1.872E-005	      9.450
     4.600	 -268.90	  3.588E-008	1.931E-005	      9.446
     4.650	 -268.68	  3.682E-008	2.034E-005	      9.441
     4.700	 -268.53	  3.744E-008	2.103E-005	      9.436
     4.750	 -268.53	  3.744E-008	2.102E-005	      9.431
     4.800	 -268.50	  3.758E-008	2.118E-005	      9.426
     4.850	 -268.47	  3.771E-008	2.134E-005	      9.421
     4.900	 -268.53	  3.744E-008	2.102E-005	      9.417
     4.950	 -268.59	  3.720E-008	2.076E-005	      9.412
     5.000	 -268.63	  3.703E-008	2.057E-005	      9.408
     5.050	 -268.82	  3.623E-008	1.969E-005	      9.403
     5.100	 -269.13	  3.497E-008	1.834E-005	      9.399
     5.150	 -269.37	  3.400E-008	1.734E-005	      9.394
     5.200	 -269.49	  3.353E-008	1.686E-005	      9.390
     5.250	 -269.33	  3.415E-008	1.749E-005	      9.385
     5.300	 -269.06	  3.523E-008	1.862E-005	      9.381
     5.350	 -268.88	  3.597E-008	1.941E-005	      9.377
     5.400	 -268.81	  3.628E-008	1.974E-005	      9.372
     5.450	 -268.95	  3.569E-008	1.910E-005	      9.368
     5.500	 -269.09	  3.513E-008	1.851E-005	      9.364
     5.550	 -269.02	  3.541E-008	1.881E-005	      9.360
     5.600	 -268.87	  3.601E-008	1.946E-005	      9.356
     5.650	 -268.76	  3.648E-008	1.996E-005	      9.352
     5.700	 -268.80	  3.632E-008	1.979E-005	      9.348
     5.750	 -268.88	  3.597E-008	1.941E-005	      9.344
     5.800	 -268.86	  3.605E-008	1.949E-005	      9.340
     5.850	 -268.75	  3.651E-008	2.000E-005	      9.336
     5.900	 -268.68	  3.681E-008	2.032E-005	      9.332
     5.950	 -268.75	  3.654E-008	2.002E-005	      9.328
     6.000	 -268.84	  3.612E-008	1.957E-005	      9.324
     6.050	 -268.97	  3.561E-008	1.902E-005	      9.320
     6.100	 -268.99	  3.552E-008	1.893E-005	      9.317
     6.150	 -269.02	  3.540E-008	1.879E-005	      9.313
     6.200	 -269.00	  3.548E-008	1.889E-005	      9.309
     6.250	 -268.89	  3.592E-008	1.935E-005	      9.306
     6.300	 -268.80	  3.629E-008	1.975E-005	      9.302
     6.350	 -268.76	  3.646E-008	1.994E-005	      9.298
     6.400	 -268.80	  3.632E-008	1.978E-005	      9.295
     6.450	 -268.86	  3.606E-008	1.950E-005	      9.291
     6.500	 -269.00	  3.547E-008	1.887E-005	      9.288
     6.550	 -269.02	  3.540E-008	1.880E-005	      9.284
     6.600	 -269.04	  3.532E-008	1.871E-005	      9.281
     6.650	 -269.08	  3.517E-008	1.856E-005	      9.277
     6.700	 -269.06	  3.524E-008	1.863E-005	      9.274
     6.750	 -268.96	  3.563E-008	1.904E-005	      9.271
     6.800	 -268.93	  3.577E-008	1.919E-005	      9.267
     6.850	 -268.96	  3.566E-008	1.908E-005	      9.264
     6.900	 -268.99	  3.553E-008	1.894E-005	      9.261
     6.950	 -269.03	  3.536E-008	1.876E-005	      9.258
     7.000	 -269.06	  3.525E-008	1.864E-005	      9.254
     7.050	 -269.06	  3.524E-008	1.863E-005	      9.251
     7.100	 -269.18	  3.476E-008	1.812E-005	      9.248
     7.150	 -269.28	  3.435E-008	1.770E-005	      9.245
     7.200	 -269.44	  3.374E-008	1.707E-005	      9.242
     7.250	 -269.57	  3.322E-008	1.655E-005	      9.238
     7.300	 -269.58	  3.319E-008	1.652E-005	      9.235
     7.350	 -269.59	  3.314E-008	1.648E-005	      9.232
     7.400	 -269.61	  3.308E-008	1.642E-005	      9.229
     7.450	 -269.62	  3.305E-008	1.639E-005	      9.226
     7.500	 -269.67	  3.287E-008	1.620E-005	      9.223
     7.550	 -269.59	  3.314E-008	1.648E-005	      9.220
     7.600	 -269.55	  3.329E-008	1.662E-005	      9.217
     7.650	 -269.51	  3.346E-008	1.680E-005	      9.214
     7.700	 -269.42	  3.381E-008	1.714E-005	      9.211
     7.750	 -269.30	  3.427E-008	1.761E-005	      9.208
     7.800	 -269.17	  3.479E-008	1.815E-005	      9.206
     7.850	 -269.19	  3.472E-008	1.808E-005	      9.203
     7.900	 -269.24	  3.451E-008	1.786E-005	      9.200
     7.950	 -269.30	  3.428E-008	1.763E-005	      9.197
     8.000	 -269.27	  3.439E-008	1.774E-005	      9.194
     8.050	 -269.29	  3.430E-008	1.765E-005	      9.191
     8.100	 -269.25	  3.448E-008	1.784E-005	      9.189
     8.150	 -269.23	  3.454E-008	1.789E-005	      9.186
     8.200	 -269.14	  3.490E-008	1.827E-005	      9.183
     8.250	 -269.08	  3.516E-008	1.854E-005	      9.180
     8.300	 -268.97	  3.559E-008	1.900E-005	      9.178
     8.350	 -268.96	  3.565E-008	1.907E-005	      9.175
     8.400	 -268.92	  3.581E-008	1.923E-005	      9.172
     8.450	 -268.83	  3.618E-008	1.963E-005	      9.170
     8.500	 -268.76	  3.648E-008	1.996E-005	      9.167
     8.550	 -268.66	  3.691E-008	2.043E-005	      9.164
     8.600	 -268.51	  3.753E-008	2.113E-005	      9.162
     8.650	 -268.47	  3.770E-008	2.132E-005	      9.159
     8.700	 -268.41	  3.799E-008	2.165E-005	      9.157
     8.750	 -268.34	  3.827E-008	2.197E-005	      9.154
     8.800	 -268.26	  3.863E-008	2.238E-005	      9.151
     8.850	 -268.21	  3.884E-008	2.263E-005	      9.149
     8.900	 -268.22	  3.882E-008	2.260E-005	      9.146
     8.950	 -268.23	  3.876E-008	2.254E-005	      9.144
     9.000	 -268.27	  3.858E-008	2.233E-005	      9.141
     9.050	 -268.35	  3.826E-008	2.195E-005	      9.139
     9.100	 -268.38	  3.809E-008	2.176E-005	      9.136
     9.150	 -268.41	  3.796E-008	2.161E-005	      9.134
     9.200	 -268.51	  3.755E-008	2.115E-005	      9.131
     9.250	 -268.58	  3.725E-008	2.081E-005	      9.129
     9.300	 -268.68	  3.681E-008	2.032E-005	      9.126
     9.350	 -268.79	  3.636E-008	1.983E-005	      9.124
     9.400	 -268.83	  3.618E-008	1.963E-005	      9.121
     9.450	 -268.88	  3.597E-008	1.941E-005	      9.119
     9.500	 -268.93	  3.579E-008	1.921E-005	      9.117
     9.550	 -268.99	  3.551E-008	1.891E-005	      9.114
     9.600	 -269.02	  3.538E-008	1.878E-005	      9.112
     9.650	 -269.03	  3.537E-008	1.877E-005	      9.109
     9.700	 -268.97	  3.560E-008	1.901E-005	      9.107
     9.750	 -268.86	  3.607E-008	1.951E-005	      9.105
     9.800	 -268.79	  3.634E-008	1.980E-005	      9.102
     9.850	 -268.81	  3.627E-008	1.973E-005	      9.100
     9.900	 -268.78	  3.638E-008	1.985E-005	      9.098
     9.950	 -268.79	  3.633E-008	1.980E-005	      9.095
    10.000	 -268.80	  3.632E-008	1.979E-005	      9.093
    10.050	 -268.80	  3.630E-008	1.976E-005	      9.091
    10.100	 -268.83	  3.616E-008	1.962E-005	      9.088
    10.150	 -268.83	  3.618E-008	1.964E-005	      9.086
    10.200	 -268.81	  3.626E-008	1.972E-005	      9.084
    10.250	 -268.88	  3.597E-008	1.941E-005	      9.081
    10.300	 -268.99	  3.551E-008	1.891E-005	      9.079
    10.350	 -269.11	  3.505E-008	1.843E-005	      9.077
    10.400	 -269.25	  3.449E-008	1.785E-005	      9.075
    10.450	 -269.38	  3.396E-008	1.730E-005	      9.072
    10.500	 -269.59	  3.314E-008	1.647E-005	      9.070
    10.550	 -269.86	  3.215E-008	1.550E-005	      9.068
    10.600	 -269.99	  3.165E-008	1.502E-005	      9.065
    10.650	 -270.09	  3.130E-008	1.469E-005	      9.063
    10.700	 -270.16	  3.104E-008	1.446E-005	      9.061
    10.750	 -270.17	  3.100E-008	1.441E-005	      9.059
    10.800	 -270.21	  3.087E-008	1.429E-005	      9.056
    10.850	 -270.22	  3.084E-008	1.427E-005	      9.054
    10.900	 -270.20	  3.091E-008	1.433E-005	      9.052
    10.950	 -270.17	  3.102E-008	1.443E-005	      9.050
    11.000	 -270.10	  3.126E-008	1.465E-005	      9.048
    11.050	 -270.04	  3.146E-008	1.485E-005	      9.045
    11.100	 -270.03	  3.150E-008	1.489E-005	      9.043
    11.150	 -270.06	  3.142E-008	1.481E-005	      9.041
    11.200	 -270.04	  3.148E-008	1.486E-005	      9.039
    11.250	 -270.00	  3.162E-008	1.500E-005	      9.037
    11.300	 -270.03	  3.151E-008	1.489E-005	      9.034
    11.350	 -270.00	  3.162E-008	1.500E-005	      9.032
    11.400	 -269.86	  3.213E-008	1.548E-005	      9.030
    11.450	 -269.75	  3.253E-008	1.588E-005	      9.028
    11.500	 -269.65	  3.291E-008	1.625E-005	      9.026
    11.550	 -269.56	  3.326E-008	1.660E-005	      9.023
    11.600	 -269.48	  3.356E-008	1.690E-005	      9.021
    11.650	 -269.33	  3.414E-008	1.749E-005	      9.019
    11.700	 -269.24	  3.450E-008	1.786E-005	      9.017
    11.750	 -269.21	  3.463E-008	1.799E-005	      9.015
    11.800	 -269.17	  3.478E-008	1.815E-005	      9.012
    11.850	 -269.20	  3.467E-008	1.803E-005	      9.010
    11.900	 -269.12	  3.499E-008	1.837E-005	      9.008
    11.950	 -268.97	  3.561E-008	1.902E-005	      9.006
    12.000	 -268.98	  3.558E-008	1.899E-005	      9.004
    12.050	 -268.97	  3.561E-008	1.902E-005	      9.002
    12.100	 -268.95	  3.570E-008	1.911E-005	      9.000
    12.150	 -268.96	  3.565E-008	1.907E-005	      8.997
    12.200	 -268.97	  3.559E-008	1.900E-005	      8.995
    12.250	 -268.96	  3.566E-008	1.908E-005	      8.993
    12.300	 -268.90	  3.589E-008	1.932E-005	      8.991
    12.350	 -268.81	  3.628E-008	1.975E-005	      8.989
    12.400	 -268.75	  3.651E-008	1.999E-005	      8.987
    12.450	 -268.80	  3.632E-008	1.978E-005	      8.984
    12.500	 -268.79	  3.637E-008	1.984E-005	      8.982
    12.550	 -268.77	  3.644E-008	1.992E-005	      8.980
    12.600	 -268.83	  3.617E-008	1.962E-005	      8.978
    12.650	 -268.88	  3.597E-008	1.940E-005	      8.976
    12.700	 -268.92	  3.580E-008	1.923E-005	      8.974
    12.750	 -268.97	  3.560E-008	1.901E-005	      8.972
    12.800	 -269.02	  3.542E-008	1.882E-005	      8.970
    12.850	 -269.01	  3.545E-008	1.885E-005	      8.967
    12.900	 -269.02	  3.540E-008	1.880E-005	      8.965
    12.950	 -269.09	  3.513E-008	1.851E-005	      8.963
    13.000	 -269.09	  3.510E-008	1.848E-005	      8.961
    13.050	 -269.01	  3.542E-008	1.882E-005	      8.959
    13.100	 -269.07	  3.521E-008	1.859E-005	      8.957
    13.150	 -269.08	  3.514E-008	1.852E-005	      8.955
    13.200	 -268.97	  3.560E-008	1.901E-005	      8.953
    13.250	 -268.92	  3.579E-008	1.922E-005	      8.950
    13.300	 -268.85	  3.611E-008	1.956E-005	      8.948
    13.350	 -268.78	  3.641E-008	1.989E-005	      8.946
    13.400	 -268.76	  3.646E-008	1.994E-005	      8.944
    13.450	 -268.71	  3.670E-008	2.021E-005	      8.942
    13.500	 -268.71	  3.667E-008	2.018E-005	      8.940
    13.550	 -268.70	  3.673E-008	2.024E-005	      8.938
    13.600	 -268.64	  3.697E-008	2.050E-005	      8.936
    13.650	 -268.64	  3.698E-008	2.051E-005	      8.934
    13.700	 -268.69	  3.677E-008	2.028E-005	      8.932
    13.750	 -268.77	  3.642E-008	1.989E-005	      8.929
    13.800	 -268.84	  3.612E-008	1.957E-005	      8.927
    13.850	 -268.88	  3.597E-008	1.941E-005	      8.925
    13.900	 -268.96	  3.564E-008	1.905E-005	      8.923
    13.950	 -269.05	  3.529E-008	1.868E-005	      8.921
    14.000	 -269.07	  3.522E-008	1.860E-005	      8.919
    14.050	 -269.05	  3.529E-008	1.868E-005	      8.917
    14.100	 -268.98	  3.558E-008	1.899E-005	      8.915
    14.150	 -268.91	  3.586E-008	1.929E-005	      8.913
    14.200	 -268.87	  3.603E-008	1.947E-005	      8.911
    14.250	 -268.86	  3.606E-008	1.951E-005	      8.909
    14.300	 -268.90	  3.590E-008	1.934E-005	      8.907
    14.350	 -268.95	  3.568E-008	1.910E-005	      8.905
    14.400	 -268.95	  3.570E-008	1.912E-005	      8.903
    14.450	 -268.95	  3.570E-008	1.911E-005	      8.900
    14.500	 -269.00	  3.547E-008	1.887E-005	      8.898
    14.550	 -269.06	  3.525E-008	1.864E-005	      8.896
    14.600	 -269.10	  3.508E-008	1.846E-005	      8.894
    14.650	 -269.12	  3.500E-008	1.838E-005	      8.892
    14.700	 -269.10	  3.509E-008	1.847E-005	      8.890
    14.750	 -269.10	  3.506E-008	1.844E-005	      8.888
    14.800	 -269.05	  3.528E-008	1.867E-005	      8.886
    14.850	 -268.86	  3.606E-008	1.950E-005	      8.884
    14.900	 -268.91	  3.586E-008	1.929E-005	      8.882
    14.950	 -268.90	  3.589E-008	1.933E-005	      8.880
    15.000	 -268.84	  3.615E-008	1.960E-005	      8.878
    15.050	 -268.82	  3.624E-008	1.970E-005	      8.876
    15.100	 -268.86	  3.606E-008	1.950E-005	      8.874
    15.150	 -268.88	  3.596E-008	1.940E-005	      8.872
    15.200	 -268.89	  3.593E-008	1.937E-005	      8.870
    15.250	 -268.87	  3.602E-008	1.947E-005	      8.868
    15.300	 -268.89	  3.594E-008	1.938E-005	      8.866
    15.350	 -268.95	  3.570E-008	1.912E-005	      8.864
    15.400	 -268.83	  3.617E-008	1.962E-005	      8.862
    15.450	 -268.78	  3.640E-008	1.987E-005	      8.860
    15.500	 -268.75	  3.652E-008	2.001E-005	      8.859
    15.550	 -268.88	  3.597E-008	1.941E-005	      8.857
    15.600	 -268.88	  3.598E-008	1.942E-005	      8.855
    15.650	 -268.82	  3.624E-008	1.970E-005	      8.853
    15.700	 -268.82	  3.622E-008	1.968E-005	      8.851
    15.750	 -268.78	  3.641E-008	1.988E-005	      8.849
    15.800	 -268.73	  3.658E-008	2.007E-005	      8.847
    15.850	 -268.69	  3.678E-008	2.030E-005	      8.845
    15.900	 -268.66	  3.688E-008	2.040E-005	      8.843
    15.950	 -268.65	  3.694E-008	2.047E-005	      8.841
    16.000	 -268.64	  3.700E-008	2.053E-005	      8.840
    16.050	 -268.59	  3.721E-008	2.076E-005	      8.838
    16.100	 -268.55	  3.737E-008	2.095E-005	      8.836
    16.150	 -268.61	  3.711E-008	2.066E-005	      8.834
    16.200	 -268.67	  3.687E-008	2.040E-005	      8.832
    16.250	 -268.65	  3.694E-008	2.047E-005	      8.830
    16.300	 -268.64	  3.697E-008	2.050E-005	      8.829
    16.350	 -268.62	  3.709E-008	2.063E-005	      8.827
    16.400	 -268.55	  3.735E-008	2.093E-005	      8.825
    16.450	 -268.50	  3.757E-008	2.117E-005	      8.823
    16.500	 -268.54	  3.742E-008	2.101E-005	      8.821
    16.550	 -268.55	  3.735E-008	2.093E-005	      8.820
    16.600	 -268.51	  3.753E-008	2.112E-005	      8.818
    16.650	 -268.52	  3.751E-008	2.111E-005	      8.816
    16.700	 -268.55	  3.738E-008	2.096E-005	      8.814
    16.750	 -268.46	  3.775E-008	2.138E-005	      8.813
    16.800	 -268.43	  3.790E-008	2.154E-005	      8.811
    16.850	 -268.46	  3.775E-008	2.138E-005	      8.809
    16.900	 -268.55	  3.736E-008	2.093E-005	      8.808
    16.950	 -268.52	  3.749E-008	2.109E-005	      8.806
    17.000	 -268.52	  3.751E-008	2.110E-005	      8.804
    17.050	 -268.55	  3.737E-008	2.095E-005	      8.803
    17.100	 -268.62	  3.706E-008	2.061E-005	      8.801
    17.150	 -268.67	  3.688E-008	2.040E-005	      8.800
    17.200	 -268.57	  3.730E-008	2.087E-005	      8.798
    17.250	 -268.61	  3.711E-008	2.065E-005	      8.796
    17.300	 -268.57	  3.728E-008	2.085E-005	      8.795
    17.350	 -268.42	  3.791E-008	2.156E-005	      8.793
    17.400	 -268.41	  3.797E-008	2.162E-005	      8.792
    17.450	 -268.39	  3.804E-008	2.171E-005	      8.790
    17.500	 -268.46	  3.776E-008	2.138E-005	      8.789
    17.550	 -268.40	  3.800E-008	2.166E-005	      8.787
    17.600	 -268.37	  3.817E-008	2.185E-005	      8.786
    17.650	 -268.35	  3.823E-008	2.192E-005	      8.785
    17.700	 -268.31	  3.843E-008	2.215E-005	      8.783
    17.750	 -268.32	  3.836E-008	2.207E-005	      8.782
    17.800	 -268.33	  3.833E-008	2.204E-005	      8.780
    17.850	 -268.32	  3.836E-008	2.207E-005	      8.779
    17.900	 -268.38	  3.811E-008	2.179E-005	      8.778
    17.950	 -268.33	  3.831E-008	2.201E-005	      8.776
    18.000	 -268.30	  3.846E-008	2.218E-005	      8.775
    18.050	 -268.33	  3.831E-008	2.201E-005	      8.774
    18.100	 -268.40	  3.802E-008	2.169E-005	      8.773
    18.150	 -268.45	  3.782E-008	2.145E-005	      8.771
    18.200	 -268.40	  3.803E-008	2.170E-005	      8.770
    18.250	 -268.56	  3.733E-008	2.091E-005	      8.769
    18.300	 -268.44	  3.786E-008	2.150E-005	      8.768
    18.350	 -268.49	  3.761E-008	2.122E-005	      8.767
    18.400	 -268.55	  3.739E-008	2.097E-005	      8.765
    18.450	 -268.64	  3.700E-008	2.054E-005	      8.764
    18.500	 -268.57	  3.726E-008	2.083E-005	      8.763
    18.550	 -268.56	  3.732E-008	2.089E-005	      8.762
    18.600	 -268.59	  3.720E-008	2.076E-005	      8.761
    18.650	 -268.62	  3.709E-008	2.063E-005	      8.760
    18.700	 -268.75	  3.654E-008	2.002E-005	      8.759
    18.750	 -268.62	  3.705E-008	2.059E-005	      8.758
    18.800	 -268.55	  3.738E-008	2.096E-005	      8.757
    18.850	 -268.50	  3.760E-008	2.121E-005	      8.756
    18.900	 -268.63	  3.702E-008	2.056E-005	      8.756
    18.950	 -268.56	  3.731E-008	2.088E-005	      8.755
    19.000	 -268.70	  3.674E-008	2.025E-005	      8.754
    19.050	 -268.68	  3.682E-008	2.033E-005	      8.753
    19.100	 -268.63	  3.703E-008	2.057E-005	      8.752
    19.150	 -268.66	  3.688E-008	2.041E-005	      8.752
    19.200	 -268.77	  3.643E-008	1.990E-005	      8.751
    19.250	 -268.58	  3.724E-008	2.080E-005	      8.750
    19.300	 -268.63	  3.703E-008	2.057E-005	      8.750
    19.350	 -268.67	  3.684E-008	2.036E-005	      8.749
    19.400	 -268.65	  3.696E-008	2.049E-005	      8.749
    19.450	 -268.58	  3.725E-008	2.081E-005	      8.748
    19.500	 -268.55	  3.737E-008	2.095E-005	      8.747
    19.550	 -268.61	  3.711E-008	2.066E-005	      8.747
    19.600	 -268.56	  3.730E-008	2.087E-005	      8.747
    19.650	 -268.53	  3.746E-008	2.105E-005	      8.746
    19.700	 -268.59	  3.721E-008	2.077E-005	      8.746
    19.750	 -268.77	  3.645E-008	1.992E-005	      8.745
    19.800	 -268.78	  3.638E-008	1.985E-005	      8.745
    19.850	 -268.74	  3.658E-008	2.007E-005	      8.745
    19.900	 -268.80	  3.630E-008	1.977E-005	      8.745
    19.950	 -268.89	  3.591E-008	1.935E-005	      8.745
    20.000	 -268.99	  3.552E-008	1.893E-005	      8.745
]


text:
using BinDeps

@BinDeps.setup

libgumbo = library_dependency("libgumbo", aliases = ["libgumbo-1"])

provides(Sources,
         URI("http://jamesporter.me/static/julia/gumbo-1.0.tar.gz"),
         libgumbo,
         unpacked_dir="gumbo-1.0")

provides(Binaries, URI("http://sourceforge.net/projects/juliadeps-win/files/gumbo-$(Sys.MACHINE).7z"),
         libgumbo, os = :Windows)

provides(BuildProcess,
         Autotools(libtarget="libgumbo.la"),
         libgumbo)

@BinDeps.install [:libgumbo => :libgumbo]


text:
function lpc2lsp(state::SpectralParamState{LinearPredictionCoef};
                 loggain=false, kargs...)
    @assert !has_loggain(state)
    lpcdef = paramdef(state)
    data = rawdata(state)
    data = SPTK.lpc2lsp(data; loggain=loggain, kargs...)
    newdef = LineSpectralPair(param_order(lpcdef))
    SpectralParamState(newdef, data,
                       has_loggain(state), gain_normalized(state))
end


text:
struct Trans2DateTimeVariate <: AbstractDateTimeVariate
    name::String
    basevariate1::AbstractDateTimeVariate
    basevariate2::AbstractDateTimeVariate
    transform::Function
end

Base.length(var::Trans2DateTimeVariate) = length(var.basevariate1)

function slice(dtvariate::Trans2DateTimeVariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    base1 = dtvariate.basevariate1
    base2 = dtvariate.basevariate2
    f = dtvariate.transform
    g = (ms1::Int64, ms2::Int64) -> begin
        dt1 = Dates.epochms2datetime(ms1)
        dt2 = Dates.epochms2datetime(ms2)
        Dates.datetime2epochms(f(dt1, dt2))
    end
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = zip(slice(base1, fromobs, toobs, slicelength), slice(base2, fromobs, toobs, slicelength))
    mapslice2(g, slices, slicelength, Int64)
end

function Base.map(dtvariate::Trans2DateTimeVariate, dataframe::AbstractDataFrame)
    Trans2DateTimeVariate(dtvariate.name, map(dtvariate.basevariate1, dataframe), map(dtvariate.basevariate1, dataframe),
                          dtvariate.transform)
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
#=
   Native Fish model.

   From
   A. Nicholson, O. Woodberry and C. Twardy
   "The 'Native Fish' Bayesian networks"
   Technical Report 2010/3, Bayesian Intelligence.
   Available from https://bayesian-intelligence.com/publications.php
   """
   A local river with tree-lined banks is known to contain native fish populations, which need
   to be conserved. Parts of the river pass through croplands, and parts are susceptible to drought
   conditions. Pesticides are known to be used on the crops. Rainfall helps native fish populations
   by maintaining water flow, which increases habitat suitability as well as connectivity between
   different habitat areas. However rain can also wash pesticides that are dangerous to fish from
   the croplands into the river. There is concern that the trees and native fish will be affected by
   drought conditions and crop pesticides.

   In short, we want to model the effect of pesticide use and rainfall on native fish abundance
   and tree condition.
   ...
   In this model, native fish abundance has two main stressors: water-related and pesticide-related.
   The model also has three variables describing the water-related stressor:
     - water flow and connectivity: More water keeps the river from fragmenting into ponds, and leads to faster
       flow, which washes out pollutants. Higher water levels are better for the fish.

     - rainfall: This is intended to be year-to-date rainfall, a relatively short-term indicator.

     - drought conditions: A long-term indicator intended to summarize historical conditions. A multi-year
       drought will leave the soil quite dry, so that rain which falls soaks into the ground before reach-
       ing the rivers. (For this reason, much of the rain in the Australian reservoir catchment areas has failed
       to reach the reservoirs.)

   Two variables describe the pesticide-related stressor:
    - Pesticide use: How much pesticide is being used in the river catchment.

    - Pesticide concentration in river: The amount of pesticide in the river itself – which for this example we
      imagine cannot easily be directly observed.
   """

    Distributions of variable drought_conditions
    no         =>    7232  (0.723200)
    yes        =>    2768  (0.276800)

    Distributions of variable annual_rainfall
    below_average =>   10000  (1.000000)

    Distributions of variable pesticide_use
    high       =>    8940  (0.894000)
    low        =>    1060  (0.106000)

    Distributions of variable river_flow
    poor       =>    7075  (0.707500)
    good       =>    2925  (0.292500)

    Distributions of variable pesticide_in_river
    low        =>    6532  (0.653200)
    high       =>    3468  (0.346800)

    Distributions of variable native_fish_abundance
    low        =>    4675  (0.467500)
    medium     =>    2745  (0.274500)
    high       =>    2580  (0.258000)

    Distributions of variable tree_condition
    good       =>    5367  (0.536700)
    damaged    =>    3584  (0.358400)
    dead       =>    1049  (0.104900)

  cf ~/webppl/native_fish.wppl

=#

using Turing, StatsPlots
include("jl_utils.jl")

@model function native_fish()
    yes = 1
    no = 2

    # Type Level;
    high = 1
    medium = 2
    low = 3
    levels = [high, medium, low]
    highLow = [high, low]

    good = 1
    damaged = 2
    dead = 3

    conditions1 = [good, damaged, dead]

    poor = 2
    conditions2 = [good, poor]

    below_average = 1
    average = 2
    above_average = 3
    rainfall = [below_average, average, above_average]

    drought_conditions ~ Categorical([0.25,0.75]) # [yes,no]
    annual_rainfall ~ Categorical([0.1,0.7,0.2])  # rainfall
    pesticide_use ~ Categorical([0.9,0.1]) # highLow

    river_flow ~
        (drought_conditions==yes && annual_rainfall==below_average) ? Categorical([0.05,0.95]) : # conditions2
        (drought_conditions==yes && annual_rainfall==average)       ? Categorical([0.15,0.85]) :
        (drought_conditions==yes && annual_rainfall==above_average) ? Categorical([0.80,0.20]) :

        (drought_conditions==no && annual_rainfall==below_average)  ? Categorical([0.40,0.60]) :
        (drought_conditions==no && annual_rainfall==average)        ? Categorical([0.60,0.40]) :
        (drought_conditions==no && annual_rainfall==above_average)  ? Categorical([0.99,0.01]) : Categorical([0.5,0.5])

    pesticide_in_river ~
        (pesticide_use==high && annual_rainfall==below_average) ? Categorical([0.3,0.7]) : # [high,low]
        (pesticide_use==high && annual_rainfall==average)       ? Categorical([0.6,0.4]) :
        (pesticide_use==high && annual_rainfall==above_average) ? Categorical([0.8,0.2]) :

        (pesticide_use==low && annual_rainfall==below_average)  ? Categorical([0.1,0.9]) :
        (pesticide_use==low && annual_rainfall==average)        ? Categorical([0.2,0.8]) :
        (pesticide_use==low && annual_rainfall==above_average)  ? Categorical([0.3,0.7]) : Categorical([0.5,0.5])

    native_fish_abundance ~
        (pesticide_in_river==high && river_flow==good) ? Categorical([0.2,0.4,0.4]) : # levels
        (pesticide_in_river==high && river_flow==poor) ? Categorical([0.01,0.1,0.89]) :
        (pesticide_in_river==low && river_flow==good)  ? Categorical([0.8,0.15,0.05]) :
        (pesticide_in_river==low && river_flow==poor)  ? Categorical([0.05,0.15,0.8]) : Categorical(ones(3) ./ 3)

    tree_condition ~
        (drought_conditions==yes && annual_rainfall==below_average) ? Categorical([0.20,0.60,0.20]) : # conditions1
        (drought_conditions==yes && annual_rainfall==average)       ? Categorical([0.25,0.60,0.15]) :
        (drought_conditions==yes && annual_rainfall==above_average) ? Categorical([0.30,0.60,0.10]) :

        (drought_conditions==no && annual_rainfall==below_average)  ? Categorical([0.70,0.25,0.05]) :
        (drought_conditions==no && annual_rainfall==average)        ? Categorical([0.80,0.18,0.02]) :
        (drought_conditions==no && annual_rainfall==above_average)  ? Categorical([0.90,0.09,0.01]) : Categorical(ones(3) ./ 3)

    # true ~ Dirac(tree_condition == dead)
    # true ~ Dirac(pesticide_use == high)

    # true ~ Dirac(annual_rainfall == above_average);
    true ~ Dirac(annual_rainfall == below_average);

end

model = native_fish()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, NUTS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :drought_conditions, ["yes","no"])
show_var_dist_pct(chns, :annual_rainfall, ["below_average","average","above_average"])
show_var_dist_pct(chns, :pesticide_use, ["high","low"])
show_var_dist_pct(chns, :river_flow, ["good","poor"])
show_var_dist_pct(chns, :pesticide_in_river, ["high","low"])
show_var_dist_pct(chns, :native_fish_abundance, ["high","medium","low"])
show_var_dist_pct(chns, :tree_condition, ["good","damaged","dead"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module TensorLabXD

# Exports
# Types:
export Sector, AbstractIrrep, Irrep
export FusionStyle, UniqueFusion, MultipleFusion, MultiplicityFreeFusion,
        SimpleFusion, GenericFusion
export BraidingStyle, SymmetricBraiding, Bosonic, Fermionic, Anyonic
export Z2Irrep, Z3Irrep, Z4Irrep, ZNIrrep, U1Irrep, SU2Irrep, CU1Irrep
export Fermion, FermionParity, FermionNumber, FermionSpin
export FibonacciAnyon
export IsingAnyon

export VectorSpace, Field, ElementarySpace, InnerProductSpace, EuclideanSpace # abstract vector spaces
export ComplexSpace, CartesianSpace, GeneralSpace, GradedSpace # concrete spaces
export ZNSpace, Z2Space, Z3Space, Z4Space, U1Space, CU1Space, SU2Space
export Vect, Rep # space constructors
export CompositeSpace, ProductSpace # composite spaces
export FusionTree
export IndexSpace, TensorSpace, AbstractTensorMap, AbstractTensor, TensorMap, Tensor # tensors and tensor properties
export TensorMapSpace, EuclideanTensorSpace, EuclideanTensorMapSpace
export AbstractEuclideanTensorMap, EuclideanTensorMap
export TruncationScheme
export SpaceMismatch, SectorMismatch, IndexError # error types

# general vector space methods
export space, field, dual, dim, dims, fuse, flip, isdual, insertunit

# partial order for vector spaces
export infimum, supremum, isisomorphic, ismonomorphic, isepimorphic

# methods for sectors and properties thereof
export sectortype, sectors, hassector, Nsymbol, Fsymbol, Rsymbol, Bsymbol,
        frobeniusschur, twist
export fusiontrees, braid, permute, transpose
export Trivial, ZNSpace, SU2Irrep, U1Irrep, CU1Irrep # Fermion
# other fusion tree manipulations, should not be exported:
# export insertat, split, merge, repartition, artin_braid,
#        bendleft, bendright, foldleft, foldright, cycleclockwise, cycleanticlockwise

# some unicode
export ⊕, ⊗, ×, ⊠, ℂ, ℝ, ℤ, ←, →, ≾, ≿, ≅, ≺, ≻
export ℤ₂, ℤ₃, ℤ₄, U₁, SU, SU₂, CU₁
export fℤ₂, fU₁, fSU₂
export ℤ₂Space, ℤ₃Space, ℤ₄Space, U₁Space, CU₁Space, SU₂Space

# tensor maps
export domain, codomain, numind, numout, numin, spacetype, storagetype, eltype
export domainind, codomainind, allind
export tensormaptype
export blocksectors, blockdim, block, blocks

# random methods for constructor
export randuniform, randnormal, randisometry, randhaar

# special purpose constructors
export zero, one, one!, id, isomorphism, unitary, isometry

# tensor algebra and factorizations
export dot, norm, normalize, normalize!, tr
export mul!, lmul!, rmul!, adjoint!, pinv, axpy!, axpby!
export leftorth, rightorth, leftnull, rightnull,
        leftorth!, rightorth!, leftnull!, rightnull!,
        tsvd!, tsvd, eigen, eigen!, eig, eig!, eigh, eigh!, exp, exp!,
        isposdef, isposdef!, ishermitian, sylvester
export braid!, permute!, transpose!, twist!
export catdomain, catcodomain

export OrthogonalFactorizationAlgorithm, QR, QRpos, QL, QLpos, LQ, LQpos, RQ, RQpos,
        SVD, SDD, Polar

# tensor operations
export @tensor, @tensoropt, @ncon, ncon, @planar, @planar2
export scalar, add!, contract!

# truncation schemes
export notrunc, truncerr, truncdim, truncspace, truncbelow

# Imports
using TupleLabXD
using TupleLabXD: StaticLength

using StridedTensorXD

using TensorContractionsXD: TensorContractionsXD, @tensor, @tensoropt, @ncon, ncon
const TO = TensorContractionsXD

using LRUCache

using HalfIntegers
using WignerSymbols

using Base: @boundscheck, @propagate_inbounds, OneTo, tail, front,
            tuple_type_head, tuple_type_tail, tuple_type_cons,
            SizeUnknown, HasLength, HasShape, IsInfinite, EltypeUnknown, HasEltype
using Base.Iterators: product, filter
import Base: min, max

import LinearAlgebra
using LinearAlgebra: norm, dot, normalize, normalize!, tr,
                        axpy!, axpby!, lmul!, rmul!, mul!,
                        adjoint, adjoint!, transpose, transpose!,
                        pinv, sylvester,
                        eigen, eigen!, svd, svd!,
                        isposdef, isposdef!, ishermitian,
                        Diagonal, Hermitian
import Base.Meta

const IndexTuple{N} = NTuple{N, Int}

# Auxiliary files
include("auxiliary/auxiliary.jl")
include("auxiliary/dicts.jl")
include("auxiliary/linalg.jl")
include("auxiliary/random.jl")

# experiment with different dictionaries
const SectorDict{K, V} = SortedVectorDict{K, V}
const FusionTreeDict{K, V} = Dict{K, V}
fusiontreedict(I) = FusionStyle(I) isa UniqueFusion ? SingletonDict : FusionTreeDict

# Exception types:
abstract type TensorException <: Exception end

# Exception type for all errors related to sector mismatch
struct SectorMismatch{S<:Union{Nothing, String}} <: TensorException
    message::S
end
SectorMismatch()=SectorMismatch{Nothing}(nothing)
Base.show(io::IO, ::SectorMismatch{Nothing}) = print(io, "SectorMismatch()")
Base.show(io::IO, e::SectorMismatch) = print(io, "SectorMismatch(", e.message, ")")

# Exception type for all errors related to vector space mismatch
struct SpaceMismatch{S<:Union{Nothing, String}} <: TensorException
    message::S
end
SpaceMismatch()=SpaceMismatch{Nothing}(nothing)
Base.show(io::IO, ::SpaceMismatch{Nothing}) = print(io, "SpaceMismatch()")
Base.show(io::IO, e::SpaceMismatch) = print(io, "SpaceMismatch(", e.message, ")")

# Exception type for all errors related to invalid tensor index specification.
struct IndexError{S<:Union{Nothing, String}} <: TensorException
    message::S
end
IndexError()=IndexError{Nothing}(nothing)
Base.show(io::IO, ::IndexError{Nothing}) = print(io, "IndexError()")
Base.show(io::IO, e::IndexError) = print(io, "IndexError(", e.message, ")")

# typerepr
type_repr(T::Type) = repr(T)

# Definitions and methods for superselection sectors (quantum numbers)
include("sectors/sectors.jl")

# Constructing and manipulating fusion trees and iterators thereof
include("fusiontrees/fusiontrees.jl")

# Definitions and methods for vector spaces
include("spaces/vectorspaces.jl")

# Definitions and methods for tensors
include("tensors/abstracttensor.jl")
include("tensors/tensortreeiterator.jl")
include("tensors/tensor.jl")
include("tensors/adjoint.jl")
include("tensors/linalg.jl")
include("tensors/tensoroperations.jl")
include("tensors/indexmanipulations.jl")
include("tensors/truncation.jl")
include("tensors/factorizations.jl")
include("tensors/planar.jl")
include("tensors/braidingtensor.jl")

# deprecations: to be removed in version 1.0 or sooner
include("auxiliary/deprecate.jl")

end


text:
ENV["JULIA_PYTHONCALL_EXE"] = "/home/zhangyong/miniconda3/bin/python"
using Pkg, Faiss, SimilaritySearch, LinearAlgebra, JLD2

function create_benchmark(dim, n, m)
    X = rand(Float32, dim, n)
    Q = rand(Float32, dim, m)
    for c in eachcol(X) normalize!(c) end
    for c in eachcol(Q) normalize!(c) end
    X, Q
end

function create_gold(X, Q, k, dist)
    db = MatrixDatabase(X)
    queries = MatrixDatabase(Q)
    ex = ExhaustiveSearch(; db, dist)
    searchbatch(ex, queries, k; parallel=true)
end

function simsearch_test(X, Q, k, dist)
     db = MatrixDatabase(X)
     k = 100
     queries = MatrixDatabase(Q)
     G = SearchGraph(; db, dist)
     index!(G; parallel_block=512)
     searchbatch(G, queries, k; parallel=true)
end

function faiss_index(X, Q, k)
	dim, n = size(X)
	str = "HNSW64,Flat"  # HNSW32
	# str = "Flat"
	metric = "L2"  # L2  IP
	gpus = ""
	idx = Faiss.Index(dim; str, metric, gpus)  # init Faiss Index
	show(idx)   # show idx info
	add(idx, permutedims(X))
	D, I = Faiss.search(idx, permutedims(Q), k)
	I, D = permutedims(I), permutedims(D)
	I, D
end

function faiss_test(X, Q, k)
     vs_gallery = permutedims(X)
     vs_query = permutedims(Q)
     D, I = local_rank(vs_query, vs_gallery, k=k, metric="IP", gpus="")
     I, D = permutedims(I), permutedims(D)
	 I, D
end

function main()
	dim = 128
	k = 100
	dbfile = "benchmark-$dim.jld2"
	goldfile = "gold-$dim.jld2"
	faissfile = "faiss-$dim.jld2"
	faissindexfile = "faissindex-$dim.jld2"
	simfile = "sim-$dim.jld2"
	dist = NormalizedCosineDistance()

	X, Q = if isfile(dbfile)
		load(dbfile, "X", "Q")
	else
		X, Q = create_benchmark(dim, 10^6, 10^3)
		jldsave(dbfile, X=X, Q=Q)
		X, Q
	end

	@info "== gold"
	gI, gD, exhaustivetime = if isfile(goldfile)
		load(goldfile, "gI", "gD", "exhaustivetime")
	else
		exhaustivetime = @elapsed gI, gD = create_gold(X, Q, k, dist)
		jldsave(goldfile; gI, gD, exhaustivetime)
		gI, gD, exhaustivetime
	end
	@show exhaustivetime
	@info "== faiss"
	fI, fD, faisstime = if isfile(faissfile)
		load(faissfile, "fI", "fD", "faisstime")
	else
		faisstime = @elapsed fI, fD = faiss_test(X, Q, k)
		jldsave(faissfile; fI, fD, faisstime)
		fI, fD, faisstime
	end
	@show faisstime

	@info "== faiss index"
	fiI, fiD, faissindextime = if isfile(faissindexfile)
		load(faissindexfile, "fiI", "fiD", "faissindextime")
	else
		faissindextime = @elapsed fiI, fiD = faiss_index(X, Q, k)
		jldsave(faissindexfile; fiI, fiD, faissindextime)
		fiI, fiD, faissindextime
	end
	@show faissindextime

	@info "== SearchGraph"
	sI, sD, simtime = if isfile(simfile)
		load(simfile, "sI", "sD", "simtime")
	else
		simtime = @elapsed sI, sD = simsearch_test(X, Q, k, dist)
		jldsave(simfile; sI, sD, simtime)
		sI, sD, simtime
	end
	@show simtime

	@show faisstime => macrorecall(gI, fI)
	@show faissindextime => macrorecall(gI, fiI)
	@show simtime => macrorecall(gI, sI)

	println(Pkg.status())
end

main()

#=
export JULIA_NUM_THREADS=40
julia --project=/home/zhangyong/codes/Faiss.jl/Project.toml "/home/zhangyong/codes/Faiss.jl/test/test_2.jl"

32
faisstime => macrorecall(gI, fI) = 28.542211967 => 1.0
faissindextime => macrorecall(gI, fiI) = 72.467151637 => 0.05765999999999996
simtime => macrorecall(gI, sI) = 1794.894965328 => 0.4336100000000001

64
faisstime => macrorecall(gI, fI) = 28.064558278 => 0.99997
faissindextime => macrorecall(gI, fiI) = 131.294182817 => 0.08665999999999986
simtime => macrorecall(gI, sI) = 88.273946672 => 0.4000899999999998

128
faisstime => macrorecall(gI, fI) = 28.228652428 => 1.0
faissindextime => macrorecall(gI, fiI) = 244.810945869 => 0.11100999999999987
simtime => macrorecall(gI, sI) = 88.660722975 => 0.40738999999999975

64   L2
faisstime => macrorecall(gI, fI) = 28.117816367 => 0.99999
faissindextime => macrorecall(gI, fiI) = 132.110294344 => 0.08639999999999981

=#


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

"""
    SamplingMethod

A method for sampling from geometric objects.
"""
abstract type SamplingMethod end

"""
    sample([rng], object, method)

Sample elements or points from geometric `object`
with `method`. Optionally, specify random number
generator `rng`.
"""
sample(object, method) =
  sample(Random.GLOBAL_RNG, object, method)

"""
    DiscreteSamplingMethod

A method for sampling from discrete representations
of geometric objects such as meshes or collections
of geometries.
"""
abstract type DiscreteSamplingMethod end

"""
    ContinuousSamplingMethod

A method for sampling from continuous representations
of geometric objects. In this case, geometric objects
are interpreted as a set of points in the embedding
space.
"""
abstract type ContinuousSamplingMethod end

function sample(rng::AbstractRNG, geometry::Union{Multi,Polygon},
                method::ContinuousSamplingMethod)
  mesh = discretize(geometry, FIST(rng))
  sample(rng, mesh, method)
end

function sample(rng::AbstractRNG, geometry::Ngon,
                method::ContinuousSamplingMethod)
  mesh = discretize(geometry, Dehn1899())
  sample(rng, mesh, method)
end

# ----------------
# IMPLEMENTATIONS
# ----------------

# discrete sampling
include("sampling/uniform.jl")
include("sampling/weighted.jl")
include("sampling/ball.jl")

# continuous sampling
include("sampling/regular.jl")
include("sampling/homogeneous.jl")
include("sampling/mindistance.jl")

# ----------
# UTILITIES
# ----------

"""
    sample([rng], object, nsamples, [weights], replace=false)

Generate `nsamples` samples from spatial `object`
uniformly or using `weights`, with or without
replacement depending on `replace` option.
"""
sample(object::DomainOrData, nsamples::Int,
       weights::AbstractVector=[]; replace=false) =
  sample(Random.GLOBAL_RNG, object, nsamples, weights, replace)

function sample(rng::AbstractRNG,
                object::DomainOrData,
                nsamples::Int,
                weights::AbstractVector,
                replace::Bool)
  method = if isempty(weights)
    UniformSampling(nsamples, replace)
  else
    WeightedSampling(nsamples, weights, replace)
  end
  sample(rng, object, method)
end


text:
using IterativeSolvers, KrylovKit, LinearAlgebra

# In this file, we provide linear solvers for the Package

abstract type AbstractLinearSolver end

# The function linsolve(J, x; kwargs...) must return whether the solve was successfull and how many steps were required for the solve.

# the following function can be used to cache some factorization, see DefaultLS() case for example
function (ls::AbstractLinearSolver)(J, rhs1, rhs2; kwargs...)
	sol1, flag1, it1 = ls(J, rhs1; kwargs...)
	sol2, flag2, it2 = ls(J, rhs2; kwargs...)
	return sol1, sol2, flag1 & flag2, (it1, it2)
end

####################################################################################################
# The following functions are used for the Continuation of Hopf points and the computation of Floquet multipliers

# this function returns a₀ * I + a₁ .* J and ensures that we don't do unnecessary computations like 0*I + 1*J
function _axpy(J, a₀, a₁)
	if a₀ == 0
		if a₁ == 1
			return J
		else
			return a₁ .* J
		end
	elseif a₀ == 1
		if a₁ == 1
			return I + J
		else
			return I + a₁ .* J
		end
	else
		return a₀ * I + a₁ .* J
	end
end

function _axpy_op(J, v::AbstractArray, a₀, a₁)
	if a₀ == 0
		if a₁ == 1
			return apply(J, v)
		else
			return a₁ .* apply(J, v)
		end
	elseif a₀ == 1
		if a₁ == 1
			return v .+ apply(J, v)
		else
			return v .+ a₁ .* apply(J, v)
		end
	else
		return a₀ .* v .+ a₁ .* apply(J, v)
	end
end

####################################################################################################
# Solvers for default \ operator (backslash)
####################################################################################################
"""
$(TYPEDEF)

This struct is used to provide the backslash operator. Can be used to solve `(a₀ * I + a₁ * J) * x = rhs`.

$(TYPEDFIELDS)
"""
@with_kw struct DefaultLS <: AbstractLinearSolver
	"Whether to catch a factorization for multiple solves. Some operators may not support LU (like ApproxFun.jl) or QR factorization so it is best to let the user decides. Some matrices do not have `factorize` like `StaticArrays.MMatrix`."
	useFactorization::Bool = true
end

# this function is used to solve (a₀ * I + a₁ * J) * x = rhs
# the options a₀, a₁ are only used for the Hopf Newton / Continuation
function (l::DefaultLS)(J, rhs; a₀ = 0, a₁ = 1, kwargs...)
	return _axpy(J, a₀, a₁) \ rhs, true, 1
end

# this function is used to solve (a₀ * I + a₁ * J) * x = rhs
# with multiple RHS. We can cache the factorization in this case
# the options a₀, a₁ are only used for the Hopf Newton / Continuation
function (l::DefaultLS)(J, rhs1, rhs2; a₀ = 0, a₁ = 1, kwargs...)
	if l.useFactorization
		Jfact = factorize(_axpy(J, a₀, a₁))
		return Jfact \ rhs1, Jfact \ rhs2, true, (1, 1)
	else
		_J = _axpy(J, a₀, a₁)
		return _J \ rhs1, _J \ rhs2, true, (1, 1)
	end
end
####################################################################################################
# Solvers for IterativeSolvers
####################################################################################################
"""
$(TYPEDEF)
Linear solver based on gmres from `IterativeSolvers.jl`. Can be used to solve `(a₀ * I + a₁ * J) * x = rhs`.
$(TYPEDFIELDS)
"""
@with_kw mutable struct GMRESIterativeSolvers{T, Tl, Tr} <: AbstractLinearSolver
	"Absolute tolerance for solver"
	abstol::T = 0.0

	"Relative tolerance for solver"
	reltol::T = 1e-8

	"Number of restarts"
	restart::Int64 = 200

	"Maximum number of iterations"
	maxiter::Int64 = 100

	"Dimension of the problem"
	N::Int64 = 0

	"Display information during iterations"
	verbose::Bool = false

	"Record information"
	log::Bool = true

	"Start with zero guess"
	initially_zero::Bool = true

	"Left preconditioner"
	Pl::Tl = IterativeSolvers.Identity()

	"Right preconditioner"
	Pr::Tr = IterativeSolvers.Identity()

	"Whether the linear operator is written inplace"
	ismutating::Bool = false
end

# this function is used to solve (a₀ * I + a₁ * J) * x = rhs
# the optional shift is only used for the Hopf Newton / Continuation
function (l::GMRESIterativeSolvers{T, Tl, Tr})(J, rhs; a₀ = 0, a₁ = 1, kwargs...) where {T, Ts, Tl, Tr}
	# no need to use fancy axpy! here because IterativeSolvers "only" handles AbstractArray
	if l.ismutating == true
		Jmap = LinearMap{T}((o, v) -> J(o, v), l.N, l.N ; ismutating = true)
		@assert ((a₀ == 0) && (a₁ == 1)) "Perturbed inplace linear problem not done yet!"
	else
		J_map = v -> _axpy_op(J, v, a₀, a₁)
		Jmap = LinearMap{T}(J_map, l.N, l.N ; ismutating = false)
	end
	res = IterativeSolvers.gmres(Jmap, rhs; abstol = l.abstol, reltol = l.reltol, log = l.log, verbose = l.verbose, restart = l.restart, maxiter = l.maxiter, initially_zero = l.initially_zero, Pl = l.Pl, Pr = l.Pr, kwargs...)
	(res[2].iters >= l.maxiter) && (@warn "IterativeSolvers.gmres iterated maxIter = $(res[2].iters) times without achieving the desired tolerance.\n")
	return res[1], length(res) > 1, res[2].iters
end
####################################################################################################
# Solvers for KrylovKit
####################################################################################################
"""
$(TYPEDEF)

Create a linear solver based on GMRES from `KrylovKit.jl`. Can be used to solve `(a₀ * I + a₁ * J) * x = rhs`.
$(TYPEDFIELDS)
!!! tip "Different linear solvers"
    By tuning the options, you can select CG, GMRES... see [here](https://jutho.github.io/KrylovKit.jl/stable/man/linear/#KrylovKit.linsolve)
"""
@with_kw mutable struct GMRESKrylovKit{T, Tl} <: AbstractLinearSolver
	"Krylov Dimension"
	dim::Int64 = KrylovDefaults.krylovdim

	"Absolute tolerance for solver"
	atol::T  = KrylovDefaults.tol

	"Relative tolerance for solver"
	rtol::T  = KrylovDefaults.tol

	"Maximum number of iterations"
	maxiter::Int64 = KrylovDefaults.maxiter

	"Verbosity ∈ {0,1,2}"
	verbose::Int64 = 0

	"If the linear map is symmetric, only meaningful if T<:Real"
	issymmetric::Bool = false

	"If the linear map is hermitian"
	ishermitian::Bool = false

	"If the linear map is positive definite"
	isposdef::Bool = false

	"Left preconditioner"
	Pl::Tl = nothing
end

# this function is used to solve (a₀ * I + a₁ * J) * x = rhs
# the optional shift is only used for the Hopf Newton / Continuation
function (l::GMRESKrylovKit{T, Tl})(J, rhs; a₀ = 0, a₁ = 1, kwargs...) where {T, Tl}
	if Tl == Nothing
		res, info = KrylovKit.linsolve(J, rhs, a₀, a₁; rtol = l.rtol, verbosity = l.verbose, krylovdim = l.dim, maxiter = l.maxiter, atol = l.atol, issymmetric = l.issymmetric, ishermitian = l.ishermitian, isposdef = l.isposdef, kwargs...)
	else # use preconditioner
		res, info = KrylovKit.linsolve(x -> (out = apply(J, x); ldiv!(l.Pl, out)), ldiv!(l.Pl, copy(rhs)), a₀, a₁; rtol = l.rtol, verbosity = l.verbose, krylovdim = l.dim, maxiter = l.maxiter, atol = l.atol, issymmetric = l.issymmetric, ishermitian = l.ishermitian, isposdef = l.isposdef, kwargs...)
	end
	info.converged == 0 && (@warn "KrylovKit.linsolve solver did not converge")
	return res, true, info.numops
end


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   BondHamiltonian{L,N}
#   --> L is the label type of bonds
#   --> N is the dimension the bond term matrix (NxN matrix)
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
abstract type AbstractBondHamiltonian{L,N} end

# export the type
export AbstractBondHamiltonian

# DOCSTRING
"""
    abstract type AbstractBondHamiltonian{L,N}

The abstract bond Hamiltonian type that describes all hamiltonian implementations on bonds. It is parametric in two types,
- `L` the label type of the bonds it covers
- `N` an integer number specifying the number of on-site degrees of freedom (e.g. number of spin components)

The abstract type is equipped with a single interface function `bondterm` which
returns the Hamiltonian matrix describing the Hamiltonian on that bond.
The function throws an adequate error when called but not implemented for concrete types.
See specific documentation for the function itself.

# Examples

```julia-REPL
julia> AbstractBondHamiltonian{Int64,3} <: AbstractBondHamiltonian{L,N} where{L,N}
true
```

"""
AbstractBondHamiltonian

################################################################################
#
#	INTERFACING / ACCESSING BOND HAMILTONIANS
#	(functions have to be overwritten by concrete types)
#
################################################################################

# Fallback if wrong label type is passed (DOES NOT HAVE TO BE OVERWRITTEN)
function bondterm(
            h :: AbstractBondHamiltonian{LH,NH},
            b :: AbstractBond{LB,NB}
        ) :: Matrix{Complex{Float64}} where {LB,NB,LH,NH}

    # print an error because implementation for concrete type is missing
    error("Passed a bond with different label type to a bond Hamiltonian\nLabel type of bond: "*string(LB)*", label type of bond Hamiltonian: "*string(LH))
end

# get bond term of Hamiltonian
function bondterm(
            h :: AbstractBondHamiltonian{L,NH},
            b :: AbstractBond{L,NB}
        ) :: Matrix{Complex{Float64}} where {L,NB,NH}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'bondterm' for bond Hamiltonian type " * string(typeof(h)))
end

export bondterm

# DOCSTRING
"""
    function bondterm(
            h :: AbstractBondHamiltonian{L,NH},
            b :: AbstractBond{L,NB}
        ) :: Matrix{Complex{Float64}} where {L,NB,NH}

The abstract bond Hamiltonian interface function. Returns a matrix with complex
entries for a passed `AbstractBond` object `b`.

E.g. for Hamiltonians of the form \$\\sum_{\\langle ij\\rangle} S_{i} J_{ij} S_{j} \\dots\$, this bond hamiltonian matrix would correspond
to the matrix \$J_{ij}\$.

# Examples

```julia-REPL
julia> bondterm(hamiltonian, b)
3×3 Array{Complex,2}:
 1.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  1.0+0.0im
```

"""
bondterm

# interface for saving and loading
function saveBondHamiltonian(
        hb :: HB,
        fn :: AbstractString,
        group :: AbstractString = "bond_hamiltonian"
        ;
        append :: Bool = false
    ) where {L,N,HB<:AbstractBondHamiltonian{L,N}}

    # print an error because implementation for concrete type is missing
    error("not implemented interface function 'saveBondHamiltonian' for bond Hamiltonian type " * string(typeof(hb)))
end

# convinience function for standard type
function loadBondHamiltonian(
        fn :: AbstractString,
        group :: AbstractString = "bond_hamiltonian"
    )

    # find out the type
    T = getBondHamiltonianType(Val(Symbol(h5readattr(fn, group)["type"])))

    # return the loaded bond hamiltonian
    return loadBondHamiltonian(T, fn, group)
end

function getBondHamiltonianType(
        ::Val{HB}
    ) where {HB}

    error("Type $(HB) could not be identified, i.e. getBondHamiltonianType(Val{:$(HB)}) is missing")
end

export saveBondHamiltonian, loadBondHamiltonian, getBondHamiltonianType


text:
using Gomah
using Gomah: L, np, @pywith, chainercv
using Flux
using Statistics

const DTYPE = Float32

using Test

function test_ch2dense()
    # get instance of Linear
    INSIZE = 10
    OUTSIZE = 20
    BSIZE = 1
    chlinear = L.Linear(in_size = INSIZE, out_size = OUTSIZE)
    dummyX = 128 * np.ones((BSIZE, INSIZE), dtype = np.float32)
    chret = reversedims(chlinear(dummyX).array)
    fldense = ch2dense(chlinear)
    flret = fldense(reversedims(dummyX))
    @test size(flret) == size(chret)
    @show maximum(abs.(flret .- chret))
    @show typeof.((fldense.W, fldense.b))
    @test all(isapprox.(flret, chret))
end

function test_ch2conv()
    # get instance of Convolution2D
    INCH = 3
    OUTCH = 8
    inH = 30
    inW = 30
    BSIZE = 1
    KSIZE = 3
    pad = 1
    s = 1
    chconv = Gomah.L.Convolution2D(
        in_channels = INCH,
        out_channels = OUTCH,
        ksize = KSIZE,
        pad = pad,
        stride = s,
    )
    # pass dummy data to create computation graph
    dummyX = np.ones((BSIZE, INCH, inH, inW), dtype = np.float32)
    chret = reversedims(chconv(dummyX).array)
    flconv = ch2conv(chconv)
    flret = flconv(ones(DTYPE, inW, inH, INCH, BSIZE))
    @test size(flret) == size(chret)
    @show maximum(abs.(flret .- chret))
    @show typeof.((flconv.weight, flconv.bias))
    @test all(isapprox.(flret, chret))
end

function test_ch2conv_nobias()
    # get instance of Convolution2D
    INCH = 3
    OUTCH = 8
    inH = 30
    inW = 30
    BSIZE = 1
    KSIZE = 3
    pad = 1
    s = 1
    chconv = Gomah.L.Convolution2D(
        ksize = KSIZE,
        pad = pad,
        stride = s,
        in_channels = INCH,
        out_channels = OUTCH,
        nobias = true,
    )
    # pass dummy data to create computation graph
    dummyX = 128 * np.ones((BSIZE, INCH, inH, inW), dtype = np.float32)
    chret = reversedims(chconv(dummyX).array)
    flconv = ch2conv(chconv)
    flret = flconv(reversedims(dummyX))
    @test size(flret) == size(chret)
    @show maximum(abs.(flret .- chret))
    @test all(isapprox.(abs.(flret - chret), 0, atol = 1e-4))
end

function test_resnetconv1ch2conv()
    # get instance of Convolution2D
    INCH = 2
    OUTCH = 3
    inH = 7
    inW = 7
    BSIZE = 1
    KSIZE = 7
    pad = 3
    s = 2
    chconv = Gomah.L.Convolution2D(
        ksize = KSIZE,
        pad = pad,
        stride = s,
        in_channels = INCH,
        out_channels = OUTCH,
        nobias = true,
    )
    # pass dummy data to create computation graph
    dummyX = 128 * np.ones((BSIZE, INCH, inH, inW), dtype = np.float32)
    chret = reversedims(chconv(dummyX).array)
    flconv = ch2conv(chconv)
    flret = flconv(reversedims(dummyX))
    @test size(flret) == size(chret)
    @show maximum(abs.(flret .- chret))
    @test all(isapprox.(abs.(flret - chret), 0, atol = 1e-4))
end

function test_ch2dwconv()
    # get instance of DepthwiseConvolution2D
    INCH = 10
    MULTIPLIER = 2
    inH = 10
    inW = 10
    BSIZE = 2
    KSIZE = 3
    pad = 1
    s = 1
    chdwconv = Gomah.L.DepthwiseConvolution2D(
        in_channels = INCH,
        channel_multiplier = MULTIPLIER,
        ksize = KSIZE,
        pad = pad,
        stride = s
    )
    # pass dummy data to create computation graph
    dummyX = 128 * np.ones((BSIZE, INCH, inH, inW), dtype = np.float32)
    chret = reversedims(chdwconv(dummyX).array)
    fldwconv = ch2dwconv(chdwconv)
    flret = fldwconv(reversedims(dummyX))
    @test size(flret) == size(chret)
    @show maximum(abs.(flret .- chret))
    @show typeof.((fldwconv.weight, fldwconv.bias))
    @test all(isapprox.(abs.(flret - chret), 0, atol = 1e-4))
end

function test_ch2bn()
    SIZE = 10
    BSIZE = 1
    dummyX = 128 * np.ones((BSIZE, SIZE), dtype = np.float32)
    chbn = L.BatchNormalization(size = SIZE)
    @pywith chainer.using_config("train", false) begin
        chret = reversedims(chbn(dummyX).array)
        flbn = ch2bn(chbn)
        Flux.testmode!(flbn, true)
        flret = flbn(reversedims(dummyX))
        @test size(flret) == size(chret)
        @show maximum(abs.(flret .- chret))
        @test all(isapprox.(abs.(flret - chret), 0, atol = 1e-4))
    end
end

@testset "bn parameter test" begin
    PyResNet = chainercv.links.model.resnet.ResNet
    num = 50
    pyresnet = PyResNet(num, pretrained_model = "imagenet")
    dummyX = rand(Float32, (1, 64, 10, 10))
    dummyX = zeros(Float32, (1, 64, 10, 10))
    pybn = pyresnet.conv1.bn
    flbn = ch2bn(pybn)
    Flux.testmode!(flbn, true)
    @pywith chainer.using_config("train", false) begin
        chret = reversedims(pybn(dummyX).array)
        flret = flbn(reversedims(dummyX))

        β = flbn.β
        γ = flbn.γ
        μ = flbn.μ
        σ² = flbn.σ²
        ϵ = flbn.ϵ
        @show typeof.((β, γ, μ, σ², ϵ))
        @test isapprox(β, pybn.beta.array)
        @test isapprox(γ, pybn.gamma.array)
        @test isapprox(μ, pybn.avg_mean)
        @test isapprox(σ², pybn.avg_var)
        @test isapprox(ϵ, pybn.eps)
    end
end

@testset "converter" begin
    test_ch2dense()
    test_ch2conv()
    test_ch2conv_nobias()
    test_resnetconv1ch2conv()
    test_ch2dwconv()
    test_ch2bn()
end


text:
#
# Data structures to represent primal and dual solutions
#
mutable struct OPFPrimalSolution <: AbstractPrimalSolution
    Pg
    Qg
    Vm
    Va
    ωt
    St
    Zt
    Sk
    Zk
    sigma_real
    sigma_imag
    sigma_lineFr
    sigma_lineTo

    function OPFPrimalSolution(opfdata::OPFData, modelinfo::ModelInfo)
        buses = opfdata.buses
        gens  = opfdata.generators
        ngen  = length(gens)
        nbus  = length(buses)
        nline = length(opfdata.lines)
        T = modelinfo.num_time_periods
        K = modelinfo.num_ctgs + 1 # base case counted separately
        @assert T >= 1
        @assert K >= 1

        Pg = zeros(ngen,K,T)
        Qg = zeros(ngen,K,T)
        Vm = zeros(nbus,K,T)
        Va = zeros(nbus,K,T)
        ωt = zeros(K,T)
        St = zeros(ngen,T)
        Zt = zeros(ngen,T)
        Sk = zeros(ngen,K,T)
        Zk = zeros(ngen,K,T)
        if modelinfo.allow_constr_infeas
            sigma_real = zeros(nbus,K,T)
            sigma_imag = zeros(nbus,K,T)
            sigma_lineFr = zeros(nline,K,T)
            sigma_lineTo = zeros(nline,K,T)
        else
            sigma_real = 0
            sigma_imag = 0
            sigma_lineFr = 0
            sigma_lineTo = 0
        end

        for i=1:ngen
            Pg[i,:,:] .= 0.5*(gens[i].Pmax + gens[i].Pmin)
            Qg[i,:,:] .= 0.5*(gens[i].Qmax + gens[i].Qmin)
        end
        for i=1:nbus
            Vm[i,:,:] .= 0.5*(buses[i].Vmax + buses[i].Vmin)
            Va[i,:,:] .= opfdata.buses[opfdata.bus_ref].Va
        end
        if T > 1
            for i=1:ngen
                St[i,2:T] .= gens[i].ramp_agc
            end
        end
        if K > 1
            for i=1:ngen
                Sk[i,2:K,:] .= gens[i].scen_agc
            end
        end

        new(Pg,Qg,Vm,Va,ωt,St,Zt,Sk,Zk,sigma_real,sigma_imag,sigma_lineFr,sigma_lineTo)
    end
end

function OPFPrimalSolution(nlp::AbstractNLPEvaluator)
    return OPFPrimalSolution(nlp.opfdata, nlp.modelinfo)
end

mutable struct OPFDualSolution <: AbstractDualSolution
    ramping
    ctgs

    function OPFDualSolution(opfdata::OPFData, modelinfo::ModelInfo)
        ngen = length(opfdata.generators)
        T = modelinfo.num_time_periods
        K = modelinfo.num_ctgs + 1 # base case counted separately
        @assert T >= 1
        @assert K >= 1

        ramping = zeros(ngen,T)
        ctgs = zeros(ngen,K,T)

        new(ramping,ctgs)
    end
end

function OPFDualSolution(nlp::AbstractNLPEvaluator)
    return OPFDualSolution(nlp.opfdata, nlp.modelinfo)
end

function get_block_view(x::OPFPrimalSolution,
                        block::CartesianIndex,
                        modelinfo::ModelInfo,
                        algparams::AlgParams)
    k = block[1]
    t = block[2]
    range_k = algparams.decompCtgs ? (k:k) : (1:(modelinfo.num_ctgs + 1))
    Pg = view(x.Pg, :, range_k, t)
    Qg = view(x.Qg, :, range_k, t)
    Vm = view(x.Vm, :, range_k, t)
    Va = view(x.Va, :, range_k, t)
    ωt = view(x.ωt, range_k, t)
    St = view(x.St, :, t)
    Zt = view(x.Zt, :, t)
    Sk = view(x.Sk, :, range_k, t)
    Zk = view(x.Zk, :, range_k, t)

    return CatView(Pg, Qg, Vm, Va, ωt, St, Zt, Sk, Zk)
end

function get_coupling_view(x::OPFPrimalSolution,
                           modelinfo::ModelInfo,
                           algparams::AlgParams)
    Pg = @view x.Pg[:]
    return Pg
end

function get_coupling_view(λ::OPFDualSolution,
                           modelinfo::ModelInfo,
                           algparams::AlgParams)
    ramping = @view λ.ramping[:]
    ctgs = @view λ.ctgs[:]

    if algparams.decompCtgs
        return CatView(ramping, ctgs)
    else
        return CatView(ramping)
    end
end

function dist(x1::OPFPrimalSolution, x2::OPFPrimalSolution,
              modelinfo::ModelInfo,
              algparams::AlgParams,
              lnorm = Inf)
    x1vec = get_coupling_view(x1, modelinfo, algparams)
    x2vec = get_coupling_view(x2, modelinfo, algparams)
    return norm(x1vec .- x2vec, lnorm) #/(1e-16 + norm(x2vec, lnorm))
end

function dist(λ1::OPFDualSolution, λ2::OPFDualSolution,
              modelinfo::ModelInfo,
              algparams::AlgParams,
              lnorm = Inf)
    λ1vec = get_coupling_view(λ1, modelinfo, algparams)
    λ2vec = get_coupling_view(λ2, modelinfo, algparams)
    return norm(λ1vec .- λ2vec, lnorm) #/(1e-16 + norm(λ2vec, lnorm))
end

max_stars_repo_name: exanauts/proxALM

max_issues_repo_name: exanauts/proxALM

max_forks_repo_name: exanauts/proxALM