text:
function rcond(A::DenseMatrix,tola::Real = 0)
   T = eltype(A)
   T1 = T <: BlasFloat ? T : T1 = promote_type(T,Float64)
   nrmA = opnorm(A,1)
   nrmA <= tola && return zero(real(T1))
   istriu(A) ? (return LinearAlgebra.LAPACK.trcon!('1','U','N',copy_oftype(A,T1))) :
       (return LinearAlgebra.LAPACK.gecon!('1', LinearAlgebra.LAPACK.getrf!(copy_oftype(A,T1))[1],real(T1)(nrmA)) )
end
"""
    lseval(A, E, B, C, D, val; atol1, atol2, rtol, fast = true) -> Gval

Evaluate `Gval`, the value of the rational matrix `G(λ) = C*inv(λE-A)*B+D` for `λ = val`.
The computed `Gval` has infinite entries if `val` is a pole (finite or infinite) of `G(λ)`.
If `val` is finite and `val*E-A` is singular or if `val = Inf` and `E` is singular,
then the entries of `Gval` are evaluated separately for minimal realizations of each input-output channel.

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 computation of minimal realizations of individual input-output channels relies on pencil manipulation algorithms,
which 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.
"""
function lseval(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling{Bool}},
                B::AbstractVecOrMat, C::AbstractMatrix, D::AbstractVecOrMat, val::Number;
                atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))),
                rtol::Real =  (size(A,1)+1)*eps(float(real(eltype(A))))*iszero(max(atol1,atol2)), fast::Bool = true)

   T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D),typeof(val))
   E == I || (T = promote_type(T,eltype(E)))
   T <: BlasFloat || (T = promote_type(Float64,T))
   # check dimensions
   n = LinearAlgebra.checksquare(A)
   if typeof(E) <: AbstractMatrix
      n == LinearAlgebra.checksquare(E) || error("A and E must have the same size")
   end
   n == 0 && (return T.(D))

   n1, m = typeof(B) <: AbstractVector ? (length(B),1) : size(B)
   n1 == n ||  error("B must have the same row size as A")
   p, n1 = size(C)
   n1 == n ||  error("C must have the same column size as A")
   p1, m1 = typeof(D) <: AbstractVector ? (length(D),1) : size(D)
   m1 == m ||  error("D must have the same column size as B")
   p1 == p ||  error("D must have the same row size as C")

   toleps = (size(A,1)+1)*eps(real(T))
   if abs(val) < Inf
      LUF = lu!(T(val)*E-A;check = false)
      tol = max(atol1,atol2)
      if rcond(LUF.U, tol) < toleps
         G = zeros(T,p,m)
         for i = 1:p
            At1, Et1, Bt1, Ct1, Dt1, = lsminreal2(A, E, B, view(C,i:i,:), view(D,i:i,:); infinite = false, noseig = false, contr = false, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
            for j = 1:m
                At11, Et11, Bt11, Ct11, Dt11, = lsminreal2(At1, Et1, view(Bt1,:,j:j), Ct1, view(Dt1,:,j:j); infinite = false, noseig = false, obs = false, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
                LUF = lu!(T(val)*Et11-At11;check = false)
                rcond(LUF.U,tol) < toleps ? G[i,j] = T(Inf) : G[i,j] = (Ct11*ldiv!(LUF,copy_oftype(Bt11,T)) + Dt11)[1,1]
            end
         end
         return G
      else
         return C*ldiv!(LUF,copy_oftype(B,T)) + D
      end
   else
      (E == I || rcond(E,atol2) > toleps) && (return T.(D))
      At, Et, Bt, Ct, Dt, = lsminreal2(A, E, B, C, D; finite = false, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
      if size(At,1) > 0 && rcond(Et,atol2) < toleps
         G = zeros(T,p,m)
         for i = 1:p
            At1, Et1, Bt1, Ct1, Dt1, = lsminreal2(At, Et, Bt, view(Ct,i:i,:), view(Dt,i:i,:); finite = false, contr = false, noseig = false, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
            for j = 1:m
                At11, Et11, Bt11, Ct11, Dt11, = lsminreal2(At1, Et1, view(Bt1,:,j:j), Ct1, view(Dt1,:,j:j); finite = false, obs = false, noseig = true, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
                size(At11,1) > 0 && rcond(Et11,atol2) < toleps ? G[i,j] = T(Inf) : G[i,j] = Dt11[1,1]
            end
         end
         return G
      else
         return Dt
      end
   end
end
"""
     lps2ls(A, E, B, F, C, G, D, H; compacted = false, atol1 = 0, atol2 = 0, atol3 = 0, rtol = min(atol1,atol2,atol3)>0 ? 0 : n*ϵ) -> (Ad,Ed,Bd,Cd,Dd)

Construct an input-output equivalent descriptor system linearizations `(Ad-λdE,Bd,Cd,Dd)` to a pencil based linearization
`(A-λE,B-λF,C-λG,D-λH)` satisfying

                -1                        -1
     Cd*(λEd-Ad)  *Bd + Dd = (C-λG)*(λE-A)  *(B-λF) + D-λH .

If `compacted = true`, a compacted linearization is determined by exploiting possible rank defficiencies of the
matrices `F`, `G`, and `H`.

The keyword arguments `atol1`, `atol2`, `atol3`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `F`, the absolute tolerance for the nonzero elements of `G`,
the absolute tolerance for the nonzero elements of `H` and the relative tolerance
for the nonzero elements of `F`, `G` and `H`. The default relative tolerance is `n*ϵ`, where `n` is the size of
 `A`, and `ϵ` is the machine epsilon of the element type of `A`.
"""
function lps2ls(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling}, B::AbstractVecOrMat, F::Union{AbstractVecOrMat,Missing},
                C::AbstractMatrix, G::Union{AbstractMatrix,Missing}, D::AbstractVecOrMat, H::Union{AbstractVecOrMat,Missing};
                compacted::Bool = false, atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))),
                atol3::Real = zero(real(eltype(A))), rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2,atol3)))

   T = promote_type(eltype(A), E == I ? Bool : eltype(E), eltype(B), eltype(C), eltype(D))
   T = promote_type(T, ismissing(F) ? T : eltype(F), ismissing(G) ? T : eltype(G), ismissing(H) ? T :  eltype(H))
   ismissing(F) && ismissing(G) && ismissing(H) && (return A, E, B, C, D)

   n = size(A,1)
   p, m = typeof(D) <: AbstractVector ? (length(D),1) : size(D)
   nm = n+m
   D1 = copy_oftype(D,T)

   (ismissing(F) || iszero(F)) ? mF = 0 : mF = m
   (ismissing(G) || iszero(G)) ? pG = 0 : pG = p
   if compacted
       if ismissing(G) || iszero(G)
           pG = 0
           C11 = zeros(T,p,pG)
           E12 = zeros(T,pG,n)
       else
          S = svd(G,full=true)
          pG = count(S.S .> max(atol2,rtol*S.S[1]))
          if pG == p
             C11 = -I
             E12 = view(G,:,:)
          else
             C11 = -S.U[:,1:pG]
             E12 = Diagonal(view(S.S,1:pG))*view(S.Vt,1:pG,:)
          end
       end
       if ismissing(F) || iszero(F)
           mF = 0
           B22 = zeros(T,mF,m)
           E23 = zeros(T,n,mF)
       else
          # svd is not working on a vector
          typeof(F) <: AbstractVector ? S = svd(reshape(F,n,1),full=true) : S = svd(F,full=true)
          mF = count(S.S .> max(atol1,rtol*S.S[1]))
          if mF == m
             B22 = -I
             E23 = view(F,:,:)
          else
             B22 = -S.Vt[1:mF,:]
             E23 = view(S.U,:,1:mF)*Diagonal(view(S.S,1:mF))
          end
       end
       A1 = [I zeros(T,pG,n+mF); zeros(T,n,pG) A zeros(T,n,mF); zeros(T,mF,pG+n) I]
       E1 = [zeros(T,pG,pG) E12 zeros(T,pG,mF); zeros(T,n,pG) E E23; zeros(T,mF,pG+n+mF)]
       B1 = [zeros(T,pG,m); B; B22]
       C1 = [C11 C zeros(T,p,mF)]
   else
       A1 = [I zeros(T,pG,n+mF); zeros(T,n,pG) A zeros(T,n,mF); zeros(T,mF,pG+n) I]
       E1 = [zeros(T,pG,pG) ismissing(G) ? zeros(T,pG,n) : G[1:pG,:] zeros(T,pG,mF); zeros(T,n,pG) E ismissing(F) ? zeros(T,n,mF) : F[:,1:mF]; zeros(T,mF,pG+n+mF)]
       B1 = [zeros(T,pG,m); B; mF > 0 ? -I : zeros(T,mF,m)]
       C1 = [pG > 0 ? -I : zeros(T,p,pG) C zeros(T,p,mF)]
   end
   (ismissing(H) || iszero(H)) &&  (return A1, E1, B1, C1, D1)

   finished = false
   if compacted
      typeof(H) <: AbstractVector ? S = svd(reshape(H,p,1),full=true) : S = svd(H,full=true)
      r = count(S.S .> max(atol3,rtol*S.S[1]))
      if r < min(p,m)
         hs2 = Diagonal(sqrt.(view(S.S,1:r)))
         A2 = Matrix{T}(I,2*r,2*r)
         E2 = [zeros(T,r,r) I; zeros(T,r,2*r)]
         B2 = [zeros(T,r,m); hs2*view(S.Vt,1:r,:)]
         C2 = [view(S.U,:,1:r)*hs2 zeros(T,p,r)]
         finished = true
      end
   end
   if !finished
      if m <= p
         A2 = Matrix{T}(I,2*m,2*m)
         E2 = [zeros(T,m,m) I; zeros(T,m,2*m)]
         B2 = [zeros(T,m,m); I]
         C2 = [H zeros(T,p,m)]
      else
         A2 = Matrix{T}(I,2*p,2*p)
         E2 = [zeros(T,p,p) I; zeros(T,p,2*p)]
         B2 = [zeros(T,p,m); -H]
         C2 = [-I zeros(T,p,p)]
      end
   end
   return blockdiag(A1,A2), blockdiag(E1,E2), [B1; B2], [C1 C2], D1
end

function blockdiag(mats::AbstractMatrix{T}...) where T
    rows = Int[size(m, 1) for m in mats]
    cols = Int[size(m, 2) for m in mats]
    res = zeros(T, sum(rows), sum(cols))
    m = 1
    n = 1
    for ind=1:length(mats)
        mat = mats[ind]
        i = rows[ind]
        j = cols[ind]
        res[m:m + i - 1, n:n + j - 1] = mat
        m += i
        n += j
    end
    return res
end
"""
    lpseval(A, E, B, F, C, G, D, H, val; atol1, atol2, rtol, fast = true) -> Gval

Evaluate `Gval`, the value of the rational matrix `G(λ) = (C-λG)*inv(λE-A)*(B-λF)+D-λH` for `λ = val`.
The computed `Gval` has infinite entries if `val` is a pole (finite or infinite) of `G(λ)`.
If `val` is finite and `val*E-A` is singular or `val` is infinite and `E` is singular,
then the entries of `Gval` are evaluated separately for each input-output chanel employing descriptor system based
minimal realizations.

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`, `F`, `G`, `H`,
and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D`, `E`, `F`, `G`, and `H`.

The computation of minimal realizations of individual input-output chanels relies on pencil manipulation algorithms,
which 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.
"""
function lpseval(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling{Bool}},
                 B::AbstractVecOrMat, F::AbstractVecOrMat, C::AbstractMatrix, G::AbstractMatrix,
                 D::AbstractVecOrMat, H::AbstractVecOrMat, val::Number; atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))),
                 rtol::Real =  (size(A,1)+1)*eps(float(real(eltype(A))))*iszero(max(atol1,atol2)), fast::Bool = true)

   # check dimensions
   n = LinearAlgebra.checksquare(A)
   if typeof(E) <: AbstractMatrix
      n == LinearAlgebra.checksquare(E) || error("A and E must have the same size")
   end

   n1, m = typeof(B) <: AbstractVector ? (length(B),1) : size(B)
   n1 == n ||  error("B must have the same row size as A")
   n1, m1 = typeof(F) <: AbstractVector ? (length(F),1) : size(F)
   (n1, m1) == (n, m) ||  error("B and F must have the same size")
   p, n1 = size(C)
   n1 == n ||  error("C must have the same column size as A")
   (p, n) == size(G) ||  error("C and G must have the same size")
   p1, m1 = typeof(D) <: AbstractVector ? (length(D),1) : size(D)
   m1 == m ||  error("D must have the same column size as B")
   p1 == p ||  error("D must have the same row size as C")
   p1, m1 = typeof(H) <: AbstractVector ? (length(H),1) : size(H)
   (p, m) == (p1, m1) || error("D and H must have the same size")

   T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D), eltype(F), eltype(G), eltype(H), typeof(val))
   E == I || (T = promote_type(T,eltype(E)))
   T <: BlasFloat || (T = promote_type(Float64,T))

   toleps = (size(A,1)+1)*eps(real(T))
   if abs(val) < Inf
      LUF = lu!(T(val)*E-A;check = false)
      tol = max(atol1,atol2)
      if rcond(LUF.U, tol) < toleps
         At, Et, Bt, Ct, Dt = lps2ls(A, E, B, F, C, G, D, H; compacted = true, atol1 = atol2, atol2 = atol2, atol3 = atol2, rtol = rtol)
         p, m = typeof(Dt) <: AbstractVector ? (length(D),1) : size(Dt)
         G = zeros(T,p,m)
         for i = 1:p
            At1, Et1, Bt1, Ct1, Dt1, = lsminreal2(At, Et, Bt, view(Ct,i:i,:), view(Dt,i:i,:); infinite = false, noseig = false, contr = false, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
            for j = 1:m
                At11, Et11, Bt11, Ct11, Dt11, = lsminreal2(At1, Et1, view(Bt1,:,j:j), Ct1, view(Dt1,:,j:j); infinite = false, noseig = false, obs = false, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
                LUF = lu!(T(val)*Et11-At11;check = false)
                rcond(LUF.U,tol) < toleps ? G[i,j] = T(Inf) : G[i,j] = (Ct11*ldiv!(LUF,copy_oftype(Bt11,T)) + Dt11)[1,1]
            end
         end
         return G
      else
         return (C-val*G)*ldiv!(LUF,copy_oftype(B-val*F,T)) + D-val*H
      end
   else
      return lseval(lps2ls(A, E, B, F, C, G, D, H; compacted = true, atol1 = atol2, atol2 = atol2, atol3 = atol2, rtol = rtol)...,val;
                    atol1 = atol1, atol2 = atol2, rtol = rtol, fast = fast)
   end
end
"""
     lsequal(A1, E1, B1, C1, D1, A2, E2, B2, C2, D2;
             fastrank = true, atol1 = 0, atol2 = 0, rtol = min(atol1,atol2)>0 ? 0 : n*ϵ) -> flag::Bool

Check if two descriptor system linearizations `(A1-λE1,B1,C1,D1)` and `(A2-λE2,B2,C2,D2)` satisfy the equivalence condition

                -1                      -1
     C1*(λE1-A1)  *B1 + D1 = C2*(λE2-A2)  *B2 + D2

The ckeck is performed by computing the normal rank `n` of the structured matrix pencil `M - λN`

              | A1-λE1   0    |   B1  |
              |   0    A2-λE2 |   B2  |
     M - λN = |---------------|-------|
              |   C1     -C2  | D1-D2 |

and verifying that `n = n1+n2`, where `n1` and `n2` are the orders of the square matrices `A1` and `A2`, respectively.

If `fastrank = true`, the rank is evaluated by counting how many singular values of `M - γ N` have magnitude
greater than `max(max(atol1,atol2), rtol*σ₁)`, where `σ₁` is the largest singular value of `M - γ N` and
`γ` is a randomly generated value [1].
If `fastrank = false`, the rank is evaluated as `nr + ni + nf + nl`, where `nr` and `nl` are the sums
of right and left Kronecker indices, respectively, while `ni` and `nf` are the number of infinite and
finite eigenvalues, respectively. The sums `nr+ni` and  `nf+nl`, are determined from an
appropriate Kronecker-like form (KLF) exhibiting the spliting of the right and left structures
of the pencil `M - λN`. For efficiency purpose, the reduction to the relevant KLF is only partially performed
using rank decisions based on rank revealing SVD-decompositions.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `M`, the absolute tolerance for the nonzero elements of `N` and the relative tolerance
for the nonzero elements of `M` and `N`. The default relative tolerance is `k*ϵ`, where `k` is the size of
the smallest dimension of `M`, and `ϵ` is the machine epsilon of the element type of `M`.

[1] A. Varga, On checking null rank conditions of rational matrices, [arXiv:1812.11396](https://arxiv.org/abs/1812.11396), 2018.
"""
function lsequal(A1::AbstractMatrix, E1::Union{AbstractMatrix,UniformScaling{Bool}},
   B1::AbstractVecOrMat, C1::AbstractMatrix, D1::AbstractVecOrMat,
   A2::AbstractMatrix, E2::Union{AbstractMatrix,UniformScaling{Bool}},
   B2::AbstractVecOrMat, C2::AbstractMatrix, D2::AbstractVecOrMat;
   atol1::Real = zero(real(eltype(A1))), atol2::Real = zero(real(eltype(A1))),
   rtol::Real =  (size(A1,1)+size(A2,1)+2)*eps(real(float(one(real(eltype(A1))))))*iszero(max(atol1,atol2)),
   fastrank::Bool = true)

   # implicit dimensional checks are performed using the try-catch scheme
   try
      T = promote_type(eltype(A1), eltype(A2))
      n1 = size(A1,1)
      n2 = size(A2,1)
      A = [A1  zeros(T,n1,n2);
           zeros(T,n2,n1) A2]
      B = [B1; B2;]
      C = [C1 -C2;]
      D = [D1-D2;]
      if E1 == I && E2 == I
         E = I
      else
         E = [E1  zeros(T,n1,n2);
              zeros(T,n2,n1) E2]
      end
      return (sprank(A,E,B,C,D,atol1 = atol1, atol2 = atol2, rtol = rtol, fastrank = fastrank) == n1+n2)
   catch err
      println("$err")
      return false
   end
end
"""
     lpsequal(A1, E1, B1, F1, C1, G1, D1, H1,
              A2, E2, B2, F2, C2, G2, D2, H2; fastrank = true, atol1 = 0, atol2 = 0, rtol = min(atol1,atol2)>0 ? 0 : n*ϵ)
              -> flag::Bool

Check if two pencil based linearizations `(A1-λE1,B1-λF1,C1-λG1,D1-λH1)` and `(A2-λE2,B2-λF2,C2-λG2,D2-λH2)`
satisfy the equivalence condition

                      -1                                      -1
     (C1-λG1)*(λE1-A1)  *(B1-λF1) + D1-λH1 = (C2-λG2)*(λE2-A2)  *(B2-λF2) + D2-λH2

The ckeck is performed by computing the normal rank `n` of the structured matrix pencil `M - λN`

              | A1-λE1    0    |    B1-λF1      |
              |   0     A2-λE2 |    B2-λF2      |
     M - λN = |----------------|----------------|
              | C1-λG1 -C2+λG2 | D1-D2-λ(H1-H2) |

and verifying that `n = n1+n2`, where `n1` and `n2` are the orders of the square matrices `A1` and `A2`, respectively.

If `fastrank = true`, the rank is evaluated by counting how many singular values of `M - γ N` have magnitude
greater than `max(max(atol1,atol2), rtol*σ₁)`, where `σ₁` is the largest singular value of `M - γ N` and
`γ` is a randomly generated value [1].
If `fastrank = false`, the rank is evaluated as `nr + ni + nf + nl`, where `nr` and `nl` are the sums
of right and left Kronecker indices, respectively, while `ni` and `nf` are the number of infinite and
finite eigenvalues, respectively. The sums `nr+ni` and  `nf+nl`, are determined from an
appropriate Kronecker-like form (KLF) exhibiting the spliting of the right and left structures
of the pencil `M - λN`. For efficiency purpose, the reduction to the relevant KLF is only partially performed
using rank decisions based on rank revealing SVD-decompositions.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `M`, the absolute tolerance for the nonzero elements of `N`,  and the relative tolerance
for the nonzero elements of `M` and `N`. The default relative tolerance is `k*ϵ`, where `k` is the size of
the smallest dimension of `M`, and `ϵ` is the machine epsilon of the element type of `M`.

[1] A. Varga, On checking null rank conditions of rational matrices, [arXiv:1812.11396](https://arxiv.org/abs/1812.11396), 2018.
"""
function lpsequal(A1::AbstractMatrix, E1::Union{AbstractMatrix,UniformScaling{Bool}},
   B1::AbstractMatrix, F1::AbstractMatrix, C1::AbstractMatrix, G1::AbstractMatrix,
   D1::AbstractMatrix, H1::AbstractMatrix,
   A2::AbstractMatrix, E2::Union{AbstractMatrix,UniformScaling{Bool}},
   B2::AbstractMatrix, F2::AbstractMatrix, C2::AbstractMatrix, G2::AbstractMatrix,
   D2::AbstractMatrix, H2::AbstractMatrix;
   atol1::Real = zero(real(eltype(A1))), atol2::Real = zero(real(eltype(A1))),
   rtol::Real =  (size(A1,1)+size(A2,1)+4)*eps(real(float(one(real(eltype(A1))))))*iszero(max(atol1,atol2)),
   fastrank::Bool = true)

   # implicit dimensional checks are performed using the try-catch scheme
   try
      T = promote_type(eltype(A1), eltype(A2))
      n1 = size(A1,1)
      n2 = size(A2,1)
      A = [A1  zeros(T,n1,n2);
           zeros(T,n2,n1) A2]
      B = [B1; B2;]
      F = [F1; F2;]
      C = [C1 -C2;]
      G = [G1 -G2;]
      D = [D1-D2;]
      H = [H1-H2;]
      if E1 == I && E2 == I
         E = I
      else
         E = [E1  zeros(T,n1,n2);
              zeros(T,n2,n1) E2]
      end
      return (prank([A B;C D],[E F;G H], atol1 = atol1, atol2 = atol2, rtol = rtol, fastrank = fastrank) == n1+n2)
   catch err
      println("$err")
      return false
   end
end
"""
    lsminreal2(A, E, B, C, D; fast = true, atol1 = 0, atol2 = 0, rtol, finite = true, infinite = true, contr = true, obs = true, noseig = true)
               -> (Ar, Er, Br, Cr, Dr, nuc, nuo, nse)

Reduce the linearization `(A-λE,B,C,D)` of a rational matrix to a reduced form `(Ar-λEr,Br,Cr,Dr)` such that

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

with the least possible order `nr` of `Ar-λEr` 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`.

The achieved dimensional reductions to fulfill conditions (1) and (2), conditions (3) and (4), and
respectively, condition (5) are returned in `nuc`, `nuo`, `nse`.

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 = false`, then condition (5) on the lack of simple infinite eigenvalues is not 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 linearization `(A-λE,B,C,D)`
to obtain an irreducible linearization 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`.

[1] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function lsminreal2(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling{Bool}},
                    B::AbstractVecOrMat, C::AbstractMatrix, D::AbstractVecOrMat;
                    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)

   emat = (typeof(E) <: AbstractMatrix)
   eident = !emat || isequal(E,I)
   n = LinearAlgebra.checksquare(A)
   emat && (n,n) != size(E) && throw(DimensionMismatch("A and E must have the same dimensions"))
   p, m = typeof(D) <: AbstractVector ? (length(D),1) : size(D)
   n1, m1 = typeof(B) <: AbstractVector ? (length(B),1) : size(B)
   (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(B), eltype(C), eltype(D))
   eident || (T = promote_type(T,eltype(E)))
   T <: BlasFloat || (T = promote_type(Float64,T))

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

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

   if eident
      A1, B1, C1, nuc, nuo  = lsminreal(A1, B1, C1, contr = contr, obs = obs, fast = fast, atol = atol1, rtol = rtol)
      return A1, emat ? Matrix{T}(I,size(A1)...) : I, B1, C1, D1, nuc, nuo, 0
   else
      # save system matrices
      Ar = copy(A1)
      Br = copy(B1)
      Cr = copy(C1)
      Dr = copy(D1)
      Er = copy(E1)
      ir = 1:n
      if finite
         if contr
            m == 0 &&  (ir = 1:0; return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr, n, 0, 0)
            _, _, _, nr, nfuc = sklf_rightfin!(Ar, Er, Br, Cr; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false)
            if nfuc > 0
               ir = 1:nr
               # save intermediary results
               A1 = Ar[ir,ir]
               E1 = Er[ir,ir]
               B1 = Br[ir,:]
               C1 = Cr[:,ir]
            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, nfuc, nr, 0)
            _, _, _, no, nfuo = sklf_leftfin!(view(Ar,ir,ir), view(Er,ir,ir), view(Cr,:,ir), view(Br,ir,:);
                                              fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false)
            if nfuo > 0
                ir = ir[end-no+1:end]
                # save intermediary results
                A1 = Ar[ir,ir]
                E1 = Er[ir,ir]
                B1 = Br[ir,:]
                C1 = Cr[:,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, n, 0, 0)
            _, _, _, nr, niuc = sklf_rightfin!(view(Er,ir,ir), view(Ar,ir,ir), view(Br,ir,:), view(Cr,:,ir);
                                              fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol,
                                              withQ = false, withZ = false)
            if niuc > 0
               ir = ir[1:nr]
               # save intermediary results
               A1 = Ar[ir,ir]
               E1 = Er[ir,ir]
               B1 = Br[ir,:]
               C1 = Cr[:,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, niuc, nr, 0)
            _, _, _, no, niuo = sklf_leftfin!(view(Er,ir,ir), view(Ar,ir,ir), view(Cr,:,ir), view(Br,ir,:);
                                              fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol,
                                              withQ = false, withZ = false)
            if niuo > 0
                ir = ir[end-no+1:end]
                # save intermediary results
                A1 = Ar[ir,ir]
                E1 = Er[ir,ir]
                B1 = Br[ir,:]
                C1 = Cr[:,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
         rE, rA22  = _svdlikeAE!(view(Ar,ir,ir), view(Er,ir,ir), nothing, nothing, view(Br,ir,:), view(Cr,:,ir),
                     fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false)
         if rA22 > 0
            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,:))
            # 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; ir[rE+rA22+1:end]]
         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, nuc, nuo, rA22
      else
         return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr, nuc, nuo, 0
      end
   end
end
"""
    lsminreal(A, B, C; fast = true, atol = 0, rtol, contr = true, obs = true, noseig = true)
              -> (Ar, Br, Cr, nuc, nuo)

Reduce the linearization `(A-λI,B,C,0)` of a strictly proper rational matrix to a reduced form `(Ar-λI,Br,Cr,0)` such that

             -1                 -1
     C*(λI-A)  *B  = Cr*(λI-Ar)  *Br

with the least possible order `nr` of `Ar` if `contr = true` and `obs = true`.
Such a realization is called `minimal` and satisfies:

     (1) rank[Br Ar-λI] = nr for all λ (controllability);

     (2) rank[Ar-λI; Cr] = nr for all λ (observability).

The achieved dimensional reductions to fulfill conditions (1) and (2) are returned in `nuc` and `nuo`, respectively.

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

If `contr = false`, then the controllability condition (1) is not enforced.

If `obs = false`, then observability condition (2) is not enforced.

To enforce conditions (1)-(2), orthogonal similarity transformations are performed on
the matrices of the original linearization `(A-λI,B,C,0)` to obtain a minimal linearization using
structured pencil reduction algorithms, as the fast versions of the reduction techniques of the
full row rank pencil [B A-λI] and full column rank pencil [A-λI;C] proposed in [1].

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 `atol` and `rtol`, specify, respectively, the absolute and relative tolerances for the
nonzero elements of matrices `A`, `B`, `C`.

[1] P. Van Dooreen, The generalized eigenstructure problem in linear system theory,
IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.
"""
function lsminreal(A::AbstractMatrix, B::AbstractVecOrMat, C::AbstractMatrix;
                   atol::Real = zero(real(eltype(A))),
                   rtol::Real =  (size(A,1)+1)*eps(real(float(one(real(eltype(A))))))*iszero(atol),
                   fast::Bool = true, contr::Bool = true, obs::Bool = true)
   n = LinearAlgebra.checksquare(A)
   n1, m = typeof(B) <: AbstractVector ? (length(B),1) : size(B)
   n == n1 || throw(DimensionMismatch("A and B must have the same number of rows"))
   (p,n1) = size(C)
   n1 == n || throw(DimensionMismatch("A and C must have the same number of columns"))
   T = promote_type(eltype(A), eltype(B), eltype(C))
   T <: BlasFloat || (T = promote_type(Float64,T))

   A1 = copy_oftype(A,T)
   B1 = copy_oftype(B,T)
   C1 = copy_oftype(C,T)

   n == 0 && (return A1, B1, C1, 0, 0)

   # save system matrices
   Ar = copy(A1)
   Br = copy(B1)
   Cr = copy(C1)
   ir = 1:n
   if contr
      m == 0 &&  (ir = 1:0; return Ar[ir,ir], Br[ir,:], Cr[:,ir], n, 0)
      _, _, nr, nuc = sklf_right!(Ar, Br, Cr; fast = fast, atol1 = atol, atol2 = atol, rtol = rtol, withQ = false)
      if nuc > 0
         ir = 1:nr
         # save intermediary results
         A1 = Ar[ir,ir]
         B1 = Br[ir,:]
         C1 = Cr[:,ir]
      else
         # restore original matrices
         Ar = copy(A1)
         Br = copy(B1)
         Cr = copy(C1)
      end
   else
      nuc = 0
      nr = n
   end
   if obs
      p == 0 &&  (ir = 1:0; return Ar[ir,ir], Br[ir,:], Cr[:,ir], nuc, nr)
      _, _, no, nuo = sklf_left!(view(Ar,ir,ir), view(Cr,:,ir), view(Br,ir,:); fast = fast, atol1 = atol, atol2 = atol, rtol = rtol, withQ = false)
      if nuo > 0
         ir = ir[end-no+1:end]
      else
         # restore saved matrices
         Ar[ir,ir] = A1
         Br[ir,:] = B1
         Cr[:,ir] = C1
      end
   else
      nuo = 0
   end
   return Ar[ir,ir], Br[ir,:], Cr[:,ir], nuc, nuo
end
"""
    lsminreal(A, E, B, C, D; fast = true, atol1 = 0, atol2, rtol, contr = true, obs = true, noseig = true)
              -> (Ar, Er, Br, Cr, Dr, nuc, nuo, nse)

Reduce the linearization `(A-λE,B,C,D)` of a rational matrix to a reduced form `(Ar-λEr,Br,Cr,Dr)` such that

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

with the least possible order `nr` of `Ar-λEr` 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`.

The achieved dimensional reductions to fulfill conditions (1) and (2), conditions (3) and (4), and
respectively, condition (5) are returned in `nuc`, `nuo`, `nse`.

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 linearization `(A-λE,B,C,D)` to obtain an irreducible linearization 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`.

[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 lsminreal(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling{Bool}},
                   B::AbstractVecOrMat, C::AbstractMatrix, D::AbstractVecOrMat;
                   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, contr::Bool = true, obs::Bool = true, noseig::Bool = true)

   emat = (typeof(E) <: AbstractMatrix)
   eident = !emat || isequal(E,I)
   n = LinearAlgebra.checksquare(A)
   emat && (n,n) != size(E) && throw(DimensionMismatch("A and E must have the same dimensions"))
   p, m = typeof(D) <: AbstractVector ? (length(D),1) : size(D)
   n1, m1 = typeof(B) <: AbstractVector ? (length(B),1) : size(B)
   (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(B), eltype(C), eltype(D))
   eident || (T = promote_type(T,eltype(E)))
   T <: BlasFloat || (T = promote_type(Float64,T))

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

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

   if eident
      A1, B1, C1, nuc, nuo  = lsminreal(A1, B1, C1, contr = contr, obs = obs, fast = fast, atol = atol1, rtol = rtol)
      return A1, emat ? Matrix{T}(I,size(A1)...) : I, B1, C1, D1, nuc, nuo, 0
   else
      # save system matrices
      Ar = copy(A1)
      Br = copy(B1)
      Cr = copy(C1)
      Dr = copy(D1)
      Er = copy(E1)
      ir = 1:n
      if contr
         m == 0 &&  (ir = 1:0; return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr, n, 0, 0)
         _, _, _, nr, nfuc, niuc = sklf_right!(Ar, Er, Br, Cr; fast = fast, atol1 = atol1, atol2 = atol2, atol3 = atol1, rtol = rtol, withQ = false, withZ = false)
         nuc = nfuc+niuc
         if nuc > 0
            ir = 1:nr
            # save intermediary results
            A1 = Ar[ir,ir]
            E1 = Er[ir,ir]
            B1 = Br[ir,:]
            C1 = Cr[:,ir]
         else
            # restore original matrices
            Ar = copy(A1)
            Er = copy(E1)
            Br = copy(B1)
            Cr = copy(C1)
         end
      else
         nuc = 0
         nr = n
      end
      if obs
         p == 0 &&  (ir = 1:0; return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr, nuc, nr, 0)
         _, _, _, no, nfuo, niuo = sklf_left!(view(Ar,ir,ir), view(Er,ir,ir), view(Cr,:,ir), view(Br,ir,:);
                                            fast = fast, atol1 = atol1, atol2 = atol2, atol3 = atol1,
                                            rtol = rtol, withQ = false, withZ = false)
          nuo = nfuo+niuo
          if nuo > 0
             ir = ir[end-no+1:end]
             # save intermediary results
             A1 = Ar[ir,ir]
             E1 = Er[ir,ir]
             B1 = Br[ir,:]
             C1 = Cr[:,ir]
          else
             # restore saved matrices
             Ar[ir,ir] = A1
             Er[ir,ir] = E1
             Br[ir,:] = B1
             Cr[:,ir] = C1
          end
      else
          nuo = 0
      end
      if noseig
         rE, rA22  = _svdlikeAE!(view(Ar,ir,ir), view(Er,ir,ir), nothing, nothing, view(Br,ir,:), view(Cr,:,ir),
                     fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false)
         if rA22 > 0
            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,:))
            # 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; ir[rE+rA22+1:end]]
         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, nuc, nuo, rA22
      else
         return Ar[ir,ir], Er[ir,ir], Br[ir,:], Cr[:,ir], Dr, nuc, nuo, 0
      end
   end
end
"""
    lpsminreal(A, E, B, F, C, G, D, H; fast = true, atol1 = 0, atol2, rtol, contr = true, obs = true)
               -> (Ar, Er, Br, Fr, Cr, Gr, Dr, Hr, V, W, nuc, nuo)

Reduce the linearization `(A-λE,B-λF,C-λG,D-λH)` of a rational matrix to a reduced form
`(Ar-λEr,Br-λFr,Cr-λGr,Dr-λHr)` such that, for appropriate
invertible upper triangular matrices `V` and `W`,

                      -1                                     -1
     V'*((C-λG)*(λE-A)  *(B-λF) + D-λH)*W = (Cr-λGr)*(λEr-Ar)  *(Br-λFr) + Dr-λHr

with the least possible order `nr` of `Ar-λEr` if `contr = true` and `obs = true`.
Such a realization is called `strongly minimal` and satisfies:

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

     (2) rank[Ar-λEr; Cr-λGr] = nr for all finite and infinite λ (strong observability)

The achieved dimensional reductions to fulfill conditions (1) and (2) are
returned in `nuc` and `nuo`, respectively.

If `contr = true`, then the strong controllability condition (1) is enforced and `W` is an invertible upper triangular matrix or
`W = I` if `nuc = 0`.
If `contr = false`, then the strong controllability condition (1) is not enforced and `W = I`.

If `obs = true`, then the strong observability condition (2) is enforced and `V` is an invertible upper triangular matrix or
`V = I` if `nuo = 0`.
If `obs = false`, then the strong observability condition (2) is not enforced and `V = I`.

To enforce conditions (1) and (2), orthogonal similarity transformations are performed on
the matrices of the original linearization `(A-λE,B-λF,C-λG,D-λH)` to obtain a strongly minimal linearization
using structured pencil reduction algorithms [1]. The resulting realization `(Ar-λEr,Br-λFr,Cr-λGr,Dr-λHr)`
fulfills the strong controllability and strong observability conditions established in [2].

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, 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 matrices `A`, `B`, `C`, `D`, the absolute tolerance for the nonzero elements of `E`, `F`, `G`, `H`
and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`, `F`, `G`, `H`.

[1] F.M. Dopico, M.C. Quintana and P. Van Dooren, Linear system matrices of rational transfer functions,
to appear in "Realization and Model Reduction of Dynamical Systems", A Festschrift to honor the 70th birthday of Thanos Antoulas",
Springer-Verlag. [arXiv:1903.05016](https://arxiv.org/pdf/1903.05016.pdf)

[2] G. Verghese, Comments on ‘Properties of the system matrix of a generalized state-space system’,
Int. J. Control, Vol.31(5) (1980) 1007–1009.
"""
function lpsminreal(A::AbstractMatrix, E::AbstractMatrix, B::AbstractMatrix, F::AbstractMatrix, C::AbstractMatrix, G::AbstractMatrix,
                   D::AbstractMatrix, H::AbstractMatrix;
                   atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))),
                   rtol::Real =  (size(A,1)+2)*eps(real(float(one(real(eltype(A))))))*iszero(min(atol1,atol2)),
                   fast::Bool = true, contr::Bool = true, obs::Bool = true, noseig::Bool = true)

    n = LinearAlgebra.checksquare(A)
    (n,n) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
    p, m = size(D)
    (n,m) == size(B) || throw(DimensionMismatch("A, B and D must have compatible dimensions"))
    (p,n) == size(C) || throw(DimensionMismatch("A, C and D must have compatible dimensions"))
    (n,m) == size(F) || throw(DimensionMismatch("B and F must have the same dimensions"))
    (p,n) == size(G) || throw(DimensionMismatch("C and G must have the same dimensions"))
    (p,m) == size(H) || throw(DimensionMismatch("D and H must have the same dimensions"))
    T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D), eltype(E), eltype(F), eltype(G), eltype(H))
    T <: BlasFloat || (T = promote_type(Float64,T))

    A1 = copy_oftype(A,T)
    B1 = copy_oftype(B,T)
    C1 = copy_oftype(C,T)
    D1 = copy_oftype(D,T)
    E1 = copy_oftype(E,T)
    F1 = copy_oftype(F,T)
    G1 = copy_oftype(G,T)
    H1 = copy_oftype(H,T)

    n == 0 && (return A1, E1, B1, F1, C1, G1, D1, H1, I, I, 0, 0)
    # save system matrices
    Ar = copy(A1)
    Br = copy(B1)
    Cr = copy(C1)
    Dr = copy(D1)
    Er = copy(E1)
    Fr = copy(F1)
    Gr = copy(G1)
    Hr = copy(H1)
    ir = 1:n
    nr = n
    if contr
      m == 0 && (ir = 1:0; return A1[ir,ir], E1[ir,ir], B1[ir,:], F1[ir,:], C1[:,ir], G1[:,ir], D1, H1, I, I, n, 0)
      _, W, nr = sklf_right!(Ar, Er, Br, Fr, Cr, Gr, Dr, Hr; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false)
      nuc = n-nr
      if nuc > 0
         ir = 1:nr
         # save intermediary results
         A1 = copy(Ar[ir,ir])
         E1 = copy(Er[ir,ir])
         B1 = copy(Br[ir,:])
         F1 = copy(Fr[ir,:])
         C1 = copy(Cr[:,ir])
         G1 = copy(Gr[:,ir])
         D1 = copy(Dr)
         H1 = copy(Hr)
         W = UpperTriangular(W)
      else
         # restore original matrices
         Ar = copy(A1)
         Br = copy(B1)
         Cr = copy(C1)
         Dr = copy(D1)
         Er = copy(E1)
         Fr = copy(F1)
         Gr = copy(G1)
         Hr = copy(H1)
         W = I
      end
   else
      W = I
      nuc = 0
   end
   if obs
      p == 0 && (ir = 1:0; return A1[ir,ir], E1[ir,ir], B1[ir,:], F1[ir,:], C1[:,ir], G1[:,ir], D1, H1, I, I, nuc, nr)
      V, _, no = sklf_left!(view(Ar,ir,ir), view(Er,ir,ir), view(Cr,:,ir), view(Gr,:,ir), view(Br,ir,:), view(Fr,ir,:), Dr, Hr;
                                         fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false)
       nuo = nr-no
       if nuo > 0
          ir = ir[end-no+1:end]
          V = UpperTriangular(V)
       else
          # restore saved matrices
          Ar[ir,ir] = A1
          Er[ir,ir] = E1
          Br[ir,:] = B1
          Fr[ir,:] = F1
          Cr[:,ir] = C1
          Gr[:,ir] = G1
          Dr = copy(D1)
          Hr = copy(H1)
          V = I
         end
   else
      V = I
      nuo = 0
   end
   return Ar[ir,ir], Er[ir,ir], Br[ir,:], Fr[ir,:], Cr[:,ir], Gr[:,ir], Dr, Hr, V, W, nuc, nuo
end


text:
## Special cases
## linear, quadratic, cubic

## return function name for special case: identity, solve_linear, solve_quadratic, solve_cubic
function special_case(ps::Vector{S}, T) where {S}
    deg = length(ps) - 1  # [p0, p1, ..., pn]
    if deg == 1
        hasmethod(solve_linear, (Vector{S}, Type{T}))    && return solve_linear
    elseif deg == 2
        hasmethod(solve_quadratic, (Vector{S}, Type{T})) && return solve_quadratic
    elseif deg == 3
        hasmethod(solve_cubic, (Vector{S}, Type{T}))     && return solve_cubic
#    elseif deg == 5
#        hasmethod(solve_quintic, (Vector{S}, Type{T}))   && return solve_quintic
    end

    return identity
end

## Assume ps = [p_0, p_1, ..., p_n] with  p_n \neq 0, but possibly p_0 = 0

## Linear
check_linear(ps) = length(ps) == 2 || throw(DomainError())
## Over C
function solve_linear(ps::Vector{Complex{T}}, ::Type{Over.CC{S}}) where {T <: Real, S}
    check_linear(ps)
    Complex{S}[-ps[1] / ps[2]]
end
solve_linear(ps::Vector{T}, ::Type{Over.CC{S}}) where {T <: Real, S} = solve_linear(complex.(ps,zeros(T,2)), Over.CC{S})

## Over R
function solve_linear(ps::Vector{T}, ::Type{Over.RR{S}}) where {T <: Real, S}
    check_linear(ps)
    S[-ps[1] / ps[2]]
end

## Over Q
function solve_linear(ps::Vector{T}, ::Type{Over.QQ{S}}) where {T <: Integer, S}
    check_linear(ps)
    Rational{S}[-ps[1] // ps[2]]
end
function solve_linear(ps::Vector{Rational{T}}, ::Type{Over.QQ{S}}) where {T <: Integer, S}
    check_linear(ps)
    Rational{S}[-ps[1] // ps[2]]
end

## Over Z
function solve_linear(ps::Vector{T}, ::Type{Over.ZZ{S}}) where {T <: Integer, S}
    check_linear(ps)
    q,r = divrem(-ps[1], ps[2])
    if r == 0
        S[q]
    else
        zeros(S, 0)
    end
end

## Quadratic
check_quadratic(ps) = length(ps) == 3 || throw(DomainError())

## Kahan quadratic equation with fma
##  https://people.eecs.berkeley.edu/~wkahan/Qdrtcs.pdf

## solve ax^2 - 2bx + c
function qdrtc(a::T, b::T, c::T) where {T <: Real}
    # z1, z2 roots of ax^2 - 2bx + c
    d = discr(a,b,c)  # (b^2 - a*c), as 2 removes 4

    if d <= 0
        r = b/a  # real
        s = sqrt(-d)/a #imag
        return (r,s,r,-s)
    else
        r = sqrt(d) * (sign(b) + iszero(b)) + b
        return (r/a, zero(T), c/r, zero(T))
    end
end

## more work could be done here.
function discr(a::T,b::T,c::T) where {T}
    pie = 3.0 # depends on 53 or 64 bit...
    d = b*b - a*c
    e = b*b + a*c

    pie*abs(d) > e && return d

    p = b*b
    dp = muladd(b,b,-p)
    q = a*c
    dq = muladd(a,c,-q)

    (p-q) + (dp - dq)
end

function solve_quadratic(ps::Vector{T}, U::Type{Over.CC{S}}) where {T <: Real, S}
    check_quadratic(ps)
    iszero(ps[1]) && return vcat(zero(S), solve_linear(ps[2:end], U))

    r1,i1,r2,i2 = qdrtc(ps[3], -(0.5) * ps[2], ps[1])
    Complex{S}[complex(r1, i1), complex(r2, i2)]
end

function solve_quadratic(ps::Vector{T}, U::Type{Over.RR{S}}) where {T <: Real, S}
    check_quadratic(ps)
    iszero(ps[1]) && return vcat(zero(S), solve_linear(ps[2:end], U))

    r1,i1,r2,i2 = qdrtc(ps[3], -(0.5) * ps[2], ps[1])
    if iszero(i1)
        S[r1, r2]
    else
        zeros(S,0)
    end
end

## Cubic
check_cubic(ps) = length(ps) == 4 || throw(DomainError())

# C -> C
function solve_cubic(ps::Vector{Complex{T}}, U::Type{Over.CC{S}}) where {T <: Real, S}
    check_cubic(ps)
    iszero(ps[1]) && return solve_quadratic(ps[2:end], U)

    [PolynomialRoots.solve_cubic_eq(ps)...]
end

# R -> C XXX Can be faster
function solve_cubic(ps::Vector{T}, U::Type{Over.CC{S}}) where {T <: Real, S}
    check_cubic(ps)
    iszero(ps[1]) && return solve_quadratic(ps[2:end], U)

    [PolynomialRoots.solve_cubic_eq(complex.(ps, zeros(T,4)))...]
end

# R -> R
function solve_cubic(ps::Vector{T}, U::Type{Over.RR{S}}) where {T <: Real, S}
    check_cubic(ps)
    iszero(ps[1]) && return solve_quadratic(ps[2:end], U)

    ## https://en.wikipedia.org/wiki/Cubic_function#Algebraic_solution

    d,c,b,a = ps
    Delta = 18 * prod(ps) - 4b^3 * d + b^2 * c^2 - 4a*c^3 - 27 * a^2 * d^2
    Delta0 = b^2 - 3a*c
    Delta1 = 2b^3 - 9*a*b*c + 27*a^2 * d

    if Delta == 0 # one or two (with a double root)

        if iszero(Delta0)
            return ones(S,3) * (-b/(3a))
        else
            return S[(9a*d - b*c)/(2Delta0), (9a*d - b*c)/(2Delta0), (4a*b*c - 9a^2*d - b^3)/(a * Delta0)]
        end

    elseif Delta < 0 # one real root

        if iszero(Delta0)
            C = cbrt((0.5) * (2*sign(Delta1)*abs(Delta1)))

        else
            Delta2 = sqrt(-27*a^2*Delta)
            C = cbrt((0.5)* ( Delta1 + Delta2))

        end
        return S[-1/(3a)*(b+C+Delta0/C)]

    elseif Delta > 0 # 3 distinct real roots

        p, q = (3a*c  - b^2)/(3a^2), (2b^3 - 9a*b*c + 27a^2*d) / (27a^3)

        if iszero(p)
            return ones(S,3) * (-b/(3a))
        else
            return S[2*sqrt(-p/3)*cos(1/3 * acos((3q)/(2p) * sqrt(-3/p)) - 2pi*k/3) for k in 0:2] - b/(3a)
        end

    end

end

check_quintic(ps::Vector) = length(ps) == 6 || throw(DomainError())
# C -> C
function solve_quintic(ps::Vector{Complex{T}}, U::Type{Over.CC{S}}) where {T <: Real, S}
    check_quintic(ps)

    PolynomialRoots.roots5(ps)
end

ext: jl


text:
mutable struct MatrixNetwork{T}
    n::Int64 # number of columns/rows
    rp::Vector{Int64} # row pointers
    ci::Vector{Int64} # column indices
    vals::Vector{T} # corresponding values
end

function MatrixNetwork(A::SparseMatrixCSC{T,Int64}) where T
    At = copy(A')
    return MatrixNetwork(size(At,2),At.colptr,At.rowval,At.nzval)
end

function _second(x)
    return x[2]
end

MatrixNetwork(edges::Vector{Tuple{Int,Int}}, n::Int) =
    MatrixNetwork(map(first,edges),map(_second,edges), n)

MatrixNetwork(ei::Vector{Int64},ej::Vector{Int64}) =
    MatrixNetwork(ei,ej,max(maximum(ei),maximum(ej)))

function MatrixNetwork(ei::Vector{Int64},ej::Vector{Int64},n::Int64)
    At = sparse(ej,ei,true,n,n);
    return MatrixNetwork(size(At,2),At.colptr,At.rowval,At.nzval)
end

function _matrix_network_direct(A::SparseMatrixCSC{T,Int64}) where T
    return MatrixNetwork(size(A,2),A.colptr,A.rowval,A.nzval)
end

function _matrix_network_direct(A::SparseMatrixCSC{T,Int64},v) where T
    nzval = ones(typeof(v),length(A.nzval))
    return MatrixNetwork(size(A,2),A.colptr,A.rowval,nzval)
end

import SparseArrays.sparse, Base.size, Base.*, LinearAlgebra.mul!, Base.convert, Base.adjoint, Base.copy

"""
Return back an adjacency matrix representation
of the transpose. This requires no work.
"""
function sparse_transpose(A::MatrixNetwork{T}) where T
    return SparseMatrixCSC(A.n,A.n,A.rp,A.ci,A.vals)
end

adjoint(A::MatrixNetwork{T}) where T = Adjoint{Any,MatrixNetwork{T}}(A)

"""
Return back an adjacency matrix representation
of the current MatrixNetwork
"""
function sparse(A::MatrixNetwork{T}) where T
    return copy(sparse_transpose(A)')
end

function size(A::MatrixNetwork)
    return (A.n,A.n)
end

import Base.ndims
ndims(op::MatrixNetwork) = 2

function size(A::MatrixNetwork, dim::Integer)
    if dim == 1 || dim == 2
        return A.n
    elseif dim > 2
        return 1
    else
        throw(DomainError(dim))
    end
end

*(M::MatrixNetwork{T}, b::AbstractVector{S}) where {T,S} =
    mul!(Array{promote_type(T,S)}(undef,size(M,2)),sparse_transpose(M)',b)

mul!(output,M::MatrixNetwork,b) = mul!(output,sparse_transpose(M)',b)

*(M::Adjoint{<:Any,<:MatrixNetwork{T}}, b::AbstractVector{S}) where {T,S} =
    mul!(Array{promote_type(T,S)}(undef,size(M.parent,2)),sparse_transpose(M.parent),b)

mul!(output,M::Adjoint{<:Any,<:MatrixNetwork},b) = mul!(output,sparse_transpose(M.parent),b)

sparse(M::Adjoint{<:Any,<:MatrixNetwork}) = sparse_transpose(M.parent)
copy(M::Adjoint{<:Any,<:MatrixNetwork}) = sparse_transpose(M.parent)

"""
`is_empty`
==========

Return true if the graph is the empty graph and
false otherwise.

Functions
---------
-`is_empty(A::MatrixNetwork)`

Example
-------
~~~~
is_empty(MatrixNetwork(Int[],Int[],0))
is_empty(erdos_renyi_undirected(0,0))
~~~~
"""
is_empty(A::MatrixNetwork) = size(A,1) == 0

"""
`is_undirected`
===============

Check the matrix associated with a matrix network
for symmetry.

Input
-----
- `A`: a matrix network

Returns
-------
- `bool` with true indicating the network is undirected
    and the matrix is symmetric
"""
function is_undirected end

function is_undirected(A::MatrixNetwork)
   M = sparse_transpose(A)
   return issymmetric(M)
end

function is_undirected(A::SparseMatrixCSC)
   return issymmetric(A)
end

"""
`is_connected`
==============

Check the matrix associated with a matrix network
for (strong) connectivity

Usage
-----
- `is_connected(A)`

Input
-----
- `A`: a `MatrixNetwork` or `SparseMatrixCSC` class

Returns
-------
- `bool` with true indicating the matrix is strongly connected
and false indicating
"""
function is_connected end

function is_connected(A::Union{MatrixNetwork,SparseMatrixCSC})
    # this is equivalent to maximum with a default value of 0
    return mapreduce(identity, max, strong_components_map(A); init=0) == 1
end

"""
`empty_graph`
================

Returns an empty graph with n vertices and zero edges

Functions
---------
* `A = empty_graph(n)` generates an empty graph on n edges.

Example
-------
~~~~
is_connected(empty_graph(0))
is_connected(empty_graph(1))
~~~~
"""
function empty_graph end

function empty_graph(n::Integer=0)
    return MatrixNetwork(n,ones(Int64,n+1),Array{Int64}(undef, 0),Array{Float64}(undef, 0))
end

"""
`random_edge`
=============
Identify a random edge of a matrix network or sparse matrix.

Functions
---------
* `random_edge(A::MatrixNetwork) -> (ei,ej,ind)`
   gets a random edge/non-zero from the matrix

* `random_edge(A::SparseMatrixCSC) -> (ei,ej,ind)`
   gets a random non-zero from the matrix

Example
-------
~~~~
G = lollipop_graph(5,3)
# count the number of edges we randomly see between the regions
C = Dict{Symbol,Int}()
M = zeros(8,8)
for i=1:1000000
  ei,ej = MatrixNetwork.random_edge(G)
  M[ei,ej] += 1
  if 1 <= ei <= 5 && 1 <= ej <= 5
    C[:Stem] = get(C, :Stem, 0) + 1
  elseif 6 <= ei <= 10 && 6 <= ej <= 10
    C[:Pop] = get(C, :Pop, 0) + 1
  else
    C[:Bridge] = get(C, :Bridge, 0) + 1
  end
end
# 4 edges in stem, 3 edges in pop, 1 edge in bridge
@show C
@show M
~~~~
"""
function random_edge(A::MatrixNetwork)
    ind = rand(1:length(A.ci)) # the index
    ej = A.ci[ind]
    ei = searchsortedlast(A.rp, ind) # uses binary search for efficiency
    @assert ei <= A.n
    return (ei,ej,ind)
end
function random_edge(A::SparseMatrixCSC)
    ind = rand(1:length(A.rowval)) # the index
    ei = A.rowval[ind]
    ej = searchsortedlast(A.colptr, ind) # uses binary search for efficiency
    @assert ej <= A.n
    return (ei,ej,ind)
end

"""
    undirected_edges(A) -> srcs,dsts

Produce lists just for the undirected edges. This assumes you want only
those edges where (target >= source).
"""
function undirected_edges(A::MatrixNetwork)
    ei = Vector{Int64}()
    ej = Vector{Int64}()
    sizehint!(ei, div(A.rp[A.n+1],2))
    sizehint!(ej, div(A.rp[A.n+1],2))
    for i=1:A.n
        for nzi=A.rp[i]:A.rp[i+1]-1
            j = A.ci[nzi]
            if j >= i
                push!(ei, i)
                push!(ej, j)
            end
        end
    end
    return ei,ej
end

"""
    directed_edges(A) -> srcs,dsts

Produce lists just for all edges of the graph, including both sides for
undirected edges. This is essentially the same as findnz for a sparse matrix,
optimized not to return the values.
"""
function directed_edges(A::MatrixNetwork)
    ei = Vector{Int64}()
    ej = Vector{Int64}()
    sizehint!(ei, A.rp[A.n+1])
    sizehint!(ej, A.rp[A.n+1])
    for i=1:A.n
        for nzi=A.rp[i]:A.rp[i+1]-1
            j = A.ci[nzi]
            push!(ei, i)
            push!(ej, j)
        end
    end
    return ei,ej
end


text:
function sector_performance(; client = GLOBAL[], datatype::Union{String, Nothing}=nothing, parser = "default")
    @argcheck in(datatype, ["json", "csv", nothing])
    params = Dict(
        "function"=>"SECTOR",
        "datatype"=>isnothing(datatype) ? "csv" : datatype,
        "apikey"=>key(client)
    )
    uri = _build_uri(client.scheme, client.host, "query", params)
    data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
    p = _parser(parser, datatype)
    return p(data)
end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))

using GLMakie, InteractiveDynamics, DynamicalSystems, OrdinaryDiffEq, LinearAlgebra
diffeq = (alg = Tsit5(), dtmax = 0.02, abstol = 1e-6, reltol = 1e-6)

# Modified version of Rozenweig-MacArthur model for predator-prey
# See https://arxiv.org/abs/2101.12107
function hassan_alkhayuon(u, p, t) # From talk of Hassan Alkhayuon
    r, c, μ, ν, α, β, χ, δ = p
    N, P = u
    common = α*N*P/(β+N)
    dN = r*N*(1 - (c/r)*N)*((N-μ)/(N+ν)) - common
    dP = χ*common - δ*P
    return SVector(dN, dP)
end

u0 = SVector(8.0, 0.01)
r1, r2 = 1.7, 2.5
# r, c, μ, ν, α, β, χ, δ = p
p = [r1, 0.19, 0.03, 0.003, 800, 1.5, 0.004, 2.2]
ds = ContinuousDynamicalSystem(hassan_alkhayuon, u0, p)

# I will re-use code from InteractiveDynamics.jl, so I'll be defining some necessary
# quantities that may not make perfect sense without looking at the source code.
# Nevertheless the source code is quite clear and is just used to initialize some
# interactive observables and stuff.
tail = 500                  # length of plotted trajectory (in units of `dt`)
lims = ((0, 15), (0, 0.04)) # limits of state space
color = to_color("#1e1528") # trajectory color

figure = Figure(resolution = (1000, 800))
pinteg = DynamicalSystems.parallel_integrator(ds, [u0]; diffeq...)
obs, finalpoints = InteractiveDynamics.init_trajectory_observables(1, pinteg, tail, SVector(1,2), identity)

main = figure[1,1] = InteractiveDynamics.init_main_trajectory_plot(
    ds, figure, SVector(1,2), lims, pinteg, [color], obs, NamedTuple(), finalpoints, 1.0
)

# Plot the stability basin and attractor
active_color = to_color("#d7800a")
passive_color = RGBAf0(0.3, 0.4, 0.4, 0.5)
r1_color = Observable(active_color)
r2_color = Observable(passive_color)

function plot_limitcycle_basin!(main, r, color)
    set_parameter!(ds, 1, r)
    # plot unstable manifold of "e3" fixed point
    zz = SVector(p[3], 0.0)
    J = ds.jacobian(zz, p, 0)
    ev = eigvecs(J)
    # system in reverse time: just negative `f`.
    pprev = (u, p, t) -> -hassan_alkhayuon(u, p, t)
    dsrev = ContinuousDynamicalSystem(pprev, SVector(1.0, 1.0), p)
    e = ev[:, 1] # eigenvector correpsonding to positive eigenvalue
    tr = trajectory(dsrev, 3.5, zz .+ 1e-3 .* e; dtmax = 0.01)
    lines!(main, columns(tr)...; linewidth = 2.0, color = color, linestyle = :dash)
    # plot limit cycle
    tr = trajectory(ds, 20.0; dt = 0.01, Ttr = 100.0)
    lines!(main, columns(tr)...; linewidth = 4.0, color = color)
    return
end

plot_limitcycle_basin!(main, r2, r2_color)
plot_limitcycle_basin!(main, r1, r1_color)

# add controls for start/stop and toggle the r value below the figure
figure[2, 1] = controllayout = GridLayout(tellwidth = false)
runbutton = controllayout[1, 1] = Button(figure; label = "run")
rbutton = controllayout[1, 2] = Button(figure; label = "toggle r")
run, rtoggle = runbutton.clicks, rbutton.clicks

# Set up interactive part of the application: start/stop time evolution
isrunning = Observable(false)
on(run) do clicks; isrunning[] = !isrunning[]; end
on(run) do clicks
    @async while isrunning[]
        DynamicalSystems.step!(pinteg)
        ob = obs[1]
        last_state = DynamicalSystems.get_state(pinteg, 1)
        ob[] = push!(ob[], last_state) # push and trigger update with `=`
        finalpoints[] = [x[][end] for x in obs]
        sleep(0.00001) # use sleep to slow down animation
        isopen(figure.scene) || break # crucial, ensures computations stop if closed window
    end
end

# Set up interactive part of the application: toggle the r value
on(rtoggle) do clicks
    r = (r1, r2)[mod(clicks, 2)+1]
    set_parameter!(ds, 1, r)
    u_modified!(pinteg, true)
    x, y = r1_color[], r2_color[]
    r1_color[], r2_color[] = y, x
end

display(figure)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Blade.jl defines the Blade type and core methods
"""

# --- Exports

# Types
export Blade

# Functions
export blade_atol

# --- Type definitions

import LinearAlgebra
using LinearAlgebra: det, diag, dot, qr

"""
    struct Blade{T<:AbstractFloat} <: AbstractBlade{T}

Blade represented with the floating-point precision of type `T`. The norm and
orientation of a Blade are encoded by its `volume`. The norm of a Blade is
equal to `abs(volume)` and the orientation of a Blade relative to its `basis`
is equal to `sign(volume)`.

Notes
=====
* The grade of a Blade type is greater than 0 and less than the dimension of
  the space that the blade is embedded in
"""
struct Blade{T<:AbstractFloat} <: AbstractBlade{T}
    #=
     Fields
      ------
      * `dim`: the dimension of the space that the blade is embedded in

      * `grade`: the dimension of the space spanned by the blade

      * `basis`: an orthonormal basis for the space spanned by the blade. Note
        that the order of the columns in `basis` defines the orientation for
        the unit blade represented by `basis`.

      * `volume`: the signed-norm (hypervolume) of the blade. The sign
         of `volume` indicates the orientation of the blade relative to the
         unit blade represented by `basis`. It is positive when the blade has
         the same orientation as `basis` and negative when the blade has the
         opposite orientation.
    =#
    dim::Int
    grade::Int
    basis::Matrix{T}
    volume::T

    """
    Construct a Blade for a geometric algebra in `dim` dimensions having the
    specified data field values. If the absolute value of `volume` is less than
    `atol`, a Scalar representing zero is returned. If the `dim` and `grade`
    are equal, a Pseudoscalar is returned.

    When `enforce_constraints` is true, constraints are enforced. When
    `copy_basis` is true, the basis of the new Blade is a copy of `basis`;
    otherwise, the basis of the new Blade is a reference to `basis`.

    Note: this inner constructor is intended primarily for use by outer
    constructors to enforce constraints.
    """
    function Blade{T}(dim::Int,
                      grade::Int,
                      basis::Matrix{T},
                      volume::Real;
                      atol::Real=blade_atol(T),
                      enforce_constraints::Bool=true,
                      copy_basis::Bool=true) where {T<:AbstractFloat}

        # --- Enforce constraints

        if enforce_constraints
            if dim < 0
                throw(ArgumentError("`dim` must be nonnegative"))
            end

            if grade < 0
                throw(ArgumentError("`grade` must be nonnegative"))
            end

            basis_size = size(basis)

            if dim != basis_size[1]
                message = string("`dim` not equal to the number of ",
                                 "rows in `basis`")
                throw(DimensionMismatch(message))
            end

            if grade != basis_size[2]
                message = string("`grade` not equal to the number of ",
                                 "columns in `basis`")
                throw(DimensionMismatch(message))
            end

            if dim > 0 && grade > 0
                for i in 1:grade
                    if LinearAlgebra.norm(basis[:, i]) ≉ 1
                        message = "columns of `basis` not normalized"
                        throw(ArgumentError(message))
                    end
                end
            end
        end

        # --- Handle edge cases

        # `volume` is effectively zero
        if abs(volume) < atol
            return zero(Blade{T})
        end

        # Return zero if dim or grade is equal to 0
        if dim == 0 || grade == 0
            return zero(Blade{T})
        end

        # Return a Pseudoscalar if the grade of the blade is equal to the
        # dimension of the embedding space.
        if grade == dim
            if det(basis) == 0
                return zero(Blade{T})
            elseif det(basis) > 0
                return Pseudoscalar{T}(dim, volume)
            else
                return Pseudoscalar{T}(dim, -volume)
            end
        end

        # Return new Blade
        copy_basis ?
            new(dim, grade, copy(basis), volume) :
            new(dim, grade, basis, volume)
    end
end

"""
    Blade(vectors::Array{T};
          volume::Union{Real, Nothing}=nothing,
          atol::Real=blade_atol(T)) where {T<:AbstractFloat}

    Blade(vectors::Array{<:Integer};
          volume::Union{Real, Nothing}=nothing,
          atol::Real=blade_atol(Float64))

Construct a Blade from a collection of vectors stored as (1) the columns of a
matrix or (2) a single vector. If the norm of the blade is less than `atol`, a
Scalar representing zero is returned. If the grade of the blade is equal to
the dimension of the space that the blade is embedded in, a Pseudoscalar is
returned.

By default, `vectors` determines the volume of the blade. However, if `volume`
is specified, `vectors` is only used to define the subspace (including
orientation) represented by the blade.

Notes
=====

Orientation
-----------

* _`vectors` contains more than one vector_. When `volume` is positive, the
  orientation of the blade is the same as the orientation of the outer product
  of the columns of `vectors` (taken in order). When `volume` is negative, the
  orientation of the blade is the opposite of the orientation implied by the
  `vectors`.

* _`vectors` contains a single vector `v`_. When `volume` is positive, the
  orientation of the blade is the same as the direction of `v`. When `volume`
  is negative, the orientation of the blade is the opposite of the direction
  of `v`.

Precision
---------

When the precision is not specified, the following rules are applied to set
the precision of the Blade.

* If `vectors` is an Array of floating-point values, the precision of the
  Blade is inferred from the precision of the elements of `vector`.

* If `vectors` is an Array of integers, the precision of the Blade
  defaults to `Float64`.
"""
function Blade{T}(vectors::Matrix{<:Real};
                  volume::Union{Real, Nothing}=nothing,
                  atol::Real=blade_atol(T)) where {T<:AbstractFloat}

    # --- Handle edge cases

    # `volume` is effectively zero
    if volume != nothing && abs(volume) < atol
        return zero(Blade{T})
    end

    # number of vectors == 0
    dims = size(vectors)
    if dims[1] == 0 || dims[2] == 0
        return zero(Blade{T})
    end

    # number of vectors > dimension of column space
    if dims[2] > dims[1]
        if dims[1] == 1
            # `vectors` is a single row vector, so convert it to a column
            # vector and call constructor for single column vector.
            return Blade{T}(reshape(vectors, dims[2]))
        else
            return zero(Blade{T})
        end
    end

    # --- Construct Blade

    # Convert `vectors` to have type `T`
    vectors = convert(Matrix{T}, vectors)

    # Compute QR factorization and signed norm
    F = qr(vectors)
    signed_norm = prod(diag(F.R))

    # Return zero if norm is below tolerance
    if abs(signed_norm) < atol
        return zero(Blade{T})
    end

    # Compute orthonormal basis for subspace
    basis = Matrix(F.Q)

    # Compute volume
    volume = isnothing(volume) ? signed_norm : volume * sign(signed_norm)

    # Return new Blade
    Blade{T}(dims[1], dims[2], basis, volume,
             atol=atol, enforce_constraints=false, copy_basis=false)
end

function Blade{T}(v::Vector{<:Real};
                  volume::Union{Real, Nothing}=nothing,
                  atol::Real=blade_atol(T)) where {T<:AbstractFloat}

    # --- Handle edge cases

    # `volume` is effectively zero
    if volume != nothing && abs(volume) < atol
        return zero(Blade{T})
    end

    # length(v) == 0
    if length(v) == 0
        return zero(Blade{T})
    end

    # --- Construct Blade

    # Convert `v` to have type `T` and compute norm
    v = convert(Vector{T}, v)
    norm_v = LinearAlgebra.norm(v)

    # Return zero if norm is below tolerance
    if abs(norm_v) < atol
        return zero(Blade{T})
    end

    # Compute basis as a column vector
    basis = reshape(v, length(v), 1) / norm_v

    # Compute volume
    volume = isnothing(volume) ? norm_v : volume

    # Return new Blade
    Blade{T}(length(v), 1, basis, volume, atol=atol,
             enforce_constraints=false, copy_basis=false)
end

Blade(vectors::Array{T};
      volume::Union{Real, Nothing}=nothing,
      atol::Real=blade_atol(T)) where {T<:AbstractFloat} =
    Blade{T}(vectors, volume=volume, atol=atol)

Blade(vectors::Array{<:Integer};
      volume::Union{Real, Nothing}=nothing,
      atol::Real=blade_atol(Float64)) =
    Blade(convert(Array{Float64}, vectors), volume=volume, atol=atol)

"""
    Blade(B::Blade;
          volume::Real=volume(B),
          atol::Real=blade_atol(typeof(volume(B))),
          copy_basis=false)

    Blade{T}(B::Blade;
          volume::Real=volume(B),
          atol::Real=blade_atol(T),
          copy_basis=false) where {T<: AbstractFloat}

Copy constructors. Construct a Blade representing the same space as `B` having
a specified oriented volume relative to `B`. A Scalar representing zero is
returned if the absolute value of `volume` is less than `atol`.

When `copy_basis` is true, the `basis` of the new Blade is a copy of the
`basis` of the original Blade; otherwise, the `basis` of the new Blade is
reference to the `basis` of the original Blade.
"""
Blade{T}(B::Blade;
         volume::Real=volume(B),
         atol::Real=blade_atol(T),
         copy_basis::Bool=false) where {T<:AbstractFloat} =

    Blade{T}(dim(B), grade(B), convert(Matrix{T}, basis(B)), volume,
             atol=atol, enforce_constraints=false, copy_basis=copy_basis)

Blade(B::Blade;
      volume::Real=B.volume,
      atol::Real=blade_atol(typeof(B.volume)),
      copy_basis=false) =
    Blade{typeof(B.volume)}(B, volume=volume, atol=atol, copy_basis=copy_basis)

"""
    Blade(x::Real)
    Blade(x::Scalar)

Convenience constructors that return a Scalar with value `x`.
"""
Blade(x::Real) = Scalar(x)
Blade(x::AbstractScalar) = Scalar(value(x))

# --- Method definitions for AbstractMultivector interface functions

dim(B::Blade) = B.dim

-(B::Blade) = Blade(B, volume=-volume(B), copy_basis=false)

Base.reverse(B::Blade) =
    mod(grade(B), 4) < 2 ?  B : Blade(B, volume=-volume(B), copy_basis=false)

function dual(B::Blade)
    # --- Extend basis(B) to an orthonormal basis for entire space.

    F = qr(basis(B))

    # --- Compute volume of dual

    # Account for orientation of Q relative to orientation of I formed from
    # standard basis
    dual_volume = volume(B) * sign(det(F.Q))

    # Account for orientation of first grade(B) columns of Q relative to
    # orientation of basis(B)
    if prod(diag(F.R)) < 0
        dual_volume = -dual_volume
    end

    # Account for sign of I^{-1} relative to I
    if mod(dim(B), 4) >= 2
        dual_volume = -dual_volume
    end

    # Account for reversals required to eliminate B
    if mod(grade(B), 4) >= 2
        dual_volume = -dual_volume
    end

    Blade{typeof(volume(B))}(dim(B), dim(B) - grade(B),
                             F.Q[:, grade(B) + 1:end], dual_volume,
                             enforce_constraints=false,
                             copy_basis=false)
end

# --- Method definitions for AbstractBlade interface functions

grade(B::Blade) = B.grade

basis(B::Blade) = B.basis

volume(B::Blade) = B.volume

inv(B::Blade) =
    mod(grade(B), 4) < 2 ?
        Blade(B, volume=1 / volume(B), copy_basis=false) :
        Blade(B, volume=-1 / volume(B), copy_basis=false)

# --- Comparison methods

==(B::AbstractBlade, C::AbstractBlade) =
    dim(B) == dim(C) && grade(B) == grade(C) &&
    volume(B) == volume(C) && basis(B) == basis(C)

function isapprox(B::Blade{T1}, C::Blade{T2};
  atol::Real=0,
  rtol::Real=atol>0 ? 0 : max(√eps(T1), √eps(T2))) where {T1<:AbstractFloat,
                                                          T2<:AbstractFloat}
    # Check dim, grade, and norm are equal
    if dim(B) != dim(C) || grade(B) != grade(C) ||
        !isapprox(norm(B), norm(C), atol=atol, rtol=rtol)

        return false
    end

    # Check that B and C represent the same space
    projection = det(transpose(basis(B)) * basis(C))
    if ≉(abs(projection), 1, atol=atol, rtol=rtol)
        return false
    end

    # Check that B and C have the same orientation
    return sign(B) * sign(C) == sign(projection)
end

function isapprox(B::Blade{T1}, C::Vector{T2};
  atol::Real=0,
  rtol::Real=atol>0 ? 0 : max(√eps(T1), √eps(T2))) where {T1<:AbstractFloat,
                                                          T2<:AbstractFloat}
    C_blade = Blade(C)
    return isapprox(B, C_blade, atol=atol, rtol=rtol)
end

function isapprox(B::Vector{T1}, C::Blade{T2};
  atol::Real=0,
  rtol::Real=atol>0 ? 0 : max(√eps(T1), √eps(T2))) where {T1<:AbstractFloat,
                                                          T2<:AbstractFloat}
    return isapprox(C, B, atol=atol, rtol=rtol)
end

# --- Utility methods

convert(::Type{T}, B::Blade) where {T<:AbstractFloat} =
    T == typeof(volume(B)) ? B : Blade{T}(B)

#=
 TODO: review numerical error in factorizations to see if a different
       tolerance would be better.
=#
"""
    blade_atol(::Type{T}) where {T<:AbstractFloat}

Return the minimum value of a nonzero blade's norm.
"""
blade_atol(::Type{T}) where {T<:AbstractFloat} = 100*eps(T)


text:
using StaticArrays

export FreeParameter, @FreeParameter
export extract_free_parameters

export ExcludeTypes, IncludeTypes
export EntireStruct, FreeParametersOnly

const FPTYPES = Real

"""
    FreeParameter{T,P,B} <: FPTYPES

Free parameter of type `T` and prior `P`.
"""
mutable struct FreeParameter{T,P,B} <: FPTYPES
  "Value used in the model"
  val::T
  "prior distribution"
  prior::P
  "bounds on value used"
  bounds::B
  function FreeParameter(val::T,
                         prior::P=nothing,
                         bounds::B=nothing
                         ) where {T,P,B}
    return new{T,P,B}(val, prior, bounds)
  end
end

#####
##### Extracting free parameters
#####

"""
    extract_free_parameters!(fp::Vector{FreeParameter}, D::Dict)

Extract a vector of `FreeParameter`'s from the dict `D`.
"""
extract_free_parameters(D::Dict) = FreeParameter[v for v in values(D) if typeof(v) <: FreeParameter]

#####
##### Annotating free parameters
#####

macro FreeParameter(s, prior=nothing, bounds=nothing)
  e = :($(s) = FreeParameter($(s), $(prior), $(bounds)))
  return :($(esc(e)))
end

"""

"""
abstract type AbstractFilterTypes end
abstract type AbstractIncludeTypes <: AbstractFilterTypes end
abstract type AbstractExcludeTypes <: AbstractFilterTypes end

"""
    IncludeTypes <: AbstractIncludeTypes

Types to include in
 - converting a struct to dict
 - exporting a dict
 - importing a dict
"""
struct IncludeTypes <: AbstractIncludeTypes
  include_types
  stop_recursion_types
end
IncludeTypes(include_types) = IncludeTypes(include_types, (FreeParameter,))

"""
    FreeParametersOnly

Convenience constructor for including only free parameters
"""
FreeParametersOnly() = IncludeTypes((FreeParameter,), (FreeParameter,))

"""
    ExcludeTypes <: AbstractExcludeTypes

Types to exclude in
 - converting a struct to dict
 - exporting a dict
 - importing a dict
"""
struct ExcludeTypes <: AbstractExcludeTypes
  exlude_types
  stop_recursion_types
end
ExcludeTypes(exlude_types) = ExcludeTypes(exlude_types, (FreeParameter,))

"""
    EntireStruct

Convenience constructor for including entire struct/dict etc.
"""
EntireStruct() = ExcludeTypes((Any,), (FreeParameter,))

struct leaf_triggers{B,AFT}
  ft::AFT
  s
  prop
  custom::B                  # C1
  emptyfieldnames::B         # C2
  isa_stop_recursion_type::B # C3
  isa_SArray::B              # C4
  isa_isbits::B              # C5
end

triggers(t::Tuple) = any(t)

function aux_leaf(lt::leaf_triggers{Bool,IncludeTypes})
  return any([lt.prop isa t for t in lt.ft.stop_recursion_types]) || any([lt.prop isa t for t in lt.ft.include_types])
end
function aux_leaf(lt::leaf_triggers{Bool,ExcludeTypes})
  return any([lt.prop isa t for t in lt.ft.stop_recursion_types]) || isbits(lt.prop)
end

triggers(lt::leaf_triggers) = triggers(lt, lt.ft)

triggers(lt::leaf_triggers, ft::IncludeTypes) = any([any([lt.prop isa t for t in ft.include_types]) && lt.emptyfieldnames,
                                                     lt.isa_stop_recursion_type
                                                     ])
triggers(lt::leaf_triggers, ft::ExcludeTypes) = any([isbits(lt.prop) && lt.emptyfieldnames,
                                                     lt.isa_stop_recursion_type,
                                                     lt.isa_SArray
                                                     ])

function __is_leaf(s, prop, ft::IncludeTypes)
  C1 = any([prop isa t for t in ft.include_types])
  C2 = isempty(fieldnames(typeof(prop)))
  C3 = any([prop isa t for t in ft.stop_recursion_types])
  C4 = (prop isa SArray)
  C5 = isbits(prop)
  return leaf_triggers(ft, s, prop, C1, C2, C3, C4, C5)
end
function __is_leaf(s, prop, ft::ExcludeTypes)
  C1 = isbits(prop)
  C2 = isempty(fieldnames(typeof(prop)))
  C3 = any([prop isa t for t in ft.stop_recursion_types])
  C4 = (prop isa SArray)
  C5 = isbits(prop)
  return leaf_triggers(ft, s, prop, C1, C2, C3, C4, C5)
end

function get_val_from_var(var::AbstractString)
  e = Meta.parse(var)
  if isbits(e) # var = "2.0" (for example)
    val = e
  else # var = "FreeParameter(...)" (for example)
    val = eval(e.args[1].args[1])(e.args[2:end]...)
  end
  return val
end

"""
    get_val(x)
Gets value of `FreeParameter`,
otherwise return identity.
"""
function get_val end
get_val(x) = x
get_val(x::FreeParameter) = x.val


text:
mutable struct CollectGemsDirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Gem}}
    agent_pos::CartesianIndex{2}
    agent_dir::AbstractDirection
    reward::T
    rng::R
    num_gem_init::Int
    num_gem_current::Int
    gem_reward::T
    gem_pos_init::Vector{CartesianIndex{2}}
    done::Bool
end

@generate_getters(CollectGemsDirected)
@generate_setters(CollectGemsDirected)

mutable struct CollectGemsUndirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Gem}}
    agent_pos::CartesianIndex{2}
    reward::T
    rng::R
    num_gem_init::Int
    num_gem_current::Int
    gem_reward::T
    gem_pos_init::Vector{CartesianIndex{2}}
    done::Bool
end

@generate_getters(CollectGemsUndirected)
@generate_setters(CollectGemsUndirected)

#####
# Directed
#####

function CollectGemsDirected(; T = Float32, height = 8, width = 8, num_gem_init = floor(Int, sqrt(height * width)), rng = Random.GLOBAL_RNG)
    objects = (AGENT, WALL, GEM)
    world = GridWorldBase(objects, height, width)

    room = Room(CartesianIndex(1, 1), height, width)
    place_room!(world, room)

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    agent_dir = RIGHT
    reward = zero(T)
    num_gem_current = num_gem_init
    gem_reward = one(T)
    gem_pos_init = CartesianIndex{2}[]
    done = false

    env = CollectGemsDirected(world, agent_pos, agent_dir, reward, rng, num_gem_init, num_gem_current, gem_reward, gem_pos_init, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::CollectGemsDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const COLLECT_GEMS_DIRECTED_LAYERS = SA.SVector(2, 3)
RLBase.state(env::CollectGemsDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_agent_dir(env), get_half_size(env), COLLECT_GEMS_DIRECTED_LAYERS)

RLBase.state_space(env::CollectGemsDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::CollectGemsDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = (copy(get_grid(env)), get_agent_dir(env))

RLBase.action_space(env::CollectGemsDirected, ::RLBase.DefaultPlayer) = DIRECTED_NAVIGATION_ACTIONS
RLBase.reward(env::CollectGemsDirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::CollectGemsDirected) = get_done(env)

function RLBase.reset!(env::CollectGemsDirected{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)
    gem_pos_init = get_gem_pos_init(env)
    num_gem_init = get_num_gem_init(env)

    world[AGENT, get_agent_pos(env)] = false
    map(pos -> world[GEM, pos] = false, gem_pos_init)
    empty!(gem_pos_init)

    for i in 1:num_gem_init
        pos = rand(rng, pos -> !any(@view world[:, pos]), env)
        world[GEM, pos] = true
        push!(gem_pos_init, pos)
    end
    set_num_gem_current!(env, num_gem_init)

    agent_start_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_agent_pos!(env, agent_start_pos)
    world[AGENT, agent_start_pos] = true

    set_agent_dir!(env, rand(rng, DIRECTIONS))

    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::CollectGemsDirected{T})(action::AbstractTurnAction) where {T}
    new_dir = turn(action, get_agent_dir(env))
    set_agent_dir!(env, new_dir)
    world = get_world(env)

    set_done!(env, get_num_gem_current(env) <= 0)
    set_reward!(env, zero(T))

    return nothing
end

function (env::CollectGemsDirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, get_agent_dir(env), agent_pos)
    if !world[WALL, dest]
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    set_reward!(env, zero(T))
    agent_pos = get_agent_pos(env)
    if world[GEM, agent_pos]
        world[GEM, agent_pos] = false
        set_num_gem_current!(env, get_num_gem_current(env) - 1)
        set_reward!(env, get_gem_reward(env))
    end

    set_done!(env, get_num_gem_current(env) <= 0)

    return nothing
end

#####
# Undirected
#####

RLBase.state_space(env::CollectGemsUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const COLLECT_GEMS_UNDIRECTED_LAYERS = SA.SVector(2, 3)
RLBase.state(env::CollectGemsUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_half_size(env), COLLECT_GEMS_UNDIRECTED_LAYERS)

RLBase.state_space(env::CollectGemsUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::CollectGemsUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = copy(get_grid(env))

RLBase.action_space(env::CollectGemsUndirected, player::RLBase.DefaultPlayer) = UNDIRECTED_NAVIGATION_ACTIONS
RLBase.reward(env::CollectGemsUndirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::CollectGemsUndirected) = get_done(env)

function CollectGemsUndirected(; T = Float32, height = 8, width = 8, num_gem_init = floor(Int, sqrt(height * width)), rng = Random.GLOBAL_RNG)
    objects = (AGENT, WALL, GEM)
    world = GridWorldBase(objects, height, width)

    room = Room(CartesianIndex(1, 1), height, width)
    place_room!(world, room)

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    reward = zero(T)
    num_gem_current = num_gem_init
    gem_reward = one(T)
    gem_pos_init = CartesianIndex{2}[]
    done = false

    env = CollectGemsUndirected(world, agent_pos, reward, rng, num_gem_init, num_gem_current, gem_reward, gem_pos_init, done)

    RLBase.reset!(env)

    return env
end

function RLBase.reset!(env::CollectGemsUndirected{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)
    gem_pos_init = get_gem_pos_init(env)
    num_gem_init = get_num_gem_init(env)

    world[AGENT, get_agent_pos(env)] = false
    map(pos -> world[GEM, pos] = false, gem_pos_init)
    empty!(gem_pos_init)

    for i in 1:num_gem_init
        pos = rand(rng, pos -> !any(@view world[:, pos]), env)
        world[GEM, pos] = true
        push!(gem_pos_init, pos)
    end
    set_num_gem_current!(env, num_gem_init)

    agent_start_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_agent_pos!(env, agent_start_pos)
    world[AGENT, agent_start_pos] = true

    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::CollectGemsUndirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, agent_pos)
    if !world[WALL, dest]
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    set_reward!(env, zero(T))
    agent_pos = get_agent_pos(env)
    if world[GEM, agent_pos]
        world[GEM, agent_pos] = false
        set_num_gem_current!(env, get_num_gem_current(env) - 1)
        set_reward!(env, get_gem_reward(env))
    end

    set_done!(env, get_num_gem_current(env) <= 0)

    return nothing
end


text:
using LinearAlgebra
using JuMP
using Hypatia
using MinimaxAdaptiveControl
using Test
##
@testset "Recreate Anders' results" begin
    A = [2 -1 1.0; 1 0 0; 0 0 0]
    B1 = Matrix{Float64}(undef, 3, 1)
    B1[:] = [0; 0; 1.0]
    B2 = -B1
    Q = I(3)
    γ = 19
    R = I(1)
    As = [A, A]
    Bs = [B1, B2]

    mac = MinimaxAdaptiveControl.MAController(As, Bs, Q, R, γ, [0.; 0; 0])
    PAnders = [
        20.628 -11.101 11.101
        -11.101 7.842 -6.842
        11.101 -6.842 7.842
    ]
    KAnders = [1.786 -1.288 1.288]

    TAnders = [155 -84.4 84.4
                -84.4 89.0 -87.5
                84.4 -87.5 89.0]
    @test norm(mac.candidates[1].P - mac.candidates[2].P) ≈ 0
    @test norm(mac.candidates[1].K + mac.candidates[2].K) ≈ 0
    @test norm(mac.candidates[1].P - PAnders) ≈ 0 atol = 1e-1
    @test norm(mac.candidates[1].K - KAnders) ≈ 0 atol = 1e-3
end

@testset "Validate Tsyn" begin
    A = [2 -1 1.0; 1 0 0; 0 0 0]
    B1 = Matrix{Float64}(undef, 3, 1)
    B1[:] = [0; 0; 1.0]
    B2 = -B1
    Q = I(3)
    γ = 20
    R = I(1)
    As = [A, A]
    Bs = [B1, B2]

    mac = MAController(As, Bs, Q, R, γ, [0.; 0; 0])
    model = Model(Hypatia.Optimizer)
    unset_silent(model) # For now broken...
    Tval, stat = Tsyn(mac, model);

    P = mac.candidates[1].P
    B = B1
    K = mac.candidates[1].K
    # Tval = Tval + 1e-4*I(3) # Might be needed
    ineq13 = P - (Q + K' * R * K + ((A - B * K)' / (inv(P) - γ^(-2) * I(3))) * (A - B * K))
    ineq14 = Tval - (Q + K' * R * K + ((A + B * K)' / (inv(P) - γ^(-2) * I(3))) * (A + B * K))
    ineq15 = Tval - (Q + K' * (R - γ^2 * B' * B) * K + (A' / (inv(Tval) - γ^(-2) * I(3)) * A))

    @test minimum(eigvals(ineq13)) > 0
    @test minimum(eigvals(ineq14)) > 0
    @test minimum(eigvals(ineq15)) > 0
    @test maximum(eigvals(Tval)) < γ^2
end

@testset "Testing X and Z" begin
    # Simple. one-dim
    A1, A2, B1, B2, Q, R = [fill(i * 1., 1, 1) for i in [1,2,3,5,7,11]]
    γ = 13
    As = [A1, A2]
    Bs = [B1, B2]
    mac = MinimaxAdaptiveControl.MAController(As, Bs, Q, R, γ, [0.; 0; 0])
    K1 = mac.candidates[1].K
    P1 = mac.candidates[1].P
    K2 = mac.candidates[2].K
    P2 = mac.candidates[2].P

    model = Model(Hypatia.Optimizer)
    @variable(model, T[1:1,1:1] in PSDCone())

    X12_ref = [(T - Q - K2' * R * K2) (A1 - B1 * K2)
            (A1 - B1 * K2)' inv(P1) - I(1) / γ^2]
    @test X12_ref == MinimaxAdaptiveControl.X(mac, T, 1, 2, 0)

    A11 = A1 - B1 * K1
    A21 = A2 - B2 * K1
    Z121_ref11 = T - Q - K1' * R * K1 + γ^2 / 2 * (A11' * A11 + A21' * A21)
    Z121_ref21 = (A11 + A21)
    Z121_ref12 = (A11 + A21)'
    Z121_ref22 = 4 / γ^4 * (γ^2 * I(1) - T)
    Z121_ref = [Z121_ref11 Z121_ref12
                Z121_ref21 Z121_ref22]

    @test Z121_ref == MinimaxAdaptiveControl.Z(mac, T, 1, 2, 1, 0)
end

@testset "test valuefunction" begin
    A = [3. -1;1 0]
    B1 = Matrix{Float64}(undef, 2, 1)
    B1[:,:] = [1. 0]
    B2 = Matrix{Float64}(undef, 2, 1)
    B2[:,:] = [5 0]

    As = [A,A]
    Bs = [B1,B2]
    Q = I(2) * 1.
    gamma = 22
    R = I(1) * 1.

    mac = MAController(As, Bs, Q, R, gamma, [0.; 0])
    model = Model(Hypatia.Optimizer)
    unset_silent(model)
    Tval, stat = Tsyn(mac, model, 0.01)

    @testset "V is monotonic" begin
        N = 100
        x = [3;-5.]
        u = [0.]
        Vs = zeros(N)
        for k = 1:N
            update!(mac, x, u)
            u = -K(mac)*x
            x = A*x + B1*u
            Vs[k] = Vbar(mac,Tval, x)
        end
        Vdiffs = Vs[2:N] - Vs[1:N-1]
        @test maximum(Vdiffs) <= 0
    end

    @testset "Reduction to ±B" begin
        A = [2 -1 1.0; 1 0 0; 0 0 0]
        B1 = Matrix{Float64}(undef, 3, 1)
        B1[:] = [0; 0; 1.0]
        B2 = -B1
        Q = I(3)
        γ = 19
        R = I(1)*1.
        As = [A, A]
        Bs = [B1, B2]

        x0 = [0.; 0; 0]
        mac = MinimaxAdaptiveControl.MAController(As, Bs, Q, R, γ,x0 )
        model = Model(Hypatia.Optimizer)
        unset_silent(model)
        Tval, stat = Tsyn(mac, model, 1e-4)

        x = [3;-5.;7]
        u = [2.]
        function Vval(mac, x,x0, u)
            P = mac.candidates[1].P
            c1 = x'*P*x - γ^2*(-x + A*x0 +B1*u)'*(-x + A*x0 +B1*u)
            c2 = x'*P*x - γ^2*(-x + A*x0 -B1*u)'*(-x + A*x0 -B1*u)
            c3 = x'*Tval*x - γ^2*((-x + A*x0)'*(-x + A*x0) + u'*B1'*B1*u)
            return maximum([c1, c2, c3])
        end

        update!(mac, x,u)
        @test Vval(mac,x,x0,u) == Vbar(mac,Tval, x)
    end
end