text:
"""
    salocd(A, C; evals, sdeg, disc = false, atol1 = 0, atol2 = 0, rtol) -> (K, Scl, blkdims)

Compute for the pair `(A,C)`, a matrix `K` such that all eigenvalues of the matrix `A+K*C` lie in the stability domain `Cs`
specified by the stability degree parameter `sdeg` and stability type parameter `disc`.
If `disc = false`, `Cs` is the set of complex numbers with real parts at most `sdeg`, while if `disc = true`,
`Cs` is the set of complex numbers with moduli at most `sdeg` (i.e., the interior of a disc of radius `sdeg` centered in the origin).
`evals` is a real or complex vector, which contains the desired eigenvalues of the matrix `A+K*C` within `Cs`.
For real data `A` and `C`, `evals` must be a self-conjugated complex set to ensure that the resulting `K` is also a real matrix.

For a pair `(A,C)` with `A` of order `n`, the number of assignable eigenvalues is `nc := n-nu`,
where `nu` is the number of fixed eigenvalues of `A`. The assignable eigenvalues are called the _observable eigenvalues_,
while the fixed eigenvalues are called the _unobservable eigenvalues_ (these are the zeros of the pencil `[A-λI; C]`).
The spectrum allocation is achieved by successively replacing the observable eigenvalues of `A` lying outside of the
stability domain `Cs` with eigenvalues provided in `evals`. All eigenvalues of `A` lying in `Cs` are kept unalterred.
If the number of specified eigenvalues in `evals` is less than the number of observable eigenvalues of `A` outside of `Cs`
(e.g., if `evals = missing`), then some of the observable eigenvalues of `A` are assigned to the nearest values
on the boundary of `Cs`. If `sdeg = missing` and `evals = missing`, the default value used for `sdeg` is -0.05
if  `disc = false` and 0.95 if `disc = true`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `C`,
and the relative tolerance for the nonzero elements of `A` and `C`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the machine epsilon of the element type of `A`.

The resulting matrix `Acl := Z'*(A+K*C)*Z` in a Schur form, where `Z` is the orthogonal/unitary matrix used
to obtain the matrix `A+K*C` in Schur form, has the form

           ( Au  *   *  )
     Acl = ( 0   Aa  *  )
           ( 0   0   Ag )

where:  `Au` contains `nu` unobservable eigenvalues of `A` lying outside `Cs`,
`Aa` contains the `na` assigned eigenvalues in `Cs` and `Ag` contains the `ng` eigenvalues of `A` in `Cs`.
The matrices `Acl` and `Z` and the vector `α` of eigenvalues of `Acl` (also of `A+K*C`)
are returned in the `Schur` object `Scl`.
The values of `nu`, `na` and `ng` are returned in the 3-dimensional vector `blkdims = [nu, na, ng]`.

Method:  The Schur method of [1] is applied to the dual pair `(A',C')` (extended to possibly unobservable pairs).

References:

[1] A. Varga.
    A Schur method for pole assignment.
    IEEE Trans. on Automatic Control, vol. 26, pp. 517-519, 1981.
"""
function salocd(A::AbstractMatrix, C::AbstractMatrix; disc::Bool = false,
                evals::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
                atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(C))),
                rtol::Real = ((size(A,1)+1)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2)))

   n = LinearAlgebra.checksquare(A)
   n == size(C,2) || throw(DimensionMismatch("A and C must have the same number of columns"))
   K, S, blkdims = saloc(copy(transpose(A)), copy(transpose(C)), evals = evals, sdeg = sdeg, disc = disc, atol1 = atol1, atol2 = atol2, rtol = rtol)
   return copy(transpose(K)), Schur(reverse(reverse(transpose(S.T),dims = 1),dims = 2), reverse(conj(S.Z),dims=2), reverse(S.values)), reverse(blkdims)
end
"""
    salocd(A, E, C; evals, sdeg, disc = false, atol1 = 0, atol2 = 0, atol3 = 0, rtol, sepinf = true, fast = true) -> (K, Scl, blkdims)

Compute for the pair `(A-λE,C)`, with `A-λE` a regular pencil, a matrix `K` such that all finite eigenvalues of the
pencil `A+K*C-λE` lie in the stability domain `Cs` specified by the stability degree parameter `sdeg` and stability type parameter `disc`.
If `disc = false`, `Cs` is the set of complex numbers with real parts at most `sdeg`, while if `disc = true`,
`Cs` is the set of complex numbers with moduli at most `sdeg` (i.e., the interior of a disc of radius `sdeg` centered in the origin).
`evals` is a real or complex vector, which contains a set of finite desired eigenvalues for the pencil `A+K*C-λE`.
For real data `A`, `E`, and `B`, `evals` must be a self-conjugated complex set to ensure that the resulting `F` is also a real matrix.

For a pair `(A-λE,C)` with `A` of order `n`, the number of assignable finite eigenvalues is `nfc := n-ninf-nfu`,
where `ninf` is the number of infinite eigenvalues of `A-λE` and `nfu` is the number of fixed finite eigenvalues of `A-λE`.
The assignable finite eigenvalues are called the _observable finite eigenvalues_,
while the fixed finite eigenvalues are called the _unobservable finite eigenvalues_ (these are the finite zeros of the pencil `[A-λE; C]`).
The spectrum allocation is achieved by successively replacing the observable finite eigenvalues of `A-λE` lying outside of the
stability domain `Cs` with eigenvalues provided in `evals`. All finite eigenvalues of `A-λE` lying in `Cs` are kept unalterred.
If the number of specified eigenvalues in `evals` is less than the number of observable finite eigenvalues of `A-λE` outside of `Cs`
(e.g., if `evals = missing`), then some of the observable finite eigenvalues of `A-λE` are assigned to the nearest values
on the boundary of `Cs`. If `sdeg = missing` and `evals = missing`, the default value used for `sdeg` is -0.05
if  `disc = false` and 0.95 if `disc = true`.

The keyword arguments `atol1`, `atol2`, , `atol3`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`, the absolute tolerance for the nonzero elements of `C`,
and the relative tolerance for the nonzero elements of `A`, `E` and `C`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the machine epsilon of the element type of `A`.

The keyword argument `sepinf` specifies the option for a preliminary separation
of the infinite eigenvalues of the pencil `A-λE` as follows: if `sepinf = false`, no separation of infinite eigenvalues is performed,
while for `sepinf = true` (the default option), a preliminary separation of the infinite eigenvalues from the finite ones is performed.
If `E` is nonsingular, then `sepinf = false` is recommended to be used. If `E` is numerically singular, then the option `sepinf = false` is used.
The separation of finite and infinite eigenvalues is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The resulting pencil `Acl-λEcl := Q'*(A+K*C-λE)*Z`,  where `Q` and `Z` are the
orthogonal/unitary matrices used to obtain the pair `(Acl,Ecl)` in a generalized Schur form (GSF), has the form

                ( Afu-λEfu    *        *       *    )
     Acl-λEcl = (   0      Afa-λEfa    *       *    )
                (   0         0     Afg-λEfg   *    )
                (   0         0        0     Ai-λEi )

where: `Afu-λEfu` contains `nfu` unobservable finite eigenvalues of `A-λE` lying outside `Cs`,
`Afa-λEfa` contains `nfa` assigned finite generalized eigenvalues in `Cs`,
`Afg-λEfg` contains `nfg` finite eigenvalues of `A-λE` in `Cs`, and
`Ai-λEi`, with `Ai` upper triangular and invertible and `Ei` upper triangular and nilpotent,
contains the `ninf` infinite eigenvalues of `A-λE`.
The matrices `Acl`, `Ecl`, `Q`, `Z` and the vectors `α` and `β` such that `α./β` are the generalized eigenvalues of
the pair `(Acl,Ecl)` are returned in `GeneralizeSchur` object `Scl`.
The values of `nfu`, `nfa`, `nfg` and `ninf` and are returned in the 4-dimensional vector `blkdims = [nfu, nfa, nfg, ninf]`.

Method:  For a pair `(A-λE,C)` with `E = I`, the dual Schur method of [1] is used, while for a general pair `(A-λE,C)` the dual generalized Schur method of [2]
is used to solve the R-stabilzation problem of [2] for the dual pair `(A'-λE',C')`.

References:

[1] A. Varga.
    A Schur method for pole assignment.
    IEEE Trans. on Automatic Control, vol. 26, pp. 517-519, 1981.

[2] A. Varga.
    On stabilization methods of descriptor systems.
    Systems & Control Letters, vol. 24, pp.133-138, 1995.
"""
function salocd(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling{Bool}}, C::AbstractMatrix; disc::Bool = false,
                evals::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
                atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))), atol3::Real = zero(real(eltype(C))),
                rtol::Real = ((size(A,1)+1)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2,atol3)),
                fast::Bool = true, sepinf::Bool = true)

   n = LinearAlgebra.checksquare(A)
   E == I || LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix"))

   n == size(C,2) || throw(DimensionMismatch("A and C must have the same number of columns"))
   K, S, blkdims = saloc(copy(transpose(A)), copy(transpose(E)), copy(transpose(C)), evals = evals, sdeg = sdeg, disc = disc, atol1 = atol1, atol2 = atol2, atol3 = atol3, sepinf = sepinf, fast = fast, rtol = rtol)
   return copy(transpose(K)), GeneralizedSchur(reverse(reverse(transpose(S.S),dims = 1),dims = 2), reverse(reverse(transpose(S.T),dims = 1),dims = 2),
                                                  reverse(S.α), reverse(S.β), reverse(conj(S.Z),dims=2), reverse(conj(S.Q),dims=2)), reverse(blkdims)
end
"""
    saloc(A, B; evals, sdeg, disc = false, atol1 = 0, atol2 = 0, rtol) -> (F, Scl, blkdims)

Compute for the pair `(A,B)`, a matrix `F` such that all eigenvalues of the matrix `A+B*F` lie in the stability domain `Cs`
specified by the stability degree parameter `sdeg` and stability type parameter `disc`.
If `disc = false`, `Cs` is the set of complex numbers with real parts at most `sdeg`, while if `disc = true`,
`Cs` is the set of complex numbers with moduli at most `sdeg` (i.e., the interior of a disc of radius `sdeg` centered in the origin).
`evals` is a real or complex vector, which contains the desired eigenvalues of the matrix `A+B*F` within `Cs`.
For real data `A` and `B`, `evals` must be a self-conjugated complex set to ensure that the resulting `F` is also a real matrix.

For a pair `(A,B)` with `A` of order `n`, the number of assignable eigenvalues is `nc := n-nu`,
where `nu` is the number of fixed eigenvalues of `A`. The assignable eigenvalues are called the controllable eigenvalues,
while the fixed eigenvalues are called the _uncontrollable eigenvalues_ (these are the zeros of the pencil `[A-λI B]`).
The spectrum allocation is achieved by successively replacing the _controllable eigenvalues_ of `A` lying outside of the
stability domain `Cs` with eigenvalues provided in `evals`. All eigenvalues of `A` lying in `Cs` are kept unalterred.
If the number of specified eigenvalues in `evals` is less than the number of controllable eigenvalues of `A` outside of `Cs`
(e.g., if `evals = missing`), then some of the controllable eigenvalues of `A` are assigned to the nearest values
on the boundary of `Cs`. If `sdeg = missing` and `evals = missing`, the default value used for `sdeg` is -0.05
if  `disc = false` and 0.95 if `disc = true`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `B`,
and the relative tolerance for the nonzero elements of `A` and `B`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the machine epsilon of the element type of `A`.

The resulting matrix `Acl := Z'*(A+B*F)*Z` in a Schur form, where `Z` is the orthogonal/unitary matrix used
to obtain the matrix `A+B*F` in Schur form, has the form

           ( Ag  *   *  )
     Acl = ( 0   Aa  *  )
           ( 0   0   Au )

where:  `Ag` contains the `ng` eigenvalues of `A` in `Cs`, `Aa` contains the `na` assigned eigenvalues in `Cs` and
`Au` contains `nu` uncontrollable eigenvalues of `A` lying outside `Cs`.
The matrices `Acl` and `Z` and the vector `α` of eigenvalues of `Acl` (also of `A+B*F`)
are returned in the `Schur` object `Scl`.
The values of `ng`, `na` and `nu` are returned in the 3-dimensional vector `blkdims = [ng, na, nu]`.

Method:  The Schur method of [1], extended to possibly uncontrollable pairs, is employed.

References:

[1] A. Varga.
    A Schur method for pole assignment.
    IEEE Trans. on Automatic Control, vol. 26, pp. 517-519, 1981.
"""
function saloc(A::AbstractMatrix, B::AbstractMatrix; disc::Bool = false,
               evals::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
               atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(B))),
               rtol::Real = ((size(A,1)+1)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2)))
   n = LinearAlgebra.checksquare(A)
   n1, m = size(B)
   n == n1 || throw(DimensionMismatch("A and B must have the same number of rows"))

   T = promote_type( eltype(A), eltype(B) )
   T <: BlasFloat || (T = promote_type(Float64,T))

   A1 = copy_oftype(A,T)
   B1 = copy_oftype(B,T)
   if ismissing(evals)
      evals1 = evals
   else
      T1 = promote_type(T,eltype(evals))
      evals1 = copy_oftype(evals,T1)
   end

   # quick exit if possible
   n == 0 && (return zeros(T,m,n), schur(zeros(T,n,n)), zeros(Int,4))

   ZERO = zero(T)

   # check for zero rows in the leading positions
   ilob = n+1
   for i = 1:n
       !iszero(view(B1,i,:)) && (ilob = i; break)
   end

   # return if B = 0
   ilob > n && (return zeros(T,m,n), schur(A1), [0,0,n] )

   # check for zero rows in the trailing positions
   ihib = ilob
   for i = n:-1:ilob+1
       !iszero(view(B1,i,:)) && (ihib = i; break)
   end

   # operate only on the nonzero rows of B
   ib = ilob:ihib
   nrmB = opnorm(view(B1,ib,:),1)

   complx = (T <: Complex)

   # sort desired eigenvalues
   if ismissing(evals1)
      evalsr = missing
      evalsc = missing
   else
      if complx
         evalsr = copy(evals1)
         evalsc = missing
      else
         evalsr = evals1[imag.(evals1) .== 0]
         isempty(evalsr) && (evalsr = missing)
         tempc = evals1[imag.(evals1) .> 0]
         if isempty(tempc)
            evalsc = missing
         else
            tempc1 = conj(evals1[imag.(evals1) .< 0])
            isequal(tempc[sortperm(real(tempc))],tempc1[sortperm(real(tempc1))]) ||
                    error("evals must be a self-conjugated complex vector")
            evalsc = [transpose(tempc[:]); transpose(conj.(tempc[:]))][:]
         end
      end
      # check that all eigenvalues are inside of the stability region
      !ismissing(sdeg) && ((disc && any(abs.(evals1) .> sdeg) )  || (!disc && any(real.(evals1) .> sdeg)))  &&
            error("The elements of evals must lie in the stability region")
   end

   # set default values of sdeg if evals = missing
   if ismissing(sdeg)
      if ismissing(evals1)
        sdeg = disc ? real(T)(0.95) : real(T)(-0.05)
        smarg = sdeg;
      else
         smarg = disc ? real(T)(0) : real(T)(-Inf)
      end
   else
      sdeg = real(T)(sdeg)
      disc && sdeg < 0 && error("sdeg must be non-negative if disc = true")
      smarg = sdeg;
   end

   nrmA = opnorm(A1,1)
   tola = max(atol1, rtol*nrmA)
   tolb = max(atol2, rtol*nrmB)

   #
   # separate stable and unstable parts with respect to sdeg
   # compute orthogonal Z  such that
   #
   #      Z^T*A*Z = [ Ag   * ]
   #                [  0  Ab ]
   #
   # where Ag has eigenvalues within the stability degree region
   # and Ab has eigenvalues outside the stability degree region.

   _, Z, α = LAPACK.gees!('V', A1)

   if disc
      select = Int.(abs.(α) .<= smarg)
   else
      select = Int.(real.(α) .<= smarg)
   end
   _, _, α = LAPACK.trsen!(select, A1, Z)
   nb = length(select[select .== 0])

   ng = n-nb
   fnrmtol = 1000*nrmA/nrmB
   fnrmtol == 0 && (fnrmtol = 1000/nrmB)

   nu = 0; na = 0
   nc = n
   F = zeros(T,m,n);
   ia = n-nb+1
   ihf = 0
   while nb > 0
      noskip = true
      if nb == 1 || complx || A1[nc,nc-1] == 0
         k = 1
      else
         k = 2
      end
      kk = nc-k+1:nc
      a2 = view(A1,kk,kk)
      evb = ordeigvals(a2)
      b2 = view(Z,ib,kk)'*view(B1,ib,:)
      if norm(b2,Inf) <= tolb
         # deflate uncontrollable block
         nb = nb-k; nc = nc-k; nu = nu+k; noskip = false
      elseif k == 1 && nb > 1 && ismissing(evalsr) && !ismissing(evalsc)
         # form a 2x2 block if there are no real eigenvalues to assign
         k = 2
         kk = nc-k+1:nc
         a2 = view(A1,kk,kk)
         if nb > 2 && A1[nc-1,nc-2] != ZERO
            # interchange last two blocks
            LAPACK.trexc!('V', nc, nc-2, A1, Z)
            evb = ordeigvals(a2)
         else
            evb = disc ? maximum(abs.(ordeigvals(a2))) : maximum(real(ordeigvals(a2)))
         end
         b2 = view(Z,ib,kk)'*view(B1,ib,:)
         if norm(b2,Inf) <= tolb
            # deflate uncontrollable block
            nb = nb-k; nc = nc-k; nu = nu+k; noskip = false
         end
      end
      if noskip
         if k == 1
            # assign a single eigenvalue
            γ, evalsr = eigselect1(evalsr, sdeg, complx ? evb[1] : real(evb[1]), disc; cflag = complx);
            if γ === nothing
                # no real eigenvalue available, adjoin a new 1x1 block if possible
                if nb == 1
                    # incompatible eigenvalues with the eigenvalue structure
                    # assign the last real pole to the default value of sdeg
                    γ = disc ? real(T)(0.95) : real(T)(-0.05)
                    f2 = -b2\(a2-I*γ)
                 else
                    # adjoin a real block or interchange the last two blocks
                    k = 2; kk = nc-1:nc;
                    a2 = view(A1,kk,kk)
                    if nb > 2 && A1[nc-1,nc-2] != ZERO
                       # interchange last two blocks
                       LAPACK.trexc!('V', nc, nc-2, A1, Z)
                       # update evb
                       evb = ordeigvals(a2)
                    else
                       # update evb
                       evb = disc ? maximum(abs.(ordeigvals(a2))) : maximum(real(ordeigvals(a2)))
                    end
                    b2 = view(Z,ib,kk)'*view(B1,ib,:)
                end
            else
               f2 = -b2\(a2-I*γ);
            end
         end
         if k == 2
            # assign a pair of eigenvalues
            γ, evalsr, evalsc = eigselect2(evalsr,evalsc,sdeg,evb[end],disc)
            f2, u = saloc2(a2,b2,γ,tola,tolb)
            if f2 === nothing  # the case b2 = 0 can not occur
               irow = 1:nc; jcol = nc-1:n;
               A1[kk,jcol] = u'*view(A1,kk,jcol); A1[irow,kk] = view(A1,irow,kk)*u;
               A1[nc,nc-1] = ZERO
               Z[:,kk] = view(Z,:,kk)*u
               nb -= 1; nc -= 1; nu += 1;
               # recover the failed selection
               imag(γ[1]) == 0 ? (ismissing(evalsr) ? evalsr = γ : evalsr = [γ; evalsr]) : (ismissing(evalsc) ? evalsc = γ : evalsc = [γ; evalsc])
            end
         end
         if f2 !== nothing
            # check for numerical stability
            norm(f2,Inf) > fnrmtol && (ihf += 1)
            X = view(B1,ib,:)*f2
            A1[1:nc,kk] += view(Z,ib,1:nc)'*X
            F += f2*view(Z,:,kk)'
            if k == 2
               # standardization step is necessary to use trexc
               i1 = 1:nc-2; lcol = nc+1:n;
               # alternative computation
               # k1 = kk[1]; k2 = kk[2]
               # RT1R, RT1I, RT2R, RT2I, CS, SN = lanv2(A1[k1,k1], A1[k1,k2], A1[k2,k1], A1[k2,k2])
               _, Z2, _ = LAPACK.gees!('V', a2)
               A1[i1,kk] = view(A1,i1,kk)*Z2
               A1[kk,lcol] = Z2'*view(A1,kk,lcol)
               Z[:,kk] = view(Z,:,kk)*Z2;
               tworeals = (A1[nc,nc-1] == ZERO)
            else
               tworeals = false
            end
            # reorder eigenvalues
            if nb > k
               try
                  LAPACK.trexc!('V', nc-k+1, ia, A1, Z)
                  tworeals && LAPACK.trexc!('V', nc, ia+1, A1, Z)
               catch
               end
            end
            nb -= k
            ia += k
            na += k
         end
      end
   end
   ihf > 0 && @warn("Possible loss of numerical reliability due to high feedback gain")
   blkdims = [ng, na, nu]
   return F, Schur(A1, Z, ordeigvals(A1)), blkdims

   # end saloc
end
"""
    saloc(A, E, B; evals, sdeg, disc = false, atol1 = 0, atol2 = 0, atol3 = 0, rtol, sepinf = true, fast = true) -> (F, Scl, blkdims)

Compute for the pair `(A-λE,B)`, with `A-λE` a regular pencil, a matrix `F` such that all finite eigenvalues of the
pencil `A+B*F-λE` lie in the stability domain `Cs` specified by the stability degree parameter `sdeg` and stability type parameter `disc`.
If `disc = false`, `Cs` is the set of complex numbers with real parts at most `sdeg`, while if `disc = true`,
`Cs` is the set of complex numbers with moduli at most `sdeg` (i.e., the interior of a disc of radius `sdeg` centered in the origin).
`evals` is a real or complex vector, which contains a set of finite desired eigenvalues for the pencil `A+B*F-λE`.
For real data `A`, `E`, and `B`, `evals` must be a self-conjugated complex set to ensure that the resulting `F` is also a real matrix.

For a pair `(A-λE,B)` with `A` of order `n`, the number of assignable finite eigenvalues is `nfc := n-ninf-nfu`,
where `ninf` is the number of infinite eigenvalues of `A-λE` and `nfu` is the number of fixed finite eigenvalues of `A-λE`.
The assignable finite eigenvalues are called the _controllable finite eigenvalues_,
while the fixed finite eigenvalues are called the _uncontrollable finite eigenvalues_ (these are the finite zeros of the pencil `[A-λE B]`).
The spectrum allocation is achieved by successively replacing the controllable finite eigenvalues of `A-λE` lying outside of the
stability domain `Cs` with eigenvalues provided in `evals`. All finite eigenvalues of `A-λE` lying in `Cs` are kept unalterred.
If the number of specified eigenvalues in `evals` is less than the number of controllable finite eigenvalues of `A-λE` outside of `Cs`
(e.g., if `evals = missing`), then some of the controllable finite eigenvalues of `A-λE` are assigned to the nearest values
on the boundary of `Cs`. If `sdeg = missing` and `evals = missing`, the default value used for `sdeg` is -0.05
if  `disc = false` and 0.95 if `disc = true`.

The keyword arguments `atol1`, `atol2`, `atol3`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`, the absolute tolerance for the nonzero elements of `B`,
and the relative tolerance for the nonzero elements of `A`, `E` and `B`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the machine epsilon of the element type of `A`.

The keyword argument `sepinf` specifies the option for a preliminary separation
of the infinite eigenvalues of the pencil `A-λE` as follows: if `sepinf = false`, no separation of infinite eigenvalues is performed,
while for `sepinf = true` (the default option), a preliminary separation of the infinite eigenvalues from the finite ones is performed.
If `E` is nonsingular, then `sepinf = false` is recommended to be used. If `E` is numerically singular, then the option `sepinf = false` is used.
The separation of finite and infinite eigenvalues is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The resulting pencil `Acl-λEcl := Q'*(A+B*F-λE)*Z`,  where `Q` and `Z` are the
orthogonal/unitary matrices used to obtain the pair `(Acl,Ecl)` in a generalized Schur form (GSF), has the form

                ( Ai-λEi    *        *       *      )
     Acl-λEcl = (   0    Afg-λEfg    *       *      )
                (   0       0     Afa-λEfa   *      )
                (   0       0        0     Afu-λEfu )

where: `Ai-λEi` with `Ai` upper triangular and invertible and `Ei` upper triangular and nilpotent, contains the `ninf` infinite eigenvalues of `A-λE`;
`Afg-λEfg` contains `nfg` finite eigenvalues of `A-λE` in `Cs`; `Afa-λEfa` contains `nfa` assigned finite generalized eigenvalues in `Cs`;
and `Afu-λEfu` contains `nfu` uncontrollable finite eigenvalues of `A-λE` lying outside `Cs`.
The matrices `Acl`, `Ecl`, `Q`, `Z` and the vectors `α` and `β` such that `α./β` are the generalized eigenvalues of
the pair `(Acl,Ecl)` are returned in the `GeneralizeSchur` object `Scl`.
The values of `ninf`, `nfg`, `nfa` and `nfu` are returned in the 4-dimensional vector `blkdims = [ninf, nfg, nfa, nfu]`.

Method:  For a pair `(A-λE,B)` with `E = I`, the Schur method of [1] is used, while for a general pair `(A-λE,B)` the generalized Schur method of [2]
is used to solve the R-stabilzation problem.

References:

[1] A. Varga.
    A Schur method for pole assignment.
    IEEE Trans. on Automatic Control, vol. 26, pp. 517-519, 1981.

[2] A. Varga.
    On stabilization methods of descriptor systems.
    Systems & Control Letters, vol. 24, pp.133-138, 1995.

"""
function saloc(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling{Bool}}, B::AbstractMatrix;
                  disc::Bool = false, evals::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
                  atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))), atol3::Real = zero(real(eltype(B))),
                  rtol::Real = ((size(A,1)+1)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2,atol3)),
                  fast::Bool = true, sepinf::Bool = true)

   n = LinearAlgebra.checksquare(A)
   if E == I
      F, S, blkdims = saloc(A, B, evals = evals, sdeg = sdeg, atol1 = atol1, atol2 = atol3, rtol = rtol)
      TS = eltype(F)
      return F, GeneralizedSchur(S.T, Matrix{TS}(I,n,n), S.values, ones(TS,n), S.Z, S.Z), [[0]; blkdims]
   end

   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix"))

   n1, m = size(B)
   n == n1 || throw(DimensionMismatch("A and B must have the same number of rows"))

   T = promote_type( eltype(A), eltype(B), eltype(E) )
   T <: BlasFloat || (T = promote_type(Float64,T))

   A1 = copy_oftype(A,T)
   E1 = copy_oftype(E,T)
   B1 = copy_oftype(B,T)
   if ismissing(evals)
      evals1 = evals
   else
      T1 = promote_type(T,eltype(evals))
      evals1 = copy_oftype(evals,T1)
   end

   # quick exit for n = 0 or B = 0

   # the following is not working for 0 dimensions
   # n == 0 || (return zeros(T,m,n), schur(zeros(T,n,n), ones(T,n,n)), zeros(T,n,m), zeros(Int,4))
   n == 0 && (return zeros(T,m,n), GeneralizedSchur(zeros(T,n,n),zeros(T,n,n),zeros(T,n),zeros(T,n),zeros(T,n,n),zeros(T,n,n)), zeros(Int,4))

   ZERO = zero(T)

   # check for zero rows in the leading positions
   ilob = n+1
   for i = 1:n
       !iszero(view(B1,i,:)) && (ilob = i; break)
   end

   # return if B = 0
   ilob > n && (return zeros(T,m,n), schur(A1,E1), [0, 0, 0, n] )

   # check for zero rows in the trailing positions
   ihib = ilob
   for i = n:-1:ilob+1
       !iszero(view(B1,i,:)) && (ihib = i; break)
   end

   # operate only on the nonzero rows of B
   ib = ilob:ihib
   nrmB = opnorm(view(B1,ib,:),1)

   complx = (T <: Complex)

   # sort desired eigenvalues
   if ismissing(evals1)
      evalsr = missing
      evalsc = missing
   else
      if complx
         evalsr = copy(evals1)
         evalsc = missing
      else
         evalsr = evals1[imag.(evals1) .== 0]
         isempty(evalsr) && (evalsr = missing)
         tempc = evals1[imag.(evals1) .> 0]
         if isempty(tempc)
            evalsc = missing
         else
            tempc1 = conj(evals1[imag.(evals1) .< 0])
            isequal(tempc[sortperm(real(tempc))],tempc1[sortperm(real(tempc1))]) ||
                    error("evals must be a self-conjugated complex vector")
            evalsc = [transpose(tempc[:]); transpose(conj.(tempc[:]))][:]
         end
      end
      # check that all eigenvalues are inside of the stability region
      !ismissing(sdeg) && ((disc && any(abs.(evals1) .> sdeg) )  || (!disc && any(real.(evals1) .> sdeg)))  &&
            error("The elements of evals must lie in the stability region")
   end

   # set default values of sdeg if evals = missing
   if ismissing(sdeg)
      if ismissing(evals1)
         sdeg = disc ? real(T)(0.95) : real(T)(-0.05)
         smarg = sdeg;
      else
         smarg = disc ? real(T)(0) : real(T)(-Inf)
      end
   else
      sdeg = real(T)(sdeg)
      disc && sdeg < 0 && error("sdeg must be non-negative if disc = true")
      smarg = sdeg;
   end

   nrmA = opnorm(A1,1)
   tola = max(atol1, rtol*nrmA)
   tolb = max(atol3, rtol*nrmB)

   Q = Matrix{T}(I,n,n)
   Z = Matrix{T}(I,n,n)

   if !sepinf
      # reduce (A,E) to generalized Schur form
      istriu(E1) ||  _qrE!(A1, E1, Q, missing; withQ = true)
      rcond = LAPACK.trcon!('1','U','N',E1)
      if rcond < eps(real(T))
         @warn("E is numerically singular: computations resumed with sepinf = true")
         sepinf = true
         A1 = copy_oftype(A,T)
         E1 = copy_oftype(E,T)
         Q = Matrix{T}(I,n,n)
      else
         gghrd!('V','V',1, n, A1, E1, Q, Z)
         _, _, α, β, _, _ = hgeqz!('V','V',1, n, A1, E1, Q, Z)
         ninf = 0
         if disc
            select = Int.((abs.(α) .<= smarg*abs.(β)))
         else
            select = Int.((real.(α ./ β) .< smarg))
         end
      end
   end

   if sepinf
      _, blkdims = fisplit!(A1, E1, Q, Z, missing, missing; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = true, withZ = true)
      ninf = blkdims[1]
      ilo = ninf+1;
      gghrd!('V','V',ilo, n, A1, E1, Q, Z)
      _, _, α, β, _, _ = hgeqz!('V','V',ilo, n, A1, E1, Q, Z)
      i2 = ilo:n
      if disc
         select2 = Int.(abs.(α[i2]) .<= smarg*abs.(β[i2]))
      else
         select2 = Int.(real.(α[i2] ./ β[i2]) .< smarg)
      end
      select = [ones(Int,ninf);select2]
   end

   # separate stable and unstable eigenvalues

   # compute orthogonal Q and Z such that
   #
   #                   [ Ai-λEi   *      *    ]
   #      Q*(A-λE)*Z = [   0    Ag-λEg   *    ]
   #                   [   0      0    Ab-λEb ]
   #
   # where Ag-λEg has eigenvalues within the stability degree region
   # and Ab-λEb has eigenvalues outside the stability degree region.

   _, _, α, β, _, _ = LAPACK.tgsen!(select, A1, E1, Q, Z)
   nb = length(select[select .== 0])

   nfg = n-nb-ninf
   fnrmtol = 1000*max(nrmA,1)/nrmB

   nfu = 0; nfa = 0
   nc = n
   F = zeros(T,m,n);
   ia = n-nb+1
   ihf = 0
   while nb > 0
      noskip = true
      if nb == 1 || complx || A1[nc,nc-1] == 0
         k = 1
      else
         k = 2
      end
      kk = nc-k+1:nc
      a2 = view(A1,kk,kk)
      e2 = view(E1,kk,kk)
      evb, = ordeigvals(a2,e2)
      b2 = view(Q,ib,kk)'*view(B1,ib,:)
      if norm(b2,Inf) <= tolb
         # deflate uncontrollable stable block
         nb = nb-k; nc = nc-k; nfu = nfu+k; noskip = false
      elseif k == 1 && nb > 1 && ismissing(evalsr) && !ismissing(evalsc)
         # form a 2x2 block if there are no real no real eigenvalues to assign
         k = 2
         kk = nc-k+1:nc
         a2 = view(A1,kk,kk)
         e2 = view(E1,kk,kk)
         if nb > 2 && A1[nc-1,nc-2] != ZERO
            # interchange last two blocks
            tgexc!(true, true, nc, nc-2, A1, E1, Q, Z)
            evb, = ordeigvals(a2,e2)
         else
            evb = disc ? maximum(abs.(ordeigvals(a2,e2)[1])) : maximum(real(ordeigvals(a2,e2)[1]))
         end
         b2 = view(Q,ib,kk)'*view(B1,ib,:)
         if norm(b2,Inf) <= tolb
            # deflate uncontrollable block
            nb = nb-k; nc = nc-k; nfu = nfu+k; noskip = false
         end
      end
      if noskip
         if k == 1
            # assign a single eigenvalue
            γ, evalsr = eigselect1(evalsr, sdeg, complx ? evb[1] : real(evb[1]), disc; cflag = complx);
            if γ === nothing
                # no real pole available, adjoin a new 1x1 block if possible
                if nb == 1
                    # incompatible eigenvalues with the eigenvalue structure
                    # assign the last real pole to the default value of sdeg
                    γ = disc ? real(T)(0.95) : real(T)(-0.05)
                    f2 = -b2\(a2-e2*γ)
                else
                    # adjoin a real block or interchange the last two blocks
                    k = 2; kk = nc-k+1:nc;
                    a2 = view(A1,kk,kk)
                    e2 = view(E1,kk,kk)
                    if nb > 2 && A1[nc-1,nc-2] != ZERO
                       # interchange last two blocks
                       tgexc!(true, true, nc, nc-2, A1, E1, Q, Z)
                    end
                    # update evb
                    evb = maximum(real(ordeigvals(a2,e2)[1]))
                    b2 = view(Q,ib,kk)'*view(B1,ib,:)
                end
            else
               f2 = -b2\(a2-e2*γ);
            end
         end
         if k == 2
            # assign a pair of eigenvalues
            γ, evalsr, evalsc = eigselect2(evalsr,evalsc,sdeg,evb[end],disc)
            f2, u, v = saloc2(a2,e2,b2,γ,tola,tolb)
            if f2 === nothing  # the case b2 = 0 can not occur
               irow = 1:nc; jcol = nc-1:n;
               A1[kk,jcol] = u'*view(A1,kk,jcol); A1[irow,kk] = view(A1,irow,kk)*v;
               A1[nc,nc-1] = ZERO
               E1[kk,jcol] = u'*view(E1,kk,jcol); E1[irow,kk] = view(E1,irow,kk)*v;
               E1[nc,nc-1] = ZERO
               Q[:,kk] = view(Q,:,kk)*u;
               Z[:,kk] = view(Z,:,kk)*v;
               nb -= 1; nc -= 1; nfu += 1
               # recover the failed selection
               imag(γ[1]) == 0 ? (ismissing(evalsr) ? evalsr = γ : evalsr = [γ; evalsr]) : (ismissing(evalsc) ? evalsc = γ : evalsc = [γ; evalsc])
            end
         end
         if f2 !== nothing
            norm(f2,Inf) > fnrmtol && (ihf += 1)
            # update matrices Acl and F
            X = view(B1,ib,:)*f2
            A1[1:nc,kk] += view(Q,ib,1:nc)'*X
            F += f2*view(Z,:,kk)'
            if k == 2
               # standardization step is necessary to use tgsen
               i1 = 1:nc-2; lcol = nc+1:n;
               _, _, _, _, Q2, Z2 = LAPACK.gges!('V','V',a2,e2)
               A1[i1,kk] = view(A1,i1,kk)*Z2
               E1[i1,kk] = view(E1,i1,kk)*Z2
               A1[kk,lcol] = Q2'*view(A1,kk,lcol)
               E1[kk,lcol] = Q2'*view(E1,kk,lcol)
               Q[:,kk] = view(Q,:,kk)*Q2;
               Z[:,kk] = view(Z,:,kk)*Z2;
               tworeals = (A1[nc,nc-1] == 0)
            else
               tworeals = false
            end
            # reorder eigenvalues
            if nb > k
               tgexc!(true, true, nc-k+1, ia, A1, E1, Q, Z)
               tworeals && tgexc!(true, true, nc, ia+1, A1, E1, Q, Z)
            end
            nb -= k
            ia += k
            nfa += k
         end
      end
   end
   ihf > 0 && @warn("Possible loss of numerical reliability due to high feedback gain")
   blkdims = [ninf, nfg, nfa, nfu]
   _, α, β = ordeigvals(A1,E1)
   return F, GeneralizedSchur(A1, E1, α, β, Q, Z), blkdims

   # end saloc
end
"""
    salocinfd(A, E, C; atol1 = 0, atol2 = 0, atol3 = 0, rtol, sepinf = true, fast = true) -> (K, L, Scl, blkdims)

Compute for the pair `(A-λE,C)`, with `A-λE` a regular pencil, the matrices `K` and `L` such that all eigenvalues of the
pencil `A+K*C-λ(E+L*C)` are infinite. For a pair `(A-λE,C)` with fixed (unobservable) finite eigenvalues, only the assignable (observable)
finite eigenvalues are moved to infinity.

For a pair `(A-λE,C)` with `A` of order `n`, the number of assignable infinite eigenvalues is `nia := n-ninf-nfu`,
where `ninf` is the number of infinite eigenvalues of `A-λE` and `nfu` is the number of fixed finite eigenvalues of `A-λE`.
The assignable finite eigenvalues are called the _observable finite eigenvalues_,
while the fixed finite eigenvalues are called the _unobservable finite eigenvalues_ (these are the finite zeros of the pencil `[A-λE; C]`).
The spectrum allocation is achieved by successively replacing the observable finite eigenvalues of `A-λE` with infinite eigenvalues.
 All infinite eigenvalues of `A-λE` are kept unalterred.

The keyword arguments `atol1`, `atol2`, , `atol3`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`, the absolute tolerance for the nonzero elements of `C`,
and the relative tolerance for the nonzero elements of `A`, `E` and `C`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the machine epsilon of the element type of `A`.

The preliminary separation of finite and infinite eigenvalues is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The resulting pencil `Acl-λEcl := Q'*(A+K*C-λ(E+L*C))*Z`,  where `Q` and `Z` are the
orthogonal/unitary matrices used to obtain the pair `(Acl,Ecl)` in a generalized Schur form (GSF), has the form

                ( Afu-λEfu    *        *     )
     Acl-λEcl = (   0      Aia-λEia    *     )
                (   0         0     Aig-λEig )

where: `Afu-λEfu`, with the pair `(Afu,Efu)` in a generalized Schur form,  contains `nfu` fixed (unobservable) finite eigenvalues of `A-λE`,
`Aia-λEia`, with `Aia` upper triangular and invertible and `Eia` upper triangular and nilpotent, contains `nia` assigned infinite generalized eigenvalues;
`Aig-λEig`, with `Aig` upper triangular and invertible and `Eig` upper triangular and nilpotent, contains the `ninf` infinite eigenvalues of `A-λE`.
The matrices `Acl`, `Ecl`, `Q`, `Z` and the vectors `α` and `β` such that `α./β` are the generalized eigenvalues of
the pair `(Acl,Ecl)` are returned in the `GeneralizeSchur` object `Scl`.
The values of `nfu`, `nia` and `ninf` and are returned in the 3-dimensional vector `blkdims = [nfu, nia, ninf]`.

Method:  For a general pair `(A-λE,C)` the dual of the modified generalized Schur method of [1] is used to determine `K` and `L` such that the pair
`(A+K*C,E+L*G)` has only infinite eigenvalues.

References:

[1] A. Varga.
    On stabilization methods of descriptor systems.
    Systems & Control Letters, vol. 24, pp.133-138, 1995.
"""
function salocinfd(A::AbstractMatrix, E::AbstractMatrix, C::AbstractMatrix; fast::Bool = true,
                atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))), atol3::Real = zero(real(eltype(C))),
                rtol::Real = ((size(A,1)+1)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2,atol3)) )

   n = LinearAlgebra.checksquare(A)
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix"))

   n == size(C,2) || throw(DimensionMismatch("A and C must have the same number of columns"))
   K, L, S, blkdims = salocinf(copy(transpose(A)), copy(transpose(E)), copy(transpose(C)), atol1 = atol1, atol2 = atol2, atol3 = atol3, fast = fast, rtol = rtol)
   #return reverse(transpose(K), dims = 2), reverse(transpose(L), dims = 2),
   #return reverse(transpose(K), dims = 1), reverse(transpose(L), dims = 1),
   return copy(transpose(K)), copy(transpose(L)),
          GeneralizedSchur(reverse(reverse(transpose(S.S),dims = 1),dims = 2), reverse(reverse(transpose(S.T),dims = 1),dims = 2),
                           reverse(S.α), reverse(S.β), reverse(conj(S.Z),dims=2), reverse(conj(S.Q),dims=2)),
          reverse(blkdims)
end
"""
    salocinf(A, E, B; atol1 = 0, atol2 = 0, atol3 = 0, rtol, fast = true) -> (F, G, Scl, blkdims)

Compute for the controllable pair `(A-λE,B)`, with `A-λE` a regular pencil, two matrices `F` and `G` such that all eigenvalues of the
pencil `A+B*F-λ(E+B*G)` are infinite. For a pair `(A-λE,B)` with fixed (uncontrollable) finite eigenvalues, only the assignable (controllable)
 finite eigenvalues are moved to infinity.

For a pair `(A-λE,B)` with `A` of order `n`, the number of assignable infinite eigenvalues is `nia := n-ninf-nfu`,
where `ninf` is the number of infinite eigenvalues of `A-λE` and `nfu` is the number of fixed finite eigenvalues of `A-λE`.
The assignable finite eigenvalues are called the _controllable finite eigenvalues_,
while the fixed finite eigenvalues are called the _uncontrollable finite eigenvalues_ (these are the finite zeros of the pencil `[A-λE B]`).
The spectrum allocation is achieved by successively replacing the controllable finite eigenvalues of `A-λE` with infinite eigenvalues.
 All infinite eigenvalues of `A-λE` are kept unalterred.

The keyword arguments `atol1`, `atol2`, `atol3`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`, the absolute tolerance for the nonzero elements of `B`,
and the relative tolerance for the nonzero elements of `A`, `E` and `B`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the machine epsilon of the element type of `A`.

The preliminary separation of finite and infinite eigenvalues is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The resulting pencil `Acl-λEcl := Q'*(A+B*F-λ(E+B*G))*Z`,  where `Q` and `Z` are the
orthogonal/unitary matrices used to obtain the pair `(Acl,Ecl)` in a generalized Schur form (GSF), has the form

                ( Aig-λEig    *        *     )
     Acl-λEcl = (   0      Aia-λEia    *     )
                (   0         0     Afu-λEfu )

where: `Aig-λEig` with `Aig` upper triangular and invertible and `Eig` upper triangular and nilpotent, contains the `ninf` infinite eigenvalues of `A-λE`;
`Aia-λEia`, with `Aia` upper triangular and invertible and `Eia` upper triangular and nilpotent, contains `nia` assigned infinite generalized eigenvalues;
and `Afu-λEfu`, with the pair `(Afu,Efu)` in a generalized Schur form,  contains `nfu` fixed (uncontrollable) finite eigenvalues of `A-λE`.
The matrices `Acl`, `Ecl`, `Q`, `Z` and the vectors `α` and `β` such that `α./β` are the generalized eigenvalues of
the pair `(Acl,Ecl)` are returned in the `GeneralizeSchur` object `Scl`.
The values of `ninf`, `nia` and `nfu` are returned in the 3-dimensional vector `blkdims = [ninf, nia, nfu]`.

Method:  For a general pair `(A-λE,B)` the modified generalized Schur method of [1] is used to determine `F` and `G` such that the pair
`(A+B*F,E+B*G)` has only infinite eigenvalues.

References:

[1] A. Varga.
    On stabilization methods of descriptor systems.
    Systems & Control Letters, vol. 24, pp.133-138, 1995.
"""
function salocinf(A::AbstractMatrix, E::AbstractMatrix, B::AbstractMatrix;
                  atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))), atol3::Real = zero(real(eltype(B))),
                  rtol::Real = ((size(A,1)+1)*eps(real(float(one(eltype(A))))))*iszero(max(atol1,atol2,atol3)),
                  fast::Bool = true)

   n = LinearAlgebra.checksquare(A)

   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix"))

   n1, m = size(B)
   n == n1 || throw(DimensionMismatch("A and B must have the same number of rows"))

   T = promote_type( eltype(A), eltype(B), eltype(E) )
   T <: BlasFloat || (T = promote_type(Float64,T))

   A1 = copy_oftype(A,T)
   E1 = copy_oftype(E,T)
   B1 = copy_oftype(B,T)

   # quick exit for n = 0 or B = 0

   # the following is not working for 0 dimensions
   # n == 0 || (return zeros(T,m,n), schur(zeros(T,n,n), ones(T,n,n)), zeros(T,n,m), zeros(Int,4))
   n == 0 && (return zeros(T,m,n), zeros(T,m,n), GeneralizedSchur(zeros(T,n,n),zeros(T,n,n),zeros(T,n),zeros(T,n),zeros(T,n,n),zeros(T,n,n)), zeros(Int,3))

   ZERO = zero(T)

   # check for zero rows in the leading positions
   ilob = n+1
   for i = 1:n
       !iszero(view(B1,i,:)) && (ilob = i; break)
   end

   # return if B = 0
   ilob > n && (return zeros(T,m,n), zeros(T,m,n), schur(A1,E1), [0, 0, n] )

   # check for zero rows in the trailing positions
   ihib = ilob
   for i = n:-1:ilob+1
       !iszero(view(B1,i,:)) && (ihib = i; break)
   end

   # operate only on the nonzero rows of B
   ib = ilob:ihib
   nrmB = opnorm(view(B1,ib,:),1)

   complx = (T <: Complex)

   nrmA = opnorm(A1,1)
   nrmE = opnorm(E1,1)
   tola = max(atol1, rtol*nrmA)
   tole = max(atol2, rtol*nrmE)
   tolb = max(atol3, rtol*nrmB)

   Q = Matrix{T}(I,n,n)
   Z = Matrix{T}(I,n,n)

   _, blkdims = fisplit!(A1, E1, Q, Z, missing, missing; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = true, withZ = true)
   ninf = blkdims[1]
   ilo = ninf+1;
   gghrd!('V','V',ilo, n, A1, E1, Q, Z)
   _, _, α, β, _, _ = hgeqz!('V','V',ilo, n, A1, E1, Q, Z)
   i2 = ilo:n

   nb = n-ninf
   fnrmtol = 1000*max(nrmA,nrmE,1)/nrmB
   scale = max(nrmA,1)/10

   nfu = 0; nia = 0
   nc = n
   F = zeros(T,m,n);
   G = zeros(T,m,n);
   ia = ninf+1
   ihf = 0

   while nb > 0
      noskip = true
      if nb == 1 || complx || A1[nc,nc-1] == 0
         k = 1
      else
         k = 2
      end
      kk = nc-k+1:nc
      a2 = view(A1,kk,kk)
      e2 = view(E1,kk,kk)
      b2 = view(Q,ib,kk)'*view(B1,ib,:)
      if norm(b2,Inf) <= tolb
         # deflate uncontrollable block
         nb = nb-k; nc = nc-k; nfu = nfu+k; noskip = false
      else
         # perturb a2 such that a2+b2*f2 is sufficiently nonzero
         if maximum(abs.(a2)) < maximum(abs.(e2)) /10000 || rank(a2,atol=tola) < k
            f2 = lmul!(scale/nrmB,rand(T,m,k).+1)
            # update A and F
            X = view(B1,ib,:)*f2
            A1[1:nc,kk] += view(Q,ib,1:nc)'*X
            F += f2*view(Z,:,kk)'
         end
         if k == 1
            # assign a single infinite eigenvalue
            g2 = -b2\e2
            X = view(B1,ib,:)*g2
            E1[1:nc,kk] += view(Q,ib,1:nc)'*X
            G += g2*view(Z,:,kk)'
            # reorder eigenvalues
            if nb > k
               tgexc!(true, true, nc-k+1, ia, A1, E1, Q, Z)
            end
            E1[ia,ia] = ZERO
         else
            # assign two infinite eigenvalues
            g2,  = saloc2(e2,a2,b2,zeros(T,2),tole,tolb)
            # update E and G
            X = view(B1,ib,:)*g2
            E1[1:nc,kk] += view(Q,ib,1:nc)'*X
            G += g2*view(Z,:,kk)'
            # perform standardization step to form two 1x1 blocks
            Q2, e2[1,1] = givens(e2[1,1],e2[2,1],1,2)
            e2[2,1] = ZERO;
            lmul!(Q2,view(E1,kk,nc:n))
            e2[2,2] = ZERO
            lmul!(Q2,view(A1,kk,nc-1:n))
            rmul!(view(Q,:,kk),Q2')
            Z2, r = givens(conj(a2[2,2]),conj(a2[2,1]),2,1)
            a2[2,2] = conj(r)
            a2[2,1] = ZERO
            rmul!(view(E1,1:nc-1,kk),Z2')
            rmul!(view(A1,1:nc-1,kk),Z2')
            rmul!(view(Z,:,kk),Z2')
            # reorder eigenvalues
            if nb > k
               tgexc!(true, true, nc-k+1, ia, A1, E1, Q, Z)
               tgexc!(true, true, nc, ia+1, A1, E1, Q, Z)
            end
            E1[ia,ia] = ZERO
            E1[ia+1,ia+1] = ZERO
         end
         nb -= k
         ia += k
         nia += k
      end
   end
   ihf > 0 && @warn("Possible loss of numerical reliability due to high feedback gain")
   blkdims = [ninf, nia, nfu]
   _, α, β = ordeigvals(A1,E1)
   return F, G, GeneralizedSchur(A1, E1, α, β, Q, Z), blkdims

   # end salocinf
end
"""
     ev = ordeigvals(A)

Compute the vector `ev` of eigenvalues of a Schur matrix `A` in their order of appearance down the diagonal of `A`.
`ordeigvals` is an order-preserving version of `eigvals` for triangular/quasitriangular matrices.
"""
function ordeigvals(A::AbstractMatrix{T}) where T
   isqtriu(A) || error("A must be in Schur form")
   n = size(A,1)
   ZERO = zero(T)
   if T <: Complex
      α = Vector{T}(undef, n)
      for i = 1:n
         α[i] = A[i,i]
      end
      return α
   else
      αr = Vector{T}(undef, n)
      αi = zeros(T,n)
      for i = 1:n
         αr[i] = A[i,i]
      end
      for i = 1:n-1
         if A[i+1,i] != ZERO
            if A[i,i] == A[i+1,i+1]
               # exploit LAPACK 2x2 block structure
               ei = sqrt(abs(A[i,i+1]))*sqrt(abs(A[i+1,i]))
               αi[i] = ei
               αi[i+1] = -ei
            else
               # an arbitrary 2x2 block
               αr[i], αi[i], αr[i+1], αi[i+1] = lanv2(A[i,i], A[i,i+1], A[i+1,i], A[i+1,i+1])
            end
            i +=1
         end
      end
      return iszero(αi) ? αr : Complex.(αr,αi)
   end
end
"""
    ordeigvals(A, B) -> (ev, α, β)

Compute the vector `ev` of generalized eigenvalues of the pair `(A,B)`, with `A` a Schur matrix and `B` an upper triangular matrix,
in their order of appearance down the diagonals of `A` and `B`.
`ordeigvals` is an order-preserving version of `eigvals` for (triangular/quasi-upper-triangular, triangular) pairs of matrices.
`α` is a complex vector and `β` is a real vector such that the generalized eigenvalues `ev` can be alternatively obtained with `α./β`.
"""
function ordeigvals(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T
   na, nb = LinearAlgebra.checksquare(A,B)
   na == nb ||  throw(DimensionMismatch("A and B must have the same size"))
   (isqtriu(A) && istriu(B)) || error("The pair (A,B) must be in a generalized Schur form")
   n = size(A,1)
   complx = T <: Complex
   if complx
      α = Vector{T}(undef, n)
      β = Vector{T}(undef, n)
      for i = 1:n
         α[i] = A[i,i]
         β[i] = B[i,i]
      end
      ev = α./β
      ev[iszero.(β)] .= Inf
      ev[iszero.(α) .& iszero.(β)] .= NaN + NaN*im
      return ev, α, β
   else
      αr = Vector{T}(undef, n)
      αi = Vector{T}(undef, n)
      β = Vector{T}(undef, n)
      ZERO = zero(T)
      small = safemin(T)
      i = 1
      while i <= n
         pair = (i < n && A[i+1,i] != ZERO)
         if pair
            ii = i:i+1
            β[i], β[i+1], αr[i], αr[i+1], wi  = lag2(A[ii,ii],B[ii,ii],small)
            iszero(wi) ? (αi[i] = wi; αi[i+1] = wi) : (αi[i] = wi; αi[i+1] = -wi; β[i+1] = β[i])
            i += 2
         else
            αr[i] = A[i,i]
            αi[i] = ZERO
            β[i] = B[i,i]
            i += 1
         end
      end
      α = iszero(αi) ? αr : Complex.(αr,αi)
      ev = α./β
      complx || (ev[imag.(ev) .> 0] = conj(ev[imag.(ev) .< 0]))
      ev[iszero.(β)] .= Inf
      ev[iszero.(α) .& iszero.(β)] .= NaN
      return ev, α, β
   end
end
"""
    saloc2(A, B, evc, tola, tolb) -> (F, U)

Compute for the real pair `(A,B)` with `A` of order two, a matrix `F` such that the eigenvalues
of the matrix `A+B*F` are equal to the complex conjugate pair contained in `evc`.
The absolute thresholds `tola` and `tolb` for the nonzero elements in `A` and `B`, respectively,
are used for controllability checks.

If the pair `(A,B)` is uncontrollable, then `F = nothing` and `U` contains an orthogonal transformation
matrix such that the transformed pair `(U'*A*U, U'*B)` is in the form

                     [ A1  X    B1 ]
    [ U'*A*U U'*B] = [             ] ,
                     [  0  A2   0  ]

where the pair `(A1,B1)` is controllable. If `norm(B) < tolb`, then `U` is the `2x2` identity matrix.
"""
function saloc2(A::AbstractMatrix{T},B::AbstractMatrix{T},evc::AbstractVector{T1},tola,tolb) where {T, T1}

   # check controllability and determine rank of B
   m = size(B,2);
   SF = svd(B; full=true); s1 = SF.S[1]
   U = SF.U
   v1 = SF.V
   if m == 1
      s2 = zero(T);
   else
      s2 = SF.S[2]
   end
   rankB = 0
   s1 > tolb && (rankB += 1)
   s2 > tolb && (rankB += 1)
   # return if norm(B) < tolb
   rankB == 0 && (return nothing, Matrix{T}(I,2,2))

   at = U'*A*U;
   if s2 <= tolb && abs(at[2,1]) <= tola
         # return if rank(B) == 1 and (A-lambda*I,B) uncontrollable
      return nothing, U
   end
   if rankB == 2
      # try with a direct solution, without inversion
      γ = [real(evc[1]) imag(evc[1]); imag(evc[2]) real(evc[2])];
      ftry = -B\(A-γ);
   end
   sc = at[1,1]+at[2,2];
   sp = real(evc[1]+evc[2]); pp = real(evc[1]*evc[2]);
   k11 = (sc-sp)/s1;
   k12 = at[2,2]/at[2,1]*k11+(pp-at[1,1]*at[2,2]+at[1,2]*at[2,1])/at[2,1]/s1;
   F = v1*[[-k11 -k12]; zeros(m-1,2)]*U';
   if rankB == 2 && norm(F) > norm(ftry)
     # choose the lower norm feedback
      F = ftry;
   end
   return F, U
   # end saloc2
end
"""
    saloc2(A, B, E, evc, tola, tolb) -> (F, U, V)

Compute for the real pair `(A-λE,B)` with `A-λE` of order two, a matrix `F` such that the
generalized eigenvalues of the matrix pair `(A+B*F,E)` are equal to the complex conjugate pair contained in `evc`.
The absolute thresholds `tola` and `tolb` for the nonzero elements in `A` and `B`, respectively,
are used for controllability checks.

If the pair `(A-λE,B)` is uncontrollable, then `F = nothing` and `U` and `V` contains orthogonal transformation
matrices such that the transformed pair `(U'*A*V-λU'*E*V, U'*B)` is in the form

                             [ A1-λE1   X    B1 ]
    [ U'*A*V-λU'*E*V U'*B] = [                  ] ,
                             [  0    A2-λE2  0  ]

where the pair `(A1-λE1,B1)` is controllable. If `norm(B) < tolb`, then `U` and `V` are `2x2` identity matrices.
"""

function saloc2(A::AbstractMatrix{T},E::AbstractMatrix{T},B::AbstractMatrix{T},evc::AbstractVector{T1},tola,tolb) where {T, T1}

    # check controllability and determine rank of B
    m = size(B,2);
    SF = svd(B; full=true); s1 = SF.S[1]
    U = SF.U
    v1 = SF.V
    if m == 1
      s2 = zero(T);
    else
      s2 = SF.S[2]
    end
    rankB = 0
    s1 > tolb && (rankB += 1)
    s2 > tolb && (rankB += 1)
    # return if norm(B) < tolb
    rankB == 0 && (return nothing, Matrix{T}(I,2,2), Matrix{T}(I,2,2))

    at = U'*A; et = U'*E;
    # determine V such that et = U'*E*V is upper triangular
    # FQ = qr([et[2,2]; et[2,1]]); V = FQ.Q*Matrix{T}(I,2,2); V = V[2:-1:1,2:-1:1]';
    G,  = LinearAlgebra.givens(et[2,2], et[2,1],1,2)
    # at = at*G;
    rmul!(at,G)
    if s2 <= tolb && abs(at[2,1]) <= tola
       # return if rank(B) == 1 and (A-lambda*E,B) uncontrollable
       # return nothing, U, V
       return nothing, U, [G.c G.s; -G.s G.c]
    end
    # reduce to standard case
    SF = svd(E\B; full=true);
    s1 = SF.S[1]; U = SF.U; v1 = SF.V
    at = U'*(E\A)*U;
    if rankB == 2
       # try with a direct solution, without inversion
       γ = [real(evc[1]) imag(evc[1]); imag(evc[2]) real(evc[2])];
       ftry = -B\(A-E*γ);
    end
    sc = at[1,1]+at[2,2];
    sp = real(evc[1]+evc[2]); pp = real(evc[1]*evc[2]);
    k11 = (sc-sp)/s1;
    k12 = at[2,2]/at[2,1]*k11+(pp-at[1,1]*at[2,2]+at[1,2]*at[2,1])/at[2,1]/s1;
    F = v1*[[-k11 -k12]; zeros(T,m-1,2)]*U'; V = Matrix{eltype(B)}(I,m,m)
    if rankB == 2 && norm(F) > norm(ftry)
      # choose the lower norm feedback
       F = ftry;
    end
    return F, U, V
    # end saloc2
end
"""
    eigselect1(ev, sdeg, evref, disc; cflag) -> (γ, evupd)

Select a real or complex eigenvalue `γ` to be assigned from a given set of eigenvalues
in the vector `ev`. The resulting `evupd` contains the remaining eigenvalues from `ev`.
The selected eigenvalue `γ` is the nearest one to the reference value `evref`.
If `ev = missing` and `disc = false`, then, if `cflag = false` (in the real case) `γ` is set equal to the
desired stability degree `sdeg`, while if `cflag = true` (in the complex case) `γ` is set
`γ = sdeg + im*imag(evref)`.
If `ev = missing` and `disc = true`, then, if `cflag = false` (in the real case) `γ` is set equal to the
desired stability degree `sdeg`, while if `cflag = true` (in the complex case) `γ` is set
`γ = sdeg/abs(evref))*evref`.
If `ev = missing` and `sdeg = missing`, then `γ = nothing`.
"""
function eigselect1(ev::Union{AbstractVector,Missing},sdeg::Union{Real,Missing},evref::Union{Real,Complex},disc::Bool = false; cflag::Bool = false )
   T = typeof(evref)
   if ismissing(ev)
      if ismissing(sdeg)
         γ = nothing
      else
         γ = disc ? (cflag ? (sdeg/abs(evref))*evref : sdeg) : (cflag ? sdeg+im*imag(evref) : sdeg)
      end
      evupd = missing
   else
      i = argmin(abs.(ev.-evref))
      γ = cflag ? ev[i] : real(ev[i]);
      evupd = [ev[1:i-1,1]; ev[i+1:end,1]]
      isempty(evupd) && (evupd = missing)
   end
   return γ, evupd

# end eigselect1
end
"""
    eigselect2(evr, evc, sdeg, evref, disc) -> (γ, evrupd, evcupd)

Select a pair of eigenvalues `γ[1]` and `γ[2]` to be assigned from given sets of real eigenvalues
in the vector `evr` and complex eigenvalues in the vector `evc`.
The selected eigenvalues are the nearest ones to the reference value `evref`.
A pair of complex eigenvalues is always selected, unless `evc = missing`, in which case a pair
of real eigenvalues is selected. The resulting `evrupd` contains the remaining real eigenvalues from `evr` and
the resulting `evcupd` contains the remaining complex eigenvalues from `evc`.
If the desired stability degree `sdeg` is missing, then a default value `sdeg = sdegdef` is used, where
`sdegdef = -0.05`, if `disc = false`, and `sdegdef = .95`, if `disc = true`.
If `evc = missing` and `evr` contains at least two eigenvalues, then
`γ[1]` and `γ[2]` are chosen the two nearest real eigenvalues to `evref`.
If `evc = missing` and `evr` contains only one value, then `γ[1] = evr` and `γ[2] = sdeg`.
If `evr = missing`, `evc = missing` and `disc = false`, then
`γ[1] = sdeg + im*imag(evref)` and `γ[2] = sdeg - im*imag(evref)`.
If `evr = missing`, `evc = missing` and `disc = true`, then
`γ[1] = sdeg/abs(evref))*evref[1]` and `γ[2] = sdeg/abs(evref))*evref[2]`.
"""
function eigselect2(evr::Union{AbstractVector,Missing},evc::Union{AbstractVector,Missing},sdeg::Union{Real,Missing},evref::Union{Real,Complex},disc::Bool)
   # corrected version to run under Julia 1.6
   evref = real(evref) + im*abs(imag(evref))
   if ismissing(evr) && ismissing(evc)
      if ismissing(sdeg)
         T = typeof(evref)
         sdegdef = disc ? real(T)(0.95) : real(T)(-0.05)
         evi = imag(evref)
         γ = [complex(sdegdef,evi), complex(sdegdef,-evi)]
      else
        if disc
           γ = [evref, conj(evref)]
           γ = (sdeg/abs(evref))*γ
        else
           evi = imag(evref)
           γ = [complex(sdeg,evi), complex(sdeg,-evi)]
        end
      end
      evrupd = missing
      evcupd = missing
   elseif ismissing(evc)
      # select two real eigenvalues
      if length(evr) < 2
         if ismissing(sdeg)
            T = typeof(evref)
            sdegdef = disc ? real(T)(0.95) : real(T)(-0.05)
            γ = [evr[1],sdegdef]
         else
            γ = [evr[1],sdeg]
         end
         evrupd = missing
         evcupd = missing
      else
         evr = evr[sortperm(abs.(evr .- evref))]
         γ = [ evr[1], evr[2]]
         evrupd = evr[3:end]
         isempty(evrupd) && (evrupd = missing)
         evcupd = missing
      end
   else
      i = argmin(abs.(evc .- evref))
      γ = [evc[i],evc[i+1]]
      evcupd = [evc[1:i-1]; evc[i+2:end]]
      isempty(evcupd) && (evcupd = missing)
      evrupd = evr
   end
   return γ, evrupd, evcupd

# end eigselect2
end

"""
    isqtriu(A::AbstractMatrix) -> Bool

Test whether `A` is a square matrix in a quasi upper triangular form
(e.g., in real or complex Schur form). In the real case, `A` may have 2x2
diagonal blocks, which however must not correspond to complex conjugate eigenvalues.
In the complex case, it is tested if `A` is upper triangular.
"""
function isqtriu(A)
   @assert !Base.has_offset_axes(A)
   m, n = size(A)
   m == n || (return false)
   m == 1 && (return true)
   eltype(A)<:Complex && (return istriu(A))
   m == 2 && (return true)
   istriu(A,-1) || (return false)
   for i = 2:m-1
       !iszero(A[i,i-1]) && !iszero(A[i+1,i]) && (return false)
   end
   return true
end
function eigvalsnosort(M; kwargs...)
   return eigvals(M; sortby=nothing, kwargs...)
end
function eigvalsnosort(M, N; kwargs...)
   ev = eigvals(M, N; sortby=nothing, kwargs...)
   eltype(M) <: Complex || (ev[imag.(ev) .> 0] = conj(ev[imag.(ev) .< 0]))
   return ev
end
function eigvalsnosort!(M; kwargs...)
   return eigvals!(M; sortby=nothing, kwargs...)
end
function eigvalsnosort!(M, N; kwargs...)
   ev = eigvals!(M, N; sortby=nothing, kwargs...)
   eltype(M) <: Complex || (ev[imag.(ev) .> 0] = conj(ev[imag.(ev) .< 0]))
   return ev
end


text:
@doc doc"""
Enumtype including `LET_*`
"""
@enum(LoopElementType,
      LET_Cut,    # [1 1; 1 1]
      LET_FMLink, # [1 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 1]
      LET_AFLink, # [0 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
      LET_Vertex, # [0 0 0 0; 0 1 1 0; 0 1 1 0; 0 0 0 0]
      LET_Cross,  # [1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1]
     )

@doc """
Loop element depicted as

```
|
o
|
```

or matrix

```
1 1 |+>
1 1 |->
```
""" LET_Cut

@doc """
Loop element depicted as

```
|  |
|--|
|  |
```

or matrix

```
1 0 0 0 |++>
0 0 0 0 |+->
0 0 0 0 |-+>
0 0 0 1 |-->
```
""" LET_FMLink

@doc """
Loop element depicted as

```
|  |
|~~|
|  |
```

or matrix

```
0 0 0 0 |++>
0 1 0 0 |+->
0 0 1 0 |-+>
0 0 0 0 |-->
```
""" LET_AFLink

@doc """
Loop element depicted as

```
|  |
~~~~
~~~~
|  |
```

or matrix

```
0 0 0 0 |++>
0 1 1 0 |+->
0 1 1 0 |-+>
0 0 0 0 |-->
```
""" LET_Vertex

@doc raw"""
Loop element depicted as

```
|  |
 \/
 /\
|  |
```

or matrix

```
1 0 0 0 |++>
0 0 1 0 |+->
0 1 0 0 |-+>
0 0 0 1 |-->
```
""" LET_Cross

@doc doc"""
(Imaginary-temporary and spatial) local operator as a perturbation with assigned loop element.
# Fields
- `let_type` : assigned loop element
- `isdiagonal` : operator is diagonal or not
    - in other words, two states connecting this perturbation are equivalent to
      each other or not.
- `time` : imaginary time ($\tau/\beta \in [0,1)$) which this perturbation acts on.
- `space` : spin or bond which this perturbation acts on.
            denotes `space` spin if space <= nspins
                or  `space - nspins` bond otherwise.
- `subspace` : subspin(s) index
- `bottom_id` :: index of node of union find assigned to a loop
- `top_id` :: index of node of union find assigned to the other loop

"""
mutable struct LocalLoopOperator
    let_type :: LoopElementType
    isdiagonal :: Bool
    time :: Float64
    space :: Int
    subspace :: Tuple{Int,Int}
    bottom_id :: Int
    top_id :: Int
end
function LocalLoopOperator(let_type::LoopElementType, time::Real, space::Int, subspace::Tuple{Int,Int})
    LocalLoopOperator(let_type, true, time, space, subspace, 0,0)
end

abstract type QuantumLocalZ2Model <: Model end

@doc doc"""
Spin-$S$ XXZ model represented as the following Hamiltonian,
\begin{equation}
\mathcal{H} = \sum_{i,j} \left[ J_{ij}^z S_i^z S_j^z
            + \frac{J_{ij}^{xy}}{2} (S_i^+ S_j^- + S_i^-S_j^+) \right]
            - \sum_i \Gamma_i S_i^x,
\end{equation}
where $S^x, S^y, S^z$ are $x, y$ and $z$ component of spin operator with length $S$,
and $S^\pm \equiv S^x \pm iS^y$ are ladder operator.
A state is represented by a product state (spins at $\tau=0$) of local $S^z$ diagonal basis and an operator string (perturbations).
A local spin with length $S$ is represented by a symmetrical summation of $2S$ sub spins with length $1/2$.
"""
mutable struct QuantumXXZ <: QuantumLocalZ2Model
    lat :: Lattice
    S2 :: Int
    spins :: Vector{Int}
    ops :: Vector{LocalLoopOperator}
    rng :: Random.MersenneTwister
    parameter :: Parameter

    function QuantumXXZ(lat::Lattice, S::Real)
        if round(2S) != 2S
            error("`S` should be integer or half-integer")
        end
        S2 = round(Int,2S)
        model = new()
        model.rng = seed!(Random.MersenneTwister(0))
        model.lat = lat
        model.S2 = S2
        model.spins = rand(model.rng,[1,-1], numsites(lat)*S2)
        model.ops = LocalLoopOperator[]
        model.parameter = Parameter()
        return model
    end
    function QuantumXXZ(lat::Lattice, S::Real, seed)
        if round(2S) != 2S
            error("`S` should be integer or half-integer")
        end
        S2 = round(Int,2S)
        model = new()
        model.rng = seed!(Random.MersenneTwister(0), seed)
        model.lat = lat
        model.S2 = S2
        model.spins = rand(model.rng,[1,-1], numsites(lat)*S2)
        model.ops = LocalLoopOperator[]
        model.parameter = Parameter()
        return model
    end
end
@doc doc"""
    QuantumXXZ(param)

Generates `QuantumXXZ` using `param["Lattice"]`, `param["S"]` and `param["Seed"]` (if defined).
Each subspin will be initialized independently and randomly.
"""
function QuantumXXZ(param::Parameter, lat::Lattice)
    S = param["S"]
    if "Seed" in keys(param)
        model = QuantumXXZ(lat,S,param["Seed"])
    else
        model = QuantumXXZ(lat,S)
    end
    merge!(model.parameter, param)
    return model
end
QuantumXXZ(param::Parameter) = QuantumXXZ(param, generatelattice(param))

@inline site2subspin(site::Integer, ss::Integer, S2::Integer) = (site-1)*S2+ss
@inline function subspin2site(subspin::Integer, S2::Integer)
    return ceil(Int, subspin/S2), mod1(subspin,S2)
end

ext: jl


text:
# A symbol is a name and the concolicly computed type
struct Sym
    name::Symbol
    _type::DataType
    Sym(base, _type::DataType) = new(Base.gensym(Symbol(base, '#', _type)), _type)
    Sym(_type::DataType) = new(Base.gensym(Symbol(_type)), _type)
end
Cassette.metadatatype(::Type{<:TraceCtx}, ::DataType) = Sym

record(args, ctx) = map(args) do arg
    if istagged(arg, ctx)
        metadata(arg, ctx)
    else
        arg
    end
end

# A callsite is a function with a concolic argument set
# children are all the callsites within the function
mutable struct Callsite
    f::Any
    args::Any
    retval::Union{Some{Any}, Nothing}
    children::Vector{Callsite}
    Callsite(f, args) = new(f, args, nothing, Callsite[])
end
Base.push!(trace::Callsite, call::Callsite) = push!(trace.children, call)

##
# We need to manually override for IntrinsicFunctions (which are the leaf-nodes we are interested in)
# Since in tagged contexts there is an automatic fallback available.
##
function Cassette.execute(ctx::TraceCtx{Callsite, <:Cassette.Tag}, f::Core.IntrinsicFunction, args...)
    call = Callsite(f, record(args, ctx))
    push!(ctx.metadata, call)
    retval = fallback(ctx, f, args...)
    if any(a -> istagged(a, ctx), args)
        retval = tag(retval, ctx, Sym(typeof(retval)))
    end
    call.retval = Some(istagged(retval, ctx) ? metadata(retval, ctx) : ctx)
    return retval
end

##
# Recursivly trace though our program.
#
# Note: This won't trace functions that are defined primitive and there are several
#       fallbacks that Cassette provides for tagged contexts.
#
# Question: Shouldn't canoverdub for primitives be false?
##
function Cassette.execute(ctx::TraceCtx, f, args...)
    # We need to push the callsite first into metadata
    # otherwise we run into issues with functions that throw errors
    call = Callsite(f, record(args, ctx))
    push!(ctx.metadata, call)

    # Recurse into the next function
    retval = if canoverdub(ctx, f, args...)
        newctx = similarcontext(ctx, metadata = call)
        Cassette.overdub(newctx, f, args...)
    else
        retval = fallback(ctx, f, args...)
        # If any of our inputs was tagged we want to tag the return value so that
        # the symbolic nature of this operation is captured.
        if any(a -> istagged(a, ctx), args) && !istagged(retval, ctx)
            retval = tag(retval, ctx, Sym(typeof(retval)))
        end
        retval
    end
    call.retval = Some(istagged(retval, ctx) ? metadata(retval, ctx) : ctx)
    return retval
end


text:
export Sphere1Sphere1, PriorSphere1, PackedSphere1Sphere1, PackedPriorSphere1

"""
$(TYPEDEF)

Factor between two Sphere1 variables.

Related

[`Sphere1`](@ref), [`PriorSphere1`](@ref), [`Polar`](@ref), [`ContinuousEuclid`](@ref)
"""
mutable struct Sphere1Sphere1{T<: SamplableBelief} <: AbstractRelativeRoots
  Z::T
  Sphere1Sphere1(z::T=Normal()) where {T <: SamplableBelief} = new{T}(z)
end
Sphere1Sphere1(::UniformScaling) = Sphere1Sphere1()

getSample(cf::CalcFactor{<:Sphere1Sphere1}, N::Int=1) = (reshape(rand(cf.factor.Z,N),:,N), )

function (cf::CalcFactor{<:Sphere1Sphere1})(meas,
                                            wxi,
                                            wxj  ) # where {M<:FactorMetadata,P<:Tuple,X<:AbstractVector}
  #
  wXjhat = addtheta(wxi[1], meas[1])
  res = difftheta(wxj[1], wXjhat)  # jXjhat =
  return res
end

"""
$(TYPEDEF)

Introduce direct observations on all dimensions of a Sphere1 variable:

Example:
--------
```julia
PriorSphere1( MvNormal([10; 10; pi/6.0], diagm([0.1;0.1;0.05].^2)) )
```

Related

[`Sphere1`](@ref), [`Prior`](@ref), [`PartialPrior`](@ref)
"""
mutable struct PriorSphere1{T<: SamplableBelief} <: AbstractPrior
  Z::T
end

# PriorSphere1(x::T) where {T <: IncrementalInference.SamplableBelief} = PriorSphere1{T}(x)
PriorSphere1(::UniformScaling) = PriorSphere1(Normal())

function getSample(cf::CalcFactor{<:PriorSphere1}, N::Int=1)
  return (reshape(rand(cf.factor.Z,N),:,N), )
end

"""
$(TYPEDEF)

Serialized object for storing PriorSphere1.
"""
mutable struct PackedPriorSphere1  <: IncrementalInference.PackedInferenceType
  datastr::String
  PackedPriorSphere1() = new()
  PackedPriorSphere1(x::String) = new(x)
end
function convert(::Type{PackedPriorSphere1}, d::PriorSphere1)
  return PackedPriorSphere1(convert(PackedSamplableBelief, d.Z))
end
function convert(::Type{PriorSphere1}, d::PackedPriorSphere1)
  distr = convert(SamplableBelief, d.datastr)
  return PriorSphere1{typeof(distr)}(distr)
end

# --------------------------------------------

"""
$(TYPEDEF)

Serialized object for storing Sphere1Sphere1.
"""
mutable struct PackedSphere1Sphere1  <: IncrementalInference.PackedInferenceType
  datastr::String
  PackedSphere1Sphere1() = new()
  PackedSphere1Sphere1(x::String) = new(x)
end
function convert(::Type{Sphere1Sphere1}, d::PackedSphere1Sphere1)
  return Sphere1Sphere1(convert(SamplableBelief, d.datastr))
end
function convert(::Type{PackedSphere1Sphere1}, d::Sphere1Sphere1)
  return PackedSphere1Sphere1(convert(PackedSamplableBelief, d.Z))
end

# --------------------------------------------


text:
# Autogenerated wrapper script for Qt5QuickControls2_jll for i686-linux-gnu-cxx03
export libqt5quickcontrols2, libqt5quicktemplates2, quickcontrols2

using Qt5Declarative_jll
JLLWrappers.@generate_wrapper_header("Qt5QuickControls2")
JLLWrappers.@declare_library_product(libqt5quickcontrols2, "libQt5QuickControls2.so.5")
JLLWrappers.@declare_library_product(libqt5quicktemplates2, "libQt5QuickTemplates2.so.5")
JLLWrappers.@declare_file_product(quickcontrols2)
function __init__()
    JLLWrappers.@generate_init_header(Qt5Declarative_jll)
    JLLWrappers.@init_library_product(
        libqt5quickcontrols2,
        "lib/libQt5QuickControls2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quicktemplates2,
        "lib/libQt5QuickTemplates2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        quickcontrols2,
        "qml/QtQuick/Controls.2",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
Data structure
"""

abstract type AbstractPHSample end

struct WeightedSample{Tv} <: AbstractPHSample
    length::Int
    maxtime::Tv
    tdat::Vector{Tv}
    wdat::Vector{Tv}
end

function _phmean(data::WeightedSample{Tv}) where Tv
    totalt = Tv(0)
    totalw = Tv(0)
    ct = Tv(0)
    for i = eachindex(data.tdat)
        ct += data.tdat[i]
        totalt += data.wdat[i] * ct
        totalw += data.wdat[i]
    end
    return totalt / totalw
end

function WeightedSample(t::Vector{Tv}, w::Vector{Tv}) where Tv
    i = sortperm(t)
    dt,maxt = itime(t[i])
    WeightedSample(length(dt), maxt, dt, w[i])
end

function PointSample(t::Vector{Tv}) where Tv
    i = sortperm(t)
    dt,maxt = itime(t[i])
    WeightedSample(length(dt), maxt, dt, ones(length(dt)))
end

function WeightedSample(f::Any, bounds::Tuple{Tv,Tv}; reltol::Tv = 1.0e-8, abstol::Tv = eps(Tv), d = 8, maxiter = 16) where Tv
    de = deint(f, bounds[1], bounds[2], reltol=reltol, abstol=abstol, d=d, maxiter=maxiter)
    WeightedSample(de.x, de.w * de.h)
end


text:
struct Coordinate{M,R1,S} <: TaskMap{M,R1,S}
    ind::Int
end
struct CoordinateT{M,R1,S} <: TaskMapT{M,R1,S}
    base_map::Coordinate{M,R1,S}
end
CoordinateT{M,R1,S}(i::Int) where {M,R1,S} = CoordinateT{M,R1,S}(Coordinate{M,R1,S}(i))

task_map_emb(::EmbRep, ::EmbRep, xm, task_map::Coordinate{ℝ{m}, R1}) where m =
    SA[xm[task_map.ind]]


text:
assert_result(result) = @assert isa(result, Fit) "`OutlierDetectionInterface.fit` must return a tuple
 of `OutlierDetectionInterface.Model` and `OutlierDetectionInterface.Scores`"

assert_scores(scores) = @assert isa(scores, Scores) "`OutlierDetectionInterface.transform` must return
`OutlierDetectionInterface.Scores`"

function MMI.fit(detector::UnsupervisedDetector, verbosity::Integer, X)
    result = fit(detector, X; verbosity=verbosity); assert_result(result);
    model, scores = result
    return model, nothing, (scores = scores,)
end

function MMI.fit(detector::SupervisedDetector, verbosity::Integer, X, y)
    result = fit(detector, X, y; verbosity=verbosity); assert_result(result);
    model, scores = result
    return model, nothing, (scores = scores,)
end

function MMI.transform(detector::Detector, model::DetectorModel, X)
    scores = transform(detector, model, X); assert_scores(scores);
    return scores
end

# specify scitypes
MMI.input_scitype(::Type{<:Detector}) = Union{AbstractMatrix{<:MMI.Continuous}, MMI.Table(MMI.Continuous)}
MMI.target_scitype(::Type{<:Detector}) = AbstractVector{<:Union{Missing,OrderedFactor{2}}}
MMI.output_scitype(::Type{<:Detector}) = AbstractVector{<:MMI.Continuous}
MMI.transform_scitype(::Type{<:Detector}) = AbstractVector{<:MMI.Continuous}

# helper function to generate colons according to the data dimensionality
ncolons(X::Data) = ntuple(_ -> Colon(), Val(ndims(X) - 1))


text:
config = MOI.Test.TestConfig()
optimize!(mock) = MOIU.mock_optimize!(mock, [0.0, 0.0],
    (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[1.0]],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[1.0]])
for mock in mocks(optimize!)
    Tests.sos_term_test(mock, config)
    Tests.sdsos_term_test(mock, config)
    Tests.dsos_term_test(mock, config)
end


text:
using Games
using Games.Generators
using BenchmarkTools

const SUITE = BenchmarkGroup()

SUITE["support_enumeration"] = BenchmarkGroup(["support_enumeration"])
SUITE["support_enumeration"]["Float"] = BenchmarkGroup()
seed = 0
rng = MersenneTwister(seed)
ns = [10, 11]
for n in ns
    sz = (n, n)
    g = random_game(rng, sz)
    SUITE["support_enumeration"]["Float"][sz] =
        @benchmarkable support_enumeration($g)
end
SUITE["support_enumeration"]["Rational"] = BenchmarkGroup()
T = Rational{Int}
ns = [7, 8]
for n in ns
    sz = (n, n)
    g = NormalFormGame(eye(T, n))
    SUITE["support_enumeration"]["Rational"][sz] =
        @benchmarkable support_enumeration($g)
end

SUITE["bimatrix_generators"] = BenchmarkGroup(["bimatrix_generators"])
seed = 0
rng = MersenneTwister(seed)

# blotto_game
hts = [(3, 62), (4, 21)]  # (h, t)
rho = 0.5
SUITE["bimatrix_generators"]["blotto_game"] = BenchmarkGroup()
for ht in hts
    SUITE["bimatrix_generators"]["blotto_game"][ht] =
        @benchmarkable blotto_game($rng, ($ht)..., $rho)
end

# ranking_game
n = 2000
SUITE["bimatrix_generators"]["ranking_game"] =
    @benchmarkable ranking_game($rng, $n)

# sgc_game
k = 500
SUITE["bimatrix_generators"]["sgc_game"] = @benchmarkable sgc_game($k)

# tournament_game
n, k = 200, 2
SUITE["bimatrix_generators"]["tournament_game"] =
    @benchmarkable tournament_game($n, $k; seed=$seed)

# unit_vector_game
n = 2000
bools = [true, false]
SUITE["bimatrix_generators"]["unit_vector_game"] = BenchmarkGroup()
for b in bools
    SUITE["bimatrix_generators"]["unit_vector_game"][b] =
        @benchmarkable unit_vector_game($rng, $n; avoid_pure_nash=$b)
end


text:
function MaxIter(arg0::jint)
    return MaxIter((jint,), arg0)
end

function get_max_iter(obj::MaxIter)
    return jcall(obj, "getMaxIter", jint, ())
end

function unlimited(::Type{MaxIter})
    return jcall(MaxIter, "unlimited", MaxIter, ())
end


text:
# -------- Computational Element --------------------------------------------- #

"""
Abstract type to encapsulate the common interface, mainly CPU-GPU transfer, for
(i) neural network layers and (ii) neural network models. We distinguish between
the two, since a Jtac `Model` contains game-specific information while a `Layer`
is agnostic of anything but its raw data input and output sizes.
"""
abstract type Element{GPU} end

"""
  swap(element)

Swap a model or layer `element` from/to gpu memory.
"""
swap(l :: Element) = error("Not implemented")

"""
  on_gpu(element)

Check whether `element` currently resides on the GPU.
"""
on_gpu(:: Element{GPU}) where {GPU} = GPU
on_gpu(:: Nothing) = false

"""
  to_cpu(element)

Move `element` to the CPU. Does nothing if the element is already there.
"""
to_cpu(el :: Element{false}) = el
to_cpu(el :: Element{true}) = swap(el)

"""
  to_gpu(element)

Try to move `element` to the GPU. Prints a warning if GPU is not available. Does
nothing if the element is already on the GPU.
"""
to_gpu(el :: Element{true})  = el

function to_gpu(el :: Element{false})
  if Knet.gpu() == -1
    @warn "No GPU was found by Knet. Element stays on CPU"
    el
  else
   swap(el)
  end
end

Base.copy(l :: Element) :: Element = error("Not implemented")

max_stars_repo_name: roSievers/jtac


text:
using Flux: onehotbatch
using Flux.Data: DataLoader, MNIST

"""
    prepare_data(x, y, [isflat=false])
"""
function prepare_data(x, y; isflat=false)
    xx = Array{Float32}(undef, size(x[1])..., 1, length(x))
    for i in 1:length(x)
        xx[:, :, :, i] = Float32.(x[i])
    end

    xx = isflat ? Flux.flatten(xx) : xx
    y = onehotbatch(y, 0:9)
    return DataLoader(xx, y, batchsize=10, shuffle=true, partial=true)
end

"""
    get_traindata(indices, [isflat=false])

Prepare the MNIST training data for the specified range.
"""
function get_traindata(indices::UnitRange; isflat=false)
    xtrain = Flux.Data.MNIST.images()[indices]
    ytrain = Flux.Data.MNIST.labels()[indices]
    return prepare_data(xtrain, ytrain, isflat=isflat)
end

"""
    get_testdata([isflat=false])

Prepare the MNIST testing data for the specified range.
"""
function get_testdata(; isflat=false)
    xtest = Flux.Data.MNIST.images(:test)
    ytest = Flux.Data.MNIST.labels(:test)
    return prepare_data(xtest, ytest, isflat=isflat)
end

"""
    getdata(isflat=false)

Prepare the MNIST dataset.
"""
function getdata(isflat=false)
    traindata = get_traindata(1:60000, isflat=isflat)
    testdata = get_testdata(isflat=isflat)
    return traindata, testdata
end


text:
##Operators

for op in (:Derivative,:Integral)
    @eval begin
        function ($op){T<:IntervalDomain}(d::Vector{T})
            n=length(d)
            R=zeros(Operator{mapreduce(eltype,promote_type,d)},n,n)
            for k=1:n
                R[k,k]=$op(d[k])
            end

            R
        end
    end
end

function Evaluation{T<:IntervalDomain}(d::Vector{T},x...)
    n=length(d)
    R=zeros(Operator{mapreduce(eltype,promote_type,d)},n,n)
    for k=1:n
        R[k,k]=Evaluation(d[k],x...)
    end

    R
end

## Construction

function Base.diagm{T<:Operator}(d::Vector{T})
    D=zeros(Operator{mapreduce(eltype,promote_type,d)},length(d),length(d))
    for k=1:length(d)
        D[k,k]=d[k]
    end
    D
end

##TODO: unify with other blkdiag
function Base.blkdiag{T<:Operator}(d1::Vector{T},d2::Vector{T})
  if isempty(d1)&&isempty(d2)
    error("Empty blkdiag")
  end
  if isempty(d1)
    TT=mapreduce(eltype,promote_type,d2)
  elseif isempty(d2)
    TT=mapreduce(eltype,promote_type,d1)
  else
     TT=promote_type(mapreduce(eltype,promote_type,d1),
                    mapreduce(eltype,promote_type,d2))
  end

  D=zeros(Operator{TT},length(d1)+length(d2),2)
  D[1:length(d1),1]=d1
  D[length(d1)+1:end,2]=d2
  D
end

Base.blkdiag(a::Operator,b::Operator) = blkdiag(Operator{promote_type(eltype(a),eltype(b))}[a],
                                                Operator{promote_type(eltype(a),eltype(b))}[b])

## broadcase

.*{N<:Number}(A::Array{N},D::Operator)=Operator{promote_type(N,eltype(D))}[A[k,j]*D for k=1:size(A,1),j=1:size(A,2)]
.*{N<:Number}(D::Operator,A::Array{N})=A.*D


text:
__precompile__()

module PisingerKnapsack

export bouknap, minknap, mulknap

const deps_file = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(deps_file)
    include(deps_file)
else
    error("PisingerKnapsack not properly installed. Please run import Pkg; Pkg.build(\"PisingerKnapsack\")")
end

include("Cwrapper.jl")
include("doubleknp.jl")

end # module


text:
module RLCore

using Random

greet() = println("Hello Reinforcement Learning!")

export
    AbstractEnvironment,
    start!,
    step!,
    get_actions,
    get_reward,
    is_terminal
include("environment.jl")

export
    AbstractAgent,
    start!,
    step!
include("agent.jl")

export
    Episode,
    run_episode!
include("episode.jl")

export
    GVF,
    Horde,
    GVFParamFuncs
include("gvf.jl")

export
    AbstractFeatureConstructors,
    create_features,
    feature_size
include("feature_constructors.jl")

# export
#     TileCoder,
#     MinMaxNormalize,
#     MeanStdNormalize
# include("features/TileCoder.jl")
# include("features/Normalize.jl")

# export
#     RandomAgent
# include("agent/random.jl")

# Not tested yet so not exported!
include("monte_carlo.jl")

end # module


text:
import .ExpressionExplorer
import .ExpressionExplorer: join_funcname_parts, FunctionNameSignaturePair
import REPL: ends_with_semicolon

"Return a copy of `old_topology`, but with recomputed results from `cells` taken into account."
function updated_topology(old_topology::NotebookTopology, notebook::Notebook, cells)

	updated_codes = Dict{Cell,ExprAnalysisCache}()
	for cell in cells
		if !(
			haskey(old_topology.codes, cell) &&
			old_topology.codes[cell].code === cell.code
		)
			updated_codes[cell] = ExprAnalysisCache(notebook, cell)
		end
	end
	new_codes = merge(old_topology.codes, updated_codes)

	updated_nodes = Dict{Cell,ReactiveNode}(cell => (
			new_codes[cell].parsedcode |>
			ExpressionExplorer.try_compute_symbolreferences |>
			ReactiveNode
		) for cell in cells)

	new_nodes = merge(old_topology.nodes, updated_nodes)

	# DONE (performance): deleted cells should not stay in the topology
	for removed_cell in setdiff(keys(old_topology.nodes), notebook.cells)
		delete!(new_nodes, removed_cell)
		delete!(new_codes, removed_cell)
	end

	NotebookTopology(nodes=new_nodes, codes=new_codes)
end


text:
"""
    module FMM3D
Wrappers for the Flatiron Institute's FMM3D
library.
# wrappers
All N-body codes return output in an
`FMMVals` structure. See documentation
of N-body codes for details.
N-body interactions with the Helmholtz kernel
- [`hfmm3d`](@ref): ``O(N)`` fast mutlipole code
- [`h3ddir`](@ref): ``O(N^2)`` direct code
N-body interactions with the Laplace kernel
- [`lfmm3d`](@ref): ``O(N)`` fast mutlipole code
- [`l3ddir`](@ref): ``O(N^2)`` direct code
# lower level routines
For a list of lower level routines see the
documentation for [`lower_level_routs`](@ref)
"""
module FMM3D

using FMM3D_jll
import Base: show, propertynames

export hfmm3d, lfmm3d, h3ddir, l3ddir, lower_level_routs
export besseljs3d

# fortran input/return types

Fd = Ref{Float64}
Fi = Ref{Int32}
Fc = Ref{ComplexF64}

# common input types

TFN = Union{Array{Float64},Nothing}
TCN = Union{Array{ComplexF64},Nothing}

"""
    FMMVals
Abstract return type for FMM and direct computation
function calls. Fields are `nothing` on return
unless requested. See individual FMM/direct
computation function documentation for specifics.
"""
abstract type FMMVals end

include("helmholtz.jl")
include("laplace.jl")
include("stokes.jl")
include("maxwell.jl")
include("lower_level_routines.jl")
include("properties.jl")

end # module


text:
using ADatasets
using Iterators
using StatsBase
using EvalCurves
using Iterators
using Lazy
using DataFrames
using SemiSupervised
using FluxExtensions
using Flux
using IPMeasures
using MLDataPattern

import ADatasets: runtest, surveydatasets
import FluxExtensions: layerbuilder, freeze

idir = filter(isdir,["/Users/tpevny/Work/Data/datasets/numerical","/mnt/output/data/datasets/numerical","/opt/output/data/datasets/numerical"])[1];
odir = filter(isdir,["/Users/tpevny/Work/Julia/results/datasets","/mnt/output/results/datasets","/opt/output/results/datasets"])[1];

function fit(data,layers,hidden,zdim,β)
  idim = size(first(data)[1],1)
  σ2 = SemiSupervised.bandwidthofkde(randn(zdim,1000),20)
  m = SemiSupervised.VAE(layerbuilder(idim,hidden,2*zdim,layers,"relu","linear","Dense"),
      layerbuilder(zdim,hidden,idim,layers,"relu","linear","Dense"),β)
  m = to32(m);
  mf = freeze(m)
  opt = Flux.Optimise.ADAM(params(m));

  FluxExtensions.learn((x) -> SemiSupervised.loss(m,getobs(x)),opt,data,() -> (),10000)

  # let's estimate the distance between the prior and posterior using MMD (zmmd)
  # and calibration distance on prior and prior
  x = IPMeasures.samplecolumns(data.data[1],1000)
  n = size(x,2)
  z = SemiSupervised.gaussiansample(SemiSupervised.hsplitsoftp(mf.q(x))...)
	kσ2 = SemiSupervised.bandwidthofkde(z,20)
	kde = x -> SemiSupervised.kde(x,z,kσ2)
  df = DataFrame(
		layers = layers,
		hidden = hidden,
		zdim = zdim,
		β = β,
		zmmd = IPMeasures.mmdg(z, randn(zdim,n), 0.5/σ2),
		ummd = IPMeasures.mmdg(randn(zdim,n), randn(zdim,n), 0.5/σ2),
		recerror =  Flux.mse(x,mf.g(z)) / n,
		likelihood =  - mean(SemiSupervised.log_normal(x, mf.g(z), β)),
		train_px = mean(SemiSupervised.px(mf,x, 100,mf.β)),
		train_pxis = mean(SemiSupervised.pxis(mf, x , 100, mf.β)),
		train_kde_pxis = mean(SemiSupervised.pxis(mf, x, kde, 100, mf.β)),
		train_pxvita = mean(SemiSupervised.pxvita(mf,x,mf.β))
		)
  (mf, kde), df
end

predicts = [(m, x) -> - SemiSupervised.px(m[1],x,100,m[1].β)[:],
	(m, x) -> - SemiSupervised.pxis(m[1],x,100, m[1].β)[:],
	(m, x) -> - SemiSupervised.pxis(m[1], x, m[2], 100, m[1].β)[:],
	(m, x) -> - SemiSupervised.pxvita(m[1],x,m[1].β)[:]]
prnames = ["px","pxis","kde_pxis","vita"]
ps = [(3,2^zdim,2^hidden,β) for zdim in 2:5 for hidden in 1:zdim for β in [0.25,0.5,1.0,2.0]]

datasets = (length(ARGS) == 0) ? ADatasets.surveydatasets : ARGS
map(d -> runtest(fit, ps, predicts, prnames, d[1],"easy",0.05,"low",idir, odir, "mvae", d[2],10000),product(datasets,1:10))


text:
function mismatchedreductions_noturbo!(
  𝛥r392,
  𝛥x923,
  𝛥ℛ,
  ℛ,
  r392,
  x923,
  𝒶𝓍k = 1:2,
  𝒶𝓍n = 1:3,
  𝒶𝓍j = 1:9,
  𝒶𝓍m = 1:9,
  𝒶𝓍i = 1:3,
)
  @inbounds @fastmath for k in 𝒶𝓍k
    for i in 𝒶𝓍i
      for m in 𝒶𝓍m
        for j in 𝒶𝓍j
          for n in 𝒶𝓍n
            ℰ𝓍1 = conj(x923[m, k, n])
            ℰ𝓍2 = 𝛥ℛ[n, j, m, i] * ℰ𝓍1
            ℰ𝓍3 = conj(r392[i, j, k])
            ℰ𝓍4 = 𝛥ℛ[n, j, m, i] * ℰ𝓍3
            𝛥r392[i, j, k] = 𝛥r392[i, j, k] + ℰ𝓍2
            𝛥x923[m, k, n] = 𝛥x923[m, k, n] + ℰ𝓍4
          end
        end
      end
    end
  end
  𝛥r392, 𝛥x923
end
function mismatchedreductions!(
  𝛥r392,
  𝛥x923,
  𝛥ℛ,
  ℛ,
  r392,
  x923,
  𝒶𝓍k = 1:2,
  𝒶𝓍n = 1:3,
  𝒶𝓍j = 1:9,
  𝒶𝓍m = 1:9,
  𝒶𝓍i = 1:3,
)
  @turbo for k in 𝒶𝓍k
    for i in 𝒶𝓍i
      for m in 𝒶𝓍m
        for j in 𝒶𝓍j
          for n in 𝒶𝓍n
            ℰ𝓍1 = conj(x923[m, k, n])
            ℰ𝓍2 = 𝛥ℛ[n, j, m, i] * ℰ𝓍1
            ℰ𝓍3 = conj(r392[i, j, k])
            ℰ𝓍4 = 𝛥ℛ[n, j, m, i] * ℰ𝓍3
            𝛥r392[i, j, k] = 𝛥r392[i, j, k] + ℰ𝓍2
            𝛥x923[m, k, n] = 𝛥x923[m, k, n] + ℰ𝓍4
          end
        end
      end
    end
  end
  𝛥r392, 𝛥x923
end

@testset "Many Loop Reductions" begin
  @show @__LINE__
  A = rand((2:6)...)
  N = ndims(A)
  T = eltype(A)
  let dims = (3, 5)
    sᵢ = size(A)
    sₒ = ntuple(Val(N)) do d
      ifelse(d ∈ dims, 1, sᵢ[d])
    end
    Tₒ = Base.promote_op(+, T, Int)
    B = similar(A, Tₒ, sₒ)

    Bᵥ = view(B, Colon(), Colon(), firstindex(B, 3), Colon(), firstindex(B, 5))
    @turbo for i_1 in indices((A, B), 1)
      for i_2 in indices((A, B), 2)
        for i_4 in indices((A, B), 4)
          Σ = zero(eltype(Bᵥ))
          for i_3 in axes(A, 3)
            for i_5 in axes(A, 5)
              Σ += A[i_1, i_2, i_3, i_4, i_5]
            end
          end
          Bᵥ[i_1, i_2, i_4] = Σ
        end
      end
    end
    @test B ≈ sum(A, dims = dims)
  end
  let dims = (1, 2, 5)

    sᵢ = size(A)
    sₒ = ntuple(Val(N)) do d
      ifelse(d ∈ dims, 1, sᵢ[d])
    end
    Tₒ = Base.promote_op(+, T, Int)
    B = similar(A, Tₒ, sₒ)

    Bᵥ = view(B, firstindex(B, 1), firstindex(B, 2), Colon(), Colon(), firstindex(B, 5))
    @turbo for i_3 in indices((A, B), 3)
      for i_4 in indices((A, B), 4)
        Σ = zero(eltype(Bᵥ))
        for i_1 in axes(A, 1)
          for i_2 in axes(A, 2)
            for i_5 in axes(A, 5)
              Σ += A[i_1, i_2, i_3, i_4, i_5]
            end
          end
        end
        Bᵥ[i_3, i_4] = Σ
      end
    end
    @test B ≈ sum(A, dims = dims)
  end

  r392 = rand(3, 9, 2)
  x923 = rand(9, 2, 3)
  K = rand(3, 9, 9, 3)
  𝛥r392_1, 𝛥x923_1, 𝛥r392_2, 𝛥x923_2, 𝛥ℛ =
    similar(r392), similar(x923), similar(r392), similar(x923), copy(K)
  𝛥r392_1 .= -1
  𝛥x923_1 .= -1
  𝛥r392_2 .= -1
  𝛥x923_2 .= -1

  mismatchedreductions_noturbo!(𝛥r392_1, 𝛥x923_1, 𝛥ℛ, K, r392, x923)
  @time mismatchedreductions!(𝛥r392_2, 𝛥x923_2, 𝛥ℛ, K, r392, x923)
  @test 𝛥r392_1 ≈ 𝛥r392_2
  @test 𝛥x923_1 ≈ 𝛥x923_2
end


text:
abstract type Move end

"""
    FastMove(moveType, stab, power, energy, cooldown)

Struct for holding fast moves that holds all information that determines
mechanics: type, STAB (same type attack bonus), power, energy, and cooldown.
Note that this is agnostic to the identity of the actual move itself.
"""
struct FastMove <: Move
    move_type::UInt8
    data::UInt16
end

"""
    FastMove(gm_move, types)

Generate a fast move from the gamemaster entry of the move, and the types of
the mon using it (to determing STAB, same type attack bonus). This should only
be used internally, as generating the move from the name is a lot cleaner.
"""
function FastMove(gm_move::Dict{String,Any}, types::Tuple{UInt8, UInt8})
    move_type = UInt8(findfirst(x -> typings[x] == gm_move["type"], 1:19))
    STAB = move_type in types ? 0x0001 : 0x0000
    power = UInt16(gm_move["power"])
    energy = UInt16(gm_move["energyGain"])
    cooldown = UInt16(gm_move["cooldown"] ÷ 500)
    return FastMove(
        move_type,
        cooldown + energy << 3 + power << 8 + STAB << 13
    )
end

"""
    FastMove(move_name, types)

Generate a fast move from the name of the move, and the types of the mon using
it (to determing STAB, same type attack bonus)
"""
function FastMove(move_name::String, types::Tuple{UInt8, UInt8})
    move_index = findfirst(isequal(move_name), map(x ->
        gamemaster["moves"][x]["moveId"], 1:length(gamemaster["moves"])))
    gm_move = gamemaster["moves"][move_index]
    return FastMove(gm_move, types)
end

get_STAB(fm::FastMove) = iszero(fm.data >> 13) ? 10 : 12
get_power(fm::FastMove) = UInt8((fm.data >> 8) & 0x001f)
get_energy(fm::FastMove) = UInt8((fm.data >> 3) & 0x001f)
get_cooldown(fm::FastMove) = UInt8(fm.data & 0x0007)

"""
    ChargedMove(moveType, stab, power, energy, buffChance, opp_buffs, self_buffs)

Struct for holding charged moves that holds all information that determines
mechanics: type, STAB (same type attack bonus), power, energy, and buff
information (chance, which buffs are applied and to whom). Note that this is
agnostic to the identity of the actual move itself.
"""
struct ChargedMove <: Move
    move_type::UInt8
    buff::StatBuffs
    data::UInt16
end

"""
    ChargedMove(gm_move, types)

Generate a charged move from the gamemaster entry of the move, and the types of
the mon using it (to determing STAB, same type attack bonus). This should only
be used internally, as generating the move from the name is a lot cleaner.
"""
function ChargedMove(gm_move::Dict{String,Any}, types::Tuple{UInt8, UInt8})
    move_type = UInt8(findfirst(x -> typings[x] == gm_move["type"], 1:19))
    STAB = move_type in types ? 0x0001 : 0x0000
    buff_target = (haskey(gm_move, "buffs") &&
        gm_move["buffTarget"] == "opponent") ? 0x0001 : 0x0000
    power = UInt16(gm_move["power"] ÷ 5)
    energy = UInt16(gm_move["energy"] ÷ 5)
    buff = !haskey(gm_move, "buffs") ? defaultBuff :
        StatBuffs(Int8(gm_move["buffs"][1]), Int8(gm_move["buffs"][2]))
    buff_chance = !haskey(gm_move, "buffs")  ? 0x0000 :
        gm_move["buffApplyChance"] == ".1"   ? 0x0001 :
        gm_move["buffApplyChance"] == ".125" ? 0x0002 :
        gm_move["buffApplyChance"] == ".2"   ? 0x0003 :
        gm_move["buffApplyChance"] == ".3"   ? 0x0004 :
        gm_move["buffApplyChance"] == ".5"   ? 0x0005 : 0x0006

    return ChargedMove(
        move_type,
        buff,
        power + energy << 6 + buff_chance << 10 + buff_target << 11 + STAB << 12
    )
end

"""
    ChargedMove(move_name, types)

Generate a charged move from the name of the move, and the types of the mon
using it (to determing STAB, same type attack bonus)
"""
function ChargedMove(move_name::String, types)
    #if move_name == "NONE"
    #    return ChargedMove(Int8(0), Int8(0), UInt8(0), Int8(0), Int8(0),
    #        defaultBuff, defaultBuff)
    #end
    move_index = findfirst(isequal(move_name), map(x ->
        gamemaster["moves"][x]["moveId"], 1:length(gamemaster["moves"])))
    gm_move = gamemaster["moves"][move_index]
    return ChargedMove(gm_move, types)
end

get_power(cm::ChargedMove) = UInt8(0x0005 * (cm.data & 0x003f))
get_energy(cm::ChargedMove) = UInt8(0x0005 * ((cm.data >> 6) & 0x000f))
get_buff_target(cm::ChargedMove) = UInt8((cm.data >> 11) & 0x0001)
get_STAB(cm::ChargedMove) = iszero(cm.data >> 12) ? 10 : 12

function get_buff_chance(cm::ChargedMove)
    buff_chance = (cm.data >> 10) & 0x03ff
    return buff_chance == 0x0000 ? 0.0   :
           buff_chance == 0x0006 ? 1.0   :
           buff_chance == 0x0001 ? 0.1   :
           buff_chance == 0x0002 ? 0.125 :
           buff_chance == 0x0003 ? 0.2   :
           buff_chance == 0x0004 ? 0.3   : 0.5
end

function buff_applies(cm::ChargedMove)
    buff_chance = (cm.data >> 10) & 0x03ff
    return buff_chance == 0x0000 ? false :
        buff_chance == 0x0006 ? true :
        buff_chance == 0x0001 ? rand(rb_rng, 1:10) == 1 :
        buff_chance == 0x0002 ? rand(rb_rng,
            (true, false, false, false, false, false, false, false)) :
        buff_chance == 0x0003 ? rand(rb_rng, 1:5) == 1 :
        buff_chance == 0x0004 ? rand(rb_rng, 1:10) < 4 :
        rand(rb_rng, (true, false))
end


text:
module vanderpol
using ReachabilityAnalysis
@taylorize function vanderpol!(dx, x, params, t)
    local μ = 1.0
    dx[1] = x[2]
    dx[2] = (μ * x[2]) * (1 - x[1]^2) - x[1]
    return dx
end

function model(X0)
    S = @system(x' = vanderpol!(x), dim:2)
    return IVP(S, X0)
end
end # module

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
# sudo apt-get install python-dev
# sudo apt-get install python-pip
# sudo pip install biopython
# sudo apt-get install numpy | unless you want to compile it..
# julia: Pkg.add("PyCall")
module PDB
export get_structure, get_chains, structure_to_matrix, export_to_pdb, get_remote_pdb

using PyCall
using Formatting
@pyimport Bio.PDB as pdb

# get a structure from a PDB File
function get_structure(filename::AbstractString)
	# using QUIET to supress discontinous chains warning
	pdbparser = pdb.PDBParser(QUIET = 1)

	# parse the file
	structure = pdbparser[:get_structure](filename, filename)
	return structure
end

# get chains from structure
function get_chains(structure::PyObject)
	chains = structure[:get_chains]()

	chainMatrices = Any[]

	for c in chains
		push!(chainMatrices,structure_to_matrix(c))
	end

	return chainMatrices
end

# Read in a PDB file and return a matrix of C_Alpha atom coordinates
function structure_to_matrix(structure::PyObject)
	atoms = structure[:get_atoms]()
	# Filter out C_Alpha atoms
	atoms = filter(x -> x[:get_name]() == "CA", atoms)
	# Get the coordinates
	atoms = map(x -> x[:get_coord](), atoms)

	n = length(atoms)

	# initialize an empty Nx3 matrix
	matrix = zeros(n,3)

	# fill it
	for i in 1:n
		matrix[i,:] = [atoms[i][1] atoms[i][2] atoms[i][3]]
	end

	matrix = reshape(matrix,n,3)

	return matrix
end

# Export a matrix of C_alpha atom coordinates to a PDB file
function export_to_pdb(residueName::AbstractString,chainID::AbstractString,matrix::Array{Float64,2}, filename::AbstractString)
	lines = AbstractString[]

	atomExpr = FormatExpr("{: <6}{: >5} {: >4}{: <1}{: >3} {: <1}{: >4}{: <1}   {: >8}{: >8}{: >8}{: >6}{: >6}      {: <4}{: >2}{: <2}\n")

	for i in 1:size(matrix)[1]
		line = format(atomExpr, "ATOM", i, "CA", " ", residueName, chainID, "1", " ", matrix[i,:][1], matrix[i,:][2], matrix[i,:][3], 1.0, 0.0, "A1", "C", " ")

		push!(lines, line)
	end

	f = open(filename, "w")
	write(f,lines)
	close(f)

	return lines
end

# Get a PDB file by ID from rcsb.org
function get_remote_pdb(id::AbstractString)
	# create cache if not present
	if !isdir(joinpath(dirname(@__FILE__), "..", "..", ".cache"))
		mkdir(joinpath(dirname(@__FILE__), "..", "..", ".cache"))
	end
	filename = joinpath(dirname(@__FILE__), "..", "..", ".cache", id)
	# check if PDB isn't already cached
	if !isreadable(filename)
		data = get("http://www.rcsb.org/pdb/files/" * id * ".pdb")
		f = open(filename, "w")
		write(f, data.data)
	end
	return filename
end

end


text:
#! /usr/bin/julia

# Rosetta Code, Multifactorial

function multifact{T<:Integer,U<:Integer}(n::T, k::U)
    -1<n && 0<k || throw(DomainError())
    1 < k || return factorial(n)
    r = one(T)
    for i in n:-k:2
        r *= i
    end
    return r
end

khi = 5
nhi = 10
println("Showing multifactorial for n in [1,", nhi, "] and k in [1,", khi, "].")
for k = 1:khi
    a = Int64[multifact(i, k) for i in 1:nhi]
    lab = "n"*"!"^k
    println(@sprintf("  %-6s =>  ", lab), a)
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code