text:
# Tensor factorization
#----------------------
const OFA = OrthogonalFactorizationAlgorithm

import LinearAlgebra: svd!, svd
const SVDAlg = Union{SVD, SDD}

Base.@deprecate(
    svd(t::AbstractTensorMap, leftind::IndexTuple, rightind::IndexTuple;
            trunc::TruncationScheme = notrunc(), p::Real = 2, alg::SVDAlg = SDD()),
    tsvd(t, leftind, rightind; trunc = trunc, p = p, alg = alg))
Base.@deprecate(
    svd(t::AbstractTensorMap;
            trunc::TruncationScheme = notrunc(), p::Real = 2, alg::SVDAlg = SDD()),
    tsvd(t; trunc = trunc, p = p, alg = alg))
Base.@deprecate(
    svd!(t::AbstractTensorMap;
            trunc::TruncationScheme = notrunc(), p::Real = 2, alg::SVDAlg = SDD()),
    tsvd(t; trunc = trunc, p = p, alg = alg))

"""
    tsvd(t::AbstractTensorMap, leftind::Tuple, rightind::Tuple;
        trunc::TruncationScheme = notrunc(), p::Real = 2, alg::Union{SVD, SDD} = SDD())
        -> U, S, V, ϵ

Compute the (possibly truncated)) singular value decomposition such that
`norm(permute(t, leftind, rightind) - U * S *V) ≈ ϵ`, where `ϵ` thus represents the truncation error.

If `leftind` and `rightind` are not specified, the current partition of left and right
indices of `t` is used. In that case, less memory is allocated if one allows the data in
`t` to be destroyed/overwritten, by using `tsvd!(t, trunc = notrunc(), p = 2)`.

A truncation parameter `trunc` can be specified for the new internal dimension, in which
case a truncated singular value decomposition will be computed. Choices are:
*   `notrunc()`: no truncation (default);
*   `truncerr(η::Real)`: truncates such that the p-norm of the truncated singular values is
    smaller than `η` times the p-norm of all singular values;
*   `truncdim(χ::Int)`: truncates such that the equivalent total dimension of the internal
    vector space is no larger than `χ`;
*   `truncspace(V)`: truncates such that the dimension of the internal vector space is
    smaller than that of `V` in any sector.
*   `truncbelow(χ::Real)`: truncates such that every singular value is larger then `χ` ;

The method `tsvd` also returns the truncation error `ϵ`, computed as the `p` norm of the
singular values that were truncated.

THe keyword `alg` can be equal to `SVD()` or `SDD()`, corresponding to the underlying LAPACK
algorithm that computes the decomposition (`_gesvd` or `_gesdd`).

Orthogonality requires `spacetype(t)<:InnerProductSpace`, and `svd(!)` is currently
only implemented for `spacetype(t)<:EuclideanSpace`.
"""
tsvd(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) =
    tsvd!(permute(t, p1, p2; copy = true); kwargs...)

"""
    leftorth(t::AbstractTensorMap, leftind::Tuple, rightind::Tuple;
                alg::OrthogonalFactorizationAlgorithm = QRpos()) -> Q, R

Create orthonormal basis `Q` for indices in `leftind`, and remainder `R` such that
`permute(t, leftind, rightind) = Q*R`.

If `leftind` and `rightind` are not specified, the current partition of left and right
indices of `t` is used. In that case, less memory is allocated if one allows the data in `t`
to be destroyed/overwritten, by using `leftorth!(t, alg = QRpos())`.

Different algorithms are available, namely `QR()`, `QRpos()`, `SVD()` and `Polar()`. `QR()`
and `QRpos()` use a standard QR decomposition, producing an upper triangular matrix `R`.
`Polar()` produces a Hermitian and positive semidefinite `R`. `QRpos()` corrects the
standard QR decomposition such that the diagonal elements of `R` are positive. Only
`QRpos()` and `Polar()` are uniqe (no residual freedom) so that they always return the same
result for the same input tensor `t`.

Orthogonality requires `spacetype(t)<:InnerProductSpace`, and `leftorth(!)` is currently
only implemented for `spacetype(t)<:EuclideanSpace`.
"""
leftorth(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) =
    leftorth!(permute(t, p1, p2; copy = true); kwargs...)

"""
    rightorth(t::AbstractTensorMap, leftind::Tuple, rightind::Tuple;
                alg::OrthogonalFactorizationAlgorithm = LQpos()) -> L, Q

Create orthonormal basis `Q` for indices in `rightind`, and remainder `L` such that
`permute(t, leftind, rightind) = L*Q`.

If `leftind` and `rightind` are not specified, the current partition of left and right
indices of `t` is used. In that case, less memory is allocated if one allows the data in `t`
to be destroyed/overwritten, by using `rightorth!(t, alg = LQpos())`.

Different algorithms are available, namely `LQ()`, `LQpos()`, `RQ()`, `RQpos()`, `SVD()` and
`Polar()`. `LQ()` and `LQpos()` produce a lower triangular matrix `L` and are computed using
a QR decomposition of the transpose. `RQ()` and `RQpos()` produce an upper triangular
remainder `L` and only works if the total left dimension is smaller than or equal to the
total right dimension. `LQpos()` and `RQpos()` add an additional correction such that the
diagonal elements of `L` are positive. `Polar()` produces a Hermitian and positive
semidefinite `L`. Only `LQpos()`, `RQpos()` and `Polar()` are uniqe (no residual freedom) so
that they always return the same result for the same input tensor `t`.

Orthogonality requires `spacetype(t)<:InnerProductSpace`, and `rightorth(!)` is currently
only implemented for `spacetype(t)<:EuclideanSpace`.
"""
rightorth(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) =
    rightorth!(permute(t, p1, p2; copy = true); kwargs...)

"""
    leftnull(t::AbstractTensor, leftind::Tuple, rightind::Tuple;
                alg::OrthogonalFactorizationAlgorithm = QRpos()) -> N

Create orthonormal basis for the orthogonal complement of the support of the indices in
`leftind`, such that `N' * permute(t, leftind, rightind) = 0`.

If `leftind` and `rightind` are not specified, the current partition of left and right
indices of `t` is used. In that case, less memory is allocated if one allows the data in `t`
to be destroyed/overwritten, by using `leftnull!(t, alg = QRpos())`.

Different algorithms are available, namely `QR()` (or equivalently, `QRpos()`), `SVD()` and
`SDD()`. The first assumes that the matrix is full rank and requires `iszero(atol)` and
`iszero(rtol)`. With `SVD()` and `SDD()`, `rightnull` will use the corresponding singular
value decomposition, and one can specify an absolute or relative tolerance for which
singular values are to be considered zero, where `max(atol, norm(t)*rtol)` is used as upper
bound.

Orthogonality requires `spacetype(t)<:InnerProductSpace`, and `leftnull(!)` is currently
only implemented for `spacetype(t)<:EuclideanSpace`.
"""
leftnull(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) =
    leftnull!(permute(t, p1, p2; copy = true); kwargs...)

"""
    rightnull(t::AbstractTensor, leftind::Tuple, rightind::Tuple;
                alg::OrthogonalFactorizationAlgorithm = LQ(),
                atol::Real = 0.0,
                rtol::Real = eps(real(float(one(eltype(t)))))*iszero(atol)) -> N

Create orthonormal basis for the orthogonal complement of the support of the indices in
`rightind`, such that `permute(t, leftind, rightind)*N' = 0`.

If `leftind` and `rightind` are not specified, the current partition of left and right
indices of `t` is used. In that case, less memory is allocated if one allows the data in `t`
to be destroyed/overwritten, by using `rightnull!(t, alg = LQpos())`.

Different algorithms are available, namely `LQ()` (or equivalently, `LQpos`), `SVD()` and
`SDD()`. The first assumes that the matrix is full rank and requires `iszero(atol)` and
`iszero(rtol)`. With `SVD()` and `SDD()`, `rightnull` will use the corresponding singular
value decomposition, and one can specify an absolute or relative tolerance for which
singular values are to be considered zero, where `max(atol, norm(t)*rtol)` is used as upper
bound.

Orthogonality requires `spacetype(t)<:InnerProductSpace`, and `rightnull(!)` is currently
only implemented for `spacetype(t)<:EuclideanSpace`.
"""
rightnull(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) =
    rightnull!(permute(t, p1, p2; copy = true); kwargs...)

"""
    eigen(t::AbstractTensor, leftind::Tuple, rightind::Tuple; kwargs...) -> D, V

Compute eigenvalue factorization of tensor `t` as linear map from `rightind` to `leftind`.

If `leftind` and `rightind` are not specified, the current partition of left and right
indices of `t` is used. In that case, less memory is allocated if one allows the data in `t`
to be destroyed/overwritten, by using `eigen!(t)`. Note that the permuted tensor on which
`eigen!` is called should have equal domain and codomain, as otherwise the eigenvalue
decomposition is meaningless and cannot satisfy
```
permute(t, leftind, rightind) * V = V * D
```

Accepts the same keyword arguments `scale`, `permute` and `sortby` as `eigen` of dense
matrices. See the corresponding documentation for more information.

See also `eig` and `eigh`
"""
LinearAlgebra.eigen(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) =
    eigen!(permute(t, p1, p2; copy = true); kwargs...)

"""
    eig(t::AbstractTensor, leftind::Tuple, rightind::Tuple; kwargs...) -> D, V

Compute eigenvalue factorization of tensor `t` as linear map from `rightind` to `leftind`.
The function `eig` assumes that the linear map is not hermitian and returns type stable
complex valued `D` and `V` tensors for both real and complex valued `t`. See `eigh` for
hermitian linear maps

If `leftind` and `rightind` are not specified, the current partition of left and right
indices of `t` is used. In that case, less memory is allocated if one allows the data in
`t` to be destroyed/overwritten, by using `eig!(t)`. Note that the permuted tensor on
which `eig!` is called should have equal domain and codomain, as otherwise the eigenvalue
decomposition is meaningless and cannot satisfy
```
permute(t, leftind, rightind) * V = V * D
```

Accepts the same keyword arguments `scale`, `permute` and `sortby` as `eigen` of dense matrices. See the corresponding documentation for more information.

See also `eigen` and `eigh`.
"""
eig(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple; kwargs...) =
    eig!(permute(t, p1, p2; copy = true); kwargs...)

"""
    eigh(t::AbstractTensorMap{<:EuclideanSpace}, leftind::Tuple, rightind::Tuple) -> D, V

Compute eigenvalue factorization of tensor `t` as linear map from `rightind` to `leftind`.
The function `eigh` assumes that the linear map is hermitian and `D` and `V` tensors with
the same `eltype` as `t`. See `eig` and `eigen` for non-hermitian tensors. Hermiticity
requires that the tensor acts on inner product spaces, and the current implementation
requires `spacetyp(t) <: EuclideanSpace`.

If `leftind` and `rightind` are not specified, the current partition of left and right
indices of `t` is used. In that case, less memory is allocated if one allows the data in
`t` to be destroyed/overwritten, by using `eigh!(t)`. Note that the permuted tensor on
which `eigh!` is called should have equal domain and codomain, as otherwise the eigenvalue
decomposition is meaningless and cannot satisfy
```
permute(t, leftind, rightind) * V = V * D
```

See also `eigen` and `eig`.
"""
eigh(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple) =
    eigh!(permute(t, p1, p2; copy = true))

"""
    isposdef(t::AbstractTensor{<:EuclideanSpace}, leftind::Tuple, rightind::Tuple) -> ::Bool

Test whether a tensor `t` is positive definite as linear map from `rightind` to `leftind`.

If `leftind` and `rightind` are not specified, the current partition of left and right
indices of `t` is used. In that case, less memory is allocated if one allows the data in
`t` to be destroyed/overwritten, by using `isposdef!(t)`. Note that the permuted tensor on
which `isposdef!` is called should have equal domain and codomain, as otherwise it is
meaningless

Accepts the same keyword arguments `scale`, `permute` and `sortby` as `eigen` of dense
matrices. See the corresponding documentation for more information.
"""
LinearAlgebra.isposdef(t::AbstractTensorMap, p1::IndexTuple, p2::IndexTuple) =
    isposdef!(permute(t, p1, p2; copy = true))

tsvd(t::AbstractTensorMap; trunc::TruncationScheme = NoTruncation(),
                            p::Real = 2, alg::Union{SVD, SDD} = SDD()) =
    tsvd!(copy(t); trunc = trunc, p = p, alg = alg)
leftorth(t::AbstractTensorMap; alg::OFA = QRpos(), kwargs...) =
    leftorth!(copy(t); alg = alg, kwargs...)
rightorth(t::AbstractTensorMap; alg::OFA = LQpos(), kwargs...) =
    rightorth!(copy(t); alg = alg, kwargs...)
leftnull(t::AbstractTensorMap; alg::OFA = QR(), kwargs...) =
    leftnull!(copy(t); alg = alg, kwargs...)
rightnull(t::AbstractTensorMap; alg::OFA = LQ(), kwargs...) =
    rightnull!(copy(t); alg = alg, kwargs...)
LinearAlgebra.eigen(t::AbstractTensorMap; kwargs...) = eigen!(copy(t); kwargs...)
eig(t::AbstractTensorMap; kwargs...) = eig!(copy(t); kwargs...)
eigh(t::AbstractTensorMap; kwargs...) = eigh!(copy(t); kwargs...)
LinearAlgebra.isposdef(t::AbstractTensorMap) = isposdef!(copy(t))

# Orthogonal factorizations (mutation for recycling memory):
# only correct if Euclidean inner product
#------------------------------------------------------------------------------------------
leftorth!(t::AdjointTensorMap{S}; alg::OFA = QRpos()) where {S<:EuclideanSpace} =
    map(adjoint, reverse(rightorth!(adjoint(t); alg = alg')))

rightorth!(t::AdjointTensorMap{S}; alg::OFA = LQpos()) where {S<:EuclideanSpace} =
    map(adjoint, reverse(leftorth!(adjoint(t); alg = alg')))

leftnull!(t::AdjointTensorMap{S}; alg::OFA = QR(), kwargs...) where {S<:EuclideanSpace} =
    adjoint(rightnull!(adjoint(t); alg = alg', kwargs...))

rightnull!(t::AdjointTensorMap{S}; alg::OFA = LQ(), kwargs...) where {S<:EuclideanSpace} =
    adjoint(leftnull!(adjoint(t); alg = alg', kwargs...))

function tsvd!(t::AdjointTensorMap{S};
                trunc::TruncationScheme = NoTruncation(),
                p::Real = 2,
                alg::Union{SVD, SDD} = SDD()) where {S<:EuclideanSpace}
    u, s, vt, err = tsvd!(adjoint(t); trunc = trunc, p = p, alg = alg)
    return adjoint(vt), adjoint(s), adjoint(u), err
end

function leftorth!(t::TensorMap{<:EuclideanSpace};
                    alg::Union{QR, QRpos, QL, QLpos, SVD, SDD, Polar} = QRpos(),
                    atol::Real = zero(float(real(eltype(t)))),
                    rtol::Real = (alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) :
                    eps(real(float(one(eltype(t)))))*iszero(atol))
    if !iszero(rtol)
        atol = max(atol, rtol*norm(t))
    end
    I = sectortype(t)
    S = spacetype(t)
    A = storagetype(t)
    Qdata = SectorDict{I, A}()
    Rdata = SectorDict{I, A}()
    dims = SectorDict{I, Int}()
    for c in blocksectors(domain(t))
        isempty(block(t,c)) && continue
        Q, R = _leftorth!(block(t, c), alg, atol)
        Qdata[c] = Q
        Rdata[c] = R
        dims[c] = size(Q, 2)
    end
    V = S(dims)
    if alg isa Polar
        @assert V ≅ domain(t)
        W = domain(t)
    elseif length(domain(t)) == 1 && domain(t) ≅ V
        W = domain(t)
    elseif length(codomain(t)) == 1 && codomain(t) ≅ V
        W = codomain(t)
    else
        W = ProductSpace(V)
    end
    return TensorMap(Qdata, codomain(t)←W), TensorMap(Rdata, W←domain(t))
end

function leftnull!(t::TensorMap{<:EuclideanSpace};
                    alg::Union{QR, QRpos, SVD, SDD} = QRpos(),
                    atol::Real = zero(float(real(eltype(t)))),
                    rtol::Real = (alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) :
                    eps(real(float(one(eltype(t)))))*iszero(atol))
    if !iszero(rtol)
        atol = max(atol, rtol*norm(t))
    end
    I = sectortype(t)
    S = spacetype(t)
    A = storagetype(t)
    V = codomain(t)
    Ndata = SectorDict{I, A}()
    dims = SectorDict{I, Int}()
    for c in blocksectors(V)
        N = _leftnull!(block(t, c), alg, atol)
        Ndata[c] = N
        dims[c] = size(N, 2)
    end
    W = S(dims)
    return TensorMap(Ndata, V←W)
end

function rightorth!(t::TensorMap{<:EuclideanSpace};
                    alg::Union{LQ, LQpos, RQ, RQpos, SVD, SDD, Polar} = LQpos(),
                    atol::Real = zero(float(real(eltype(t)))),
                    rtol::Real = (alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) :
                    eps(real(float(one(eltype(t)))))*iszero(atol))
    if !iszero(rtol)
        atol = max(atol, rtol*norm(t))
    end
    I = sectortype(t)
    S = spacetype(t)
    A = storagetype(t)
    Ldata = SectorDict{I, A}()
    Qdata = SectorDict{I, A}()
    dims = SectorDict{I, Int}()
    for c in blocksectors(codomain(t))
        isempty(block(t,c)) && continue
        L, Q = _rightorth!(block(t, c), alg, atol)
        Ldata[c] = L
        Qdata[c] = Q
        dims[c] = size(Q, 1)
    end
    V = S(dims)
    if alg isa Polar
        @assert V ≅ codomain(t)
        W = codomain(t)
    elseif length(codomain(t)) == 1 && codomain(t) ≅ V
        W = codomain(t)
    elseif length(domain(t)) == 1 && domain(t) ≅ V
        W = domain(t)
    else
        W = ProductSpace(V)
    end
    return TensorMap(Ldata, codomain(t)←W), TensorMap(Qdata, W←domain(t))
end

function rightnull!(t::TensorMap{<:EuclideanSpace};
                    alg::Union{LQ, LQpos, SVD, SDD} = LQpos(),
                    atol::Real = zero(float(real(eltype(t)))),
                    rtol::Real = (alg ∉ (SVD(), SDD())) ? zero(float(real(eltype(t)))) :
                    eps(real(float(one(eltype(t)))))*iszero(atol))
    if !iszero(rtol)
        atol = max(atol, rtol*norm(t))
    end
    I = sectortype(t)
    S = spacetype(t)
    A = storagetype(t)
    V = domain(t)
    Ndata = SectorDict{I, A}()
    dims = SectorDict{I, Int}()
    for c in blocksectors(V)
        N = _rightnull!(block(t, c), alg, atol)
        Ndata[c] = N
        dims[c] = size(N, 1)
    end
    W = S(dims)
    return TensorMap(Ndata, W←V)
end

function tsvd!(t::TensorMap{<:EuclideanSpace};
                trunc::TruncationScheme = NoTruncation(),
                p::Real = 2,
                alg::Union{SVD, SDD} = SDD())
    S = spacetype(t)
    I = sectortype(t)
    A = storagetype(t)
    Ar = similarstoragetype(t, real(eltype(t)))
    Udata = SectorDict{I, A}()
    Σmdata = SectorDict{I, Ar}() # this will contain the singular values as matrix
    Vdata = SectorDict{I, A}()
    dims = SectorDict{sectortype(t), Int}()
    if isempty(blocksectors(t))
        W = S(dims)
        truncerr = zero(real(eltype(t)))
        return TensorMap(Udata, codomain(t)←W), TensorMap(Σmdata, W←W),
                    TensorMap(Vdata, W←domain(t)), truncerr
    end
    for (c, b) in blocks(t)
        U, Σ, V = _svd!(b, alg)
        Udata[c] = U
        Vdata[c] = V
        if @isdefined Σdata # cannot easily infer the type of Σ, so use this construction
            Σdata[c] = Σ
        else
            Σdata = SectorDict(c=>Σ)
        end
        dims[c] = length(Σ)
    end
    if !isa(trunc, NoTruncation)
        Σdata, truncerr = _truncate!(Σdata, trunc, p)
        truncdims = SectorDict{I, Int}()
        for c in blocksectors(t)
            truncdim = length(Σdata[c])
            if truncdim != 0
                truncdims[c] = truncdim
                if truncdim != dims[c]
                    Udata[c] = Udata[c][:, 1:truncdim]
                    Vdata[c] = Vdata[c][1:truncdim, :]
                end
            else
                delete!(Udata, c)
                delete!(Vdata, c)
                delete!(Σdata, c)
            end
        end
        dims = truncdims
        W = S(dims)
    else
        W = S(dims)
        if length(domain(t)) == 1 && domain(t)[1] ≅ W
            W = domain(t)[1]
        elseif length(codomain(t)) == 1 && codomain(t)[1] ≅ W
            W = codomain(t)[1]
        end
        truncerr = abs(zero(eltype(t)))
    end
    for (c, Σ) in Σdata
        Σmdata[c] = copyto!(similar(Σ, length(Σ), length(Σ)), Diagonal(Σ))
    end
    return TensorMap(Udata, codomain(t)←W), TensorMap(Σmdata, W←W),
            TensorMap(Vdata, W←domain(t)), truncerr
end

function LinearAlgebra.ishermitian(t::TensorMap)
    domain(t) == codomain(t) || return false
    spacetype(t) <: EuclideanSpace || return false # hermiticity only defined for euclidean
    for (c, b) in blocks(t)
        ishermitian(b) || return false
    end
    return true
end

LinearAlgebra.eigen!(t::TensorMap) = ishermitian(t) ? eigh!(t) : eig!(t)

function eigh!(t::TensorMap{<:EuclideanSpace}; kwargs...)
    domain(t) == codomain(t) ||
        throw(SpaceMismatch("`eigh!` requires domain and codomain to be the same"))
    S = spacetype(t)
    I = sectortype(t)
    A = storagetype(t)
    Ar = similarstoragetype(t, real(eltype(t)))
    Ddata = SectorDict{I, Ar}()
    Vdata = SectorDict{I, A}()
    dims = SectorDict{I, Int}()
    for (c, b) in blocks(t)
        values, vectors = eigen!(Hermitian(b); kwargs...)
        d = length(values)
        Ddata[c] = copyto!(similar(values, (d, d)), Diagonal(values))
        Vdata[c] = vectors
        dims[c] = d
    end
    if length(domain(t)) == 1
        W = domain(t)[1]
    else
        W = S(dims)
    end
    return TensorMap(Ddata, W←W), TensorMap(Vdata, domain(t)←W)
end

function eig!(t::TensorMap; kwargs...)
    domain(t) == codomain(t) ||
        throw(SpaceMismatch("`eig!` requires domain and codomain to be the same"))
    S = spacetype(t)
    I = sectortype(t)
    T = complex(eltype(t))
    Ac = similarstoragetype(t, T)
    Ddata = SectorDict{I, Ac}()
    Vdata = SectorDict{I, Ac}()
    dims = SectorDict{I, Int}()
    for (c, b) in blocks(t)
        values, vectors = eigen!(b; kwargs...)
        d = length(values)
        Ddata[c] = copyto!(similar(values, T, (d, d)), Diagonal(values))
        if eltype(vectors) == T
            Vdata[c] = vectors
        else
            Vdata[c] = copyto!(similar(vectors, T), vectors)
        end
        dims[c] = d
    end
    if length(domain(t)) == 1
        W = domain(t)[1]
    else
        W = S(dims)
    end
    return TensorMap(Ddata, W←W), TensorMap(Vdata, domain(t)←W)
end

function LinearAlgebra.isposdef!(t::TensorMap)
    domain(t) == codomain(t) ||
        throw(SpaceMismatch("`isposdef` requires domain and codomain to be the same"))
    spacetype(t) <: EuclideanSpace || return false
    for (c, b) in blocks(t)
        isposdef!(b) || return false
    end
    return true
end


text:
### Exports ###

export GenericNConfig, GenericNState, GenericNeurons

### Types ###

GenericNConfig = NConfig{Int}
const GenericNState = NState{Vector{Float32}, Vector{Bool}}
const GenericNeurons = Neurons{GenericNConfig, GenericNState}
GenericNeurons(size::Int; a=0.02, b=0.2, c=-65, d=8, thresh=30, θRate=0.5, traceRate=0.5, boostRate=0.0) =
  GenericNeurons(
    GenericNConfig(
      size,
      a,
      b,
      c,
      d,
      thresh,
      θRate,
      traceRate,
      boostRate),
    GenericNState(
      fill(c, size),
      zeros(Float32, size),
      zeros(Float32, size),
      zeros(Bool, size),
      zeros(Float32, size),
      zeros(Float32, size),
      zeros(Float32, size)))

### Methods ###

function Base.show(io::IO, neurons::GenericNeurons)
  print(io, "GenericNeurons ($(neurons.conf.size))")
end

Base.getindex(neurons::GenericNeurons, idx) =
  neurons.state.F[idx]

# TODO: Make this actually random?
function randnode(::Type{GenericNeurons}, size::Int=rand(1:20))
  GenericNeurons(size)
end

ext: jl


text:
using Test

@testset "AlgebraicNets" begin
  include("AlgebraicNets.jl")
end

@testset "MathFormulas" begin
  include("MathFormulas.jl")
end

@testset "MarkovCategories" begin
  include("Markov.jl")
end


text:
Erm = [5.2091e-8 + 4.7554e-10im 5.1973e-8 + 1.3643e-9im 5.1613e-8 + 2.9788e-9im
4.2874e-9 + 1.7963e-10im 4.1918e-9 + 4.7889e-10im 3.9374e-9 + 9.3153e-10im
1.0775e-9 + 9.9565e-11im 1.0018e-9 + 2.4289e-10im 8.2961e-10 + 4.0922e-10im
4.0849e-10 + 6.3458e-11im 3.4952e-10 + 1.3979e-10im 2.3699e-10 + 1.9941e-10im
1.9052e-10 + 4.3411e-11im 1.4511e-10 + 8.5273e-11im 7.4032e-11 + 1.0108e-10im
1.0036e-10 + 3.1001e-11im 6.576e-11 + 5.3654e-11im 2.2352e-11 + 5.2092e-11im
5.7189e-11 + 2.2803e-11im 3.1081e-11 + 3.4386e-11im 5.4567e-12 + 2.7148e-11im
3.4415e-11 + 1.7153e-11im 1.4873e-11 + 2.2317e-11im 2.4467e-13 + 1.4389e-11im
2.1557e-11 + 1.3141e-11im 7.0349e-12 + 1.4637e-11im -1.0647e-12 + 7.8713e-12im
1.3929e-11 + 1.0226e-11im 3.2008e-12 + 9.7079e-12im -1.1879e-12 + 4.5351e-12im
9.2269e-12 + 8.0675e-12im 1.341e-12 + 6.5276e-12im -1.039e-12 + 2.7974e-12im
6.2413e-12 + 6.443e-12im 4.6648e-13 + 4.468e-12im -8.8411e-13 + 1.8526e-12im
4.2983e-12 + 5.2025e-12im 8.0958e-14 + 3.1295e-12im -7.7447e-13 + 1.3001e-12im
3.0059e-12 + 4.2424e-12im -6.9763e-14 + 2.2547e-12im -7.0025e-13 + 9.4357e-13im
2.1314e-12 + 3.4901e-12im -1.1167e-13 + 1.6785e-12im -6.4217e-13 + 6.9387e-13im
1.5301e-12 + 2.8938e-12im -1.0869e-13 + 1.2942e-12im -5.8859e-13 + 5.0816e-13im
1.1108e-12 + 2.4164e-12im -9.1134e-14 + 1.0329e-12im -5.3408e-13 + 3.6698e-13im
8.1514e-13 + 2.0307e-12im -7.2337e-14 + 8.5115e-13im -4.7832e-13 + 2.5928e-13im
6.0454e-13 + 1.7165e-12im -5.7319e-14 + 7.2048e-13im -4.2302e-13 + 1.7835e-13im
4.5317e-13 + 1.4587e-12im -4.7451e-14 + 6.2302e-13im -3.7003e-13 + 1.1862e-13im
3.439e-13 + 1.2459e-12im -4.198e-14 + 5.4771e-13im -3.2117e-13 + 7.539e-14im
2.6426e-13 + 1.0693e-12im -4.0044e-14 + 4.8692e-13im -2.775e-13 + 4.4841e-14im
2.0623e-13 + 9.2216e-13im -4.0497e-14 + 4.3646e-13im -2.3919e-13 + 2.3657e-14im
1.6372e-13 + 7.9896e-13im -4.2373e-14 + 3.9333e-13im -2.0623e-13 + 9.1698e-15im
1.3239e-13 + 6.9545e-13im -4.4978e-14 + 3.5557e-13im -1.7824e-13 - 5.9266e-16im
1.0936e-13 + 6.0824e-13im -4.7833e-14 + 3.2214e-13im -1.5451e-13 - 7.0435e-15im
9.229e-14 + 5.3452e-13im -5.0562e-14 + 2.9218e-13im -1.3444e-13 - 1.1312e-14im
7.951e-14 + 4.7201e-13im -5.2943e-14 + 2.6511e-13im -1.1747e-13 - 1.416e-14im
6.9914e-14 + 4.1888e-13im -5.4899e-14 + 2.4057e-13im -1.0303e-13 - 1.6004e-14im
6.261e-14 + 3.7358e-13im -5.6378e-14 + 2.1828e-13im -9.0672e-14 - 1.717e-14im
5.6883e-14 + 3.3483e-13im -5.7353e-14 + 1.9799e-13im -7.9995e-14 - 1.7912e-14im
5.2295e-14 + 3.0158e-13im -5.7858e-14 + 1.795e-13im -7.0718e-14 - 1.8344e-14im
4.8529e-14 + 2.7296e-13im -5.7942e-14 + 1.6265e-13im -6.2613e-14 - 1.853e-14im
4.5319e-14 + 2.4824e-13im -5.765e-14 + 1.4729e-13im -5.5482e-14 - 1.8525e-14im
4.2455e-14 + 2.2676e-13im -5.7023e-14 + 1.3331e-13im -4.9166e-14 - 1.8372e-14im
3.9831e-14 + 2.0805e-13im -5.6113e-14 + 1.2057e-13im -4.3554e-14 - 1.8096e-14im
3.7375e-14 + 1.9167e-13im -5.497e-14 + 1.0896e-13im -3.856e-14 - 1.7717e-14im
3.5028e-14 + 1.7726e-13im -5.3634e-14 + 9.8398e-14im -3.411e-14 - 1.7251e-14im
3.274e-14 + 1.6451e-13im -5.2145e-14 + 8.8783e-14im -3.0135e-14 - 1.6711e-14im
3.0487e-14 + 1.5316e-13im -5.0534e-14 + 8.0035e-14im -2.6583e-14 - 1.6112e-14im
2.8268e-14 + 1.4302e-13im -4.8831e-14 + 7.2079e-14im -2.3409e-14 - 1.5468e-14im
2.6086e-14 + 1.339e-13im -4.7062e-14 + 6.4843e-14im -2.0575e-14 - 1.4792e-14im
2.3939e-14 + 1.2566e-13im -4.525e-14 + 5.8266e-14im -1.8047e-14 - 1.4095e-14im
2.183e-14 + 1.1816e-13im -4.3414e-14 + 5.2291e-14im -1.5794e-14 - 1.3386e-14im
1.9763e-14 + 1.113e-13im -4.157e-14 + 4.6864e-14im -1.3788e-14 - 1.2673e-14im
1.7748e-14 + 1.05e-13im -3.9732e-14 + 4.1939e-14im -1.2004e-14 - 1.1964e-14im
1.5793e-14 + 9.9185e-14im -3.791e-14 + 3.747e-14im -1.0419e-14 - 1.1263e-14im
1.3906e-14 + 9.3799e-14im -3.6113e-14 + 3.3418e-14im -9.0132e-15 - 1.0577e-14im
1.2094e-14 + 8.8785e-14im -3.4351e-14 + 2.9746e-14im -7.7682e-15 - 9.9101e-15im
1.036e-14 + 8.4095e-14im -3.2631e-14 + 2.6421e-14im -6.6671e-15 - 9.2655e-15im
8.71e-15 + 7.9695e-14im -3.0957e-14 + 2.3412e-14im -5.695e-15 - 8.6458e-15im
7.145e-15 + 7.5558e-14im -2.9332e-14 + 2.0692e-14im -4.8384e-15 - 8.0518e-15im
5.6671e-15 + 7.1661e-14im -2.776e-14 + 1.8234e-14im -4.0853e-15 - 7.4845e-15im
4.2776e-15 + 6.7983e-14im -2.6243e-14 + 1.6017e-14im -3.4245e-15 - 6.9447e-15im
2.9778e-15 + 6.4503e-14im -2.4782e-14 + 1.4017e-14im -2.8457e-15 - 6.4328e-15im
1.7683e-15 + 6.1205e-14im -2.3379e-14 + 1.2216e-14im -2.3399e-15 - 5.9494e-15im
6.4854e-16 + 5.8072e-14im -2.2036e-14 + 1.0596e-14im -1.8987e-15 - 5.4943e-15im
-3.8461e-16 + 5.5095e-14im -2.075e-14 + 9.1407e-15im -1.5153e-15 - 5.0667e-15im
-1.3344e-15 + 5.2266e-14im -1.9522e-14 + 7.8351e-15im -1.1835e-15 - 4.6655e-15im
-2.2038e-15 + 4.9575e-14im -1.8351e-14 + 6.6658e-15im -8.9736e-16 - 4.2897e-15im
-2.9958e-15 + 4.7015e-14im -1.7235e-14 + 5.6203e-15im -6.5176e-16 - 3.9385e-15im
-3.7131e-15 + 4.4579e-14im -1.6174e-14 + 4.6871e-15im -4.4179e-16 - 3.6109e-15im
-4.3585e-15 + 4.2259e-14im -1.5166e-14 + 3.8555e-15im -2.6302e-16 - 3.3061e-15im
-4.9347e-15 + 4.005e-14im -1.4211e-14 + 3.1157e-15im -1.1143e-16 - 3.0232e-15im
-5.4459e-15 + 3.7947e-14im -1.3307e-14 + 2.4594e-15im 1.612e-17 - 2.7611e-15im
-5.8964e-15 + 3.5945e-14im -1.2451e-14 + 1.8787e-15im 1.2237e-16 - 2.5185e-15im
-6.2903e-15 + 3.4038e-14im -1.1642e-14 + 1.3668e-15im 2.0984e-16 - 2.2941e-15im
-6.6314e-15 + 3.2224e-14im -1.0878e-14 + 9.1697e-16im 2.8086e-16 - 2.0869e-15im
-6.9234e-15 + 3.0498e-14im -1.0157e-14 + 5.2318e-16im 3.3758e-16 - 1.8958e-15im
-7.1695e-15 + 2.8856e-14im -9.4778e-15 + 1.7975e-16im 3.82e-16 - 1.7198e-15im
-7.3728e-15 + 2.7295e-14im -8.838e-15 - 1.1857e-16im 4.1596e-16 - 1.5581e-15im
-7.5364e-15 + 2.5812e-14im -8.2364e-15 - 3.7663e-16im 4.4115e-16 - 1.4097e-15im
-7.6632e-15 + 2.4402e-14im -7.6714e-15 - 5.9852e-16im 4.5883e-16 - 1.2738e-15im
-7.7565e-15 + 2.3064e-14im -7.1408e-15 - 7.8772e-16im 4.6997e-16 - 1.1494e-15im
-7.8192e-15 + 2.1792e-14im -6.643e-15 - 9.4746e-16im 4.7545e-16 - 1.0357e-15im
-7.8541e-15 + 2.0585e-14im -6.1762e-15 - 1.0808e-15im 4.7608e-16 - 9.3199e-16im
-7.8638e-15 + 1.944e-14im -5.7388e-15 - 1.1904e-15im 4.7263e-16 - 8.374e-16im
-7.8507e-15 + 1.8353e-14im -5.3291e-15 - 1.2791e-15im 4.6581e-16 - 7.5127e-16im
-7.8171e-15 + 1.7323e-14im -4.9459e-15 - 1.3493e-15im 4.5628e-16 - 6.7294e-16im
-7.765e-15 + 1.6346e-14im -4.5877e-15 - 1.4032e-15im 4.4465e-16 - 6.018e-16im
-7.6966e-15 + 1.5421e-14im -4.2533e-15 - 1.443e-15im 4.3149e-16 - 5.3727e-16im
-7.6135e-15 + 1.4545e-14im -3.9412e-15 - 1.4705e-15im 4.1725e-16 - 4.7881e-16im
-7.5174e-15 + 1.3716e-14im -3.6503e-15 - 1.4871e-15im 4.0218e-16 - 4.2594e-16im
-7.4098e-15 + 1.2931e-14im -3.3792e-15 - 1.4941e-15im 3.8647e-16 - 3.7821e-16im
-7.2922e-15 + 1.2188e-14im -3.1268e-15 - 1.4925e-15im 3.7031e-16 - 3.352e-16im
-7.1659e-15 + 1.1486e-14im -2.8918e-15 - 1.4835e-15im 3.5388e-16 - 2.9651e-16im
-7.032e-15 + 1.0821e-14im -2.6733e-15 - 1.4681e-15im 3.3734e-16 - 2.6176e-16im
-6.8918e-15 + 1.0193e-14im -2.4703e-15 - 1.4471e-15im 3.2084e-16 - 2.3061e-16im
-6.7464e-15 + 9.5997e-15im -2.2816e-15 - 1.4216e-15im 3.0452e-16 - 2.0272e-16im
-6.5968e-15 + 9.0389e-15im -2.1066e-15 - 1.3924e-15im 2.8852e-16 - 1.7777e-16im
-6.4439e-15 + 8.5094e-15im -1.9442e-15 - 1.3601e-15im 2.7295e-16 - 1.5547e-16im
-6.2887e-15 + 8.0096e-15im -1.7937e-15 - 1.3256e-15im 2.5794e-16 - 1.3555e-16im
-6.1318e-15 + 7.5379e-15im -1.6543e-15 - 1.2894e-15im 2.4352e-16 - 1.1777e-16im
-5.9739e-15 + 7.093e-15im -1.5252e-15 - 1.2518e-15im 2.297e-16 - 1.0196e-16im
-5.8153e-15 + 6.6733e-15im -1.4059e-15 - 1.2131e-15im 2.1647e-16 - 8.7953e-17im
-5.6565e-15 + 6.2777e-15im -1.2957e-15 - 1.1736e-15im 2.0383e-16 - 7.5578e-17im
-5.498e-15 + 5.9048e-15im -1.194e-15 - 1.1335e-15im 1.9176e-16 - 6.4689e-17im
-5.3402e-15 + 5.5534e-15im -1.1001e-15 - 1.0932e-15im 1.8028e-16 - 5.514e-17im
-5.1835e-15 + 5.2224e-15im -1.0136e-15 - 1.0528e-15im 1.6936e-16 - 4.6795e-17im
-5.0283e-15 + 4.9106e-15im -9.3383e-16 - 1.0126e-15im 1.5901e-16 - 3.9526e-17im
-4.8748e-15 + 4.6172e-15im -8.6043e-16 - 9.7283e-16im 1.4922e-16 - 3.3209e-17im
-4.7234e-15 + 4.3409e-15im -7.9288e-16 - 9.3365e-16im 1.3999e-16 - 2.773e-17im
-4.5744e-15 + 4.081e-15im -7.3072e-16 - 8.9526e-16im 1.313e-16 - 2.2978e-17im
-4.4281e-15 + 3.8365e-15im -6.7355e-16 - 8.5783e-16im 1.2316e-16 - 1.8851e-17im
-4.2847e-15 + 3.6065e-15im -6.2098e-16 - 8.2144e-16im 1.1554e-16 - 1.5274e-17im
-4.1441e-15 + 3.3904e-15im -5.727e-16 - 7.8613e-16im 1.0842e-16 - 1.2196e-17im
-4.0066e-15 + 3.1872e-15im -5.2843e-16 - 7.5188e-16im 1.0176e-16 - 9.572e-18im
-3.872e-15 + 2.9964e-15im -4.8786e-16 - 7.1872e-16im 9.5549e-17 - 7.3563e-18im
-3.7406e-15 + 2.8172e-15im -4.5073e-16 - 6.8664e-16im 8.975e-17 - 5.5068e-18im
-3.6123e-15 + 2.649e-15im -4.1679e-16 - 6.5566e-16im 8.4342e-17 - 3.9834e-18im
-3.4872e-15 + 2.4911e-15im -3.8578e-16 - 6.2578e-16im 7.93e-17 - 2.7478e-18im
-3.3652e-15 + 2.343e-15im -3.5747e-16 - 5.9699e-16im 7.4604e-17 - 1.7636e-18im
-3.2466e-15 + 2.204e-15im -3.3163e-16 - 5.6931e-16im 7.023e-17 - 9.9619e-19im
-3.1312e-15 + 2.0737e-15im -3.0806e-16 - 5.4274e-16im 6.616e-17 - 4.1275e-19im
-3.0191e-15 + 1.9515e-15im -2.8655e-16 - 5.1728e-16im 6.2372e-17 + 1.8045e-20im
-2.9103e-15 + 1.8369e-15im -2.669e-16 - 4.9293e-16im 5.8849e-17 + 3.2597e-19im
-2.8049e-15 + 1.7295e-15im -2.4893e-16 - 4.6969e-16im 5.5572e-17 + 5.3934e-19im
-2.7029e-15 + 1.6287e-15im -2.3247e-16 - 4.4756e-16im 5.2526e-17 + 6.8275e-19im
-2.6041e-15 + 1.5343e-15im -2.1741e-16 - 4.2649e-16im 4.9695e-17 + 7.6676e-19im
-2.5086e-15 + 1.4458e-15im -2.0366e-16 - 4.0645e-16im 4.7067e-17 + 7.9708e-19im
-2.4162e-15 + 1.363e-15im -1.9113e-16 - 3.8739e-16im 4.463e-17 + 7.7914e-19im
-2.3269e-15 + 1.2856e-15im -1.7974e-16 - 3.6927e-16im 4.2371e-17 + 7.1814e-19im
-2.2405e-15 + 1.2133e-15im -1.6941e-16 - 3.5205e-16im 4.0281e-17 + 6.1903e-19im
-2.157e-15 + 1.1457e-15im -1.6007e-16 - 3.3569e-16im 3.8347e-17 + 4.8657e-19im
-2.0763e-15 + 1.0827e-15im -1.5164e-16 - 3.2016e-16im 3.656e-17 + 3.2528e-19im
-1.9984e-15 + 1.0239e-15im -1.4405e-16 - 3.0541e-16im 3.491e-17 + 1.3949e-19im
-1.9232e-15 + 9.6918e-16im -1.3724e-16 - 2.9142e-16im 3.3387e-17 - 6.6654e-20im
-1.8506e-15 + 9.182e-16im -1.3114e-16 - 2.7815e-16im 3.1981e-17 - 2.8921e-19im
-1.7805e-15 + 8.7077e-16im -1.2569e-16 - 2.6557e-16im 3.0685e-17 - 5.244e-19im
-1.713e-15 + 8.2666e-16im -1.2083e-16 - 2.5364e-16im 2.949e-17 - 7.6861e-19im
-1.6478e-15 + 7.8567e-16im -1.1651e-16 - 2.4235e-16im 2.8387e-17 - 1.0184e-18im
-1.585e-15 + 7.4758e-16im -1.1267e-16 - 2.3166e-16im 2.737e-17 - 1.2704e-18im
-1.5245e-15 + 7.1221e-16im -1.0926e-16 - 2.2154e-16im 2.643e-17 - 1.5216e-18im
-1.4663e-15 + 6.7937e-16im -1.0623e-16 - 2.1197e-16im 2.5561e-17 - 1.7689e-18im
-1.4103e-15 + 6.4888e-16im -1.0354e-16 - 2.0293e-16im 2.4755e-17 - 2.0095e-18im
-1.3564e-15 + 6.2056e-16im -1.0113e-16 - 1.9438e-16im 2.4006e-17 - 2.2405e-18im
-1.3047e-15 + 5.9424e-16im -9.8963e-17 - 1.863e-16im 2.3309e-17 - 2.4595e-18im
-1.255e-15 + 5.6977e-16im -9.6999e-17 - 1.7867e-16im 2.2656e-17 - 2.6638e-18im
-1.2072e-15 + 5.47e-16im -9.5198e-17 - 1.7148e-16im 2.2042e-17 - 2.8511e-18im
-1.1615e-15 + 5.2577e-16im -9.3519e-17 - 1.6469e-16im 2.1461e-17 - 3.019e-18im
-1.1176e-15 + 5.0594e-16im -9.1926e-17 - 1.5829e-16im 2.0909e-17 - 3.1655e-18im
-1.0756e-15 + 4.8738e-16im -9.0384e-17 - 1.5226e-16im 2.0379e-17 - 3.2885e-18im
-1.0354e-15 + 4.6996e-16im -8.8858e-17 - 1.4658e-16im 1.9868e-17 - 3.3859e-18im
-9.9699e-16 + 4.5354e-16im -8.7316e-17 - 1.4123e-16im 1.937e-17 - 3.4558e-18im
-9.6032e-16 + 4.3801e-16im -8.5727e-17 - 1.3619e-16im 1.888e-17 - 3.4965e-18im
-9.2534e-16 + 4.2325e-16im -8.4061e-17 - 1.3145e-16im 1.8396e-17 - 3.5062e-18im
-8.9202e-16 + 4.0915e-16im -8.2289e-17 - 1.2698e-16im 1.7911e-17 - 3.4833e-18im
-8.6033e-16 + 3.9559e-16im -8.0383e-17 - 1.2278e-16im 1.7423e-17 - 3.4262e-18im
-8.3022e-16 + 3.8248e-16im -7.8319e-17 - 1.1883e-16im 1.6927e-17 - 3.3335e-18im
-8.0165e-16 + 3.6972e-16im -7.6069e-17 - 1.1512e-16im 1.642e-17 - 3.2036e-18im]

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
@testset "metric" begin
    p = MeterProvider()
    e = InMemoryExporter()
    r = MetricReader(p, e)

    m = Meter("test"; provider = p)

    c = Counter{Int}("fruit_counter", m)
    c(; name = "apple", color = "red")
    c(2; name = "lemon", color = "yellow")
    c(1; name = "lemon", color = "yellow")
    c(2; name = "apple", color = "green")
    c(5; name = "apple", color = "red")
    c(4; name = "lemon", color = "yellow")

    @test length(metrics()) == 0  # the default meter provider is a dummy one

    @test length(metrics(p)) == 1
    mt = first(metrics(p))
    @test length(collect(mt)) == 3

    @test resource(mt) == resource(p)

    udc = UpDownCounter{Int}("stock", m)
    udc(200)
    udc(-50)
    udc(150)

    function fake_observer(x)
        i = x
        function ()
            res = i * 2
            i += 1
            res
        end
    end

    oc = ObservableCounter{Int}(fake_observer(0), "observable_counter", m)

    oudc = ObservableUpDownCounter{Int}(fake_observer(-2), "observable_up_down_counter", m)

    og = ObservableGauge{Int}(fake_observer(-3), "observable_gauge", m)

    h = Histogram{Int}("n_requests", m)
    h(3)
    h(5)
    h(8)
    h(12)
    h(2)
    h(1_000_000)

    r()

    c_points = first(m for m in e if m.name == "fruit_counter")
    @test length(c_points) == 3
    @test c_points[StaticAttrs((color = "red", name = "apple"))].value == 6
    @test c_points[StaticAttrs((color = "yellow", name = "lemon"))].value == 7
    @test c_points[StaticAttrs((color = "green", name = "apple"))].value == 2

    h_points = first(m for m in e if m.name == "n_requests")
    @test length(h_points) == 1
    hist = h_points[StaticAttrs()].value.counts
    @test hist[1] == 0
    @test hist[2] == 3
    @test hist[3] == 1
    @test hist[4] == 1
    @test hist[end] == 1
    @test sum(hist) == 6

    udc_points = first(m for m in e if m.name == "stock")
    @test udc_points[StaticAttrs()].value == 300

    oc_points = first(m for m in e if m.name == "observable_counter")
    @test oc_points[StaticAttrs()].value == 0

    oudc_points = first(m for m in e if m.name == "observable_up_down_counter")
    @test oudc_points[StaticAttrs()].value == -4

    og_points = first(m for m in e if m.name == "observable_gauge")
    @test og_points[StaticAttrs()].value == -6

    # observe again!
    # sync instruments should have no changes
    # async instruments should have been changed
    empty!(e)
    r()

    c_points = first(m for m in e if m.name == "fruit_counter")
    @test length(c_points) == 3
    @test c_points[StaticAttrs((color = "red", name = "apple"))].value == 6
    @test c_points[StaticAttrs((color = "yellow", name = "lemon"))].value == 7
    @test c_points[StaticAttrs((color = "green", name = "apple"))].value == 2

    h_points = first(m for m in e if m.name == "n_requests")
    @test length(h_points) == 1
    hist = h_points[StaticAttrs()].value.counts
    @test hist[1] == 0
    @test hist[2] == 3
    @test hist[3] == 1
    @test hist[4] == 1
    @test hist[end] == 1
    @test sum(hist) == 6

    udc_points = first(m for m in e if m.name == "stock")
    @test udc_points[StaticAttrs()].value == 300

    oc_points = first(m for m in e if m.name == "observable_counter")
    @test oc_points[StaticAttrs()].value == 0 + 2

    oudc_points = first(m for m in e if m.name == "observable_up_down_counter")
    @test oudc_points[StaticAttrs()].value == -4 + (-2)

    og_points = first(m for m in e if m.name == "observable_gauge")
    @test og_points[StaticAttrs()].value == -4 # the last value

    # observe again
    empty!(e)
    r()

    oc_points = first(m for m in e if m.name == "observable_counter")
    @test oc_points[StaticAttrs()].value == 0 + 2 + 4

    oudc_points = first(m for m in e if m.name == "observable_up_down_counter")
    @test oudc_points[StaticAttrs()].value == -4 + (-2) + 0

    og_points = first(m for m in e if m.name == "observable_gauge")
    @test og_points[StaticAttrs()].value == -2 # the last value

    periodic_reader = PeriodicMetricReader(
        CompositMetricReader(r, MetricReader(p, ConsoleExporter()));
        export_interval_seconds = 1,
    )

    periodic_reader()

    sleep(3)
    shut_down!(periodic_reader)

    global_meter_provider!(p)
    mr = MetricReader()  # read to console by default
end

@testset "AggregationStore" begin
    store = OpenTelemetrySDK.AggregationStore{OpenTelemetrySDK.DataPoint{Int,Nothing}}(;
        n_max_points = 3,
        n_max_attrs = 4,
    )
    init_data = () -> OpenTelemetrySDK.DataPoint{Int}()

    p1 = get!(init_data, store, StaticAttrs())
    p2 = get!(init_data, store, StaticAttrs((a = 1, b = 2, c = 3)))
    p3 = get!(init_data, store, StaticAttrs((b = 2, a = 1, c = 3)))
    @test p2 === p3

    p4 = get!(init_data, store, StaticAttrs((b = 2, c = 3, a = 1)))
    @test p2 === p4

    @test (@test_logs (
        :warn,
        "maximum cached keys in agg store reached, please consider increase `n_max_points`",
    ) get!(init_data, store, StaticAttrs((c = 3, b = 2, a = 1)))) === p2

    p5 = get!(init_data, store, StaticAttrs((; x = 123)))
    p6 = get!(init_data, store, StaticAttrs((; y = 123)))
    @test isnothing(p6)
end

@testset "Views" begin
    b = (5.0, 10.0, 25.0, 50.0, 100.0)

    p = MeterProvider(
        views = [
            View("X"),
            View("Foo"; instrument_name = "Y"),
            View(
                "Bar";
                instrument_name = "Y",
                aggregation = HistogramAgg{Int}(boundaries = b),
            ),
        ],
    )

    e = InMemoryExporter()
    r = MetricReader(p, e)
    m = Meter("test"; provider = p)

    X = Counter{Int}("X", m)
    X(2)
    X(3)
    X(1)

    Y = Histogram{Int}("Y", m)
    Y(-1)
    Y(1)
    Y(2000)
    Y(200)

    Z = UpDownCounter{Int}("Z", m)
    Z(-10)
    Z(10)

    r()

    @test length(e.pool) == 3

    @test first(m for m in e if m.name == "X")[StaticAttrs()].value == 2 + 3 + 1

    foo = first(m for m in e if m.name == "Foo")[StaticAttrs()].value.counts
    # DEFAULT_HISTOGRAM_BOUNDARIES = (0.0, 5.0, 10.0, 25.0, 50.0, 75.0, 100.0, 250.0, 500.0, 1000.0)
    @test foo[1] == 1
    @test foo[2] == 1
    @test foo[8] == 1
    @test foo[end] == 1
    @test sum(foo) == 4

    bar = first(m for m in e if m.name == "Bar")[StaticAttrs()].value.counts
    # b = (5.0, 10.0, 25.0, 50.0, 100.0)
    @test bar[1] == 2
    @test bar[end] == 2
    @test sum(bar) == 4
end

@testset "Drop Agg" begin
    p = MeterProvider(
        views = [
            View(
                "*";
                instrument_type = Counter,
                meter_name = "test",
                meter_version = v"0.0.1-dev",
                meter_schema_url = "",
            ),
            View("X"; aggregation = DROP),
        ],
    )

    e = InMemoryExporter()
    r = MetricReader(p, e)

    m = Meter("test"; provider = p)
    x = Counter{Int}("X", m)
    x(1)

    y = Counter{Int}("Y", m)
    tracer = Tracer()
    with_span("test", tracer) do
        y(2)
    end

    r()

    @test length(e) == 1
    @test e.pool[1].name == "Y"
end


text:
include("/home/chad/jlTools/Slice.jl")

using StatsFuns
using Plots
using Distributions
pyplot()

function randSphere(d)
    x = randn(d)
    return x / vecnorm(x)
end

function hitrunStep(logf, x0, y0)
    v = randSphere(size(x0))
    logf1 = k ->  logf((x0 .+ k .* v))
    (k, y) = sliceStep(logf1, 0.0, y0, 1.0)
    x = x0 .+ k .* v
    (x, y)
end

function test(n)
    xs = zeros(n,2)
    logf = x -> logpdf(MvNormal(2,1.0), x)
    x = zeros(2)
    y = logf(x)
    for i in 1:n
        (x,y) = hitrunStep(logf, x, y)
        xs[i,:] = x
    end
    xs
end

x = test(1000)


text:
export SeisData

# This is type-stable for S = SeisData() but not for keyword args
@doc """
    SeisData

A custom structure designed to contain the minimum necessary information for
processing univariate geophysical data.

    SeisChannel

A single channel designed to contain the minimum necessary information for
processing univariate geophysical data.

## Fields

| **Field** | **Description**                                               |
|:-------|:------                                                           |
| :n     | Number of channels [^1]                                          |
| :c     | TCP connections feeding data to this object [^1]                 |
| :id    | Channel id. Uses NET.STA.LOC.CHA format when possible            |
| :name  | Freeform channel name                                            |
| :loc   | Location (position) vector; any subtype of InstrumentPosition    |
| :fs    | Sampling frequency in Hz; fs=0.0 for irregularly-sampled data.   |
| :gain  | Scalar gain                                                      |
| :resp  | Instrument response; any subtype of InstrumentResponse           |
| :units | String describing data units. UCUM standards are assumed.        |
| :src   | Freeform string describing data source.                          |
| :misc  | Dictionary for non-critical information.                         |
| :notes | Timestamped notes; includes automatically-logged information.    |
| :t     | Matrix of time gaps in integer μs, formatted [Sample# Length]    |
| :x     | Time-series data                                                 |

[^1]: Not present in SeisChannel objects.

See also: `InstrumentPosition`, `PZResp`
""" SeisData
mutable struct SeisData <: GphysData
  n::Int64
  id::Array{String,1}                 # id
  name::Array{String,1}               # name
  loc::Array{InstrumentPosition,1}    # loc
  fs::Array{Float64,1}                # fs
  gain::Array{Float64,1}              # gain
  resp::Array{InstrumentResponse,1}   # resp
  units::Array{String,1}              # units
  src::Array{String,1}                # src
  misc::Array{Dict{String,Any},1}     # misc
  notes::Array{Array{String,1},1}     # notes
  t::Array{Array{Int64,2},1}          # time
  x::Array{FloatArray,1}              # data
  c::Array{TCPSocket,1}               # connections

  function SeisData()
    return new(0,
                Array{String,1}(undef,0),
                Array{String,1}(undef,0),
                Array{InstrumentPosition,1}(undef,0),
                Array{Float64,1}(undef,0),
                Array{Float64,1}(undef,0),
                Array{InstrumentResponse,1}(undef,0),
                Array{String,1}(undef,0),
                Array{String,1}(undef,0),
                Array{Dict{String,Any},1}(undef,0),
                Array{Array{String,1},1}(undef,0),
                Array{Array{Int64,2},1}(undef,0),
                Array{FloatArray,1}(undef,0),
                Array{TCPSocket,1}(undef,0)
              )
  end

  function SeisData( n::Int64,
            id::Array{String,1}                 , # id
            name::Array{String,1}               , # name
            loc::Array{InstrumentPosition,1}    , # loc
            fs::Array{Float64,1}                , # fs
            gain::Array{Float64,1}              , # gain
            resp::Array{InstrumentResponse,1}   , # resp
            units::Array{String,1}              , # units
            src::Array{String,1}                , # src
            misc::Array{Dict{String,Any},1}     , # misc
            notes::Array{Array{String,1},1}     , # notes
            t::Array{Array{Int64,2},1}          , # time
            x::Array{FloatArray,1})

    return new(n,
      id, name, loc, fs, gain, resp, units, src, misc, notes, t, x,
      Array{TCPSocket,1}(undef,0)
      )
  end

  function SeisData(n::UInt)
    S = new(n,
                Array{String,1}(undef,n),
                Array{String,1}(undef,n),
                Array{InstrumentPosition,1}(undef,n),
                Array{Float64,1}(undef,n),
                Array{Float64,1}(undef,n),
                Array{InstrumentResponse,1}(undef,n),
                Array{String,1}(undef,n),
                Array{String,1}(undef,n),
                Array{Dict{String,Any},1}(undef,n),
                Array{Array{String,1},1}(undef,n),
                Array{Array{Int64,2},1}(undef,n),
                Array{FloatArray,1}(undef,n),
                Array{TCPSocket,1}(undef,0)
              )

    # Fill these fields with something to prevent undefined reference errors
    fill!(S.id, "")                                         #  id
    fill!(S.name, "")                                       # name
    fill!(S.src, "")                                        # src
    fill!(S.units, "")                                      # units
    fill!(S.fs, 0.0)                                        # fs
    fill!(S.gain, 1.0)                                      # gain
    for i = 1:n
      S.notes[i]  = Array{String,1}(undef,0)                # notes
      S.misc[i]   = Dict{String,Any}()                      # misc
      S.t[i]      = Array{Int64,2}(undef,0,2)               # t
      S.x[i]      = Array{Float32,1}(undef,0)               # x
      S.loc[i]    = GeoLoc()                                # loc
      S.resp[i]   = PZResp()                                # resp
    end
    return S
  end
  SeisData(n::Int) = n > 0 ? SeisData(UInt(n)) : SeisData()
end

# This intentionally undercounts exotic objects in :misc (e.g. a nested Dict)
# because those objects aren't written to disk or created by SeisIO
function sizeof(S::SeisData)
  s = sizeof(S.c) + 120
  for f in datafields
    V = getfield(S,f)
    s += sizeof(V)
    for i = 1:S.n
      v = getindex(V, i)
      s += sizeof(v)
      if f == :notes
        if !isempty(v)
          s += sum([sizeof(j) for j in v])
        end
      elseif f == :misc
        k = collect(keys(v))
        s += sizeof(k) + 64 + sum([sizeof(j) for j in k])
        for p in values(v)
          s += sizeof(p)
          if typeof(p) == Array{String,1}
            s += sum([sizeof(j) for j in p])
          end
        end
      end
    end
  end
  return s
end

function read(io::IO, ::Type{SeisData})
  Z = getfield(BUF, :buf)
  L = getfield(BUF, :int64_buf)

  # read begins ------------------------------------------------------
  N     = fastread(io, Int64)
  checkbuf_strict!(L, 2*N)
  fast_readbytes!(io, Z, 3*N)
  c1    = copy(Z[1:N])
  c2    = copy(Z[N+1:2*N])
  y     = code2typ.(getindex(Z, 2*N+1:3*N))
  cmp   = fastread(io)
  read!(io, L)
  nx    = getindex(L, N+1:2*N)

  if cmp == 0x01
    checkbuf_8!(Z, maximum(nx))
  end

  return SeisData(N,
    read_string_vec(io, Z),
    read_string_vec(io, Z),
    InstrumentPosition[read(io, code2loctyp(getindex(c1, i))) for i = 1:N],
    fastread(io, Float64, N),
    fastread(io, Float64, N),
    InstrumentResponse[read(io, code2resptyp(getindex(c2, i))) for i = 1:N],
    read_string_vec(io, Z),
    read_string_vec(io, Z),
    [read_misc(io, Z) for i = 1:N],
    [read_string_vec(io, Z) for i = 1:N],
    [read!(io, Array{Int64, 2}(undef, getindex(L, i), 2)) for i = 1:N],
    FloatArray[cmp == 0x01 ?
      (fast_readbytes!(io, Z, getindex(nx, i)); Blosc.decompress(getindex(y,i), Z)) :
      read!(io, Array{getindex(y,i), 1}(undef, getindex(nx, i)))
      for i = 1:N])
end

function write(io::IO, S::SeisData)
  N     = getfield(S, :n)
  LOC   = getfield(S, :loc)
  RESP  = getfield(S, :resp)
  T     = getfield(S, :t)
  X     = getfield(S, :x)
  MISC  = getfield(S, :misc)
  NOTES = getfield(S, :notes)

  cmp = false
  if KW.comp != 0x00
    nx_max = maximum([sizeof(getindex(X, i)) for i = 1:S.n])
    if (nx_max > KW.n_zip) || (KW.comp == 0x02)
      cmp = true
      Z = getfield(BUF, :buf)
      checkbuf_8!(Z, nx_max)
    end
  end

  codes = Array{UInt8,1}(undef, 3*N)
  L = Array{Int64,1}(undef, 2*N)

  # write begins ------------------------------------------------------
  write(io, N)
  p = fastpos(io)
  fastskip(io, 19*N+1)

  write_string_vec(io, S.id)                                          # id
  write_string_vec(io, S.name)                                        # name
  i = 0                                                               # loc
  while i < N
    i = i + 1
    loc = getindex(LOC, i)
    setindex!(codes, loctyp2code(loc), i)
    write(io, loc)
  end
  write(io, S.fs)                                                     # fs
  write(io, S.gain)                                                   # gain
  i = 0                                                               # resp
  while i < N
    i = i + 1
    resp = getindex(RESP, i)
    setindex!(codes, resptyp2code(resp), N+i)
    write(io, resp)
  end
  write_string_vec(io, S.units)                                       # units
  write_string_vec(io, S.src)                                         # src

  for i = 1:N; write_misc(io, getindex(MISC, i)); end                 # misc
  for i = 1:N; write_string_vec(io, getindex(NOTES, i)); end          # notes
  i = 0                                                               # t
  while i < N
    i = i + 1
    t = getindex(T, i)
    setindex!(L, size(t,1), i)
    write(io, t)
  end
  i = 0                                                               # x
  while i < N
    i = i + 1
    x = getindex(X, i)
    nx = lastindex(x)
    if cmp
      l = zero(Int64)
      while l == 0
        l = Blosc.compress!(Z, x, level=5)
        (l > zero(Int64)) && break
        nx_max = nextpow(2, nx_max)
        checkbuf_8!(Z, nx_max)
        @warn(string("Compression ratio > 1.0 for channel ", i, "; are data OK?"))
      end
      xc = view(Z, 1:l)
      write(io, xc)
      setindex!(L, l, N+i)
    else
      write(io, x)
      setindex!(L, nx, N+i)
    end
    setindex!(codes, typ2code(eltype(x)), 2*N+i)
  end
  q = fastpos(io)

  fastseek(io, p)
  write(io, codes)
  write(io, cmp)
  write(io, L)
  fastseek(io, q)
  # write ends ------------------------------------------------------
  return nothing
end


text:
function solve()
    ns = [parse(Int, x) for x in readline()]
    l = length(ns)
    m = -1
    for i=1:2^l-1
        s = 0
        for j=0:l-1
            if i & (1<<j) != 0
                s += ns[j+1]
            end
        end
        if s % 3 == 0
            m = max(m, count_ones(i))
        end
    end
    m == -1 ? -1 : l-m
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
export rule

@rule Dirichlet(:out, Marginalisation) (m_a::PointMass{ <: AbstractVector }, ) = Dirichlet(mean(m_a))

@rule Dirichlet(:out, Marginalisation) (q_a::PointMass{ <: AbstractVector }, ) = Dirichlet(mean(q_a))


text:
function test_disc_log_picard(P::GrpAbFinGen, mP::Hecke.MapPicardGrp, O::Hecke.AlgAssAbsOrd)
  for i = 1:5
    I = ideal(O, rand(O, 10))
    while !isinvertible(I)[1]
      I = ideal(O, rand(O, 10))
    end
    if !iszero(mP\I)
      return false
    end
  end

  if ngens(P) == 0
    return true
  end

  for i = 1:3
    c = rand(1:10, ngens(P))
    p = P(c)
    I1 = mP(p)
    I2 = mP(P[1])^c[1]
    for j = 2:ngens(P)
      I2 *= mP(P[j])^c[j]
    end
    if mP\I1 != mP\I2 || mP\I1 != p
      return false
    end
  end
  return true
end

@testset "Picard group of maximal orders of algebras" begin

  Qx, x = FlintQQ["x"]
  f = x^3 - 10x^2 - 3x - 2
  g = x^2 - 9x + 1
  A = AlgAss(f*g)
  O = maximal_order(A)
  P, mP = picard_group(O)
  @test issnf(P)
  @test P.snf == fmpz[ 2 ]
  @test test_disc_log_picard(P, mP, O)
  I = mP(P[1])
  @test_throws ErrorException Hecke.principal_generator(I)
  I2 = I^2
  a = Hecke.principal_generator(I2)
  @test I2 == ideal(O, a)

end

@testset "Picard group of non maximal orders of algebras" begin

  G = abelian_group([2, 3])
  A = AlgGrp(FlintQQ, G)
  O = Order(A, basis(A))
  P, mP = picard_group(O)
  @test issnf(P)
  @test ngens(P) == 0
  @test mP\ideal(O, one(O)) in P

  # To make sure it also works with AlgAss
  B, mB = AlgAss(A)
  O = Order(B, basis(B))
  P, mP = picard_group(O)
  @test issnf(P)
  @test ngens(P) == 0
  @test mP\ideal(O, one(O)) in P

  G = abelian_group([3, 3])
  A = AlgGrp(FlintQQ, G)
  O = Order(A, basis(A))
  P, mP = picard_group(O)
  @test issnf(P)
  @test P.snf == fmpz[ 3 ]
  @test test_disc_log_picard(P, mP, O)

  Qx, x = FlintQQ["x"]
  f = x^3 - 10x^2 - 3x - 2
  g = x^2 - 9x + 1
  A = AlgAss(f*g)
  O = Order(A, basis(A))
  P, mP = picard_group(O, true)
  @test issnf(P)
  @test P.snf == fmpz[ 2 ]
  @test test_disc_log_picard(P, mP, O)
  I = mP(P[1])
  @test_throws ErrorException principal_generator(I)
  I2 = I^2
  a = principal_generator(I2)
  @test I2 == ideal(O, a)
  I = ideal(O, 7*one(O)) # not coprime to the conductor of O in the maximal order
  a = principal_generator(I)
  @test I == ideal(O, a)
  # Test the refined discrete logarithm
  @test isdefined(mP, :betas)
  I = rand(O, -5:5)*mP(P[1])
  r, p = Hecke.refined_disc_log_picard_group(I, mP)
  @test evaluate(r)*mP(P[1]) == I
  @test p == P[1]
end

@testset "Unit groups of orders of algebras" begin

  Qx, x = FlintQQ["x"]
  f = x^3 - 12*x^2 - 6324*x + 459510
  K, a = number_field(f, "a", cached = false)
  OK = equation_order(K)
  UK, mUK = unit_group(OK)

  A = AlgAss(f)
  OA = Order(A, basis(A))
  UA, mUA = unit_group(OA)
  @test issnf(UA)
  @test UA.snf == fmpz[ 2, 0 ]
  G, GtoUK = sub(UK, [ mUK\OK(K(coefficients(elem_in_algebra(mUA(UA[i]), copy = false), copy = false))) for i = 1:ngens(UA) ])
  for i = 1:ngens(UK)
    @test haspreimage(GtoUK, UK[i])[1]
  end

end


text:
using Behavior
using Niancat.Languages

@given("that the word {String} is in the Swedish dictionary") do context, word
    context[:dictionary] = SwedishDictionary([word])
end

@when("a look-up of the word with spelling {String} is made") do context, word
    dictionary = context[:dictionary]
    context[:isindictionary]= word in dictionary
end

@then("it is found in the dictionary") do context
    isindictionary = context[:isindictionary]
    @expect isindictionary == true
end

@then("it is not found in the dictionary") do context
    isindictionary = context[:isindictionary]
    @expect isindictionary == false
end


text:
# SPDX-License-Identifier: MIT

# Joint Normal distribution with fixed identity covariance matrix.

struct JointNormalFixedCov{T1<:AbstractMatrix} <: AbstractNormal
    M::T1 # QxD mean
end

η(p::JointNormalFixedCov) = vec(p.M)
ξ(p::JointNormalFixedCov, η) = η

function unpack(p::JointNormalFixedCov, μ)
	Q, D = size(p.M)
    X = reshape(μ, Q, D)
	(X=X,)
end

A(p::JointNormalFixedCov, η) = (1/2) * dot(η, η)

function sample(p::JointNormalFixedCov)
	Q, D = size(p.M)
    p.M + randn(Q, D)
end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
using OffsetArrays

using ..Ellipticals: Weierstrass, term_expansion_normalized!

struct PoleTerm
    term ::OffsetVector{Dict{Variable, Coefficient}}
end

struct Cohesive
    solution ::Dict{RationalComplex, PoleTerm}
end

struct Pole
    taylor_rows::UnitRange{Int64}
    laurent_rows::IntOffsetVector
    coords::RationalComplex
    r_outer::Float64
end

function Pole(upper_term::Int64,
              first_row::Int64,
              coords::RationalComplex,
              r_outer::Float64)

    last_taylor_row = first_row + 2*(upper_term+1)-1
    taylor_rows = first_row : last_taylor_row

    laurent_rows = OffsetArray{Int64}(undef, -upper_term : -1)
    laur_row = last_taylor_row + 1
    for i in eachindex(laurent_rows)
        laurent_rows[i] = laur_row
        laur_row += 2
    end

    return Pole(taylor_rows, laurent_rows, coords, r_outer)
end

function last_row(pole::Pole)
    return maximum(pole.laurent_rows)+1
end

struct EllipticTerm
    re_column::Int64
    # im_column = re_column + 1
    derivative_n::Int64
end

struct SinglePoleElliptic
    terms::Vector{EllipticTerm}
    pole::Pole
    upper_power::Int64
end

function SinglePoleElliptic(upper_term::Int64,
                            first_column::Int64,
                            pole::Pole)

    terms = Vector{EllipticTerm}(undef, 0)

    re_column = first_column
    for i in 1 : upper_term
        push!(terms, EllipticTerm(re_column, i-2))
        re_column += 2
    end

    return SinglePoleElliptic(terms, pole, upper_term)
end

function last_column(elliptic::SinglePoleElliptic)
    return elliptic.terms[end].re_column + 1
end

struct GeneralElliptic
    sp_terms::Vector{SinglePoleElliptic}
end

function fill_slae_coupling_taylor!(slae::Array{Float64, 2},
                                    rows::UnitRange{Int64},
                                    term::EllipticTerm,
                                    delta::RationalComplex,
                                    series_dict::RComplex2OffsetMatrix,
                                    r_source::Float64,
                                    r_dest::Float64,
                                    temp_vector::ComplexOffsetVector)
    upper_term = div(size(rows)[1], 2) - 1

    term_expansion_normalized!(temp_vector, series_dict, delta,
                               term.derivative_n, upper_term, r_source, r_dest)

    first_row = first(rows)
    re_col = term.re_column
    im_col = re_col + 1

    for i in 0 : upper_term
        c, d = reim(temp_vector[i])

        slae[first_row + 2*i,   re_col] =  c
        slae[first_row + 2*i+1, re_col] =  d
    end
    for i in 0 : upper_term
        c, d = reim(temp_vector[i])

        slae[first_row + 2*i,   im_col] = -d
        slae[first_row + 2*i+1, im_col] =  c
    end
end

function fill_slae_coupling_laurent!(slae::Array{Float64, 2},
                                     rows::OffsetArray{Int64, 1, Array{Int64, 1}},
                                     term::EllipticTerm)
    row = rows[term.derivative_n]
    val = (-1.0)^term.derivative_n
    slae[row,   term.re_column]   = val
    slae[row+1, term.re_column]   = 0.0
    slae[row,   term.re_column+1] = 0.0
    slae[row+1, term.re_column+1] = val
end

function doubly_periodic_distance(z1::RationalComplex, z2::RationalComplex)
    delta = z1 - z2

    if real(delta) > 1//2
        delta -= 1
    end
    if real(delta) < -1//2
        delta += 1
    end
    if imag(delta) > 1//2
        delta -= 1im
    end
    if imag(delta) < -1//2
        delta += 1im
    end

    return delta
end

function rational_complex_distance(z1::RationalComplex, z2::RationalComplex,
                                   wei::Weierstrass)
        a, b = reim(z1-z2)
        return abs(a*wei.ω1 + b*wei.ω3)
end

function fill_slae_coupling!(slae::Array{Float64, 2},
                             elliptic::GeneralElliptic,
                             series_dict::RComplex2OffsetMatrix)

    max_n_terms = 0
    for el in elliptic.sp_terms
        # В системе уравнений на каждое комплексное равенство отводится
        # две строки для действительной и мнимой частей.
        n_terms = el.upper_power
        if n_terms > max_n_terms
            max_n_terms = n_terms
        end
    end
    temp_vector = OffsetVector{ComplexF64}(undef, 0 : max_n_terms-1)

    for spe in elliptic.sp_terms
        for term in spe.terms
            fill_slae_coupling_laurent!(slae, spe.pole.laurent_rows, term)
        end
    end

    for source in elliptic.sp_terms
        r_source = source.pole.r_outer
        for term in source.terms
            for dest in elliptic.sp_terms
                delta = doubly_periodic_distance(dest.pole.coords, source.pole.coords)
                r_dest= dest.pole.r_outer
                fill_slae_coupling_taylor!(slae, dest.pole.taylor_rows, term,
                                           delta, series_dict,
                                           r_source, r_dest, temp_vector)
            end
        end
    end
end

function push_to_wei_deltas!(deltas_dict::RComplex2IntDict, delta::RationalComplex,
                              new_value::Int64)
    if real(delta) < 0
        delta = -delta
    end

    old_value = get(deltas_dict, delta, 0)
    if new_value > old_value
        deltas_dict[delta] = new_value
    end
end

#####################################################################
##  TESTING SECTION
#####################################################################

# function test_construct_dict_of_wei_series()
#     ω1 = complex(1.0)
#     ω3 = exp(1im)
#     wei = Weierstrass(ω1, ω3)

#     test_dict = Dict{RationalComplex, Tuple{Int, Bool}}()
#     n = 1000
#     for i in 1 : n
#         a, b, c, d = rand(Int16, 4)
#         r = abs(a) < abs(b) ? a//b : b//a
#         i = abs(c) < abs(d) ? c//d : d//c
#         z = complex(r, i)

#         test_dict[z] = (15, rand(Bool))
#     end

#     return construct_dict_of_wei_series(test_dict, wei)
# end

# function test_term_expansion_normalized!()
#     ω1 = complex(1.0)
#     ω3 = exp(1im)
#     wei = Weierstrass(ω1, ω3)

#     test_dict = Dict{RationalComplex, Tuple{Int, Bool}}()
#     n = 1000
#     n_terms = 15
#     for i in 1 : n
#         a, b, c, d = rand(Int16, 4)
#         r = abs(a) < abs(b) ? a//b : b//a
#         i = abs(c) < abs(d) ? c//d : c//d
#         z = complex(r, i)
#         if r < 0
#             z = -z
#         end

#         test_dict[z] = (n_terms, rand(Bool))
#     end
#     test_dict[complex(0//1)] = (n_terms, rand(Bool))

#     wei_series_dict = construct_dict_of_wei_series(test_dict, wei)

#     output = OffsetVector{ComplexF64}(undef, 0:n_terms)
#     delta = complex(0//1)
#     r_s = 1.0
#     r_d = 1.0
#     @time begin
#         for n_d in -1:13
#             term_expansion_normalized!(output, wei_series_dict, delta, n_d, 12, r_s, r_d)
#         end
#     end
# end

# function test_fill_slae_coupling()

#     ω1 = complex(1.0)
#     ω3 = exp(1im)
#     wei = Weierstrass(ω1, ω3)

#     poles = Vector{Pole}(undef, 0)
#     elliptics = Vector{SinglePoleElliptic}(undef, 0)

#     upper_term = 15
#     first_row = 1
#     first_column = 1
#     central_coords = complex(1//2, 1//2)
#     r_central = 0.4
#     central_fiber = Pole(upper_term, first_row, central_coords, r_central)
#     elliptic = SinglePoleElliptic(upper_term, first_column, central_fiber)
#     push!(elliptics, elliptic)

#     upper_term = 6
#     nx = 20
#     ny = 20
#     r_subfiber = min(abs(ω1)/nx, abs(ω3)/ny) / 2.0

#     for i in 1 : nx
#         for j in 1 : ny
#             coords = complex(i//nx, j//ny)
#             dist = rational_complex_distance(coords, central_coords, wei)

#             if dist > r_central + r_subfiber
#                 first_row = last_row(elliptics[end].pole) + 1
#                 subfiber_pole = Pole(upper_term, first_row, coords, r_subfiber)

#                 first_column = last_column(elliptics[end]) + 1
#                 elliptic = SinglePoleElliptic(upper_term, first_column, subfiber_pole)

#                 push!(elliptics, elliptic)
#             end
#         end
#     end

#     general_elliptic = GeneralElliptic(elliptics)

#     wei_deltas = Dict{RationalComplex, Tuple{Int64, Bool}}()
#     n_poles = size(elliptics)[1]
#     for i in 1 : n_poles
#         for j in i : n_poles
#             delta = elliptics[i].pole.coords - elliptics[j].pole.coords
#             upper_power = max(elliptics[i].upper_power, elliptics[j].upper_power)
#             sigma = false
#             push_to_wei_deltas!(wei_deltas, delta, (upper_power, sigma))
#         end
#     end
# end

# function run_tests()
#     #test_theta()
#     #test_weierstrass()
#     #test_construct_dict_of_wei_series()
#     #test_term_expansion_normalized!()
#     test_fill_slae_coupling()

#     println("Tests passed")
# end

# function bench_weierstrass()
#     ω1 = complex(1.0)
#     ω3 = exp(1im)
#     wei = Weierstrass(ω1, ω3)

#     z_rand = rand(ComplexF64, 10000)
#     output = OffsetVector{ComplexF64}(undef, -2:10)

#     z = z_rand[1]
#     @show @allocated complute_normalized!(z, abs(z), wei, 10, output)

#     println("Benchmark: complute_normalized!()")
#     @time begin
#         for z in z_rand
#             complute_normalized!(z, abs(z), wei, 10, output)
#         end
#     end
# end

# function run_benchmarks()
#     bench_weierstrass()
# end

# run_tests()
# # run_benchmarks()