text:
"""
    gss2ss(sys; Eshape = "ident", atol = 0, atol1 = atol, atol2 = atol, rtol = nϵ) -> (sysr, r)

Convert the descriptor system `sys = (A-λE,B,C,D)` to an input-output equivalent descriptor system realization
`sysr = (Ar-λEr,Br,Cr,Dr)` without non-dynamic modes and having the same transfer function matrix.
The resulting `Er` is in the SVD-like form `Er = blockdiag(E1,0)`, with `E1` an `r × r` nonsingular matrix,
where `r` is the rank of `E`.

The keyword argument `Eshape` specifies the shape of `E1` as follows:

if `Eshape = "ident"` (the default option), `E1` is an identity matrix of order `r` and if `E` is nonsingular,
then the resulting system `sysr` is a standard state-space system;

if `Eshape = "diag"`, `E1` is a diagonal matrix of order `r`, where the diagonal elements are the
decreasingly ordered nonzero singular values of `E`;

if `Eshape = "triu"`, `E1` is an upper triangular nonsingular matrix of order `r`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.

If `Eshape = "triu"`, the reductions of `E` and `A` are performed using rank decisions based on rank revealing
QR-decompositions with column pivoting.  If `Eshape = "ident"` or `Eshape = "diag"` the reductions are performed
using the more reliable SVD-decompositions.
"""
function gss2ss(sys::DescriptorStateSpace{T}; Eshape = "ident", atol::Real = zero(real(T)), atol1::Real = atol,
                atol2::Real = atol, rtol::Real =  sys.nx*eps(real(float(one(real(T)))))*iszero(max(atol1,atol2))) where T

    # finish for a standard state space system or pure gain
    n = sys.nx
    (n == 0 || sys.E == I) && (return sys, n)

    T1 = T <: BlasFloat ? T : promote_type(Float64,T)

    Ar, Er, Br, Cr, Dr = dssdata(T1,sys)

    # exploit the upper triangular form of a nonsingular E
    epsm = eps(real(T1))
    if  istriu(Er) && rcond(UpperTriangular(Er),atol2) > n*epsm
        Eshape == "triu" && (return sys, n)
        if Eshape == "ident"
            eltype(Er) <: Complex &&
               (return dss(ldiv!(UpperTriangular(Er),Ar), ldiv!(UpperTriangular(Er),Br), Cr, Dr, Ts = sys.Ts), n)
            # make diagonal elements of E positive
            indneg = (diag(Er) .< 0)
            if any(indneg)
                Ar[indneg,:] = -Ar[indneg,:]
                Er[indneg,:] = -Er[indneg,:]
                Br[indneg,:] = -Br[indneg,:]
            end
            e2 = UpperTriangular(sqrt(Er))
            ldiv!(e2,Ar)
            rdiv!(Ar,e2)
            ldiv!(e2,Br)
            rdiv!(Cr,e2)
            return dss(Ar, Br, Cr, Dr, Ts = sys.Ts), n
        end
    end

    # Using orthogonal/unitary transformation matrices Q and Z, reduce the
    # matrices A, E, B and C to the forms
    #
    #               [ At11  At12 At13 ]                  [ Et11  0  0 ]
    # Ar = Q'*A*Z = [ At21  At22  0   ] ,  Er = Q'*E*Z = [  0    0  0 ] ,
    #               [ At31   0    0   ]                  [  0    0  0 ]
    #
    #             [ Bt1 ]
    # Br = Q'*B = [ Bt2 ] ,  Cr = C*Z = [ Ct1  Ct2  Ct3 ]
    #             [ Bt3 ]
    #
    fast = (Eshape == "triu")
    rE, rA22  = _svdlikeAE!(Ar, Er, nothing, nothing, Br, Cr,
                     fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false)
    # employ state residualzation formulas to remove non-dynamic modes
    if rA22 > 0
        i1 = 1:rE
        i2 = rE+1:rE+rA22
        # make At22 = I
        fast ? (A22 = UpperTriangular(Ar[i2,i2])) : (A22 = Diagonal(Ar[i2,i2]))
        ldiv!(A22,view(Ar,i2,i1))
        ldiv!(A22,view(Br,i2,:))
        # apply simplified residualization formulas
        Dr -= Cr[:,i2]*Br[i2,:]
        Br[i1,:] -= Ar[i1,i2]*Br[i2,:]
        Cr[:,i1] -= Cr[:,i2]*Ar[i2,i1]
        Ar[i1,i1] -= Ar[i1,i2]*Ar[i2,i1]
        ir = [i1; rE+rA22+1:n]
    else
        i1 = 1:rE
        ir = 1:n
    end
    # bring E to the required shape
    if fast || Eshape == "diag"
       # Eshape == "triu" or Eshape == "diag": we are already done
       return dss(view(Ar,ir,ir),view(Er,ir,ir),view(Br,ir,:),view(Cr,:,ir),Dr, Ts = sys.Ts), rE
    elseif Eshape == "ident"
        tid = Diagonal(1 ./sqrt.(diag(Er[i1,i1])))
        Er[i1,i1] = eye(T1,rE)
        lmul!(tid,view(Ar,i1,ir))
        rmul!(view(Ar,ir,i1),tid)
        lmul!(tid,view(Br,i1,:))
        rmul!(view(Cr,:,i1),tid)
        return rE+rA22 == n ?  (dss(view(Ar,ir,ir),view(Br,ir,:),view(Cr,:,ir),Dr, Ts = sys.Ts), rE) :
                               (dss(view(Ar,ir,ir),view(Er,ir,ir),view(Br,ir,:),view(Cr,:,ir),Dr, Ts = sys.Ts), rE)
    else
        error("improper shape option for E")
    end

    # end GSS2SS
end
"""

     dss2ss(sys, x0 = 0; state-mapping = false, simple_infeigs = true, fast = true, atol1, atol2, rtol)
               -> (sysr, xr0, Mx, Mu)

Return for a proper descriptor system `sys = (A-λE,B,C,D)` and initial state `x0`,
the equivalent reduced order standard system `sysr = (Ar-λI,Br,Cr,Dr)` and
the corresponding reduced consistent initial state `xr0`.

If `state_mapping = true`, the state mapping matrices `Mx` and `Mu` are also determined such that
the values `x(t)` and `xr(t)` of the state vectors of the systems `sys` and `sysr`, respectively,
and the input vector `u(t)` are related as `x(t) = Mx*xr(t)+Mu*u(t)`.
In this case, higher order uncontrollable infinite eigenvalues can be eliminated if `simple_infeigs = false`.

By default, `state_mapping = false` and `Mx = nothing` and `Mu = nothing`.
In this case, higher order uncontrollable or unobservable infinite eigenvalues
can be eliminated if `simple_infeigs = false`.

By default, `simple_infeigs = true`, and simple infinite eigenvalues for the pair `(A,E)` are assumed and eliminated.

The underlying pencil manipulation algorithms employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true` (default), or the SVD-decomposition,
if `fast = false`. The rank decision based on the SVD-decomposition is generally more reliable,
but the involved computational effort is higher.

The keyword arguments `atol1`, `atol2` and `rtol` specify, respectively,
the absolute tolerance for the nonzero elements of `A`, `B`, `C`, `D`, the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`.
The default relative tolerance is `n*ϵ`, where `n` is the order of the square matrices `A` and `E`, and  `ϵ` is the working machine epsilon.
"""
function dss2ss(sys::DescriptorStateSpace{T}, x0::Vector = zeros(T,sys.nx); state_mapping::Bool = false, simple_infeigs::Bool = true,
                   fast::Bool = true, atol::Real = zero(float(real(T))), atol1::Real = atol, atol2::Real = atol,
                   rtol::Real = sys.nx*eps(real(float(one(T))))*iszero(min(atol1,atol2))) where T
   T1 = T <: BlasFloat ? T : promote_type(Float64,T)
   n = sys.nx
   m = sys.nu
   nullx0 = iszero(x0)
   sys.E == I && (return sys, copy_oftype(x0,T1), state_mapping ? I : nothing, state_mapping ? zeros(T1,n,m) : nothing)

   # eliminate uncontrollable/unobservable infinite and nonzero eigenvalues
   simple_infeigs || state_mapping ||
      ((sys, L, Z) = gir_lrtran(sys; fast, rtran = !nullx0, finite = false, noseig = false, atol1, atol2, rtol) )

   A, E, B, C, D = dssdata(T1,sys)
   n1, m = size(B)
   if n1 < n
      state_mapping && (@warn "state mapping disabled being not feasible: try simple_infeigs = true"; state_mapping = false)
      # adjust initial condition
      nullx0 ? x0 = zeros(T1,n1) : x0 = Z[:,1:n1]'*x0
   end

   ONE = one(T1)
   state_mapping || (Mx = nothing; Mu = nothing)
   if simple_infeigs
      Z = Matrix{T1}(I,n1,n1)
      n, rA22  = _svdlikeAE!(A, E, nothing, Z, B, C; fast, atol1, atol2, rtol, withQ = false)
      n+rA22 == n1 || error("The system is possibly improper: try with simple_infeigs = false")
      if rA22 == 0
         # get rid of E matrix, exploit upper triangular or diagonal shape
         fast ? (F = UpperTriangular(E)) : (F = Diagonal(E))
         ldiv!(F,A); ldiv!(F,B)
         state_mapping && (Mx = Z; Mu = zeros(T1,n,m))
         return dss(A, B, C, D, Ts = sys.Ts), Z'*x0, Mx, Mu
      else
         i1 = 1:n
         i2 = n+1:n1
         # adjust initial condition
         Z1 = view(Z,:,i1)
         Z2 = view(Z,:,i2)
         A11 = view(A,i1,i1)
         E11 = view(E,i1,i1)
         A12 = view(A,i1,i2)
         A21 = view(A,i2,i1)
         B1 = view(B,i1,:)
         B2 = view(B,i2,:)
         C1 = view(C,:,i1)
         C2 = view(C,:,i2)
         xt = Z1'*x0
         # make A22 = I, exploit upper triangular or diagonal shape
         fast ? (A22 = UpperTriangular(A[i2,i2])) : (A22 = Diagonal(A[i2,i2]))
         ldiv!(A22,A21)  # inv(A22)*A21 -> A21
         ldiv!(A22,B2)   # inv(A22)*B2 -> B2
         # apply simplified residualization formulas
         mul!(D, C2, B2, -ONE, ONE)     # D -> D - C2*B2
         mul!(B1, A12, B2, -ONE, ONE)   # B1 -> B1 - A12*B2
         mul!(C1, C2, A21, -ONE, ONE)   # C1 -> C1 - C2*A21
         mul!(A11, A12, A21, -ONE, ONE) # A11 -> A11 - A12*A21
         # state_mapping && rank(B2; atol = atol1, rtol) < rA22 &&
         #    (@warn "state mapping disabled being not feasible"; state_mapping = false; Mx = nothing; Mu = nothing)
         state_mapping && (mul!(Z1, Z2, A21, ONE, ONE); Mx = copy(Z1); Mu = -Z2*B2)
         # quick exit in constant case
         n == 0 && (return dss(A11, B1, C1, D, Ts = sys.Ts), xt, Mx, Mu)
         # get rid of E matrix, exploit upper triangular or diagonal shape
         fast ? (F = UpperTriangular(E11)) : (F = Diagonal(E11))
         ldiv!(F,A11); ldiv!(F,B1)
         return dss(A11, B1, C1, D, Ts = sys.Ts), xt, Mx, Mu
      end
   else
      # separate infinite eigenvalues in the trailing part; leading E11 is upper triangular
      Z = Matrix{T1}(I,n1,n1)
      _, blkdims = MatrixPencils.fisplit!(A, E, nothing, Z, B, C; fast, finite_infinite = true, atol1, atol2, rtol, withQ = false)
      nf, ni = blkdims
      if ni == 0
         # get rid of E matrix, exploit upper triangular shape
         F = UpperTriangular(E)
         ldiv!(F,A); ldiv!(F,B)
         state_mapping && (Mx = Z; Mu = zeros(T1,n,m))
         return dss(A, B, C, D, Ts = sys.Ts), Z'*x0, Mx, Mu
      end
      i1 = 1:nf
      # separate uncontrollable infinite eigenvalues
      i2 = nf+1:nf+ni
      # save trailing matrices of A, E, B and C
      Ai = A[i2,i2]
      Ei = E[i2,i2]
      Bi = B[i2,:]
      Ci = C[:,i2]
      Qi = nothing
      Zi = Matrix{T}(I,ni,ni)
      Qi, Zi, _, _, niuc = sklf_rightfin!(view(E,i2,i2), view(A,i2,i2), view(B,i2,:), view(C,:,i2);
                                          fast, atol1, atol2, rtol, withQ = false, withZ = true)
      if niuc == 0
         # restore original matrices if all infinite eigenvalues are controllable
         A[i2,i2] = copy(Ai)
         E[i2,i2] = copy(Ei)
         B[i2,:] = copy(Bi)
         C[:,i2] = copy(Ci)
      else
         # apply the left transformation Zi to Z2, A12 and E12
         Z[:,i2] = view(Z,:,i2)*Zi
         A[i1,i2] = view(A,i1,i2)*Zi
         E[i1,i2] = view(E,i1,i2)*Zi
      end

      # refine structure information
      nic = ni-niuc
      i2 = nf+1:nf+nic
      if norm(view(E,i2,i2),Inf) <= atol2 + rtol*norm(E,Inf)
         # only simple controllable infinite eigenvalues
         Z1 = view(Z,:,i1)
         Z2 = view(Z,:,i2)
         A11 = view(A,i1,i1)
         A12 = view(A,i1,i2)
         B1 = view(B,i1,:)
         B2 = view(B,i2,:)
         C1 = view(C,:,i1)
         C2 = view(C,:,i2)
         E11 = view(E,i1,i1)
         E12 = view(E,i1,i2)
         # xt = Z1'*x0
         # make A22 = I
         A22 = UpperTriangular(A[i2,i2])
         ldiv!(A22,B2)   # inv(A22)*B2 -> B2
         # make E11 = I
         F = UpperTriangular(E11)
         # inv(E1)[A11 A12] -> [A11 A12]; inv(E1)*E12 -> E12; inv(E1)*B1 -> B1
         ldiv!(F,view(A,i1,:)); ldiv!(F,E12); ldiv!(F,B1)
         # make leading part of E block diagonal (i.e., diag(I,0))
         mul!(C2, C1, E12, -ONE, ONE)   # C2 -> C2 - C1*E12
         mul!(A12, A11, E12, -ONE, ONE)   # A12 -> A12 - A11*E12
         # apply simplified residualization formulas
         mul!(D, C2, B2, -ONE, ONE)     # D -> D - C2*B2
         mul!(B1, A12, B2, -ONE, ONE)   # B1 -> B1 - A12*B2
         xt = Z1'*x0 + E12*(Z2'*x0)     # determine consistent reduced initial state
         mul!(Z2, Z1, E12, -ONE, ONE)   # Z2 -> Z2 - Z1*E12
         state_mapping && (Mx = copy(Z1); Mu = -Z2*B2)
         return dss(A11, B1, C1, D, Ts = sys.Ts), xt, Mx, Mu
      else
         error("controllable higher order infinite eigenvalues present")
      end
   end
   # end DSS2SS
end
"""
    sysr = gir(sys; finite = true, infinite = true, contr = true, obs = true, noseig = false,
               fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol = nϵ)

Compute for a descriptor system `sys = (A-λE,B,C,D)` of order `n` a reduced order descriptor system
`sysr = (Ar-λEr,Br,Cr,Dr)` of order `nr ≤ n` such that `sys` and `sysr` have the same transfer function matrix, i.e.,

             -1                    -1
     C*(λE-A)  *B + D = Cr*(λEr-Ar)  *Br + Dr .

The least possible order `nr` is achieved if `finite = true`, `infinite = true`, `contr = true`,
`obs = true` and `nseig = true`. Such a realization is called `minimal` and satisfies:

     (1) rank[Br Ar-λEr] = nr for all finite λ (finite controllability)

     (2) rank[Br Er] = nr (infinite controllability)

     (3) rank[Ar-λEr; Cr] = nr for all finite λ (finite observability)

     (4) rank[Er; Cr] = nr (infinite observability)

     (5) Ar-λEr has no simple infinite eigenvalues

A realization satisfying only conditions (1)-(4) is called `irreducible` and is computed by default.

Some reduction steps can be skipped by appropriately selecting the keyword arguments
`contr`, `obs`, `finite`, `infinite` and `nseig`.

If `contr = false`, then the controllability conditions (1) and (2) are not enforced.
If `contr = true` and `finite = true`, then the finite controllability condition (1) is enforced.
If `contr = true` and `finite = false`, then the finite controllability condition (1) is not enforced.
If `contr = true` and `infinite = true`, then the infinite controllability condition (2) is enforced.
If `contr = true` and `infinite = false`, then the infinite controllability condition (2) is not enforced.

If `obs = false`, then observability condition (3) and (4) are not enforced.
If `obs = true` and `finite = true`, then the finite observability condition (3) is enforced.
If `obs = true` and `finite = false`, then the finite observability condition (3) is not enforced.
If `obs = true` and `infinite = true`, then the infinite observability condition (4) is enforced.
If `obs = true` and `infinite = false`, then the infinite observability condition (4) is not enforced.

If `nseig = true`, then condition (5) on the lack of simple infinite eigenvalues is also enforced.

To enforce conditions (1)-(4), the `Procedure GIR` in `[1, page 328]` is employed, which performs
orthogonal similarity transformations on the matrices of the original realization `(A-λE,B,C,D)`
to obtain an irreducible realization using structured pencil reduction algorithms.
To enforce condition (5), residualization formulas (see, e.g., `[1, page 329]`) are employed which
involves matrix inversions.

The underlying pencil manipulation algorithms employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of matrices `A`, `B`, `C`, `D`, the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`.
The default relative tolerance is `nϵ`, where `ϵ` is the working _machine epsilon_
and `n` is the order of the system `sys`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.

[1] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function gir(SYS::DescriptorStateSpace{T}; atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
             rtol::Real =  SYS.nx*eps(real(float(one(real(T)))))*iszero(max(atol1,atol2)),
             fast::Bool = true, finite::Bool = true, infinite::Bool = true,
             contr::Bool = true, obs::Bool = true, noseig::Bool = false) where T
    if SYS.E == I
        A, B, C = lsminreal(SYS.A, SYS.B, SYS.C; fast = fast, atol = atol1, rtol = rtol, contr = contr, obs = obs)
        T1 = eltype(A)
        return DescriptorStateSpace{T1}(A, I, B, C, copy_oftype(SYS.D,T1), SYS.Ts)
    else
        A, E, B, C, D = lsminreal2(SYS.A, SYS.E, SYS.B, SYS.C, SYS.D;
            fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, finite = finite, infinite = infinite,
            contr = contr, obs = obs, noseig = noseig)
        return DescriptorStateSpace{eltype(A)}(A, E, B, C, D, SYS.Ts)
    end
end
"""
    gir_lrtran(sys; ltran = false, rtran = false, finite = true, infinite = true, contr = true, obs = true,
             noseig = false, fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol = nϵ) -> (sysr, Q, Z)

This is a special version of the function `gir` to additionally determine the left and right
transformation matrices `Q = [Q1 Q2]` and `Z = [Z1 Z2]`, respectively, such that the matrices `Ar`, `Er`, `Br` and `Cr` of the resulting
descriptor system `sysr = (Ar-λEr,Br,Cr,Dr)` are given by `Ar = Q1'*A*Z1`, `Er = Q1'*E*Z1`, `Br = Q1'*B`, `Cr = C*Z1`,
where the number of columns of `Q1` and `Z1` is equal to the order of matrix `Ar`. `Q` and `Z` result orthogonal if `noseig = false`.
`Q = nothing` if `ltran = false` and `Z = nothing` if `rtran = false`.
See [`gir`](@ref) for details on the rest of keyword parameters.
"""
function gir_lrtran(SYS::DescriptorStateSpace{T}; ltran::Bool = false, rtran::Bool = false, atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
                    rtol::Real =  SYS.nx*eps(real(float(one(real(T)))))*iszero(max(atol1,atol2)),
                    fast::Bool = true, finite::Bool = true, infinite::Bool = true,
                    contr::Bool = true, obs::Bool = true, noseig::Bool = false) where T
    T1 = T <: BlasFloat ? T : promote_type(Float64,T)
    A, E, B, C, D = dssdata(T1,SYS)

    E == I && (E = eye(T,size(A,1)))
    A, E, B, C, D, Q, Z, = lsminreal2_lrtran(A, E, B, C, D; withQ = ltran, withZ = rtran,
                                             fast, atol1, atol2, rtol, finite, infinite, contr, obs, noseig)
    return dss(A, E, B, C, D, Ts = SYS.Ts), Q, Z
end
"""
    sysr = gminreal(sys; contr = true, obs = true, noseig = true, fast = true,
                    atol = 0, atol1 = atol, atol2 = atol, rtol = nϵ)

Compute for a descriptor system `sys = (A-λE,B,C,D)` of order `n` a reduced order descriptor system
`sysr = (Ar-λEr,Br,Cr,Dr)` of order `nr ≤ n` such that `sys` and `sysr` have the same transfer function matrix, i.e.,

             -1                    -1
     C*(λE-A)  *B + D = Cr*(λEr-Ar)  *Br + Dr .

The least possible order `nr` is achieved if `contr = true`, `obs = true` and `nseig = true`.
Such a realization is called `minimal` and satisfies:

     (1) rank[Br Ar-λEr] = nr for all finite λ (finite controllability)

     (2) rank[Br Er] = nr (infinite controllability)

     (3) rank[Ar-λEr; Cr] = nr for all finite λ (finite observability)

     (4) rank[Er; Cr] = nr (infinite observability)

     (5) Ar-λEr has no simple infinite eigenvalues

A realization satisfying only conditions (1)-(4) is called `irreducible`.

Some reduction steps can be skipped by appropriately selecting the keyword arguments
`contr`, `obs` and `nseig`.

If `contr = false`, then the controllability conditions (1) and (2) are not enforced.

If `obs = false`, then observability condition (3) and (4) are not enforced.

If `nseig = false`, then condition (5) on the lack of simple infinite eigenvalues is not enforced.

To enforce conditions (1)-(4), orthogonal similarity transformations are performed on
the matrices of the original realization `(A-λE,B,C,D)` to obtain an irreducible realization using
structured pencil reduction algorithms, as the fast versions of the reduction techniques of the
full row rank pencil [B A-λE] and full column rank pencil [A-λE;C] proposed in [1].
To enforce condition (5), residualization formulas (see, e.g., `[2, page 329]`) are employed which
involves matrix inversions.

The underlying pencil manipulation algorithms employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of matrices `A`, `B`, `C`, `D`, the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`.
The default relative tolerance is `nϵ`, where `ϵ` is the working _machine epsilon_
and `n` is the order of the system `sys`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.

[1] P. Van Dooreen, The generalized eigenstructure problem in linear system theory,
IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.

[2] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function gminreal(SYS::DescriptorStateSpace{T}; atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
    rtol::Real =  SYS.nx*eps(real(float(one(real(T)))))*iszero(max(atol1,atol2)),
    fast::Bool = true, contr::Bool = true, obs::Bool = true, noseig::Bool = true) where T
    if SYS.E == I
        A, B, C = lsminreal(SYS.A, SYS.B, SYS.C; fast = fast, atol = atol1, rtol = rtol, contr = contr, obs = obs)
        T1 = eltype(A)
        return DescriptorStateSpace{T1}(A, I, B, C, copy_oftype(SYS.D,T1), SYS.Ts)
    else
        A, E, B, C, D = lsminreal(SYS.A, SYS.E, SYS.B, SYS.C, SYS.D;
            fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol,
            contr = contr, obs = obs, noseig = noseig)
        return DescriptorStateSpace{eltype(A)}(A, E, B, C, D, SYS.Ts)
    end
end
"""
    gbalmr(sys, balance = false, matchdc = false, ord = missing, atolhsv = 0, rtolhsv = nϵ,
           atolmin = atolhsv, rtolmin = rtolhsv,
           atol = 0, atol1 = atol, atol2 = atol, rtol, fast = true) -> (sysr, hs)

Compute for a proper and stable descriptor system `sys = (A-λE,B,C,D)` with the transfer function
matrix `G(λ)`, a reduced order realization `sysr = (Ar-λEr,Br,Cr,Dr)` and the vector `hs` of decreasingly
ordered Hankel singular values of the system `sys`. If `balance = true`, a balancing-based approach
is used to determine a reduced order minimal realization
of the form `sysr = (Ar-λI,Br,Cr,Dr)`. For a continuous-time system `sys`, the resulting realization `sysr`
is balanced, i.e., the controllability and observability grammians are equal and diagonal.
If additonally `matchdc = true`, the resulting `sysr` is computed using state rezidualization formulas
(also known as _singular perturbation approximation_) which additionally preserves the DC-gain of `sys`.
In this case, the resulting realization `sysr` is balanced (for both continuous- and discrete-time systems).
If `balance = false`, an enhanced accuracy balancing-free approach is used to determine the
reduced order system `sysr`.

If `ord = nr`, the resulting order of `sysr` is `min(nr,nrmin)`, where `nrmin` is the order of a minimal
realization of `sys` determined as the number of Hankel singular values exceeding `max(atolmin,rtolmin*HN)`, with
`HN`, the Hankel norm of `G(λ)`. If `ord = missing`, the resulting order is chosen as the number of Hankel
singular values exceeding `max(atolhsv,rtolhsv*HN)`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of `A`, `B`, `C`, `D`,
the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`.
The default relative tolerance is `nϵ`, where `ϵ` is the working _machine epsilon_
and `n` is the order of the system `sys`. The keyword argument `atol` can be used
to simultaneously set `atol1 = atol` and `atol2 = atol`.

If `E` is singular, the uncontrollable infinite eigenvalues of the pair `(A,E)` and
the non-dynamic modes are elliminated using minimal realization techniques.
The rank determinations in the performed reductions
are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

Method:  For the order reduction of a standard system, the balancing-free method of [1] or
the balancing-based method of [2] are used. For a descriptor system the balancing related order reduction
methods of [3] are used. To preserve the DC-gain of the original system, the singular perturbation
approximation method of [4] is used in conjunction with the balancing-based or balancing-free
approach of [5].

References

[1] A. Varga.
    Efficient minimal realization procedure based on balancing.
    In A. El Moudni, P. Borne, and S.G. Tzafestas (Eds.),
    Prepr. of the IMACS Symp. on Modelling and Control of Technological
    Systems, Lille, France, vol. 2, pp.42-47, 1991.

[2] M. S. Tombs and I. Postlethwaite.
    Truncated balanced realization of a stable non-minimal state-space
    system. Int. J. Control, vol. 46, pp. 1319–1330, 1987.

[3] T. Stykel.
    Gramian based model reduction for descriptor systems.
    Mathematics of Control, Signals, and Systems, 16:297–319, 2004.

[4] Y. Liu Y. and B.D.O. Anderson
    Singular Perturbation Approximation of Balanced Systems,
    Int. J. Control, Vol. 50, pp. 1379-1405, 1989.

[5] Varga A.
    Balancing-free square-root algorithm for computing singular perturbation approximations.
    Proc. 30-th IEEE CDC,  Brighton, Dec. 11-13, 1991, Vol. 2, pp. 1062-1065.
"""
function gbalmr(sys::DescriptorStateSpace{T}; balance::Bool = false, matchdc::Bool = false, fast::Bool = true,
    ord::Union{Int,Missing} = missing, atolhsv::Real = zero(real(T)),
    rtolhsv::Real = sqrt(eps(real(float(one(T)))))*iszero(atolhsv), atolmin::Real = atolhsv, rtolmin::Real = rtolhsv,
    offset::Real = sqrt(eps(float(real(T)))), atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
    rtol::Real = (size(sys.A,1)*eps(real(float(one(T)))))*iszero(min(atol1,atol2))) where T

    T1 = T <: BlasFloat ? T : promote_type(Float64,T)
    s2eps = offset
    disc = !iszero(sys.Ts)
    standsys = sys.E == I
    ONE = one(T1)

    if  standsys
        n = order(sys)
        # for a non-dynamic system, we set the Hankel norm to zero,
        # but the Hankel singular values are empty
        n == 0 && (return sys, zeros(real(T1),0))
        # reduce the system to Schur coordinate form
        SF = schur(sys.A)
        # check stability
        ((disc && maximum(abs.(SF.values)) >= 1-s2eps) || (!disc && maximum(real(SF.values)) >= -s2eps)) &&
              error("The system sys is unstable")
        bs = SF.Z'*sys.B
        cs = sys.C*SF.Z
        S = plyaps(SF.T, bs; disc = disc)
        R = plyaps(SF.T', cs'; disc = disc)
        SV = svd!(R*S); hs = SV.S
    else
        # eliminate uncontrollable infinite eigenvalues and non-dynamic modes if possible
        if LinearAlgebra.LAPACK.gecon!('1', LinearAlgebra.LAPACK.getrf!(copy(sys.E))[1],opnorm(sys.E,1)) < s2eps
           sys = gir(sys, finite = false, noseig = true, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
           # the resulting sys may still be improper
           LinearAlgebra.LAPACK.gecon!('1', LinearAlgebra.LAPACK.getrf!(copy(sys.E))[1],opnorm(sys.E,1)) < s2eps &&
               error("The system sys is not proper")
        end
        n = order(sys)
        # for a non-dynamic system, we set the Hankel norm to zero,
        # but the Hankel singular values are empty
        n == 0 && (return sys, zeros(real(T1),0))
        # reduce the system to generalized Schur coordinate form
        SF = schur(sys.A,sys.E)
        ((disc && maximum(abs.(SF.values)) >= 1-s2eps) || (!disc && maximum(real(SF.values)) >= -s2eps)) &&
              error("The system sys is unstable")
        bs = SF.Q'*sys.B
        cs = sys.C*SF.Z
        S = plyaps(SF.S, SF.T, bs; disc = disc)
        R = plyaps(SF.S', SF.T', cs'; disc = disc)
        SV = svd!(R*UpperTriangular(SF.T)*S); hs = SV.S
    end
    # determine the order nrmin of a minimal realization
    atolhsv >= 0 || error("The threshold atolhsv must be non-negative")
    rtolhsv >= 0 || error("The threshold rtolhsv must be non-negative")
    atolmin >= 0 || error("The threshold atolmin must be non-negative")
    rtolmin >= 0 || error("The threshold rtolmin must be non-negative")
    nrmin = count(hs .> max(atolmin,rtolmin*hs[1]))
    # determine the order nr of the reduced model
    if ismissing(ord)
       nr = min(nrmin, count(hs .> max(atolhsv,rtolhsv*hs[1])))
    else
       ord >= 0 || error("Desired order ord must be non-negative, got instead $ord")
       nr = min(nrmin,ord)
    end
    nr == 0 && (return dss(sys.D,Ts = sys.Ts), hs)
    i1 = 1:nr
    if balance
       # apply balancing-based formulas: resulting Er = I
       if matchdc && nr < nrmin
          indmin = 1:nrmin
          hsi2 = Diagonal(1 ./sqrt.(view(hs,indmin)))
          Tl = lmul!(R',view(SV.U,:,indmin))*hsi2
          Tr = lmul!(S,SV.V[:,indmin])*hsi2
          standsys ? amin = Tl'*SF.T*Tr : amin = Tl'*SF.S*Tr
          bmin = Tl'*bs
          cmin = cs*Tr
          i2 = nr+1:nrmin
          Ar = view(amin,i1,i1)
          Br = view(bmin,i1,:)
          Cr = view(cmin,:,i1)
          Dr = copy_oftype(sys.D,T1)
          disc && (amin[i2,i2] -= I)
          LUF = lu!(view(amin,i2,i2))
          ldiv!(LUF,view(amin,i2,i1))
          ldiv!(LUF,view(bmin,i2,:))
          # apply state residualization formulas
          mul!(Dr,view(cmin,:,i2),view(bmin,i2,:),-ONE, ONE)
          mul!(Br,view(amin,i1,i2),view(bmin,i2,:),-ONE, ONE)
          mul!(Cr,view(cmin,:,i2),view(amin,i2,i1),-ONE, ONE)
          mul!(Ar,view(amin,i1,i2),view(amin,i2,i1),-ONE, ONE)
          # return the minimal balanced system
          return dss(Ar, I, Br, Cr, Dr, Ts = sys.Ts), hs
       else
          hsi2 = Diagonal(1 ./sqrt.(view(hs,i1)))
          Tl = lmul!(R',view(SV.U,:,i1))*hsi2
          Tr = lmul!(S,SV.V[:,i1])*hsi2
          # return the minimal balanced system
          return dss(standsys ? Tl'*SF.T*Tr : Tl'*SF.S*Tr, I, Tl'*bs, cs*Tr, sys.D, Ts = sys.Ts), hs
       end
    else
        # apply balancing-free formulas
        if nr < n
            if matchdc
                i2 = nr+1:nrmin
                Tl = [Matrix(qr!(lmul!(R',view(SV.U,:,i1))).Q) Matrix(qr!(lmul!(R',view(SV.U,:,i2))).Q)]
                Tr = [Matrix(qr!(lmul!(S,SV.V[:,i1])).Q) Matrix(qr!(lmul!(S,SV.V[:,i2])).Q)]
                #standsys ? (amin = Tl'*SF.T*Tr; emin = Tl'*Tr) : (amin = Tl'*SF.S*Tr; emin = Tl'*SF.T*Tr)
                standsys ? amin = Tl'*SF.T*Tr : amin = Tl'*SF.S*Tr
                bmin = Tl'*bs
                cmin = cs*Tr
                Ar = view(amin,i1,i1)
                Er = standsys ? view(Tl,:,i1)'*view(Tr,:,i1) : view(Tl,:,i1)'*UpperTriangular(SF.T)*view(Tr,:,i1)
                Br = view(bmin,i1,:)
                Cr = view(cmin,:,i1)
                Dr = copy_oftype(sys.D,T1)
                disc && (amin[i2,i2] -= standsys ? view(Tl,:,i2)'*view(Tr,:,i2) : view(Tl,:,i2)'*UpperTriangular(SF.T)*view(Tr,:,i2))
                LUF = lu!(view(amin,i2,i2))
                ldiv!(LUF,view(amin,i2,i1))
                ldiv!(LUF,view(bmin,i2,:))
                # apply state residualization formulas
                mul!(Dr,view(cmin,:,i2),view(bmin,i2,:),-ONE, ONE)
                mul!(Br,view(amin,i1,i2),view(bmin,i2,:),-ONE, ONE)
                mul!(Cr,view(cmin,:,i2),view(amin,i2,i1),-ONE, ONE)
                mul!(Ar,view(amin,i1,i2),view(amin,i2,i1),-ONE, ONE)
                standsys || (return dss(Ar, Er, Br, Cr, Dr, Ts = sys.Ts), hs)
                # determine a standard reduced system if sys is a standard system
                SV = svd!(Er)
                di2 = Diagonal(1 ./sqrt.(SV.S))
                return dss(di2*SV.U'*Ar*SV.Vt'*di2, I, di2*(SV.U'*Br), (Cr*SV.Vt')*di2, Dr, Ts = sys.Ts), hs
            else
                Tl = Matrix(qr!(lmul!(R',view(SV.U,:,i1))).Q)
                Tr = Matrix(qr!(lmul!(S,SV.V[:,i1])).Q)
                # build the minimal system
                standsys || (return dss(Tl'*SF.S*Tr, Tl'*SF.T*Tr, Tl'*bs, cs*Tr, sys.D, Ts = sys.Ts), hs)
                # determine a standard reduced system if sys is a standard system
                SV = svd!(Tl'*Tr)
                di2 = Diagonal(1 ./sqrt.(SV.S))
                return dss(di2*SV.U'*Tl'*SF.T*Tr*SV.Vt'*di2, I, di2*SV.U'*(Tl'*bs), (cs*Tr)*SV.Vt'*di2, sys.D, Ts = sys.Ts), hs
            end
         else
            return sys, hs  # keep original system if order is preserved
         end
    end
    # end GBALMR
end
function lsminreal2_lrtran(A::AbstractMatrix, E::AbstractMatrix,
                    B::AbstractVecOrMat, C::AbstractMatrix, D::AbstractVecOrMat;
                    withQ::Bool = false, withZ::Bool = false,
                    atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))),
                    rtol::Real =  (size(A,1)+1)*eps(real(float(one(real(eltype(A))))))*iszero(max(atol1,atol2)),
                    fast::Bool = true, finite::Bool = true, infinite::Bool = true,
                    contr::Bool = true, obs::Bool = true, noseig::Bool = true)
   #
   # lsminreal2_ltran(A, E, B, C, D; withQ = false, withZ = false, fast = true, atol1 = 0, atol2 = 0, rtol,
   #                  finite = true, infinite = true, contr = true, obs = true, noseig = true)
   #                  -> (Ar, Er, Br, Cr, Dr, Q, Z, nuc, nuo, nse)

   # This is a special version of lsminreal2 to also determine the left and right transformation matrices
   # Q = [Q1 Q2] and Z = [Z1 Z2], respectively, such that the matrices Ar, Er, Br, and Cr of the resulting descriptor system
   # (Ar-λEr,Br,Cr,Dr) are given by Ar = Q1'*A*Z1, Er = Q1'*E*Z1, Br = Q1'*B, Cr = C*Z1, where the number of columns of Q1 and Z1
   # is equal to the order of matrix Ar. Q and Z result orthogonal if noseig = false.
   # Q = nothing if withQ = false and Z = nothing if withZ = false.
   n = LinearAlgebra.checksquare(A)
   (n,n) != size(E) && throw(DimensionMismatch("A and E must have the same dimensions"))
   p, m = size(D,1), size(D,2)
   n1, m1 = size(B,1), size(B,2)
   (n,m) == (n1, m1) || throw(DimensionMismatch("A, B and D must have compatible dimensions"))
   (p,n) == size(C) || throw(DimensionMismatch("A, C and D must have compatible dimensions"))
   T = promote_type(eltype(A), eltype(E), eltype(B), eltype(C), eltype(D))
   T <: BlasFloat || (T = promote_type(Float64,T))
   ONE = one(T)

   A1 = copy_oftype(A,T)
   E1 = copy_oftype(E,T)
   B1 = copy_oftype(B,T)
   C1 = copy_oftype(C,T)
   D1 = copy_oftype(D,T)

   withQ ? Q = Matrix{T}(I,n,n) : Q = nothing
   withZ ? Z = Matrix{T}(I,n,n) : Z = nothing

   n == 0 && (return A1, E1, B1, C1, D1, Q, Z, 0, 0, 0)

   # save system matrices
   Ar = copy(A1)
   Br = copy(B1)
   Cr = copy(C1)
   Dr = copy(D1)
   Er = copy(E1)
   ir = 1:n
   iz1 = ir; iz2 = n+1:n;
   if finite
      if contr
         m == 0 &&  (ir = 1:0; return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr, Q, Z, n, 0, 0)
         Q1, Z1, _, nr, nfuc = sklf_rightfin!(Ar, Er, Br, Cr; fast, atol1, atol2, rtol, withQ, withZ)
         if nfuc > 0
            ir = 1:nr
            # save intermediary results
            A1 = Ar[ir,ir]
            E1 = Er[ir,ir]
            B1 = Br[ir,:]
            C1 = Cr[:,ir]
            iz1 = ir; iz2 = nr+1:n
            withQ && (Q = copy(Q1))
            withZ && (Z = copy(Z1))
         else
            # restore original matrices
            Ar = copy(A1)
            Er = copy(E1)
            Br = copy(B1)
            Cr = copy(C1)
         end
      else
         nfuc = 0
         nr = n
      end
      if obs
         p == 0 &&  (ir = 1:0; return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr,
                                      withQ ? Q[:,[iz1;iz2]] : Q, withZ ? Z[:,[iz1;iz2]] : Z, nfuc, nr, 0)
         Q1, Z1, _, no, nfuo = sklf_leftfin!(view(Ar,ir,ir), view(Er,ir,ir), view(Cr,:,ir), view(Br,ir,:);
                                           fast, atol1, atol2, rtol, withQ, withZ)
         if nfuo > 0
             iz2 = [ir[1:end-no];iz2]
             ir = ir[end-no+1:end]
             # save intermediary results
             A1 = Ar[ir,ir]
             E1 = Er[ir,ir]
             B1 = Br[ir,:]
             C1 = Cr[:,ir]
             withQ && (Q[:,iz1] = Q[:,iz1]*Q1)
             withZ && (Z[:,iz1] = Z[:,iz1]*Z1)
             iz1 = ir
          else
             # restore saved matrices
             Ar[ir,ir] = A1
             Er[ir,ir] = E1
             Br[ir,:] = B1
             Cr[:,ir] = C1
         end
      else
         nfuo = 0
      end
   else
      nfuc = 0
      nfuo = 0
   end
   if infinite
      if contr
         m == 0 &&  (ir = 1:0; return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr,
                                      withQ ? Q[:,[iz1;iz2]] : Q, withZ ? Z[:,[iz1;iz2]] : Z, n, 0, 0)
         Q1, Z1, _, nr, niuc = sklf_rightfin!(view(Er,ir,ir), view(Ar,ir,ir), view(Br,ir,:), view(Cr,:,ir);
                                           fast, atol1, atol2, rtol, withQ, withZ)
         if niuc > 0
            iz2 = [ir[nr+1:end];iz2]
            ir = ir[1:nr]
            # save intermediary results
            A1 = Ar[ir,ir]
            E1 = Er[ir,ir]
            B1 = Br[ir,:]
            C1 = Cr[:,ir]
            withQ && (Q[:,iz1] = Q[:,iz1]*Q1)
            withZ && (Z[:,iz1] = Z[:,iz1]*Z1)
            iz1 = ir
        else
            # restore original matrices
            Ar[ir,ir] = A1
            Er[ir,ir] = E1
            Br[ir,:] = B1
            Cr[:,ir] = C1
         end
      else
         niuc = 0
      end
      if obs
        p == 0 &&  (ir = 1:0; return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr,
                                      withQ ? Q[:,[iz1;iz2]] : Q, withZ ? Z[:,[iz1;iz2]] : Z, niuc, nr, 0)
         Q1, Z1, _, no, niuo = sklf_leftfin!(view(Er,ir,ir), view(Ar,ir,ir), view(Cr,:,ir), view(Br,ir,:);
                                           fast, atol1, atol2, rtol, withQ, withZ)
         if niuo > 0
             iz2 = [ir[1:end-no];iz2]
             ir = ir[end-no+1:end]
             # save intermediary results
             A1 = Ar[ir,ir]
             E1 = Er[ir,ir]
             B1 = Br[ir,:]
             C1 = Cr[:,ir]
             withQ && (Q[:,iz1] = Q[:,iz1]*Q1)
             withZ && (Z[:,iz1] = Z[:,iz1]*Z1)
             iz1 = ir
         else
             # restore saved matrices
             Ar[ir,ir] = A1
             Er[ir,ir] = E1
             Br[ir,:] = B1
             Cr[:,ir] = C1
         end
      else
          niuo = 0
      end
   else
      niuc = 0
      niuo = 0
   end
   nuc = nfuc+niuc
   nuo = nfuo+niuo
   if noseig
      nm = length(ir)
      withQ ? Q1 = Matrix{T}(I,nm,nm) : Q1 = nothing
      withZ ? Z1 = Matrix{T}(I,nm,nm) : Z1 = nothing
      rE, rA22  = _svdlikeAE!(view(Ar,ir,ir), view(Er,ir,ir), Q1, Z1, view(Br,ir,:), view(Cr,:,ir);
                              fast, atol1, atol2, rtol, withQ, withZ)
      if rA22 > 0
         withQ && (Q[:,iz1] = Q[:,iz1]*Q1)
         withZ && (Z[:,iz1] = Z[:,iz1]*Z1)
         i1 = ir[1:rE]
         i2 = ir[rE+1:rE+rA22]
         # make A22 = I
         fast ? (A22 = UpperTriangular(Ar[i2,i2])) : (A22 = Diagonal(Ar[i2,i2]))
         ldiv!(A22,view(Ar,i2,i1))
         ldiv!(A22,view(Br,i2,:))
         withQ && rdiv!(view(Q,:,i2),A22')
         # apply simplified residualization formulas
         mul!(Dr, view(Cr,:,i2), view(Br,i2,:), -ONE, ONE)               # Dr -= Cr[:,i2]*Br[i2,:]
         mul!(view(Br,i1,:), view(Ar,i1,i2), view(Br,i2,:), -ONE, ONE)   # Br[i1,:] -= Ar[i1,i2]*Br[i2,:]
         mul!(view(Cr,:,i1), view(Cr,:,i2), view(Ar,i2,i1), -ONE, ONE)   # Cr[:,i1] -= Cr[:,i2]*Ar[i2,i1]
         mul!(view(Ar,i1,i1), view(Ar,i1,i2), view(Ar,i2,i1), -ONE, ONE) # Ar[i1,i1] -= Ar[i1,i2]*Ar[i2,i1]
         withQ &&  mul!(view(Q,:,i1), view(Q,:,i2), Ar[i1,i2]', -ONE, ONE) # (Q[:,i1] -= Q[:,i2]*Ar[i1,i2]')
         withZ &&  mul!(view(Z,:,i1), view(Z,:,i2), view(Ar,i2,i1), -ONE, ONE) # (Z[:,i1] -= Z[:,i2]*Ar[i2,i1])
         iz2 = [ir[rE+1:rE+rA22];iz2]
         ir =  [i1; ir[rE+rA22+1:end]]
         iz1 = ir
      else
         # restore saved matrices
         Ar[ir,ir] = A1
         Er[ir,ir] = E1
         Br[ir,:] = B1
         Cr[:,ir] = C1
      end
      return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr, withQ ? Q[:,[iz1;iz2]] : Q, withZ ? Z[:,[iz1;iz2]] : Z, nuc, nuo, rA22
   else
      return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr, withQ ? Q[:,[iz1;iz2]] : Q, withZ ? Z[:,[iz1;iz2]] : Z, nuc, nuo, 0
   end
end


text:
using Compose, DataFrames, Gadfly

set_default_plot_size(9inch, 6inch)

theme1 = Theme(point_size=6pt)
coord1 = Coord.cartesian(xmin=0.0, xmax=6.0)
D = DataFrame(x=1:5,
              y=[0.768448, 0.940515, 0.673959, 0.395453, 0.313244],
              V1=["A","A","B","B","D"],
              V2 = string.([1,2,2,3,3])  )

pa = plot(x=1:5, y=[0.77, 0.94, 0.67, 0.39, 0.31], shape=["A","A","B","B","D"], theme1, coord1,
    Guide.shapekey(title="Key",labels=["α","β","δ"]),
    Guide.title("Guide.shapekey") )

pb = plot(D, x=:x, y=:y, shape=:V1, color=:V1,
        Scale.shape_discrete(levels=["D","A","B"]),
        theme1, coord1, Guide.title("Shape==Color") )

pc = plot(D, x=:x, y=:y, shape=:V1, color=:V2,  coord1,
        Guide.colorkey(title="Color"),
        Guide.shapekey(title="Shape "),
        Theme(point_size=6pt, key_swatch_color="slategrey"),
        Guide.title("Shape!=Color") )

pd = plot(D, x=:x, y=:y, shape=[Shape.square], size=[8pt])
pe = plot()
pf = plot(Theme(point_size=6pt, alphas=[0.6]),
            layer(x=0.5:4, y=D.y[1:4], shape=["Gr1"]),
             layer(x=1:4, y=D.y[1:4], shape=["Gr2"]) )

gridstack([pa pb pc; pd pe pf])

ext: jl


text:
using RandomMatrixDistributions
using Test

n = 1000
p = 400
γ = p/n

# test white Wishart
for beta in [1, 4, 8]
    println("Testing white Wishart with β = ", beta)

    λs = randeigvals(SpikedWishart(beta, n, p))

    @test length(λs) == p
    @test eltype(λs) <: Real

    @test minimum(λs) > 0

    println("Testing Jacobi with β = ", beta)
    λs = randeigvals(Jacobi(beta, n, n, p))

    @test length(λs) == p
    @test eltype(λs) <: Real

    @test minimum(λs) > 0
end

println("Testing Spiked Wishart with β = 1")
s = 8
λs = randeigvals(SpikedWishart(1, n, p, [2, 4, s]))/n

@test length(λs) == p
@test eltype(λs) <: Real

# Check if the maximum spike is reasonable
# This tests can fail in principle, but only with negligible probability
μ = (1 + s) * (1 + γ/s)
σ = (1 + s) * sqrt(γ * 2 * (1 - γ/s^2))

@test abs(maximum(λs) - μ)/σ < 4


text:
@testset "misc" begin
    # finite.jl
    a = ["aa", "bb", "aa", "bb"] |> categorical
    @test scitype(a[1]) == Multiclass{2}

    # schema show
    df = DataFrame(x=[1.0,2.0,3.0],y=["a","b","c"])
    s = schema(df)
    io = IOBuffer()
    show(io, MIME("text/plain"), MLJScientificTypes.schema(df))
    @test String(take!(io)) == "┌─────────┬─────────┬────────────┐\n│ _.names │ _.types │ _.scitypes │\n├─────────┼─────────┼────────────┤\n│ x       │ Float64 │ Continuous │\n│ y       │ String  │ Textual    │\n└─────────┴─────────┴────────────┘\n_.nrows = 3\n"

    # coerce
    x = Any['a', 5]
    @test (@test_logs (:info, "Char value encountered, such value will be coerced according to the corresponding numeric value (e.g. 'A' to 65).") coerce(x, Count)) == [97, 5]
    x = categorical(['a','b','a','b'])
    @test coerce(x, Continuous) == [1.0,2.0,1.0,2.0]
    y = [missing, 1, 2]
    x = y[2:end]
    c = MLJScientificTypes._tighten_if_needed(x, eltype(x), true)
    @test c == [1, 2]
    @test eltype(c) == Int
    y = categorical([missing,1,2])
    x = y[2:end]
    @test eltype(x) >: Missing
    c = MLJScientificTypes._tighten_if_needed(x, eltype(x), true)
    @test c == [1,2]
    @test eltype(c) == Int
    c = coerce(x, OrderedFactor, tight=true)
    @test c == categorical([1,2])
    @test !(eltype(c) >: Missing)

    # increase autotype coverage
    M = MLJScientificTypes
    @test M.string_to_multiclass(String, ["a","b"], 0) == String
    @test M.string_to_multiclass(Textual, ["a","b"], 0) == Multiclass
    @test M.string_to_multiclass(Textual, ["a","b", missing], 0) == Union{Missing,Multiclass}

    # explicit scitype test
    S = ScientificTypes
    @test S.Scitype(Int, M.MLJ()) == Count
    @test S.Scitype(Float64, M.MLJ()) == Continuous
    @test S.Scitype(SubString, M.MLJ()) == Textual

    X = [1,2,3]
    @test elscitype(X) == Count
    S.Scitype(::Type{Float16}, ::MLJScientificTypes.MLJ) = Count
    Xf = Float16[1,2,3]
    @test elscitype(Xf) == Count
end

@testset "Schema" begin
    M = MLJScientificTypes
    sch = M.Schema((:a, :b), (Int, Int), (Count, Count), 5)
    @test sch isa M.Schema{(:a, :b),Tuple{Int64,Int64},Tuple{Count,Count},5}
    @test sch.names == (:a, :b)
    @test sch.types == (Int, Int)
    @test sch.scitypes == (Count, Count)
    @test sch.nrows == 5

    @test_throws ArgumentError sch.something
    @test propertynames(sch) == (:names, :types, :scitypes, :nrows)

    X = [1,2,3]
    @test_throws ArgumentError schema(X)
end

@testset "Err" begin
    x = [1,2,3,4]
    @test_throws MLJScientificTypes.CoercionError coerce!(x, Continuous)
    x = (1,2,3,4)
    @test_throws MLJScientificTypes.CoercionError coerce(x, Continuous)
end


text:
# <editor-fold desc="reporting of calculation progress and error handling"

# XXX return elapsed time since Start_date
function getElapsed(start::DateTime)
    elapSec_per = Dates.value(floor(now() - start,Dates.Second(1)))
    if elapSec_per < 3600*24
        elap_str = Dates.format(DateTime(2015,01,01,Int(floor(elapSec_per / 3600)),Int(floor(elapSec_per % 3600/ 60)),elapSec_per % 60), "HH:MM:SS")
    else
        elap_str = Dates.format(DateTime(2015,01,Int(floor(elapSec_per / (3600*24))),Int(floor(elapSec_per % (3600*24) / 3600)),Int(floor(elapSec_per % 3600/ 60)),elapSec_per % 60), "dd:HH:MM:SS")
    end
    return elap_str
end

# XXX teste for errors so far and optional writes report file, even if no serious errrors occured yet
function errorTest(report::Array{Tuple,1},options::modOptions;write::Bool = false, inCode::Bool = false)
    errStatus_dic = Dict(1 => :green, 2 => :yellow,3 => :red)
    if any(getindex.(report,1) .== 3)
		output_df = DataFrame(map(idx -> getindex.(report, idx), eachindex(first(report))), [:type, :section, :location, :message])
        CSV.write("$(options.outDir)/reporting_$(options.outStamp).csv",  insertcols!(output_df[!,2:end], 1, :errStatus => map(x -> errStatus_dic[x],output_df[!,:type])))
		printstyled("$(inCode ? "" : " - " )Errors encountered! Wrote reporting_$(options.outStamp).csv for details!"; color = :light_red)
        error()
    else
		numWarn = length(findall(getindex.(report,1) .== 2))
        if write && length(report) > 0
			output_df = DataFrame(map(idx -> getindex.(report, idx), eachindex(first(report))), [:type, :section, :location, :message])
            CSV.write("$(options.outDir)/reporting_$(options.outStamp).csv",  insertcols!(output_df[!,2:end], 1, :errStatus => map(x -> errStatus_dic[x],output_df[!,:type])))
			printstyled("$(inCode ? "" : " - " )No errors and $numWarn warning(s) encountered. Wrote reporting_$(options.outStamp).csv for details! \n"; color = numWarn > 0 ? :light_yellow : :light_green)
        else
			printstyled("$(inCode ? "" : " - " )No errors and $numWarn warning(s) encountered. \n"; color = numWarn > 0 ? :light_yellow : :light_green)
        end
    end
end

# XXX produces a output message and tests for errors accordingly to globally set reporting values
function produceMessage(options::modOptions,report::Array{Tuple,1},currentLvl::Int64,fixedString::String,dynamicString::Any="")
	sty_dic = Dict(1 => :bold, 2 => :normal, 3 => :light_black)

	sty_dic[currentLvl]
    if options.reportLvl >= currentLvl
		if options.errCheckLvl >= currentLvl
			printstyled(options.objName; color = :underline); printstyled(" ", getElapsed(options.startTime), fixedString, dynamicString; color = sty_dic[currentLvl])
		else
			printstyled(options.objName; color = :underline); printstyled(" ",getElapsed(options.startTime), fixedString, dynamicString, "\n"; color = sty_dic[currentLvl])
		end
	end
    if options.errCheckLvl >= currentLvl errorTest(report,options,write = options.errWrtLvl >= currentLvl) end
end

# </editor-fold>

# <editor-fold desc="miscellaneous data processing"

# XXX new plus function to avoid error when one element being added up is nothing
plus(a::Int,b::Int) = a + b
plus(a::Int,b::Nothing) = a
plus(a::Nothing,b::Int) = b

# XXX provides names of columns as array of symbols ('names' function itself was changed from strings to symbols)
namesSym(df::DataFrame) = map(x -> Symbol(x),names(df))
namesSym(df::DataFrameRow) = map(x -> Symbol(x),names(df))

# XXX returns dataframe columns without value column
removeVal(input_df::DataFrame) = filter(x -> !(x in (:val,:ratio)),namesSym(input_df))
removeVal(col_arr::Array{Symbol,1}) = filter(x -> !(x in (:val,:ratio)),col_arr)

# XXX return an empty integer array instead of an error, if a key is not in a dictionary
getDicEmpty(dic::Dict,key::Any) = key in keys(dic) ? dic[key] : Int[]

# XXX get names of column of type integer
intCol(in_df::DataFrame) = getindex.(filter(x -> eltype(x[2]) <: Int, collect(pairs(eachcol(in_df)))),1)
intCol(in_df::DataFrame,add_sym::Symbol) = union(intCol(in_df),intersect(namesSym(in_df),[add_sym]))

# XXX puts relevant dimensions in consistent order and adds remaining entries at the end
orderDim(inDim_arr::Array{Symbol,1},intCol_arr::Array{Symbol,1}) = intersect([:Ts_exp, :Ts_expSup, :Ts_disSup, :Ts_dis, :R_exp, :R_dis, :R_from, :R_to, :C, :Te], intersect(inDim_arr,intCol_arr)) |> (x -> [x...,setdiff(inDim_arr,x)...])
orderDim(inDim_arr::Array{Symbol,1}) = intersect([:Ts_exp, :Ts_expSup, :Ts_disSup, :Ts_dis, :R_exp, :R_dis, :R_from, :R_to, :R_a, :R_b, :C, :Te], inDim_arr) |> (x -> [x...,setdiff(inDim_arr,x)...])

# XXX puts dataframes columns in consistent order
orderDf(in_df::DataFrame) = select(in_df,orderDim(namesSym(in_df),intCol(in_df) |> (z -> isempty(z) ? Symbol[] : z)))

# XXX writes all tuples occuring in a tuple of pairs and tuples
mixedTupToTup(x) = typeof(x) <: Pair ? map(y -> mixedTupToTup(y),collect(x)) :  x

# XXX check if dataframe should be considered, if energy balance is created for carriers in array
filterCarrier(var_df::DataFrame,c_arr::Array{Int,1}) = :C in namesSym(var_df) ? filter(r -> r.C in c_arr,var_df) : var_df

# XXX create dataframe with all potential dimensions for carrier provided
function createPotDisp(c_arr::Array{Int,1},anyM::anyModel)

	lvl_arr = map(x -> anyM.cInfo[x], c_arr) |> (y -> map(z -> getfield.(y,z),[:tsDis, :rDis]))
	allLvl_df = DataFrame(C = c_arr, lvlTs = lvl_arr[1], lvlR = lvl_arr[2])
	tsDis_dic, rDis_dic = [Dict(x => getfield.(getNodesLvl(anyM.sets[z[2]],x),:idx) for x in unique(lvl_arr[z[1]])) for z in enumerate([:Ts,:R])]

	allLvl_df[!,:Ts_dis] = map(x -> tsDis_dic[x],allLvl_df[!,:lvlTs])
	allLvl_df[!,:R_dis] = map(x -> rDis_dic[x],allLvl_df[!,:lvlR])

	var_df = flatten(flatten(select(allLvl_df,Not([:lvlTs,:lvlR])),:Ts_dis),:R_dis)

	# add column for superordinate dispatch timestep
	supTs_dic =  Dict(x => getAncestors(x,anyM.sets[:Ts],:int,anyM.supTs.lvl)[end] for x in unique(var_df[!,:Ts_dis]))
	var_df[!,:Ts_disSup] = map(x -> supTs_dic[x], var_df[!,:Ts_dis])

	return var_df
end

# XXX gets technology name as symbol from id and the other way around
techSym(tInt::Int,tech_tree::Tree) = Symbol(getUniName(tInt,tech_tree)[end])
techInt(tSym::Symbol,tech_tree::Tree) = filter(x -> x.val == string(tSym),collect(values(tech_tree.nodes)))[1].idx

# </editor-fold>

# <editor-fold desc="data frame based manipulations"

# XXX finds entries where expansion or capacity would be fixed to zero
function filterZero(src_df::DataFrame,par_obj::ParElement,anyM::anyModel)
	if isdefined(par_obj,:name)
	# copies parameter obj and adds ":up" to inheritance for any dimensions, otherwise variables would be created, but fixed to zero due to a zero limit on a higher level in the tree
		modPar_obj = par_obj
		modPar_obj.herit = modPar_obj.herit |> (y -> tuple(vcat(y..., map(x -> x => :up,getindex.(y,1))...)...))
		# filter zero cases
		zero_df = select!(filter(r -> r.val == 0, matchSetParameter(src_df, modPar_obj, anyM.sets)),Not(:val))
	else
		zero_df = src_df[[],:]
	end
	return zero_df
end

# XXX removes all entries occuring in remove array from input table
function removeEntries(remove_arr::Array{DataFrame,1},input_df::DataFrame)
    if !isempty(remove_arr)
        remove_df = length(remove_arr) == 1 ? remove_arr[1] : vcat(remove_arr...)
        colRemove_arr = namesSym(remove_df)
		out_df = antijoin(input_df,remove_df; on = colRemove_arr)
		return out_df
    else
        return input_df
    end
end

# XXX merge provided dataframe into prep_dic
function mergePrepDic!(key_sym::Symbol,prep_dic::Dict{Symbol,NamedTuple},capaResi_df::DataFrame,capaRatio_df::DataFrame = DataFrame())
	if key_sym in keys(prep_dic)
		prep_dic[key_sym]= (var = prep_dic[key_sym].var, ratio = capaRatio_df, resi = capaResi_df)
	else
		prep_dic[key_sym] = (var = intCol(capaResi_df) |> (x -> DataFrame(Pair.(x,fill(Int[],length(x))))), ratio = capaRatio_df, resi = capaResi_df)
	end
end

# XXX performs a left or outer join operation and replaces any missing values
function joinMissing(leftData_df::DataFrame, rightData_df::DataFrame, key_arr::Union{Array{Symbol,1},Array{Pair{Symbol,Symbol},1}}, how_sym::Symbol, missVal_dic::Dict, uni_boo::Bool = false)

	# perform join operation
	if how_sym == :left
		joinData_df = leftjoin(leftData_df,rightData_df; on = key_arr, makeunique = uni_boo)
	elseif how_sym == :outer
		joinData_df = outerjoin(leftData_df,rightData_df; on = key_arr, makeunique = uni_boo)
	end

	miss_col = filter(x -> any(ismissing.(x[2])), collect(pairs(eachcol(joinData_df))))
    # check, if any column contains missing values
    if isempty(miss_col) return dropmissing(joinData_df) end

    # replace missing value, cases differ depending if data type needs to be adjusted
    for col in miss_col
        joinData_df[!,col[1]] = map(x -> coalesce(x,missVal_dic[col[1]]),col[2])
    end

    return dropmissing(joinData_df)
end

# XXX get array of scaling factors for add_df
function getResize(add_df::DataFrame,time_obj::Tree,supDis::NamedTuple{(:lvl,:step,:sca),Tuple{Int,Tuple{Vararg{Int,N} where N},Dict{Tuple{Int,Int},Float64}}})
    tsDisLvl_dic = Dict(x => x == 0 ? 1 : getfield(time_obj.nodes[x],:lvl) for x in unique(add_df[!,:Ts_dis]))
	lvl_arr = map(x -> tsDisLvl_dic[x],add_df[!,:Ts_dis])
	aboveSupResize_fl = maximum(values(supDis.sca)) * length(supDis.step) # scaling value used for variables above the superordinate dispatch level
	sca_arr = map(x -> supDis.lvl > x[1] ? aboveSupResize_fl : supDis.sca[(x[2],x[1])] ,zip(lvl_arr,add_df[!,:Ts_disSup]))
    return sca_arr
end

# XXX gets the upper bound used for dispatch variables
function getUpBound(in_df::DataFrame,dispBound_fl::Float64,supTs::NamedTuple{(:lvl,:step,:sca),Tuple{Int,Tuple{Vararg{Int,N} where N},Dict{Tuple{Int,Int},Float64}}},treeTs::Tree)
	if !isnan(dispBound_fl)
		upBound_arr = dispBound_fl * getResize(in_df,treeTs,supTs)
	else
		upBound_arr = fill(NaN,size(in_df,1))
	end
	return upBound_arr
end

# </editor-fold>

# <editor-fold desc="functions and sub-functions to aggregate variables"

# XXX aggregates variables in aggEtr_df to rows in srcEtr_df, function used, if all entries of search have the same resolution (all entries in a relevant column are on the same level)
function aggUniVar(aggEtr_df::DataFrame, srcEtr_df::DataFrame, agg_arr::Array{Symbol,1},srcRes_tup::NamedTuple,sets_dic::Dict{Symbol,Tree})
	if isempty(aggEtr_df) return fill(AffExpr(),size(srcEtr_df,1)) end

	# only selects relevant columns
    aggEtr_df = select(aggEtr_df,vcat(:var,agg_arr...))
	srcEtr_df = select(srcEtr_df,agg_arr)

	# adjusts entries in aggregation dataframe to comply with resolution of search dataframe
	for dim in intersect(keys(srcRes_tup),agg_arr)
		set_sym = Symbol(split(string(dim),"_")[1])
		dim_dic = Dict(x => getAncestors(x,sets_dic[set_sym],:int,getfield(srcRes_tup,dim))[end] for x in unique(aggEtr_df[!,dim]))
		aggEtr_df[!,dim] .= map(x -> dim_dic[x],aggEtr_df[!,dim])
	end

	aggEtrGrp_df = combine(groupby(aggEtr_df,agg_arr), :var => (x -> sum(x)) => :var)
	var_arr = joinMissing(srcEtr_df,aggEtrGrp_df,agg_arr,:left,Dict(:var => AffExpr()))[!,:var]
	return var_arr
end

# XXX aggregates variables in aggEtr_df to rows in srcEtr_df, function used, if entries of search can have different resolutions (not all entries in a relevant column are on the same level)
function aggDivVar(aggEtr_df::DataFrame, srcEtr_df::DataFrame, agg_tup::Tuple, sets_dic::Dict{Symbol,Tree}; aggFilt::Tuple = ())

	# XXX sanity checks regarding columns
	if all(namesSym(aggEtr_df) |> (y -> map(x -> !(x in y),agg_tup))) error("tried to perform aggregation on column not existing in dataframe to be aggregated") end
	if all(namesSym(srcEtr_df) |> (y -> map(x -> !(x in y),agg_tup))) error("tried to perform aggregation on column not existing in dataframe to aggregate") end

	select!(aggEtr_df,intCol(aggEtr_df,:var))
	# XXX filter entries from aggEtr_df, that based on isolated analysis of columns will not be aggregated
	for dim in intersect(aggFilt,agg_tup)
		set_sym = Symbol(split(string(dim),"_")[1])
		allSrc_set = unique(srcEtr_df[!,dim]) |> (z -> union(BitSet(z),map(x -> BitSet(getDescendants(x,sets_dic[set_sym],true)),z)...))
		aggEtr_df = aggEtr_df[findall(map(x -> (x in allSrc_set),aggEtr_df[!,dim])),:]
	end

	if isempty(aggEtr_df) return fill(AffExpr(),size(srcEtr_df,1)) end

	# XXX filter entries from srcEtr_df, that based on isolated anlysis of columns will not have any values aggregated to
	idxRel_set = BitSet(1:size(srcEtr_df,1))
	for dim in agg_tup
		set_sym = Symbol(split(string(dim),"_")[1])
		allAgg_set = unique(aggEtr_df[!,dim]) |> (z -> union(BitSet(z),map(y -> BitSet(getAncestors(y,sets_dic[set_sym],:int,0)), z)...))
		idxRel_set = intersect(idxRel_set,BitSet(findall(map(x -> x in allAgg_set, srcEtr_df[!,dim]))))
	end
	srcEtrAct_df = srcEtr_df[collect(idxRel_set),:]
	# group aggregation dataframe to relevant columns and removes unrequired columns
	aggEtrGrp_df = combine(groupby(aggEtr_df,collect(agg_tup)), :var => (x -> sum(x)) => :var)

	# XXX create dictionaries in each dimension that assign rows suited for aggregation for each value
	chldRows = Dict{Symbol,Dict{Int,BitSet}}()
	for col in agg_tup
		# row that are potentially aggregated
		findCol_arr = aggEtrGrp_df[!,col]
		findCol_set = BitSet(findCol_arr)

		# entries that other entries can be aggregated to
		searchVal_set = BitSet(unique(srcEtrAct_df[!,col]))

		# to every unique value in column the value itself and its children are assigned
		set_sym = Symbol(split(string(col),"_")[1])
		idxChild_dic = Dict(x => intersect(findCol_set,[x,getDescendants(x,sets_dic[set_sym],true)...]) for x in searchVal_set)

		# for each unique value in column the rows with children are assigned
		grp_df = groupby(DataFrame(val = findCol_arr, id = 1:length(findCol_arr)),:val)
		dicVal_dic = Dict(x.val[1] => BitSet(sort(x[!,:id])) for x in grp_df) |> (dic -> Dict(x => union(map(y -> dic[y],collect(idxChild_dic[x]))...) for x in keys(idxChild_dic)))
		# excludes column from search, if based on it, every entry in find could be aggregated to every row in search
		# (if this holds true for all columns, make an exception for the last one and dont remove it to, because otherwise correct aggregation cannot happen )
		if all(length.(values(dicVal_dic)) .== length(findCol_arr)) && !(col == agg_tup[end] && length(chldRows) < 1)
			select!(srcEtrAct_df,Not(col)); continue
		else
			chldRows[col] = dicVal_dic
		end
	end

	# XXX finds aggregation by intersecting suited rows in each dimension
	if isempty(chldRows)
		aggRow_arr = fill(BitSet(),size(srcEtrAct_df,1))
	else
		aggRow_arr = collectKeys(keys(chldRows)) |> (y -> map(x -> intersect(map(y -> chldRows[y][x[y]],y)...) ,eachrow(srcEtrAct_df)))
	end

	# XXX aggregates values according to lookup
	out_arr = Array{AffExpr}(undef,size(srcEtr_df,1))
	out_arr[collect(idxRel_set)] =  map(x -> sum(aggEtrGrp_df[x,:var]), collect.(aggRow_arr))
	out_arr[setdiff(1:size(srcEtr_df,1),idxRel_set)] .= AffExpr()

	return out_arr
end

# </editor-fold>

# <editor-fold desc="manipulate model related data frames"

# XXX add superordinate dispatch timestep to expansion dataframe
function addSupTsToExp(expMap_df::DataFrame,para_obj::Dict{Symbol,ParElement},type_sym::Symbol,tsYear_dic::Dict{Int,Int},anyM::anyModel)
	if !isempty(expMap_df)
		lftm_df = matchSetParameter(flatten(expMap_df,:Ts_expSup),para_obj[Symbol(:life,type_sym)],anyM.sets,newCol = :life)
		lftmDel_df = matchSetParameter(lftm_df,para_obj[Symbol(:del,type_sym)],anyM.sets,newCol = :del)
		lftmDel_df[!,:Ts_disSup] = map(x -> filter(y -> (tsYear_dic[y] >= tsYear_dic[x.Ts_expSup] + x.del) && (tsYear_dic[y] <= tsYear_dic[x.Ts_expSup] + x.life + x.del),collect(anyM.supTs.step)), eachrow(lftmDel_df))
		select!(lftmDel_df,Not([:life,:del]))
		grpCol_arr = intCol(expMap_df) |> (x -> :ratio in namesSym(expMap_df) ? vcat(:ratio,x...) : x)
		expMap_df = combine(groupby(lftmDel_df,grpCol_arr), [:Ts_expSup,:Ts_disSup] .=> (x -> [x]) .=> [:Ts_expSup,:Ts_disSup])
	else
		expMap_df[!,:Ts_disSup] = Array{Array{Int,1},1}()
	end
	return expMap_df
end

# XXX expand expansion dataframe to capacity dataframe
function expandExpToCapa(in_df::DataFrame)

	noExpCol_arr = intCol(in_df)

	allDf_arr = map(eachrow(in_df)) do x
		l_arr = length.(x.Ts_disSup)
		rem_df = repeat(DataFrame(x[noExpCol_arr]), inner = sum(l_arr), outer = 1)
		ext_df = DataFrame(Ts_expSup = vcat(map(y -> fill(x.Ts_expSup[y],l_arr[y]),1:length(l_arr))...), Ts_disSup = vcat(x.Ts_disSup...))
		return hcat(rem_df,ext_df)
	end

	if !isempty(allDf_arr)
		capa_df = select(vcat(allDf_arr...),orderDim(namesSym(allDf_arr[1])))[!,Not(:Ts_exp)]
	else
		 capa_df = select(in_df,Not(:Ts_exp)); capa_df[!,:Ts_disSup] = Int[];
	end

	return orderDf(capa_df)
end

# XXX expands any table including columns with temporal and spatial dispatch levels and the corresponding expansion regions and superordinate dispatch steps to full dispatch table
function expandExpToDisp(inData_df::DataFrame,ts_dic::Dict{Tuple{Int,Int},Array{Int,1}},r_dic::Dict{Tuple{Int64,Int64},Array{Int64,1}},preserveTsSupTs::Bool = false)
    # adds regional timesteps and check if this causes non-unique values (because spatial expansion level can be below dispatch level)
	expR_df = unique(combine(x -> (R_dis = r_dic[(x.R_exp[1],x.lvlR[1])],), groupby(inData_df,namesSym(inData_df)))[!,Not([:R_exp,:lvlR])])
	expTs_df = combine(x -> (Ts_dis = ts_dic[(x.Ts_disSup[1],x.lvlTs[1])],), groupby(expR_df,namesSym(expR_df)))[!,Not(:lvlTs)]

    # adds dispatch timesteps to table and returns
	if !preserveTsSupTs select!(expTs_df,Not(:Ts_disSup)) end
	return expTs_df
end

# XXX obtains residual capacities for technologies
function checkResiCapa(var_sym::Symbol, stockCapa_df::DataFrame, part::AbstractModelPart, anyM::anyModel, addSym::Symbol = Symbol())
  resiPar_sym = Symbol(var_sym,:Resi,addSym)
   if resiPar_sym in tuple(keys(part.par)...)
	   # search for defined residual values
	  stock_df = filter(r -> r.val != 0.0, matchSetParameter(stockCapa_df, part.par[resiPar_sym], anyM.sets))
   else
	   stock_df = filter(x -> false,stockCapa_df)
	   stock_df[!,:val] = Float64[]
   end

	# convers returned value to affine expression
   stock_df[!,:var] =  AffExpr.(stock_df[!,:val])
   select!(stock_df,Not(:val))

   return stock_df
end

# XXX get a dataframe with all variable of the specified type
function getAllVariables(va::Symbol,anyM::anyModel; reflectRed::Bool = true, filterFunc::Function = x -> true)

	varToPart_dic = Dict(:exc => :exc, :capaExc => :exc, :oprCapaExc => :exc, :expExc => :exc, :crt => :bal, :lss => :bal, :trdSell => :trd, :trdBuy => :trd, :emission => Symbol())
	techSym_arr = collect(keys(anyM.parts.tech))

	if va == :excDir
		va = :exc
	end

	if !(va in keys(varToPart_dic)) # get all variables for technologies
		va_dic = Dict(:stIn => (:stExtIn, :stIntIn), :stOut => (:stExtOut, :stIntOut))
		techType_arr = filter(x -> !isempty(x[2]),[(vaSpec,filter(y -> vaSpec in keys(anyM.parts.tech[y].var), techSym_arr)) for vaSpec in (va in keys(va_dic) ? va_dic[va] : (va,))])
		if !isempty(techType_arr)
			allVar_df = vcat(map(x -> anyM.parts.tech[x[2]].var[x[1]], vcat(map(x -> collect(zip(fill(x[1],length(x[2])),x[2])),techType_arr)...))...)
		else
			allVar_df = DataFrame()
		end
	elseif va != :emission # get variables from other parts
		if va in keys(getfield(anyM.parts,varToPart_dic[va]).var)
			allVar_df = getfield(anyM.parts,varToPart_dic[va]).var[va]
		else
			allVar_df = DataFrame()
		end
	else va == :emission # for emission all use variables are obtained and then already matched with emission factors

		if !(:emissionFac in keys(anyM.parts.lim.par))
			lock(anyM.lock)
			push!(anyM.report,(2,"limits","emissionUp","upper emission limits but no emission factors provided"))
			unlock(anyM.lock)
			allVar_df = DataFrame()
		else
			# get all carriers and technologies that might be relevant to compute emissions
			if :Te in namesSym(anyM.parts.lim.par[:emissionFac].data)
				emC_arr = unique(vcat(map(x -> [x,getDescendants(x,anyM.sets[:C],true)...],unique(filter(x -> x.Te == 0, anyM.parts.lim.par[:emissionFac].data)[!,:C]))...))
				emTe_arr = unique(vcat(map(x -> [x,getDescendants(x,anyM.sets[:Te],true)...],unique(filter(x -> x.Te != 0, anyM.parts.lim.par[:emissionFac].data)[!,:Te]))...))
			else
				emC_arr = unique(vcat(map(x -> [x,getDescendants(x,anyM.sets[:C],true)...],unique(anyM.parts.lim.par[:emissionFac].data[!,:C]))...))
				emTe_arr = Array{Int64,1}()
			end

			# get use variables
			allVar_df = getAllVariables(:use,anyM, filterFunc = x -> x.C in emC_arr || x.Te in emTe_arr)

			# get expressions for storage and exchange losses, if this is enabled
			if anyM.options.emissionLoss

				# get all carriers being stored
				allSt_arr = unique(vcat(vcat(map(x -> map(y -> collect(x.carrier[y]),intersect(keys(x.carrier),(:stExtIn,:stExtOut,:stIntIn,:stIntOut))),values(anyM.parts.tech))...)...))
				if !isempty(intersect(emC_arr,vcat(map(x -> [x,getDescendants(x,anyM.sets[:C],true)...],allSt_arr)...)))
					# get all storage variables where storage losses can lead to emissions
					stVar_dic = Dict((string(st) |> (y -> Symbol(uppercase(y[1]),y[2:end]))) => getAllVariables(st,anyM, filterFunc = x -> x.C in emC_arr || x.Te in emTe_arr) for st in (:stIn,:stOut))
					stLvl_df = getAllVariables(:stLvl,anyM, filterFunc = x -> x.C in emC_arr)

					# loop over relevant storage technologies to obtain loss vallues
					tSt_arr = unique(stLvl_df[!,:Te])
					for tInt in tSt_arr
						part = anyM.parts.tech[techSym(tInt,anyM.sets[:Te])]
						# add expression quantifying storage losses for storage in- and and output
						for st in keys(stVar_dic)
							stVar_df = stVar_dic[st]
							stVar_df = matchSetParameter(filter(x -> x.Te == tInt,stVar_df),part.par[Symbol(:eff,st)],anyM.sets)
							stVar_df[!,:var] = stVar_df[!,:var] .* (1 .- stVar_df[!,:val])
							select!(stVar_df,Not(:val))
							allVar_df = vcat(allVar_df,stVar_df)
						end

						# add expression quantifying storage losses for storage discharge
						if :stDis in keys(part.par)
							sca_arr = getResize(stLvl_df,anyM.sets[:Ts],anyM.supTs)
							stLvl_df = matchSetParameter(filter(x -> x.Te == tInt,stLvl_df),part.par[:stDis],anyM.sets)
							stLvl_df[!,:var] = stLvl_df[!,:var] .* (1 .- (1 .- stLvl_df[!,:val]) .^ sca_arr)
							select!(stLvl_df,Not(:val))
							allVar_df = vcat(allVar_df,stLvl_df)
						end
					end
				end

				# add expressions for exchange losses
                if :exc in keys(anyM.parts.exc.var)
					exc_df = getAllVariables(:exc,anyM, filterFunc = x -> x.C in emC_arr)
					exc_df = getExcLosses(convertExcCol(exc_df),anyM.parts.exc.par,anyM.sets)
					# exchange losses are equally split between import and export region
					filter!(x -> x.loss != 0.0,exc_df)
					if !isempty(exc_df)
		                exc_df[!,:var] = exc_df[!,:var] .* exc_df[!,:loss] .* 0.5
						exc_df = rename(combine(groupby(vcat(exc_df,rename(exc_df,:R_a => :R_b,:R_b => :R_a)),filter(x -> x != :R_b,intCol(exc_df))),:var => (x -> sum(x)) => :var),:R_a => :R_dis)
						# dimensions not relevant for exchange are set to 0
						exc_df[!,:Te] .= 0; exc_df[!,:Ts_expSup] .= 0; exc_df[!,:M] .= 0
						allVar_df = vcat(allVar_df,exc_df)
					end
				end
			end

			allVar_df = matchSetParameter(allVar_df,anyM.parts.lim.par[:emissionFac],anyM.sets)
		end

		if !isempty(allVar_df)
			allVar_df[!,:var] = allVar_df[!,:val]  ./ 1e6 .* allVar_df[!,:var]
			select!(allVar_df,Not(:val))
		end
	end

	if !(va in (:capaConv,:capaStIn,:capaStOut,:capaStSize,:oprCapaConv,:oprCapaStIn,:oprCapaStOut,:oprCapaStSize,:expConv,:expStIn,:expStOut,:expStSize)) && !isempty(allVar_df) && reflectRed
		allVar_df[!,:var] .= allVar_df[!,:var] .* anyM.options.redStep
	end

	return filter(filterFunc,allVar_df)
end

# XXX replaces orginal carriers in var_df with all leafes connected to respective carrier (and itself) and flattens it
function replCarLeafs(var_df::DataFrame,c_tree::Tree;cCol::Symbol=:C,noLeaf::Array{Int,1} = Int[])

	cToLeafs_dic = Dict(x => filter(y -> isempty(c_tree.nodes[y].down) || y in noLeaf,[x,getDescendants(x,c_tree,true)...]) for x in unique(var_df[!,cCol]))
	var_df[!,:C] = map(x -> cToLeafs_dic[x],var_df[!,cCol])
	var_df = flatten(var_df,:C)

	return var_df
end

# XXX returns array of technologies and respective dispatch variables relevant for input carrier
function getRelTech(c::Int,tech_dic::Dict{Symbol,TechPart},c_tree::Tree)

	techSym_arr = collect(keys(tech_dic))
	relTech_arr = Array{Tuple{Symbol,Symbol},1}()
	for tSym in techSym_arr
		addConvTech_arr = intersect((:use,:gen),filter(y -> c in tech_dic[tSym].carrier[y], collect(keys(tech_dic[tSym].carrier))))
		if isempty(c_tree.nodes[c].down) # actual dispatch variables for storage only exists for carriers that are leaves
			addStTech_arr = intersect((:stExtIn,:stExtOut),filter(y -> c in union(map(z -> union([z],getDescendants(z,c_tree,true)),tech_dic[tSym].carrier[y])...), collect(keys(tech_dic[tSym].carrier))))
		else
			addStTech_arr = Array{Tuple{Int,Symbol},1}()
		end
		union(addConvTech_arr,addStTech_arr) |> (y -> append!(relTech_arr,collect(zip(fill(tSym,length(y)),y))))
	end

	return relTech_arr
end

# </editor-fold>


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""SecretReference represents a Secret Reference. It has enough information to retrieve secret in any namespace

    IoK8sApiCoreV1SecretReference(;
        name=nothing,
        namespace=nothing,
    )

    - name::String : Name is unique within a namespace to reference a secret resource.
    - namespace::String : Namespace defines the space within which the secret name must be unique.
"""
mutable struct IoK8sApiCoreV1SecretReference <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    namespace::Any # spec type: Union{ Nothing, String } # spec name: namespace

    function IoK8sApiCoreV1SecretReference(;name=nothing, namespace=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1SecretReference, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1SecretReference, Symbol("namespace"), namespace)
        setfield!(o, Symbol("namespace"), namespace)
        o
    end
end # type IoK8sApiCoreV1SecretReference

const _property_map_IoK8sApiCoreV1SecretReference = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("namespace")=>Symbol("namespace"))
const _property_types_IoK8sApiCoreV1SecretReference = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("namespace")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1SecretReference }) = collect(keys(_property_map_IoK8sApiCoreV1SecretReference))
Swagger.property_type(::Type{ IoK8sApiCoreV1SecretReference }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1SecretReference[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1SecretReference }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1SecretReference[property_name]

function check_required(o::IoK8sApiCoreV1SecretReference)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1SecretReference }, name::Symbol, val)
end


text:
using DrWatson
using BSON

includet("core.jl")
includet("bacterium_step.jl")
includet("food_step.jl")
includet("runner.jl")
includet("plotter.jl")

max_stars_repo_name: jiaangou/SimpleEvolution

max_issues_repo_name: jiaangou/SimpleEvolution

max_forks_repo_name: jiaangou/SimpleEvolution


text:
abstract type StochasticDiffEqCache <: DECache end
abstract type StochasticDiffEqConstantCache <: StochasticDiffEqCache end
abstract type StochasticDiffEqMutableCache <: StochasticDiffEqCache end

mutable struct StochasticCompositeCache{T,F} <: StochasticDiffEqCache
  caches::T
  choice_function::F
  current::Int
end

function alg_cache(alg::algType,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{T}}) where {T,algType<:StochasticCompositeAlgorithm}
  caches = map((x)->alg_cache(x,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,Val{T}),alg.algs)
  StochasticCompositeCache(caches,alg.choice_function,1)
end


text:
using PortHamiltonian
using PyPlot
using LinearAlgebra
using ForwardDiff

include("dataexperiment.jl")

# load structural parameters
pslosh, pbeam, prigid = dataexperiment(0.4);

L = pslosh["a"];
N = 20;
disc = discrete_phs_closed(N, 0, L)

Br = [Matrix(1I, N, N); zeros(N,N);zeros(1,N)];
B = blkdiag(zeros(2*N,1),[1]);
D = blkdiag(disc.D,[0]);
Jrb = [0];
invM = disc.Q[N+1:end,N+1:end]
R = 000000.00*Br*norm(invM)*Br';
J = blkdiag(disc.J,Jrb)-R;

# defining Q for saint-venant equation (no rotation)
rho = pslosh["rho"];
g = pslosh["g"];
b = pslosh["b"];
hbar = pslosh["h"];
xi,w = PortHamiltonian.lgwt(N,0,L);
Q = blkdiag(disc.Q,[0])
pfluid = PortHamiltonian.Phs(J, B, D, Q);
#PortHamiltonian.set_constraint!(pfluid, B[:,1],[0])

function Hamiltonian(X)
	alpha1 = X[1:N]
	alpha2 = X[N+1:2*N]
	p = X[end]
	M = w' * (alpha1 .* alpha2)
	H = (w'*((alpha1 .* (alpha2.^2))/rho .+ rho*g*((alpha1).^2)/b)/2 .+ ((p.-M).^2)/(2*prigid["mass"]))[1]
	return H
end

hess = x->ForwardDiff.hessian(Hamiltonian, x)
xeq = [ones(N,1)*b*hbar;zeros(N,1);0.001*prigid["mass"]*0][:]
Qlin = hess(xeq[:])
pfluid.Q = Qlin
pfluid.Hamiltonian = x-> Hamiltonian(x)
pfluid.GradHam = x-> ForwardDiff.gradient(Hamiltonian, x)
pfluid.hessian = x-> hess(x)

#pfluid.GradHam = x -> Q*x
function dynamics(t, X, Xd)
	Xd[1:end]=  pfluid.J * pfluid.GradHam(X[1:end]) + 0.2*pfluid.B[:,2] * sin(2*pi*0.5*t)
	#r[end:end] = pfluid.G' * pfluid.hessian(X[1:end-1]) * (pfluid.J * pfluid.GradHam(X[1:end-1]) + pfluid.B[:,3] * sin(2*pi*t) *0.5+ pfluid.G * X[end])
end

#
 using Sundials
#
 t = range(0,stop = 8,length = 1000)
 yout = Sundials.cvode(dynamics,xeq, collect(t))
 #yout, ypout = Sundials.idasol(dynamics2, [xeq;l0;0], [xeq;0;0]*0, [t])
 pfluid.GradHam = x-> Qlin*x
 youtlin = Sundials.cvode(dynamics,xeq, collect(t))
 #youtlin, ypoutlin = Sundials.idasol(dynamics2, [xeq;l0; 0], [xeq;0; 0]*0, [t])
 #pn = PortHamiltonian.constraint_elimination(pfluid)
 using PyPlot

 figure(1);
 surf(xi,t,yout[:,N+1:end-1])
 #surf(xi,t,yout[:,1:N])
 #include("anima.jl")

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
const Mul{Style, Factors<:Tuple} = Applied{Style, typeof(*), Factors}

const MulArray{T, N, Args} = ApplyArray{T, N, typeof(*), Args}

const MulVector{T, Args} = MulArray{T, 1, Args}
const MulMatrix{T, Args} = MulArray{T, 2, Args}

Mul(A...) = applied(*, A...)

check_applied_axes(A::Mul) = check_mul_axes(A.args...)

size(M::Mul, p::Int) = size(M)[p]
axes(M::Mul, p::Int) = axes(M)[p]
ndims(M::Mul) = ndims(last(M.args))

_mul_ndims(::Type{Tuple{A}}) where A = ndims(A)
_mul_ndims(::Type{Tuple{A,B}}) where {A,B} = ndims(B)
ndims(::Type{<:Mul{<:Any,Args}}) where Args = _mul_ndims(Args)

length(M::Mul) = prod(size(M))
size(M::Mul) = length.(axes(M))

@inline _eltypes() = tuple()
@inline _eltypes(A, B...) = tuple(eltype(A), _eltypes(B...)...)

@inline eltype(M::Mul) = _mul_eltype(_eltypes(M.args...)...)

@inline mulaxes1(::Tuple{}) = ()
@inline mulaxes1(::Tuple{}, B, C...) = mulaxes1(B, C...)
@inline mulaxes1(A::Tuple, C...) = first(A)
@inline mulaxes2(::Tuple{}) = ()
@inline mulaxes2(::Tuple{}, B, C...) = mulaxes2(B, C...)
@inline mulaxes2(A::Tuple{<:Any}, C...) = ()
@inline mulaxes2(A::Tuple{<:Any,<:Any}, C...) = last(A)

@inline _combine_axes(::Tuple{}, ::Tuple{}) = ()
@inline _combine_axes(a, ::Tuple{}) = (a,)
@inline _combine_axes(a, b) = (a,b)
@inline mulaxes(ax...) = _combine_axes(mulaxes1(ax...), mulaxes2(reverse(ax)...))

@inline axes(M::Mul) = mulaxes(map(axes,M.args)...)
@inline axes(M::Mul{<:Any, Tuple{}}) = ()

# *(A::Mul, B::Mul) = apply(*,A.args..., B.args...)
# *(A::Mul, B) = apply(*,A.args..., B)
# *(A, B::Mul) = apply(*,A, B.args...)
⋆(A...) = Mul(A...)

function show(io::IO, A::Mul)
    if length(A.args) == 0
        print(io, "⋆()")
        return
    end
    print(io, first(A.args))
    for a in A.args[2:end]
        print(io, '⋆', a)
    end
end

####
# Matrix * Array
####

ApplyStyle(::typeof(*), args::Type{<:AbstractArray}...) = mulapplystyle(MemoryLayout.(args)...)

"""
   lmaterialize(M::Mul)

materializes arrays iteratively, left-to-right.
"""
lmaterialize(M::Mul) = _lmaterialize(M.args...)

_lmaterialize(A, B) = apply(*,A,B)
_lmaterialize(A, B, C, D...) = _lmaterialize(apply(*,A,B), C, D...)

# arguments for something that is a *
@inline _arguments(::ApplyLayout{typeof(*)}, A) = arguments(A)
@inline _arguments(_, A) = (A,)
@inline _arguments(A) = _arguments(MemoryLayout(typeof(A)), A)

@inline __flatten(A::Tuple{<:Any}, B::Tuple) = (A..., _flatten(B...)...)
@inline __flatten(A::Tuple, B::Tuple) = _flatten(A..., B...)

@inline _flatten() = ()
@inline _flatten(A, B...) = __flatten(_arguments(A), B)
@inline flatten(A) = _mul(_flatten(_arguments(A)...)...)

copy(M::Mul{DefaultArrayApplyStyle,<:Tuple{<:Any,<:Any}}) = copyto!(similar(M), M)
copy(A::Mul{DefaultArrayApplyStyle}) = flatten(lmaterialize(A))

struct FlattenMulStyle <: ApplyStyle end

copy(A::Mul{FlattenMulStyle}) = materialize(flatten(A))

####
# MulArray
#####

_mul(A) = A
_mul(A,B,C...) = Mul(A,B,C...)

_mul_colsupport(j, Z) = colsupport(Z,j)
_mul_colsupport(j, Z::AbstractArray) = colsupport(Z,j)
_mul_colsupport(j, Z, Y...) = axes(Z,1) # default is return all
function _mul_colsupport(j, Z::AbstractArray, Y...)
    rws = colsupport(Z,j)
    a = size(Z,1)+1
    b = 0
    for k in rws
        cs = _mul_colsupport(k, Y...)
        a = min(a,first(cs))
        b = max(b,last(cs))
    end
    a:b
end

colsupport(B::Mul, j) = _mul_colsupport(j, reverse(B.args)...)
colsupport(B::MulArray, j) = _mul_colsupport(j, reverse(B.args)...)

function _getindex(M::Mul, ::Tuple{<:Any}, k::Integer)
    A,Bs = first(M.args), tail(M.args)
    B = _mul(Bs...)
    ret = zero(eltype(M))
    for j = rowsupport(A, k) ∩ colsupport(B,1)
        ret += A[k,j] * B[j]
    end
    ret
end

_getindex(M::Mul, ax, k::Integer) = M[Base._ind2sub(ax, k)...]
getindex(M::Mul, k::Integer) = _getindex(M, axes(M), k)
@propagate_inbounds getindex(A::Mul{LazyArrayApplyStyle}, k::Integer) = Applied{DefaultArrayApplyStyle}(A)[k]

getindex(M::Mul, k::CartesianIndex{1}) = M[convert(Int, k)]
getindex(M::Mul, kj::CartesianIndex{2}) = M[kj[1], kj[2]]

function getindex(M::Mul, k::Integer, j::Integer)
    A,Bs = first(M.args), tail(M.args)
    B = _mul(Bs...)
    ret = zero(eltype(M))
    @inbounds for ℓ in (rowsupport(A,k) ∩ colsupport(B,j))
        ret += A[k,ℓ] * B[ℓ,j]
    end
    ret
end

@propagate_inbounds getindex(A::Mul{LazyArrayApplyStyle}, k::Integer, j::Integer) =
    Applied{DefaultArrayApplyStyle}(A)[k,j]

_flatten(A::MulArray, B...) = _flatten(Applied(A), B...)
flatten(A::MulArray) = ApplyArray(flatten(Applied(A)))

adjoint(A::MulArray) = ApplyArray(*, reverse(map(adjoint,A.args))...)
transpose(A::MulArray) = ApplyArray(*, reverse(map(transpose,A.args))...)

###
# sub materialize
###

# determine rows/cols of multiplication
__mul_args_rows(kr, a) = (kr,)
__mul_args_rows(kr, a, b...) =
    (kr, __mul_args_rows(rowsupport(a,kr), b...)...)
_mul_args_rows(kr, a, b...) = __mul_args_rows(rowsupport(a,kr), b...)
__mul_args_cols(jr, z) = (jr,)
__mul_args_cols(jr, z, y...) =
    (__mul_args_cols(colsupport(z,jr), y...)..., jr)
_mul_args_cols(jr, z, y...) = __mul_args_cols(colsupport(z,jr), y...)

sublayout(::ApplyLayout{typeof(*)}, _...) = ApplyLayout{typeof(*)}()

call(::ApplyLayout{typeof(*)}, V::SubArray) = *

function _mat_mul_arguments(args, (kr,jr))
    kjr = intersect.(_mul_args_rows(kr, args...), _mul_args_cols(jr, reverse(args)...))
    map(view, args, (kr, kjr...), (kjr..., jr))
end

_vec_mul_view(a...) = view(a...)
_vec_mul_view(a::AbstractVector, kr, ::Colon) = view(a, kr)

# this is a vector view of a MulVector
function _vec_mul_arguments(args, (kr,))
    kjr = intersect.(_mul_args_rows(kr, args...), _mul_args_cols(Base.OneTo(1), reverse(args)...))
    _vec_mul_view.(args, (kr, kjr...), (kjr..., :))
end

# this is a vector view of a MulMatrix
_vec_mul_arguments(args, (kr,jr)::Tuple{AbstractVector,Number}) =
    _mat_mul_arguments(args, (kr,jr))

# this is a row-vector view
_vec_mul_arguments(args, (kr,jr)::Tuple{Number,AbstractVector}) =
    _vec_mul_arguments(reverse(map(transpose, args)), (jr,kr))

_mat_mul_arguments(V) = _mat_mul_arguments(arguments(parent(V)), parentindices(V))
_vec_mul_arguments(V) = _vec_mul_arguments(arguments(parent(V)), parentindices(V))

arguments(::ApplyLayout{typeof(*)}, V::SubArray{<:Any,2}) = _mat_mul_arguments(V)
arguments(::ApplyLayout{typeof(*)}, V::SubArray{<:Any,1}) = _vec_mul_arguments(V)

@inline sub_materialize(::ApplyLayout{typeof(*)}, V) = apply(*, arguments(V)...)
@inline copyto!(dest::AbstractArray{T,N}, src::SubArray{T,N,<:ApplyArray{T,N,typeof(*)}}) where {T,N} =
    copyto!(dest, Applied(src))

##
# adoint Mul
##

adjointlayout(::Type, ::ApplyLayout{typeof(*)}) = ApplyLayout{typeof(*)}()
transposelayout(::ApplyLayout{typeof(*)}) = ApplyLayout{typeof(*)}()

call(::ApplyLayout{typeof(*)}, V::Adjoint) = *
call(::ApplyLayout{typeof(*)}, V::Transpose) = *

arguments(::ApplyLayout{typeof(*)}, V::Adjoint) = reverse(adjoint.(arguments(V')))
arguments(::ApplyLayout{typeof(*)}, V::Transpose) = reverse(transpose.(arguments(V')))

##
# * specialcase
##

for op in (:*, :\)
    @eval broadcasted(::DefaultArrayStyle{N}, ::typeof($op), a::Number, b::ApplyArray{<:Number,N,typeof(*)}) where N =
        ApplyArray(*, broadcast($op,a,first(b.args)), tail(b.args)...)
end

broadcasted(::DefaultArrayStyle{N}, ::typeof(/), b::ApplyArray{<:Number,N,typeof(*)}, a::Number) where N =
        ApplyArray(*, most(b.args)..., broadcast(/,last(b.args),a))

for Typ in (:Lmul, :Rmul)
    @eval $Typ(M::Mul) = $Typ(M.args...)
end

##
# L/Rmul
##

struct LmulStyle <: AbstractArrayApplyStyle end
struct RmulStyle <: AbstractArrayApplyStyle end

similar(M::Applied{LmulStyle}, ::Type{T}) where T = similar(Lmul(M), T)
copy(M::Applied{LmulStyle}) = copy(Lmul(M))

similar(M::Applied{RmulStyle}, ::Type{T}) where T = similar(Rmul(M), T)
copy(M::Applied{RmulStyle}) = copy(Rmul(M))

@inline copyto!(dest::AbstractVecOrMat, M::Mul{LmulStyle}) = copyto!(dest, Lmul(M.args...))
@inline copyto!(dest::AbstractVecOrMat, M::Mul{RmulStyle}) = copyto!(dest, Rmul(M.args...))

@inline materialize!(M::Mul{LmulStyle}) = materialize!(Lmul(M))
@inline materialize!(M::Mul{RmulStyle}) = materialize!(Rmul(M))

mulapplystyle(::QLayout, _) = LmulStyle()
mulapplystyle(::QLayout, ::LazyLayout) = LazyArrayApplyStyle()

factorizestyle(_) = DefaultArrayApplyStyle()

for op in (:factorize, :qr, :lu, :cholesky)
    @eval begin
        $op(B::LazyMatrix) = apply($op, B)
        ApplyStyle(::typeof($op), B::Type{<:AbstractMatrix}) = factorizestyle(MemoryLayout(B))
        materialize(A::Applied{DefaultArrayApplyStyle,typeof($op),<:Tuple{<:AbstractMatrix{T}}}) where T =
            Base.invoke($op, Tuple{AbstractMatrix{T}}, A.args...)

        eltype(A::Applied{<:Any,typeof($op)}) = float(eltype(first(A.args)))
    end
end


text:
module Interact

using Reexport

@reexport using InteractBase
import InteractBase: notifications
import Widgets: Widget, @layout, @nodeps
import Observables: @on, @map!, @map

@reexport using OrderedCollections
@reexport using Observables
@reexport using Knockout
@reexport using CSSUtil
@reexport using WebIO
@reexport using Widgets
using RelocatableFolders, Memoize

struct Bulma<:InteractBase.WidgetTheme; end

folder_dir = @path joinpath(@__DIR__, "..")
notebookdir = joinpath(folder_dir, "doc", "notebooks")

@memoize bulma_css() = joinpath(folder_dir, "assets", "bulma.min.css")
@memoize bulma_confined_css() = joinpath(folder_dir, "assets", "bulma_confined.min.css")

function InteractBase.libraries(::Bulma)
    bulmalib = InteractBase.isijulia() ? bulma_confined_css() : bulma_css()
    vcat(InteractBase.font_awesome(), InteractBase.style_css(), bulmalib)
end

function __init__()
    InteractBase.registertheme!(:bulma, Bulma())
    settheme!(Bulma())
    nothing
end

end


text:
export moments

abstract type PolynomialSet end

function JuMP.in_set_string(print_mode, set::PolynomialSet)
    return string(JuMP._math_symbol(print_mode, :in), ' ', set)
end

struct ZeroPoly <: PolynomialSet end
struct NonNegPoly <: PolynomialSet end
struct PosDefPolyMatrix <: PolynomialSet end

function JuMP.function_string(::Type{JuMP.REPLMode},
                              p::MultivariatePolynomials.APL)
    return sprint(show, MIME"text/plain"(), p)
end
function JuMP.function_string(::Type{JuMP.IJuliaMode},
                              p::MultivariatePolynomials.APL)
    # `show` prints `$$` around what `_show` prints.
    return sprint(MultivariatePolynomials._show, MIME"text/latex"(), p)
end

### Shapes for polynomial/moments primal-dual pair ###

# Inspired from `JuMP.dual_shape` docstring example
struct PolynomialShape{BT <: MB.AbstractPolynomialBasis} <: JuMP.AbstractShape
    basis::BT
end
function JuMP.reshape_vector(x::Vector, shape::PolynomialShape)
    return polynomial(x, shape.basis)
end
struct MomentsShape{BT <: MB.AbstractPolynomialBasis} <: JuMP.AbstractShape
    basis::BT
end
function JuMP.reshape_vector(x::Vector, shape::MomentsShape)
    return measure(x, shape.basis)
end
JuMP.dual_shape(shape::PolynomialShape) = MomentsShape(shape.basis)
JuMP.dual_shape(shape::MomentsShape) = PolynomialShape(shape.basis)

JuMP.reshape_set(::ZeroPolynomialSet, ::PolynomialShape) = ZeroPoly()
function JuMP.moi_set(::ZeroPoly, basis::MB.AbstractPolynomialBasis;
                      domain::AbstractSemialgebraicSet=FullSpace()
                     )
    return ZeroPolynomialSet(domain, basis)
end

"""
    moments(cref::JuMP.ConstraintRef)

Return the [`MomentsAttribute`](@ref) of `cref`.
"""
function MultivariateMoments.moments(cref::JuMP.ConstraintRef)
    return MOI.get(cref.model, MomentsAttribute(), cref)
end

"""
    bridges(F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})

Return a list of bridges that may be needed to bridge `F`-in-`S` constraints but
not the bridges that may be needed by constraints added by the bridges.
"""
bridges(::Type{<:MOI.AbstractFunction}, ::Type{<:MOI.AbstractSet}) = []

"""
    bridges(S::Type{<:MOI.AbstractSet})

Return a list of bridges that may be needed to bridge constrained variables
in `S` but not the bridges that may be needed by constraints added by the
bridges.
"""
bridges(S::Type{<:MOI.AbstractSet}) = bridges(MOI.VectorOfVariables, S)

function bridges(::Type{<:MOI.AbstractVectorFunction},
                 ::Type{<:ZeroPolynomialSet{FullSpace}})
    return [ZeroPolynomialBridge]
end

function bridges(::Type{<:MOI.AbstractVectorFunction},
                 ::Type{<:ZeroPolynomialSet{<:AbstractAlgebraicSet}})
    return [ZeroPolynomialInAlgebraicSetBridge]
end

function bridges(::Type{<:MOI.AbstractVectorFunction},
                 ::Type{<:PlusMinusSet})
    return [PlusMinusBridge]
end

"""
    bridgeable(c::JuMP.AbstractConstraint, F::Type{<:MOI.AbstractFunction},
               S::Type{<:MOI.AbstractSet})

Wrap the constraint `c` in `JuMP.BridgeableConstraint`s that may be needed to
bridge `F`-in-`S` constraints.
"""
function bridgeable end

function _concrete(bridge_type::Type{<:MOIB.Variable.AbstractBridge},
                   S::Type{<:MOI.AbstractSet})
    MOIB.Variable.concrete_bridge_type(bridge_type, S)
end
function _concrete(bridge_type::Type{<:MOIB.Constraint.AbstractBridge},
                   S::Type{<:MOI.AbstractSet})
    MOIB.Constraint.concrete_bridge_type(bridge_type, MOI.VectorOfVariables, S)
end
function bridgeable(c::JuMP.AbstractConstraint,
                    S::Type{<:MOI.AbstractSet})
    bridge_types = bridges(S)
    for bridge_type in bridge_types
        c = BridgeableConstraint(c, bridge_type)
        concrete_bridge_type = _concrete(bridge_type{Float64}, S)
        for (ST,) in MOIB.added_constrained_variable_types(concrete_bridge_type)
            c = bridgeable(c, ST)
        end
        for (FT, ST) in MOIB.added_constraint_types(concrete_bridge_type)
            c = bridgeable(c, FT, ST)
        end
    end
    return c
end

function bridgeable(c::JuMP.AbstractConstraint,
                    F::Type{<:MOI.AbstractFunction},
                    S::Type{<:MOI.AbstractSet})
    bridge_types = bridges(F, S)
    for bridge_type in bridge_types
        c = BridgeableConstraint(c, bridge_type)
        concrete_bridge_type = MOIB.Constraint.concrete_bridge_type(
            bridge_type{Float64}, F, S)
        for (ST,) in MOIB.added_constrained_variable_types(concrete_bridge_type)
            c = bridgeable(c, ST)
        end
        for (FT, ST) in MOIB.added_constraint_types(concrete_bridge_type)
            c = bridgeable(c, FT, ST)
        end
    end
    return c
end

### @constraint/@SDconstraint macros ###

non_constant(a::Vector{<:Number}) = convert(Vector{AffExpr}, a)
non_constant(a::Vector{<:JuMP.AbstractJuMPScalar}) = a
non_constant_coefficients(p) = non_constant(coefficients(p))

## ZeroPoly
function JuMP.build_constraint(_error::Function, p::AbstractPolynomialLike,
                               s::ZeroPoly;
                               domain::AbstractSemialgebraicSet=FullSpace(),
                               basis_type = MB.MonomialBasis,
                               kws...)
    coefs = non_constant_coefficients(p)
    basis = MB.basis_covering_monomials(basis_type, monomials(p))
    if domain isa BasicSemialgebraicSet
        # p(x) = 0 for all x in a basic semialgebraic set. We replace it by
        # p(x) ≤ 0 and p(x) ≥ 0 for all x in the basic semialgebraic set.
        # We need to determine the cone two use for `NonNegPoly` which is stored
        # in the `ext` field of the `model` so we need the `model` hence we need
        # to go to `JuMP.add_constraint`.

        # We need to recreate the full list of keyword arguments. `kws` is a
        # `Iterators.Pairs, `kws.data` is a `NamedTuple` and `kws.itr` is a
        # `Tuple`.
        all_kws = Iterators.Pairs(merge((domain=domain,), kws.data),
                                  (:domain, kws.itr...))
        return Constraint(_error, p, s, all_kws)
    else
        set = JuMP.moi_set(s, basis; domain=domain, kws...)
        constraint = JuMP.VectorConstraint(coefs, set, PolynomialShape(basis))
        return bridgeable(constraint, JuMP.moi_function_type(typeof(coefs)),
                          typeof(set))
    end
end
function JuMP.build_constraint(_error::Function, p::AbstractPolynomialLike,
                               s::MOI.EqualTo; kws...)
    return JuMP.build_constraint(_error, p-s.value, ZeroPoly(); kws...)
end

## NonNegPoly and PosDefPolyMatrix
# The `model` is not given in `JuMP.build_constraint` so we create a custom
# `Constraint` object and transform the `set` in `JuMP.add_constraint`.
struct Constraint{FT, PT, ST<:PolynomialSet, KWT<:Iterators.Pairs} <: JuMP.AbstractConstraint
    _error::FT
    polynomial_or_matrix::PT
    set::ST
    kws::KWT
end
function JuMP.build_constraint(_error::Function, polynomial_or_matrix,
                               set::Union{NonNegPoly, PosDefPolyMatrix};
                               kws...)
    return Constraint(_error, polynomial_or_matrix, set, kws)
end

# FIXME the domain will not appear in the printing, it should be a field of
#       `ZeroPoly` which is `FullSpace` by default
JuMP.reshape_set(::PlusMinusSet, ::PolynomialShape) = ZeroPoly()
function JuMP.add_constraint(model::JuMP.Model,
                             constraint::Constraint{<:Any, <:Any, ZeroPoly},
                             name::String = "")
    cone = getdefault(model, NonNegPoly())
    coefs = non_constant_coefficients(constraint.polynomial_or_matrix)
    basis = MB.MonomialBasis(monomials(constraint.polynomial_or_matrix))
    set = PlusMinusSet(JuMP.moi_set(cone, basis; constraint.kws...))
    new_constraint = JuMP.VectorConstraint(coefs, set, PolynomialShape(basis))
    bridgeable_con = bridgeable(
        new_constraint, JuMP.moi_function_type(typeof(coefs)), typeof(set))
    return JuMP.add_constraint(model, bridgeable_con, name)
end

function JuMP.add_constraint(model::JuMP.Model, constraint::Constraint,
                             name::String = "")
    set = getdefault(model, constraint.set)
    new_constraint = JuMP.build_constraint(constraint._error,
                                           constraint.polynomial_or_matrix, set;
                                           constraint.kws...)
    return JuMP.add_constraint(model, new_constraint, name)
end

# NonNegPoly
function JuMP.build_constraint(_error::Function, p::AbstractPolynomialLike,
                               s::MOI.GreaterThan; kws...)
    return JuMP.build_constraint(_error, p-s.lower, NonNegPoly(); kws...)
end
function JuMP.build_constraint(_error::Function, p::AbstractPolynomialLike,
                               s::MOI.LessThan; kws...)
    return JuMP.build_constraint(_error, s.upper-p, NonNegPoly(); kws...)
end

# PosDefPolyMatrix
# there is already a method for AbstractMatrix in PSDCone in JuMP so we need a
# more specific here to avoid ambiguity
function JuMP.build_constraint(_error::Function,
                               p::AbstractMatrix{<:AbstractPolynomialLike},
                               s::PSDCone; kws...)
    return JuMP.build_constraint(_error, p, PosDefPolyMatrix(); kws...)
end


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

@userplot HScatter

@recipe function f(hs::HScatter; lag=0, tol=1e-1, distance=Euclidean())
  # get inputs
  sdata = hs.args[1]
  var₁ = hs.args[2]
  var₂ = length(hs.args) == 3 ? hs.args[3] : var₁

  # lookup valid data
  locs₁ = findall(!ismissing, sdata[var₁])
  locs₂ = findall(!ismissing, sdata[var₂])
  𝒟₁ = view(sdata, locs₁)
  𝒟₂ = view(sdata, locs₂)
  X₁ = [coordinates(centroid(𝒟₁, i)) for i in 1:nelements(𝒟₁)]
  X₂ = [coordinates(centroid(𝒟₂, i)) for i in 1:nelements(𝒟₂)]
  z₁ = 𝒟₁[var₁]
  z₂ = 𝒟₂[var₂]

  # compute pairwise distance
  m, n = length(z₁), length(z₂)
  pairs = [(i,j) for j in 1:n for i in j:m]
  ds = [evaluate(distance, X₁[i], X₂[j]) for (i,j) in pairs]

  # find indices with given lag
  match = findall(abs.(ds .- lag) .< tol)

  if isempty(match)
    @warn "no points were found with lag = $lag, skipping..."
    return nothing
  end

  # scatter plot coordinates
  mpairs = view(pairs, match)
  x = z₁[first.(mpairs)]
  y = z₂[last.(mpairs)]

  x̄ = mean(x)
  ȳ = mean(y)

  xguide --> var₁
  yguide --> var₂

  # plot h-scatter
  @series begin
    seriestype --> :scatter
    seriescolor --> :black
    label --> "samples"

    x, y
  end

  # plot regression line
  @series begin
    seriestype --> :path
    seriescolor --> :red
    label --> "regression"

    X = [x ones(length(x))]
    ŷ = X * (X \ y)

    x, ŷ
  end

  # plot identity line
  @series begin
    seriestype --> :path
    seriescolor --> :black
    linestyle --> :dash
    label --> "identity"

    xmin, xmax = extrema(x)
    ymin, ymax = extrema(y)
    vmin = min(xmin, ymin)
    vmax = max(xmax, ymax)

    [vmin, vmax], [vmin, vmax]
  end

  # plot mean lines
  @series begin
    primary --> false
    seriestype --> :vline
    seriescolor --> :green
    [x̄]
  end
  @series begin
    primary --> false
    seriestype --> :hline
    seriescolor --> :green
    [ȳ]
  end
  @series begin
    primary --> false
    seriestype --> :scatter
    seriescolor --> :green
    marker --> :square
    markersize --> 4
    [(x̄, ȳ)]
  end
end