text:
"""
    poly2pm(PM; grade = k) -> P

Build a grade `k` matrix polynomial representation `P(λ)` from a polynomial matrix, polynomial vector or scalar polynomial `PM(λ)`.

`PM(λ)` is a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

`P(λ)` is a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`, for which
the coefficient matrices `P_l`, `l = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,l]` contains the `l`-th coefficient matrix `P_l` (multiplying `λ**(l-1)`).
If `grade = missing`, then `k` is chosen the largest degree of the elements of `PM`.
The coefficients of the degree `d` element `(i,j)` of `PM(λ)` result in `P[i,j,1:d+1]`.
"""
function poly2pm(PM::Matrix{<:Polynomial}; grade::Union{Int,Missing} = missing)
   p, m = size(PM)
   degs = degree.(PM)
   d = maximum(degs)
   ismissing(grade) ? k = d+1 : k = max(d,grade)+1
   T = eltype(eltype(PM))
   k == 0 && (return zeros(T,p,m,1))
   P = zeros(T,p,m,k)
   for j = 1:m
      for i = 1:p
         degs[i,j] < 0 || (P[i,j,1:degs[i,j]+1] = coeffs(PM[i,j]))
      end
   end
   return P
end
function poly2pm(PM::Matrix{T}; grade::Union{Int,Missing} = missing) where T <: Number
   p, m = size(PM)
   d = 0
   ismissing(grade) ? k = d+1 : k = max(d,grade)+1
   k == 0 && (return zeros(T,p,m,1))
   if k == 1
      P = reshape(PM,p,m,1)
   else
      P = zeros(T,p,m,k)
      P[:,:,1] = PM
   end
   return P
end
function poly2pm(PM::Vector{<:Polynomial}; grade::Union{Int,Missing} = missing)
   m = length(PM)
   degs = degree.(PM)
   d = maximum(degs)
   ismissing(grade) ? k = d+1 : k = max(d,grade)+1
   T = eltype(eltype(PM))
   k == 0 && (return zeros(T,m,1,1))
   P = zeros(T,m,1,k)
   for i = 1:m
      degs[i] < 0 || (P[i,1,1:degs[i]+1] = coeffs(PM[i]))
   end
   return P
end
function poly2pm(PM::Vector{T}; grade::Union{Int,Missing} = missing) where T <: Number
   m = length(PM)
   d = 0
   ismissing(grade) ? k = d+1 : k = max(d,grade)+1
   k == 0 && (return zeros(T,m,1,1))
   if k == 1
      P = reshape(PM,m,1,1)
   else
      P = zeros(T,m,1,k)
      P[:,:,1] = PM
   end
   return P
end
function poly2pm(PM::Union{Adjoint{T,Vector{T}},Transpose{T,Vector{T}}}; grade::Union{Int,Missing} = missing) where T <: Number
   m = length(PM)
   d = 0
   ismissing(grade) ? k = d+1 : k = max(d,grade)+1
   k == 0 && (return zeros(T,1,m,1))
   adj = typeof(PM) <: Adjoint
   if k == 1
      adj ? P = reshape(conj(PM.parent),1,m,1) : P = reshape(PM.parent,1,m,1)
   else
      P = zeros(T,1,m,k)
      adj ? P[:,:,1] = reshape(conj(PM.parent),1,m,1) : P[:,:,1] = reshape(PM.parent,1,m,1)
   end
   return P
end
function poly2pm(PM::Polynomial{T}; grade::Union{Int,Missing} = missing) where T
   d = degree(PM)
   ismissing(grade) ? k = d+1 : k = max(d,grade)+1
   k == 0 && (return zeros(T,1,1,1))
   P = zeros(T,1,1,k)
   d < 0 || (P[1,1,1:d+1] = coeffs(PM))
   return P
end
function poly2pm(PM::Number; grade::Union{Int,Missing} = missing)
   T = typeof(PM)
   PM == zero(T) ? d = -1 : d = 0
   ismissing(grade) ? k = d+1 : k = max(d,grade)+1
   k == 0 && (return zeros(T,1,1,1))
   P = zeros(T,1,1,k)
   d < 0 || (P[1,1,1] = PM)
   return P
end
"""
    pm2poly(P[,var = 'x']) -> PM

Build the polynomial matrix `PM(λ)` from its matrix polynomial representation `P(λ)`.

`P(λ)` is a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`, for which
the coefficient matrices `P_l`, `l = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,l]` contains the `l`-th coefficient matrix `P_l` (multiplying `λ**(l-1)`).

`PM(λ)` is a matrix of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
The element `(i,j)` of `PM(λ)` is built from the coefficients contained in `P[i,j,1:k+1]`.
The symbol to be used for the indeterminate `λ` can be specified in the optional input variable `var`.
"""
function pm2poly(PM::AbstractArray{T,3},var::Union{Char, AbstractString, Symbol}='x') where T
   m, n, k = size(PM)
   P = zeros(Polynomial{T},m,n)
   for i = 1:m
      for j = 1:n
          P[i,j] = Polynomial(PM[i,j,1:k],var)
      end
   end
   return P
end
"""
    pmdeg(P) -> deg

Determine the degree `deg` of a polynomial matrix `P(λ)`.

`P(λ)` is a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`, for which
the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).
The degree of `P(λ)` is `deg = j-1`, where `j` is the largest index for which `P[:,:,j]` is nonzero. The degree of
the zero polynomial matrix is defined to be `deg = -1`.

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
The degree of `P(λ)` is the largest degree of the elements of `P(λ)`.
The degree of the zero polynomial matrix is defined to be `-1`.
"""
function pmdeg(P::AbstractArray{T,3}) where T
   for j = size(P,3):-1:1
       norm(P[:,:,j],Inf) > 0 && return j-1
   end
   return -1
end
function pmdeg(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number})
   typeof(P) <: Number && (P == 0 ? (return -1) : (return 0) )
   return maximum(degree.(P))
end
"""
    pmeval(P,val) -> R

Evaluate `R = P(val)` for a polynomial matrix `P(λ)`, using Horner's scheme.

`P(λ)` is a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`, for which
the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
"""
function pmeval(P::AbstractArray{T,3},val::Number) where {T}
   # Horner's scheme
   p, m, k1 = size(P)
   nd = pmdeg(P)+1
   S = typeof(val)
   nd == 0 && return zeros(promote_type(T,S),p,m)
   R = P[:,:,nd]*one(S)
   for k = nd-1:-1:1
      R = R*val+ P[:,:,k]
   end
   return R
end
pmeval(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial},val::Number) = pmeval(poly2pm(P),val)
pmeval(P::Union{Number,AbstractVecOrMat{<:Number}},val::Number = 0) = P
"""
    pmreverse(P[,j]) -> Q

Build `Q(λ) = λ^j*P(1/λ)`, the `j`-reversal of a polynomial matrix `P(λ)` for `j ≥ deg(P(λ))`.
If `j` is not specified, the default value `j = deg(P(λ))` is used.

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

If deg(P(λ)), then `Q(λ)` is a grade `j` polynomial matrix of the form
`Q(λ) = Q_1 + λ Q_2 + ... + λ**j Q_(j+1)`, for which
the coefficient matrices `Q_i`, `i = 1, ..., j+1`, are stored in the 3-dimensional matrix `Q`,
where `Q[:,:,i]` contains the `i`-th coefficient matrix `Q_i` (multiplying `λ**(i-1)`).
The coefficient matrix `Q_i` is either `0` if `i ≤ j-d` or `Q_(j-d+i) = P_(d-i+1)` for `i = 1, ..., d`.
"""
function pmreverse(P::AbstractArray{T,3}, j::Int = pmdeg(P)) where T
   d = pmdeg(P)
   j < d && error("j must be at least $d")
   m, n, k1 = size(P)
   Q = zeros(eltype(P),m,n,j+1)
   Q[:,:,j-d+1:j+1] = reverse(P[:,:,1:d+1],dims=3)
   return Q
end
pmreverse(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number}; kwargs...) =
      pmreverse(poly2pm(P); kwargs...)
"""
    pmdivrem(N,D) -> (Q, R)

Compute the quotients in `Q(λ)` and remainders in `R(λ)` of the elementwise polynomial divisions `N(λ)./D(λ)`.

`N(λ)` is a polynomial matrix of the form `N(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`, for which
the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `N`,
where `N[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

`D(λ)` is a polynomial matrix of the form `D(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`, for which
the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `D`,
where `D[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

Alternatively, `N(λ)` and `D(λ)` can be specified as matrices of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The polynomial matrices of quotients `Q(λ)` and remainders `R(λ)` are stored in the 3-dimensional
matrices `Q` and  `R`, respectively, where `Q[:,:,i]` and `R[:,:,i]` contain the `i`-th
coefficient matrix multiplying `λ**(i-1)`.
"""
function pmdivrem(N::AbstractArray{T1,3},D::AbstractArray{T2,3}) where {T1,T2}
    p, m, k = size(N)
    p1, m1, k1 = size(D)
    (p,m) == (p1, m1) || error("Numerator and denominator polynomial matrices must have the same size")

    degQ1 = 0
    degD1 = 0
    for j = 1:m
        for i = 1:p
            n1 = poldeg1(N[i,j,1:k])
            d1 = poldeg1(D[i,j,1:k1])
            d1 == 0 && error("DivideError: zero denominator polynomial")
            n1 < d1 || (degQ1 = max(degQ1,n1-d1+1))
            degD1 = max(degD1,d1)
        end
    end
    T = eltype(one(T1)/one(T2))
    R = zeros(T,p,m,max(degD1-1,1))
    Q = zeros(T,p,m,max(degQ1,1))
    for j = 1:m
        for i = 1:p
            q, r = poldivrem(N[i,j,1:k],D[i,j,1:k1])
            nq1 = poldeg1(q)
            nr1 = poldeg1(r)
            nq1 == 0 || (Q[i,j,1:nq1] = copy_oftype(q[1:nq1],T))
            nr1 == 0 || (R[i,j,1:nr1] = copy_oftype(r[1:nr1],T))
        end
    end
    return Q[:,:,1:pmdeg(Q)+1], R[:,:,1:pmdeg(R)+1]
end
pmdivrem(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
         D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}) =
         pmdivrem(poly2pm(N),poly2pm(D))
"""
     pm2lpCF1(P; grade = l) -> (M, N)

Build a strong linearization `M - λN` of a polynomial matrix `P(λ)` in the first companion Frobenius form.

`P(λ)` is a grade `k` polynomial matrix assumed of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`, with
the coefficient matrices `P_i`, `i = 1, ..., k+1` stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).
The effective grade `l` to be used for linearization can be specified via the keyword argument `grade` as
`grade = l`, where `l` must be chosen equal to or greater than the degree of `P(λ)`.
The default value used for `l` is `l = deg(P(λ))`.

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

If `P(λ)` is a `m x n` polynomial matrix of effective grade `l` and degree `d`,
then the resulting matrix pencil `M - λN` satisfies the following conditions [1]:

(1) `M - λN` has dimension `(m+n*(l-1)) x n*l` and `M - λN` is regular if `P(λ)` is regular;

(2) `M - λN` and `P(λ)` have the same finite eigenvalues;

(3) the partial multiplicities of infinite eigenvalues of `M - λN` are in excess with `l-d` to the
partial multiplicities of the infinite eigenvalues of `P(λ)`;

(4) `M - λN` and `P(λ)` have the same number of right Kronecker indices and the right
Kronecker indices of `M - λN` are in excess with `l-1` to the right Kronecker indices of `P(λ)`;

(5) `M - λN` and `P(λ)` have the same left Kronecker structure (i.e., the same left Kronecker indices).

[1] F. De Terán, F. M. Dopico, D. S. Mackey, Spectral equivalence of polynomial matrices and
the Index Sum Theorem, Linear Algebra and Its Applications, vol. 459, pp. 264-333, 2014.
"""
function pm2lpCF1(P::AbstractArray{T,3}; grade::Int = pmdeg(P)) where T
   m, n, k1 = size(P)
   deg = pmdeg(P)
   grade < deg && error("The selected grade must be at least $deg")
   grade == -1 && (return zeros(T,m,n), nothing )
   grade == 0 && (return P[:,:,1], nothing )
   grade == 1 && (grade > deg ? (return P[:,:,1], zeros(T,m,n)) : (return P[:,:,1], -P[:,:,2]) )
   nd = n*grade
   nd1 = nd-n
   grade > deg ? (N = [ zeros(T,nd1+m,n)  [zeros(T,m,nd1); I] ]) :
                 (N = [ [P[:,:,grade+1]; zeros(T,nd1,n)]  [zeros(T,m,nd1); I] ])
   M = [ zeros(T,m,nd);  [I zeros(T,nd1,n)] ]
   deg == -1 && (return M, N)
   k = nd
   it = 1:m
   grade > deg ? ne = deg+1 : ne = max(1,deg)
   for i = 1:ne
      M[it,k-n+1:k] = -P[:,:,i]
      k -= n
   end
   return M, N
end
pm2lpCF1(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,AbstractVecOrMat{<:Number}, Number}; kwargs...) =
       pm2lpCF1(poly2pm(P); kwargs...)
"""
    pm2lpCF2(P; grade = l) -> (M, N)

Build a strong linearization `M - λN` of a polynomial matrix `P(λ)` in the second companion Frobenius form.

`P(λ)` is a grade `k` polynomial matrix assumed of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`, with
the coefficient matrices `P_i`, `i = 1, ..., k+1` stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).
The effective grade `l` to be used for linearization can be specified via the keyword argument `grade` as
`grade = l`, where `l` must be chosen equal to or greater than the degree of `P(λ)`.
The default value used for `l` is `l = deg(P(λ))`.

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

If `P(λ)` is a `m x n` polynomial matrix of effective grade `l` and degree `d`, then the resulting matrix pencil
`M - λN` satisfies the following conditions [1]:

(1) `M - λN` has dimension `l*m x (n+(l-1)*m)` and `M - λN` is regular if `P(λ)` is regular;

(2) `M - λN` and `P(λ)` have the same finite eigenvalues;

(3) the partial multiplicities of infinite eigenvalues of `M - λN` are in excess with `l-d` to the
partial multiplicities of the infinite eigenvalues of `P(λ)`;

(4) `M - λN` and `P(λ)` have the same right Kronecker structure (i.e., the same right Kronecker indices);

(5) `M - λN` and `P(λ)` have the same number of left Kronecker indices and the left
Kronecker indices of `M - λN` are in excess with `l-1` to the left Kronecker indices of `P(λ)`.

[1] F. De Terán, F. M. Dopico, D. S. Mackey, Spectral equivalence of polynomial matrices and
the Index Sum Theorem, Linear Algebra and Its Applications, vol. 459, pp. 264-333, 2014.
"""
function pm2lpCF2(P::AbstractArray{T,3}; grade::Int = pmdeg(P)) where T
   m, n, k1 = size(P)
   deg = pmdeg(P)
   grade < deg && error("The selected grade must be at least $deg")
   grade == -1 && (return zeros(T,m,n), nothing )
   grade == 0 && (return P[:,:,1], nothing  )
   grade == 1 && (grade > deg ? (return P[:,:,1], zeros(T,m,n)) : (return P[:,:,1], -P[:,:,2]) )
   md = m*grade
   md1 = md-m
   grade > deg ? (N = [ zeros(T,md,n)  [zeros(T,m,md1); I] ]) :
                 (N = [ [P[:,:,grade+1]; zeros(T,md1,n)]  [zeros(T,m,md1); I] ])
   M = [ zeros(T,md,n)  [I; zeros(T,m,md1)] ]
   deg == -1 && (return M, N)
   k = md
   it = 1:n
   grade > deg ? ne = deg+1 : ne = max(1,deg)
   for i = 1:ne
       M[k-m+1:k,it] = -P[:,:,i]
       k -= m
   end
   return M, N
end
pm2lpCF2(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,AbstractVecOrMat{<:Number}, Number}; kwargs...) =
       pm2lpCF2(poly2pm(P); kwargs...)
"""
     pm2ls(P; contr = false, obs = false, noseig = false, minimal = false,
              fast = true, atol = 0, rtol) -> (A, E, B, C, D)

Build a structured linearization as a system matrix `S(λ)` of the form

            | A-λE | B |
     S(λ) = |------|---|
            |  C   | D |

of the polynomial matrix `P(λ)` which preserves a part of the Kronecker structure of `P(λ)`.

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

If `d` is the degree of `P(λ)` and `n` is the order of `A-λE`, then the computed linearization satisfies:

(1) `A-λE` is regular and `P(λ) = C*inv(λE-A)*B+D`;

(2) `rank[B A-λE] = n` (controllability) if `minimal = true` or `contr = true`, in which case
the right Kronecker structure is preserved;

(3) `rank[A-λE; C] = n` (observability)  if `minimal = true` or `contr = true`, in which case
the left Kronecker structure is preserved;

(4) `A-λE` has no non-dynamic modes if `minimal = true` or `noseig = true`.

If conditions (1)-(4) are satisfied, the linearization is called `minimal` and the resulting order `n`
is the least achievable order. If conditions (1)-(3) are satisfied, the linearization is called `irreducible`
and the resulting order `n` is the least achievable order using orthogonal similarity transformations.
For an irreducible linearization `S(λ)` preserves the pole-zero structure (finite and infinite) and the
left and right Kronecker structures of `P(λ)`.

The underlying pencil manipulation algorithms [1] and [2] to compute reduced order linearizations
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition, if `fast = false`.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The keyword arguments `atol` and `rtol`, specify the absolute and relative tolerances for the
nonzero coefficients of `P(λ)`, respectively.

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

[2] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function pm2ls(P::AbstractArray{T,3}; minimal::Bool = false, contr::Bool = false, obs::Bool = false, noseig::Bool = false,
               fast::Bool = true, atol::Real = zero(real(T)),
               rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
   minimal && (contr = true; obs = true; noseig = true)
   p, m, k1 = size(P)
   nd = pmdeg(P)+1
   nd == 0 && (return zeros(T,0,0), zeros(T,0,0), zeros(T,0,m), zeros(T,p,0), zeros(T,p,m))
   D = P[:,:,1]
   nd == 1 && (return zeros(T,0,0), zeros(T,0,0), zeros(T,0,m), zeros(T,p,0), D)
   if xor(contr,obs)
      if obs
         # build an observable linearization
         n = p*nd
         E = [zeros(T,n,p) [I; zeros(T,p,p*(nd-1))]]
         B = zeros(T,n,m)
         k = p
         for i = 2:nd
             B[k+1:k+p,:] = P[:,:,i]
             k += p
         end
         C = [ -I  zeros(T,p,p*(nd-1))  ]
      else
         # build a controllable linearization
         n = m*nd
         E = [zeros(T,n,m) [I; zeros(T,m,m*(nd-1))]]
         B = [zeros(T,m*(nd-1),m); -I ]
         C = zeros(T,p,n)
         k = 0
         for i = 1:nd-1
             C[:,k+1:k+m] = P[:,:,nd-i+1]
             k += m
         end
      end
      A = Matrix{T}(I,n,n)
   else
      if p <= m
         # build an observable linearization
         n = p*nd
         E = [zeros(T,n,p) [I; zeros(T,p,p*(nd-1))]]
         B = zeros(T,n,m)
         k = p
         for i = 2:nd
             B[k+1:k+p,:] = P[:,:,i]
             k += p
         end
         C = [ -I  zeros(T,p,p*(nd-1))  ]
         if contr
            # remove uncontrollable part
            T <: BlasFloat ? T1 = T : T1 = promote_type(Float64,T)
            Er = copy_oftype(E,T1)
            Br = copy_oftype(B,T1)
            Cr = copy_oftype(C,T1)
            _, _, nr, nuc = sklf_right!(Er, Br, Cr; fast = fast, atol1 = atol, atol2 = atol, rtol = rtol, withQ = false)
            if nuc > 0
               ir = 1:nr
               # save intermediary results
               E = Er[ir,ir]
               B = Br[ir,:]
               C = Cr[:,ir]
            end
            A = Matrix{T1}(I,nr,nr)
         else
            A = Matrix{T}(I,n,n)
         end
      else
         # build a controllable linearization
         n = m*nd
         E = [zeros(T,n,m) [I; zeros(T,m,m*(nd-1))]]
         B = [zeros(T,m*(nd-1),m); -I ]
         C = zeros(T,p,n)
         k = 0
         for i = 1:nd-1
             C[:,k+1:k+m] = P[:,:,nd-i+1]
             k += m
         end
         if obs
            # remove unobservable part
            T <: BlasFloat ? T1 = T : T1 = promote_type(Float64,T)
            Er = copy_oftype(E,T1)
            Br = copy_oftype(B,T1)
            Cr = copy_oftype(C,T1)
            _, _, nr, nuo = sklf_left!(Er, Cr, Br; fast = fast, atol1 = atol, atol2 = atol, rtol = rtol, withQ = false)
            if nuo > 0
               ir = n-nr+1:n
               # save intermediary results
               E = Er[ir,ir]
               B = Br[ir,:]
               C = Cr[:,ir]
            end
            A = Matrix{T1}(I,nr,nr)
         else
            A = Matrix{T}(I,n,n)
         end
      end
   end
   if noseig
      A, E, B, C, D  = lsminreal(A,E,B,C,D,contr = false, obs = false, fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
   end
   return A, E, B, C, D
end
pm2ls(P::Union{AbstractVecOrMat{Polynomial{T}},Polynomial{T},AbstractVecOrMat{Polynomial{T,X}},Polynomial{T,X}}; kwargs...) where {T,X} =
         pm2ls(poly2pm(P); kwargs...)
pm2ls(P::Union{AbstractVecOrMat{T},Number}; kwargs...) where {T <: Number} =
         pm2ls(poly2pm(P); kwargs...)
# pm2ls(P::Union{AbstractVecOrMat{Polynomial{T,X}},Polynomial{T,X}}; kwargs...) where {T,X} =
#          pm2ls(poly2pm(P); kwargs...)
"""
     pm2lps(P; contr = false, obs = false) -> (A, E, B, F, C, G, D, H)

Build a structured linearization as a system matrix `S(λ)` of the form

            | A-λE | B-λF |
     S(λ) = |------|------|
            | C-λG | D-λH |

of the polynomial matrix `P(λ)` which preserves a part of the Kronecker structure of `P(λ)`.

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

If `d` is the degree of the `p x m` polynomial matrix `P(λ)`, then the computed linearization satisfies:

(1) `A-λE` is a `n x n` regular pencil, where `n = p(d-1)` if `contr = false` and `p <= m`
and `n = m(d-1)` otherwise;

(2) `P(λ) = (C-λG)*inv(λE-A)*(B-λF)+D-λH`;

(3) `rank[B-λF A-λE] = n` for any finite and infinite `λ` (strong controllability) if `contr = true`, in which case
the right Kronecker structure is preserved;

(4) `rank[A-λE; C-λG] = n` for any finite and infinite `λ` (strong observability)  if `obs = true`, in which case
the left Kronecker structure is preserved.

If conditions (1)-(4) are satisfied, the linearization is called `strongly minimal`, the resulting order `n`
is the least achievable order and `S(λ)` preserves the pole-zero structure (finite and infinite) and the
left and right Kronecker structures of `P(λ)`.

The pencil based linearization is built using the methods described in [1].

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function pm2lps(P::AbstractArray{T,3}; contr::Bool = false, obs::Bool = false) where T
   p, m, k1 = size(P)
   d = pmdeg(P)
   nd = d+1
   H = zeros(T,p,m)
   d == -1 && (return zeros(T,0,0), zeros(T,0,0), zeros(T,0,m), zeros(T,0,m), zeros(T,p,0), zeros(T,p,0), H, H)
   D = P[:,:,1]
   d == 0 && (return zeros(T,0,0), zeros(T,0,0), zeros(T,0,m), zeros(T,0,m), zeros(T,p,0), zeros(T,p,0), D, H)
   d == 1 && (return zeros(T,0,0), zeros(T,0,0), zeros(T,0,m), zeros(T,0,m), zeros(T,p,0), zeros(T,p,0), D, -P[:,:,2])
   if obs || (!contr && p <= m)
      # build a strongly observable linearization
      n = p*(d-1)
      E = [zeros(T,n,p) [I; zeros(T,p,n-p)]]
      B = zeros(T,n,m)
      F = zeros(T,n,m)
      k = 0
      for i = 2:d
          B[k+1:k+p,:] = P[:,:,i]
          k += p
      end
      F[n-p+1:n,:] = -P[:,:,nd]
      C = zeros(T,p,n)
      G = [I zeros(T,p,n-p) ]
      A = Matrix{T}(I,n,n)
   else
     # build a strongly controllable linearization
     n = m*(d-1)
     E = [zeros(T,n,m) [I; zeros(T,m,n-m)]]
     B = zeros(T,n,m)
     F = [zeros(T,n-m,m); I]
     C = zeros(T,p,n)
     G = [ -P[:,:,nd] zeros(T,p,n-m)]
     k = n
     for i = 2:d
         #C[:,k-m+1:k] = P[:,:,nd-i]
         C[:,k-m+1:k] = P[:,:,i]
         k -= m
     end
     A = Matrix{T}(I,n,n)
   end
   return A, E, B, F, C, G, D, H
end
pm2lps(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number}; kwargs...) where {T} =
         pm2lps(poly2pm(P); kwargs...)
"""
    ls2pm(A, E, B, C, D; fast = true, atol1 = 0, atol2 = 0, gaintol = 0, rtol = min(atol1,atol2) > 0 ? 0 : n*ϵ, val) -> P

Build the polynomial matrix `P(λ) = C*inv(λE-A)*B+D` corresponding to its structured linearization

     | A-λE | B |
     |------|---|
     |  C   | D |

by explicitly determining for each polynomial entry, its coefficients from its roots and a corresponding gain.

The keyword arguments `atol1` and `atol2` specify the absolute tolerances for the elements of `A`, `B`, `C`, `D`, and,
respectively, of `E`, and `rtol` specifies the relative tolerances for the nonzero elements of `A`, `B`, `C`, `D` and `E`.
The default relative tolerance is `(n+1)*ϵ`, where `n` is the size of the size dimension of `A`, and `ϵ` is the
machine epsilon of the element type of coefficients of `A`.

The keyword argument `gaintol` specifies the threshold for the magnitude of the nonzero elements of the gain matrix
`C*inv(γE-A)*B+D`, where `γ = val` if `val` is a number or `γ` is a randomly chosen complex value of unit magnitude,
if `val = missing`. Generally, `val` should not be a root of any of entries of `P`.
"""
function ls2pm(A::AbstractMatrix, E::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix;
                fast::Bool = true, atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))),
                gaintol::Real = zero(real(eltype(A))), val::Union{Number,Missing} = missing,
                rtol::Real = ((size(A,1)+1)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)))
   n = LinearAlgebra.checksquare(A)
   (n,n) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
   p, m = size(D)
   (n,m) == size(B) || throw(DimensionMismatch("A, B and D must have compatible dimensions"))
   (p,n) == size(C) || throw(DimensionMismatch("A, C and D must have compatible dimensions"))

   T = promote_type(eltype(A),eltype(E), eltype(B),eltype(C),eltype(D))
   T <: BlasFloat || (T = promote_type(Float64,T))
   compl = T <: Complex

   P = zeros(T,p,m,n+1)
   ismissing(val) && (val = exp(rand()*im))
   A1, E1, B1, C1, D1 = lsminreal2(A, E, B, C, D, infinite = false, noseig = false, atol1 = atol1, atol2 = atol2)
   Pval = lseval(A1, E1, B1, C1, D1, val)
   isunimodular(A1, E1, atol1 = atol1, atol2 = atol2, rtol = rtol) ||
                error("The given linearization cannot be converted to a polynomial form")
   for i = 1:p
       for j = 1:m
           if abs(Pval[i,j]) > gaintol
              zer, iz, = spzeros(A1, E1, B1[:,j:j], C1[i:i,:], D1[i:i,j:j];
                                 fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
              c, pval = polcoeffval(zer[1:(length(zer)-sum(iz))],val)
              P[i,j,1:length(c)] = compl ? c*(Pval[i,j]/pval) : real(c*(Pval[i,j]/pval))
           end
       end
   end
   return P[:,:,1:pmdeg(P)+1]
end
"""
    lps2pm(A, E, B, F, C, G, D, H; fast = true, atol1 = 0, atol2 = 0, gaintol = 0, rtol = min(atol1,atol2) > 0 ? 0 : n*ϵ, val) -> P

Build the polynomial matrix `P(λ) = (C-λG)*inv(λE-A)*(B-λF)+D-λH` corresponding to its structured linearization

     | A-λE | B-λF |
     |------|------|
     | C-λG | D-λH |

by explicitly determining for each polynomial entry, its coefficients from its roots and corresponding gain.

The keyword arguments `atol1` and `atol2` specify the absolute tolerances for the elements of `A`, `B`, `C`, `D`, and of
`E`, `F`, `G`, `H`, respectively,  and `rtol` specifies the relative tolerances for the nonzero elements of
`A`, `B`, `C`, `D`, `E`, F`, `G`, `H`.
The default relative tolerance is `(n+2)*ϵ`, where `n` is the size of the size dimension of `A`, and `ϵ` is the
machine epsilon of the element type of coefficients of `A`.

The keyword argument `gaintol` specifies the threshold for the magnitude of the nonzero elements of the gain matrix
`C*inv(γE-A)*B+D`, where `γ = val` if `val` is a number or `γ` is a randomly chosen complex value of unit magnitude,
if `val = missing`. Generally, `val` should not be a root of any of entries of `P`.
"""
function lps2pm(A::AbstractMatrix, E::AbstractMatrix,B::AbstractMatrix, F::AbstractMatrix,
   C::AbstractMatrix, G::AbstractMatrix, D::AbstractMatrix, H::AbstractMatrix;
   fast::Bool = true, atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))),
   gaintol::Real = zero(real(eltype(A))), val::Union{Number,Missing} = missing,
   rtol::Real = ((size(A,1)+2)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)))
   n = LinearAlgebra.checksquare(A)
   (n,n) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
   p, m = size(D)
   (n,m) == size(B) || throw(DimensionMismatch("A, B and D must have compatible dimensions"))
   (p,n) == size(C) || throw(DimensionMismatch("A, C and D must have compatible dimensions"))
   (n,m) == size(F) || throw(DimensionMismatch("B and F must have the same dimensions"))
   (p,n) == size(G) || throw(DimensionMismatch("C and G must have the same dimensions"))
   (p,m) == size(H) || throw(DimensionMismatch("D and H must have the same dimensions"))
   T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D), eltype(E), eltype(F), eltype(G), eltype(H))
   T <: BlasFloat || (T = promote_type(Float64,T))

   compl = T <: Complex
   P = zeros(T,p,m,n+2)
   ismissing(val) && (val = exp(rand()*im))
   A1, E1, B1, F1, C1, G1, D1, H1, V1, W1 = lpsminreal(A, E, B, F, C, G, D, H, atol1 = atol1, atol2 = atol2, rtol = rtol)
   Pval = V1'\lpseval(A1, E1, B1, F1, C1, G1, D1, H1, val)/W1
   isunimodular(A1, E1, atol1 = atol1, atol2 = atol2, rtol = rtol) ||
                error("The given linearization cannot be converted to a polynomial form")
   M1 = [A1 B1/W1; V1'\C1 V1'\D1/W1]
   N1 = [E1 F1/W1; V1'\G1 V1'\H1/W1]
   n = size(A1,1)
   indi = [1:n;[1]]
   indj = [1:n;[1]]
   n1 = n+1
   for i = 1:p
      indi[n1] = n+i
      for j = 1:m
         if abs(Pval[i,j]) > gaintol
            indj[n1] = n+j
            zer, iz, = pzeros(view(M1,indi,indj),view(N1,indi,indj);
                              fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
            c, pval = polcoeffval(zer[1:(length(zer)-sum(iz))],val)
            P[i,j,1:length(c)] = compl ? c*(Pval[i,j]/pval) : real(c*(Pval[i,j]/pval))
         end
      end
   end
   return P[:,:,1:pmdeg(P)+1]
end
"""
     spm2ls(T, U, V, W; fast = true, contr = false, obs = false, minimal = false, atol = 0, rtol) -> (A, E, B, C, D)

Build a structured linearization as a system matrix `S(λ)` of the form

            | A-λE | B |
     S(λ) = |------|---|
            |  C   | D |

of the structured polynomial matrix

              | -T(λ) | U(λ) |
       P(λ) = |-------|------|
              | V(λ)  | W(λ) |

such that `V(λ)*inv(T(λ))*U(λ)+W(λ) = C*inv(λE-A)*B+D`. The resulting linearization `S(λ)` preserves a part,
if `minimal = false`, or the complete Kronecker structure, if `minimal = true`, of `P(λ)`. In the latter case,
the order `n` of `A-λE` is the least possible one and `S(λ)` is a strong linearization of `P(λ)`.

`T(λ)`, `U(λ)`, `V(λ)`, and `W(λ)` can be specified as polynomial matrices of the form `X(λ) = X_1 + λ X_2 + ... + λ**k X_(k+1)`,
for `X = T`, `U`, `V`, and `W`, for which the coefficient matrices `X_i`, `i = 1, ..., k+1`, are stored in
the 3-dimensional matrices `X`, where `X[:,:,i]` contains the `i`-th coefficient matrix `X_i` (multiplying `λ**(i-1)`).

`T(λ)`, `U(λ)`, `V(λ)`, and `W(λ)` can also be specified as matrices, vectors or scalars of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The computed structured linearization satisfies:

(1) `A-λE` is regular;

(2) `rank[B A-λE] = n` (controllability) if `minimal = true` or `contr = true`, in which case
the finite and right Kronecker structures are preserved;

(3) `rank[A-λE; C] = n` (observability)  if `minimal = true` or `obs = true`, in which case
the finite and left Kronecker structures are preserved;

(4) `A-λE` has no simple infinite eigenvalues if `minimal = true`, in which case the complete Kronecker structure is preserved.

The keyword arguments `atol` and `rtol`, specify, respectively, the absolute and relative tolerance for the
nonzero coefficients of the matrices `T(λ)`, `U(λ)`, `V(λ)` and `W(λ)`. The default relative tolerance is `nt*ϵ`,
where `nt` is the size of the square matrix `T(λ)` and `ϵ` is the machine epsilon of the element type of its coefficients.

The structured linearization is built using the methods described in [1].

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function spm2ls(T::Union{AbstractArray{T1,3},AbstractArray{T1,2}},U::Union{AbstractArray{T2,3},AbstractArray{T2,2}},
                V::Union{AbstractArray{T3,3},AbstractArray{T3,2}},W::Union{AbstractArray{T4,3},AbstractArray{T4,2}};
                contr::Bool = false, obs::Bool = false, minimal::Bool = false,
                fast::Bool = true, atol::Real = zero(real(T1)),
                rtol::Real = size(T,1)*eps(real(float(one(T1))))*iszero(atol)) where {T1, T2, T3, T4}
   if ndims(T) == 2
      n, nt = size(T)
      n == nt || throw(DimensionMismatch("T(λ) must be a square polynomial matrix"))
      n == rank(T, atol=atol, rtol=rtol) || error("T(λ) must be a regular square polynomial matrix")
      ndT = 1
      T = reshape(T,n,nt,ndT)
   else
      n, nt, ndT = size(T)
      n == nt || throw(DimensionMismatch("T(λ) must be a square polynomial matrix"))
      n == pmrank(T, atol=atol, rtol=rtol) || error("T(λ) must be a regular square polynomial matrix")
      ndT = max(pmdeg(T)+1,1)
   end
   if ndims(U) == 2
      nt, m = size(U)
      ndU = 1
      U = reshape(U,nt,m,ndU)
   else
      nt, m, ndU = size(U)
      ndU = max(pmdeg(U)+1,1)
   end
   n == nt || throw(DimensionMismatch("T(λ) and U(λ) must have the same number of rows"))
   if ndims(V) == 2
      p, nt = size(V)
      ndV = 1
      V = reshape(V,p,nt,ndV)
   else
      p, nt, ndV = size(V)
      ndV = max(pmdeg(V)+1,1)
   end
   n == nt || throw(DimensionMismatch("T(λ) and V(λ) must have the same number of columns"))
   if ndims(W) == 2
      pt, mt = size(W)
      ndW = 1
      W = reshape(W,pt, mt,ndW)
   else
      pt, mt, ndW = size(W)
      ndW = max(pmdeg(W)+1,1)
   end
   pt == p || throw(DimensionMismatch("W(λ) and V(λ) must have the same number of rows"))
   mt == m || throw(DimensionMismatch("W(λ) and U(λ) must have the same number of columns"))
   nd = max(ndT,ndU,ndV,ndW)
   TT = promote_type(T1, T2, T3, T4)
   if nd == 1
      if minimal
         Ar,Er,Br,Cr,Dr = lsminreal(T[:,:,1], zeros(TT,n,n), U[:,:,1], V[:,:,1], W[:,:,1],
                                    contr = false, obs = false, atol1 = atol, atol2 = atol, rtol = rtol)
         return Ar,Er,Br,Cr,Dr
      else
        return T[:,:,1], zeros(TT,n,n), U[:,:,1], V[:,:,1], W[:,:,1]
      end
   end

   # build the compound polynomial matrix [-T(λ) U(λ); V(λ) W(λ)]
   P = zeros(TT,n+p,n+m,nd)
   ia = 1:n
   jb = n+1:n+m
   ic = n+1:n+p
   P[ia,ia,1:ndT] = -T[:,:,1:ndT]
   P[ia,jb,1:ndU] = U[:,:,1:ndU]
   P[ic,ia,1:ndV] = V[:,:,1:ndV]
   P[ic,jb,1:ndW] = W[:,:,1:ndW]

   # build a linearization of the compound polynomial matrix [-T(λ) U(λ); V(λ) W(λ)]
   A, E, B, C, D = pm2ls(P, contr = contr, obs = obs, noseig = false,
                           fast = fast, atol = atol, rtol = rtol)

   # form the linearization of P(λ) = V(λ)*inv(T(λ))*U(λ)+W(λ)
   nr = size(A,1)
   Ar = [A B[:,ia]; C[ia,:] D[ia,ia]]
   Er = [E zeros(TT,nr,n); zeros(TT,n,n+nr)]
   Br = [B[:,jb]; D[ia,jb]]
   Cr = [C[ic,:] D[ic,ia]]
   Dr = D[ic,jb]
   if minimal
      Ar,Er,Br,Cr,Dr = lsminreal(Ar,Er,Br,Cr,Dr, fast=fast, atol1 = atol, atol2 = atol, rtol = rtol)
   end
   return Ar, Er, Br, Cr, Dr
end
spm2ls(T::Union{AbstractVecOrMat{Polynomial{T1}},Polynomial{T1},Number,AbstractMatrix{T1}}, U::Union{AbstractVecOrMat{Polynomial{T2}},Polynomial{T2},Number,AbstractVecOrMat{T2}},
       V::Union{AbstractVecOrMat{Polynomial{T3}},Polynomial{T3},Number,AbstractVecOrMat{T3}}, W::Union{AbstractVecOrMat{Polynomial{T4}},Polynomial{T4},Number,AbstractVecOrMat{T4}}; kwargs...) where {T1, T2, T3, T4} =
       spm2ls(poly2pm(T),poly2pm(U),poly2pm(V),poly2pm(W); kwargs...)
"""
     spm2lps(T, U, V, W; fast = true, contr = false, obs = false, minimal = false, atol = 0, rtol) -> (A, E, B, F, C, G, D, H)

Build a structured linearization

              | A-λE | B-λF |
     M - λN = |------|------|
              | C-λG | D-λH |

of the structured polynomial matrix

              | -T(λ) | U(λ) |
       P(λ) = |-------|------|
              |  V(λ) | W(λ) |

such that `V(λ)*inv(T(λ))*U(λ)+W(λ) = (C-λG))*inv(λE-A)*(B-λF)+D-λH`. The resulting linearization `M - λN` preserves a part,
if `minimal = false`, or the complete Kronecker structure, if `minimal = true`, of `P(λ)`. In the latter case,
the order `n` of `A-λE` is the least possible one and `M - λN` is a strong linearization of `P(λ)`.

`T(λ)`, `U(λ)`, `V(λ)`, and `W(λ)` can be specified as polynomial matrices of the form `X(λ) = X_1 + λ X_2 + ... + λ**k X_(k+1)`,
for `X = T`, `U`, `V`, and `W`, for which the coefficient matrices `X_i`, `i = 1, ..., k+1`, are stored in
the 3-dimensional matrices `X`, where `X[:,:,i]` contains the `i`-th coefficient matrix `X_i` (multiplying `λ**(i-1)`).

`T(λ)`, `U(λ)`, `V(λ)`, and `W(λ)` can also be specified as matrices, vectors or scalars of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The computed structured linearization satisfies:

(1) `A-λE` is regular;

(2) `rank[B-λF A-λE] = n` (strong controllability) if `minimal = true` or `contr = true`, in which case
the finite and right Kronecker structures are preserved;

(3) `rank[A-λE; C-λG] = n` (strong observability)  if `minimal = true` or `obs = true`, in which case
the finite and left Kronecker structures are preserved.

The keyword arguments `atol` and `rtol`, specify, respectively, the absolute and relative tolerance for the
nonzero coefficients of the matrices `T(λ)`, `U(λ)`, `V(λ)` and `W(λ)`. The default relative tolerance is `nt*ϵ`,
where `nt` is the size of the square matrix `T(λ)` and `ϵ` is the machine epsilon of the element type of its coefficients.
"""
function spm2lps(T::Union{AbstractArray{T1,3},AbstractArray{T1,2}},U::Union{AbstractArray{T2,3},AbstractArray{T2,2}},
                V::Union{AbstractArray{T3,3},AbstractArray{T3,2}},W::Union{AbstractArray{T4,3},AbstractArray{T4,2}};
                contr::Bool = false, obs::Bool = false, minimal::Bool = false,
                fast::Bool = true, atol::Real = zero(real(T1)),
                rtol::Real = size(T,1)*eps(real(float(one(T1))))*iszero(atol)) where {T1, T2, T3, T4}
   if ndims(T) == 2
      n, nt = size(T)
      n == nt || throw(DimensionMismatch("T(λ) must be a square polynomial matrix"))
      n == rank(T, atol=atol, rtol=rtol) || error("T(λ) must be a regular square polynomial matrix")
      ndT = 1
      T = reshape(T,n,nt,ndT)
   else
      n, nt, ndT = size(T)
      n == nt || throw(DimensionMismatch("T(λ) must be a square polynomial matrix"))
      n == pmrank(T, atol=atol, rtol=rtol) || error("T(λ) must be a regular square polynomial matrix")
      ndT = max(pmdeg(T)+1,1)
   end
   if ndims(U) == 2
      nt, m = size(U)
      ndU = 1
      U = reshape(U,nt,m,ndU)
   else
      nt, m, ndU = size(U)
      ndU = max(pmdeg(U)+1,1)
   end
   n == nt || throw(DimensionMismatch("T(λ) and U(λ) must have the same number of rows"))
   if ndims(V) == 2
      p, nt = size(V)
      ndV = 1
      V = reshape(V,p,nt,ndV)
   else
      p, nt, ndV = size(V)
      ndV = max(pmdeg(V)+1,1)
   end
   n == nt || throw(DimensionMismatch("T(λ) and V(λ) must have the same number of columns"))
   if ndims(W) == 2
      pt, mt = size(W)
      ndW = 1
      W = reshape(W,pt, mt,ndW)
   else
      pt, mt, ndW = size(W)
      ndW = max(pmdeg(W)+1,1)
   end
   pt == p || throw(DimensionMismatch("W(λ) and V(λ) must have the same number of rows"))
   mt == m || throw(DimensionMismatch("W(λ) and U(λ) must have the same number of columns"))
   nd = max(ndT,ndU,ndV,ndW)
   TT = promote_type(T1, T2, T3, T4)
   if nd == 1
      if minimal
         Ar,Er,Br,Cr,Dr = lsminreal(T[:,:,1], zeros(TT,n,n), U[:,:,1], V[:,:,1], W[:,:,1],
                                    contr = false, obs = false, atol1 = atol, atol2 = atol, rtol = rtol)
         nr = size(Ar,1)
         TW = eltype(Ar)
         return Ar,Er,Br,zeros(TW,nr,m),Cr,zeros(TW,p,nr),Dr,zeros(TW,p,m)
      else
        return T[:,:,1], zeros(TT,n,n), U[:,:,1], zeros(TT,n,m), V[:,:,1], zeros(TT,p,n), W[:,:,1], zeros(TT,p,m)
      end
   end

   # build the compound polynomial matrix [-T(λ) U(λ); V(λ) W(λ)]
   P = zeros(TT,n+p,n+m,nd)
   ia = 1:n
   jb = n+1:n+m
   ic = n+1:n+p
   P[ia,ia,1:ndT] = -T[:,:,1:ndT]
   P[ia,jb,1:ndU] = U[:,:,1:ndU]
   P[ic,ia,1:ndV] = V[:,:,1:ndV]
   P[ic,jb,1:ndW] = W[:,:,1:ndW]

   # build a linearization of the compound polynomial matrix [-T(λ) U(λ); V(λ) W(λ)]
   A, E, B, F, C, G, D, H = pm2lps(P, contr = contr, obs = obs)

   # form the linearization of P(λ) = V(λ)*inv(T(λ))*U(λ)+W(λ)
   nr = size(A,1)
   Ar = [A B[:,ia]; C[ia,:] D[ia,ia]]
   Er = [E F[:,ia]; G[ia,:] H[ia,ia]]
   Br = [B[:,jb]; D[ia,jb]]
   Fr = [F[:,jb]; H[ia,jb]]
   Cr = [C[ic,:] D[ic,ia]]
   Gr = [G[ic,:] H[ic,ia]]
   Dr = D[ic,jb]
   Hr = H[ic,jb]
   if minimal
      Ar,Er,Br,Fr,Cr,Gr,Dr,Hr,Vr,Wr = lpsminreal(Ar,Er,Br,Fr,Cr,Gr,Dr,Hr, fast=fast, atol1 = atol, atol2 = atol, rtol = rtol)
      return Ar,Er,Br/Wr,Fr/Wr,Vr'\Cr,Vr'\Gr,Vr'\Dr/Wr,Vr'\Hr/Wr
   else
      return Ar,Er,Br,Fr,Cr,Gr,Dr,Hr
   end
end
function spm2lps(T::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
        U::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
        V::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
        W::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
   # TODO: checking that all entries have the same variable
    return spm2lps(poly2pm(T),poly2pm(U),poly2pm(V),poly2pm(W); kwargs...)
end


text:
# modified example from GradDescent.jl documentation
using GradDescent
using Yota
using Statistics
using Random

obj(Y, X, b) = mean((Y .- X * b) .^ 2.0) # objective to minimize

srand(1) # set seed
n = 1000 # number of observations
d = 10   # number of covariates
X = randn(n, d) # simulated covariates
b = randn(d)    # generated coefficients
ϵ = randn(n) * 0.1 # noise
Y = X * b + ϵ # observed outcome

epochs = 100 # number of epochs

θ = randn(d) # initialize model parameters
opt = Adam(α=1.0)  # initalize optimizer with learning rate 1.0

for i in 1:epochs
    val, g = grad(obj, Y, X, θ)
    println("Epoch: $i; loss = $val")

    δ = update(opt, g[3])
    θ -= δ
end

ext: jl


text:
import Blosc
import CodecZlib
import JSON

abstract type Compressor end
getCompressor(compdict::Dict) = getCompressor(compressortypes[compdict["id"]],compdict)
getCompressor(::Nothing) = NoCompressor()

struct BloscCompressor <: Compressor
    blocksize::Int
    clevel::Int
    cname::String
    shuffle::Bool
end

"""
    BloscCompressor(;blocksize=0, clevel=5, cname="lz4", shuffle=true)

Returns a `BloscCompressor` struct that can serve as a Zarr array compressor. Keyword arguments are:

* `clevel=5` the compression level, number between 0 (no compression) and 9 (max compression)
* `cname="lz4"` compressor name, can be one of `"blosclz"`, `"lz4"`, and `"lz4hc"`
* `shuffle=true` enables/disables bit-shuffling
"""
BloscCompressor(;blocksize=0, clevel=5, cname="lz4", shuffle=true) =
    BloscCompressor(blocksize, clevel, cname, shuffle)

function getCompressor(::Type{BloscCompressor}, d::Dict)
    BloscCompressor(d["blocksize"], d["clevel"], d["cname"], d["shuffle"] > 0)
end

zuncompress(a, r::AbstractArray, ::BloscCompressor) = copyto!(a, Blosc.decompress(Base.nonmissingtype(eltype(a)), r))

function zcompress(a, f::AbstractArray, c::BloscCompressor)
    Blosc.set_compressor(c.cname)
    r = Blosc.compress(a, level=c.clevel, shuffle=c.shuffle)
    empty!(f)
    append!(f, r)
end

JSON.lower(c::BloscCompressor) = Dict("id"=>"blosc", "cname"=>c.cname,
    "clevel"=>c.clevel, "shuffle"=>c.shuffle ? 1 : 0, "blocksize"=>c.blocksize)

"""
    NoCompressor()

Creates an object that can be passed to ZArray constructors without compression.
"""
struct NoCompressor <: Compressor end

function zuncompress(a, r::AbstractArray, ::NoCompressor)
  copyto!(a, reinterpret(eltype(a),r))
end

function zcompress(a, f::AbstractArray, ::NoCompressor)
  a2 = reinterpret(UInt8,a)
  empty!(f)
  append!(f, a2)
end

JSON.lower(::NoCompressor) = nothing

compressortypes = Dict("blosc"=>BloscCompressor, nothing=>NoCompressor)

"""
    ZlibCompressor(clevel=-1)
Returns a `ZlibCompressor` struct that can serve as a Zarr array compressor. Keyword arguments are:
* `clevel=-1` the compression level, number between -1 (Default), 0 (no compression) and 9 (max compression)
*  default is -1 compromise between speed and compression (currently equivalent to level 6).
"""
struct ZlibCompressor <: Compressor
    clevel::Int
end

ZlibCompressor(;clevel=-1) = ZlibCompressor(clevel)

function getCompressor(::Type{ZlibCompressor}, d::Dict)
    ZlibCompressor(d["level"])
end

function zuncompress(a, r::AbstractArray, ::ZlibCompressor)
    result = transcode(CodecZlib.ZlibDecompressor,r)
    copyto!(a, reinterpret(Base.nonmissingtype(eltype(a)),result))
end

function zcompress(a, f::AbstractArray, ::ZlibCompressor)
    a_uint8 = reinterpret(UInt8,a)
    r = transcode(CodecZlib.ZlibCompressor,a_uint8[:])
    empty!(f)
    append!(f,r)
end

JSON.lower(z::ZlibCompressor) = Dict("id"=>"zlib", "level" => z.clevel)

Zarr.compressortypes["zlib"] = ZlibCompressor


text:
push!(LOAD_PATH, "src")

using RegERMs, Base.Test

tests = [
    "loss",
    "model",
    "mercer_map",
    "svm",
    "logistic_regression",
    "ridge_regression"
]

println("Running tests:")

for t in tests
    test_fn = "$t.jl"
    println("* $test_fn")
    include(test_fn)
end


text:
#!/usr/bin/env julia

using DandelionSlack
using DandelionSlack.Util

using DocOpt

doc = """Find the Slack channel id given a channel name.

Usage:
    find_channel_id.jl <team> <name>
    find_channel_id.jl --version

Options:
    --version   Show version
"""

arguments = docopt(doc, version=v"0.0.1")

team = arguments["<team>"]
name = arguments["<name>"]

token = find_token(team)

channels_list = ChannelsList(token, Nullable{Int64}())
try
    status, response = makerequest(channels_list, real_requests)
    channel_index = findfirst(c -> c.name == name, response.channels)
    if channel_index == 0
        println("No channel named $(name) in team $(team)")
    else
        channel = response.channels[channel_index]
        println("Channel id for $(name) in team $(team): $(channel.id)")
    end
catch ex
    dump(ex)
end


text:
@compat abstract type StepSizeRule end
mutable struct constantStep <: StepSizeRule
    α::Float64
end
mutable struct constantLength <: StepSizeRule
    α::Float64
end
mutable struct poliakRule <: StepSizeRule
    η::Float64
end
mutable struct diminishingAffine <: StepSizeRule
    a::Float64
    b::Float64
end
mutable struct diminishingSqrt <: StepSizeRule
    a::Float64
end

function computeStepSize(rule::constantStep, epoch, ∇q, UB, LB)
  return rule.α
end
function computeStepSize(rule::constantLength, epoch, ∇q, UB, LB)
  return rule.α / norm(∇q)
end
function computeStepSize(rule::poliakRule, epoch, ∇q, UB, LB)
  return rule.η*(UB-LB)/ dot(∇q, ∇q)
end
function computeStepSize(rule::diminishingAffine, epoch, ∇q, UB, LB)
  return rule.a / (rule.b + epoch)
end
function computeStepSize(rule::diminishingSqrt, epoch, ∇q, UB, LB)
  return rule.a / sqrt(epoch)
end

function slashStep!(rule::constantStep, factor)
  rule.α /= factor
end
function slashStep!(rule::constantLength, factor)
  rule.α /= factor
end
function slashStep!(rule::poliakRule, factor)
  rule.η /= factor
end
function slashStep!(rule::diminishingAffine, factor)
  rule.a /= factor
end
function slashStep!(rule::diminishingSqrt, factor)
  rule.a /= factor
end


text:
using FredData
using Test

include("test_without_key.jl")

# Normal usage - API key must be present in ENV
if haskey(ENV, "FRED_API_KEY")
    include("test_with_key.jl")
end


text:
#unum-oddsandends.jl
#mathematical odds and ends

################################################################################

doc"""
  `Unums.frac_val(::UInt64)`
  `Unums.frac_val(::ArrayNum)`
  returns the integral value for a number used as a fraction.
"""
frac_val(x::UInt64) = x
function frac_val{FSS}(v::ArrayNum{FSS})
  sum = big(0)
  for i = 1:length(v.a)
    sum <<= 64
    sum += v.a[i]
  end
  sum
end

doc"""`Unums.calculate(x::Unum)` returns a bigfloat equivalent of the unum.  NB:
currently doesn't work so well for FSS > 9"""
@universal function calculate(x::Unum)
  sign = (x.flags & UNUM_SIGN_MASK != 0) ? -1 : 1
  #the sub`normal case
  if (x.exponent == 0)
    big(2.0)^(decode_exp(x) + 1) * sign * (frac_val(x.fraction)) / big(2.0)^(64 * __ulength(x.fraction))
  else #the normalcase
    big(2.0)^(decode_exp(x)) * sign * (1 + frac_val(x.fraction) / big(2.0)^(64 * __ulength(x.fraction)))
  end
end
export calculate

__ulength(::UInt64) = 1
__ulength{FSS}(::ArrayNum{FSS}) = __cell_length(FSS)

################################################################################

doc"""
  `additiveinverse!` creates the additive inverse value of a unum, by flipping
  the sign.  This can be better than the `-` operator because it doesn't copy
  the unum.  A reference to the unum is returned.
"""
@universal additiveinverse!(x::Unum) = (x.flags $= UNUM_SIGN_MASK; return x)
export additiveinverse!

doc"""
  `coerce_sign!(a::Unum, b)` maps the sign bit from b onto a.  b can either be a
  UInt16 or a Unum.
"""
@universal coerce_sign!(a::Unum, b::Unum) = coerce_sign!(a, b.flags)
@universal coerce_sign!(a::Unum, sgn::UInt16) = (a.flags = (a.flags & ~UNUM_SIGN_MASK) | (sgn & UNUM_SIGN_MASK); return a)

doc"""
  `abs!(::Unum)` forces the value of the unum to be positive.  Returns the
  unum for chaining purposes.
"""
@universal abs!(x::Unum) = ((x.flags &= ~UNUM_SIGN_MASK); return x)
@universal Base.abs(x::Unum) = abs!(copy(x))
export abs!

@universal function Base.copy!(dest::Unum, src::Unum)
  dest.fsize = src.fsize
  dest.esize = src.esize
  dest.flags = src.flags & UNUM_FLAG_MASK
  dest.exponent = src.exponent

  (FSS < 7) ? (dest.fraction = src.fraction) : (copy!(dest.fraction, src.fraction))

  return dest  #for chaining purposes
end

doc"""
  `Unums.normalize!(::Unum)` takes a unum that is purportedly subnormal form and
  normalizes it.  This entails shifting just past the top bit.  this function
  returns the number of places shifted.

  This function should not be run on a fraction that is all zero, nor on a
  function which is not subnormal.
"""
@universal function normalize!(x::Unum)
  leftshift = clz(x.fraction) + o16
  frac_lsh!(x, leftshift)
  x.fsize = (leftshift > x.fsize) ? z16 : x.fsize - leftshift
  return leftshift
end

doc"""
  `Unums.resolve_degenerates!(::Unum)` checks for degeneracy in unum values,
  and resolves to "canonical" form - which means all nonzero subnormals are
  converted to normal form if possible, and exact zeros are the smallest zero.
  fsize is maximally trimmed for exact values.
"""
@universal function resolve_degenerates!(x::Unum)
  (x.exponent != 0) && return x   #kick out if our exponent is not zero.

  (x.esize == max_esize(ESS)) && return x #kick out if we're not a strange subnormal.

  if is_all_zero(x.fraction)
    is_exact(x) && return zero(typeof(x))
    return x #if we're actually zero or a zero+ulp subnormal we can't shift.
  end

  true_exponent = decode_exp(x)
  #now, count leading zeros, be prepared to shift left.
  leftshift = normalize!(x)
  true_exponent -= leftshift - o16
  (x.esize, x.exponent) = encode_exp(true_exponent)
  exact_trim!(x)

  return x
end

################################################################################
## variadic macros that trigger exactitude checking.

if options[:devmode]
  macro ensure_exact(x)
    esc(options[:devmode] ? :(is_exact($x) || throw(ArgumentError("passed parameter must be exact"))) : :())
  end
else
  macro ensure_exact(x); :(); end
end

if options[:devmode]
  macro ensure_ulp(x)
    esc(options[:devmode] ? :(is_ulp($x) || throw(ArgumentError("passed parameter must be an ulp"))) : :())
  end
else
  macro ensure_ulp(x); :(); end
end

################################################################################
## sophisticated exactitude functions.

doc"""
  `Unums.outer_ulp!(::Unum)` returns the smallest-width ulp immediately above the
  current (exact) unum.
"""
@universal function outer_ulp!(x::Unum)
  @ensure_exact(x)

  resolve_degenerates!(x)

  x.fsize = max_fsize(FSS)
  make_ulp!(x)
end
@universal outer_ulp(x::Unum) = outer_ulp!(copy(x))

doc"""
  `Unums.inner_ulp!(::Unum)` returns the smallest-width ulp immediately below the
  current unum.
"""
@universal function inner_ulp!(x::Unum)
  @ensure_exact(x)
  resolve_degenerates!(x)
  make_ulp!(x)

  borrowed = frac_sub_ubit!(x, max_fsize(FSS))

  x.fsize = max_fsize(FSS)

  #if we borrowed, then fraction must have been zero.
  if borrowed
    #there is no ulp inward of x.  Consider replacing this with "Throw an error"
    (x.exponent == 0) && return nan!(x)
    _xexp = decode_exp(x)
    (_xexp == min_exponent(ESS)) && (x.exponent = 0; return x)
    (x.esize, x.exponent) = encode_exp(_xexp - 1)
  end

  return x
end
@universal inner_ulp(x::Unum) = inner_ulp!(copy(x))

doc"""
"""
@universal function outer_exact!(x::Unum)
  @ensure_ulp(x)
  resolve_degenerates!(x)
  carry = frac_add_ubit!(x, x.fsize)
  if carry
    if (is_subnormal(x))
      (x.exponent = o64)
    else
      exp = decode_exp(x) + 1
      if exp > max_exponent(ESS)
        inf!(x, @signof(x))
      else
        (x.esize, x.exponent) = encode_exp(exp)
      end
    end
  end
  exact_trim!(make_exact!(x))
end
@universal outer_exact(x::Unum) = outer_exact!(copy(x))

@universal function inner_exact!(x::Unum)
  @ensure_ulp(x)
  make_exact!(x)
end
@universal inner_exact(x::Unum) = inner_exact!(copy(x))

#real number/dedekind cut formulas.
@universal function lub(x::Unum)
  is_exact(x) && return x
  is_positive(x) && return outer_exact(x)
  return inner_exact(x)
end
@universal function glb(x::Unum)
  is_exact(x) && return x
  is_positive(x) && return inner_exact(x)
  return outer_exact(x)
end

################################################################################
## recast as upper and lower versions.

@universal upper_exact!(x::Unum) = is_positive(x) ? outer_exact!(x) : inner_exact!(x)
@universal lower_exact!(x::Unum) = is_positive(x) ? inner_exact!(x) : outer_exact!(x)
@universal upper_ulp!(x::Unum) = is_positive(x) ? outer_ulp!(x) : inner_ulp!(x)
@universal lower_ulp!(x::Unum) = is_positive(x) ? inner_ulp!(x) : outer_ulp!(x)

@universal upper_exact(x::Unum) = upper_exact!(copy(x))
@universal lower_exact(x::Unum) = lower_exact!(copy(x))
@universal upper_ulp(x::Unum) = upper_ulp!(copy(x))
@universal lower_ulp(x::Unum) = lower_ulp!(copy(x))

################################################################################

@universal function next_unum!(x::Unum)
  @ensure_exact(x)
  resolve_degenerates!(x)
  carried = frac_add_ubit!(x, max_fsize(FSS))
  if carried
    exponent = decode_exp(x)
    (exponent > max_exponent(ESS)) && return U(Inf)
    (x.esize, x.exponent) = encode_exp(exponent + 1)
  end
  exact_trim!(x)
  return x
end

frac_ctz{ESS,FSS}(x::UnumSmall{ESS,FSS}) = ctz(x.fraction >> (64 - 1 << FSS))
frac_ctz{ESS,FSS}(x::UnumLarge{ESS,FSS}) = ctz(x.fraction)

frac_cto{ESS,FSS}(x::UnumSmall{ESS,FSS}) = cto(x.fraction >> (64 - 1 << FSS))
frac_cto{ESS,FSS}(x::UnumLarge{ESS,FSS}) = cto(x.fraction)

@universal next_unum(x::Unum) = next_unum!(copy(x))
################################################################################
## dumb exactitude functions.

doc"""`Unums.make_exact(::Unum)` forces the ubit of a unum to be 0."""
@universal make_exact!(x::Unum) = (x.flags &= ~UNUM_UBIT_MASK; x)
@universal make_exact(x::Unum) = make_exact!(copy(x))

doc"""`Unums.make_ulp(::Unum)` forces the ubit of a unum to be 1."""
@universal make_ulp!(x::Unum) = (x.flags |= UNUM_UBIT_MASK; x)
@universal make_ulp(x::Unum) = make_exact!(copy(x))

################################################################################
## carry resolution

doc"""
  `Unums.resolve_carry!(carry::UInt64, ::Unum, exponent::Int64)` resolves a
  carry (invisible bit) value that may exceed one after calculation events.
  You should pass this function an exponent value that will be returned,
  appropriately modified.
"""
@universal function resolve_carry!(carry::UInt64, x::Unum, exponent::Int64)
  leftzeroes = clz(carry)
  if (leftzeroes < 0x003F) #less than 63 zeroes
    shift = 0x003F - leftzeroes
    rsh_and_set_ubit!(x, shift)
    #now copy the bits over from the carried segment.
    frac_copy_top!(x, (((o64 << shift) - o64) & carry) << (leftzeroes + o16))
    exponent += shift
  end
  (exponent > max_exponent(ESS)) && mmr!(x)  #set it to mmr, if the exponent is too large.
  (x.esize, x.exponent) = encode_exp(exponent)
  return exponent
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using Documenter, Plots

makedocs(modules=[Plots], doctest=false)

deploydocs(
    deps = Deps.pip("mkdocs==0.17.5",
        "mkdocs-material==2.9.4" ,"python-markdown-math", "pygments", "pymdown-extensions"),
    repo   = "github.com/JuliaPlots/PlotDocs.jl.git",
    julia  = "1.0",
    osname = "linux"
)


text:
# ------------------------------------------------------------
# Does it work with multiple arguments?
# It should becasue the arguments are slurped into a tuple
# which becomes the key of the cache.

function memoize(f)
    memo = Dict()
    (args...; kwargs...) -> begin
        x = (args, kwargs)
        if haskey(memo, x)
            return memo[x]
        else
            value = f(args...; kwargs...)
            memo[x] = value
            return value
        end
    end
end

# Simulate a slow function with positional arguments and keyword arguments
slow_op = (a, b = 2; c = 3, d) -> begin
    sleep(2)
    a + b + c + d
end

op = memoize(slow_op);

@time op(2, d = 5);
@time op(2, d = 5);
@time op(1, c = 4, d = 5);
@time op(1, c = 4, d = 5);
#=
julia> @time op(2, d = 5);
  2.006106 seconds (22 allocations: 864 bytes)

julia> @time op(2, d = 5);
  0.000024 seconds (17 allocations: 672 bytes)

julia> @time op(1, c = 4, d = 5);
  2.525788 seconds (255.28 k allocations: 12.640 MiB)

julia> @time op(1, c = 4, d = 5);
  0.000023 seconds (16 allocations: 624 bytes)
=#


text:
## ------------------------------------------------------
## Interface
## ------------------------------------------------------
invpseudohessian(x) = (t = typeof(x); error("Not defined for type $t", ))
momentmatrix(x) = (t = typeof(x); error("Not defined for type $t"))

numobs(x) = StatsBase.nobs(x)
dof_resid(x) = numobs(x) - length(coef(x))

## ------------------------------------------------------
## HAC
## ------------------------------------------------------
function setupkernelweights!(k::T, m) where T<:HAC
    cf = coef(m)
    resize!(k.weights, length(cf))
    fill!(k.weights, 1)
end

setupkernelweights!(k::T, m::F) where {T<:Smoothed, F} = return nothing

function StatsBase.vcov(k::T, m;
    prewhite=false, dof_adjustment::Bool=true, scale::Real=1) where T<:Union{HAC, Smoothed}
    return _vcov(k, m, prewhite, dof_adjustment, scale)
end

function _vcov(k::T, m, prewhite::Bool, dof_adjustment::Bool, scale::Real) where T<:Union{HAC, Smoothed}
    setupkernelweights!(k, m)
    B  = invpseudohessian(m)
    mm = momentmatrix(m)
    A = lrvar(k, mm; prewhite=prewhite, demean=false)
    scale *= (dof_adjustment ? numobs(m)^2/dof_resid(m) : numobs(m))
    V = Symmetric((B*A*B).*scale)
    return V
end

vcovmatrix(
    k::T,
    m,
    factorization=Cholesky;
    prewhite=false,
    dof_adjustment::Bool=true,
    scale::Real=1,
) where T<:Union{HAC, Smoothed} = _vcovmatrix(k, m, prewhite, dof_adjustment, scale, factorization)

function _vcovmatrix(
    k::T,
    m,
    prewhite::Bool,
    dof_adjustment::Bool,
    scale::Real,
    ::Type{Cholesky},
) where T<:Union{HAC, Smoothed}
    V = _vcov(k, m, prewhite, dof_adjustment, scale)
    return CovarianceMatrix(cholesky(V, check=true), k, V)
end

function _vcovmatrix(
    k::T,
    m,
    prewhite::Bool,
    dof_adjustment::Bool,
    scale::Real,
    ::Type{SVD},
) where T<:Union{HAC, Smoothed}
    V = _vcov(k, m, prewhite, dof_adjustment, scale)
    return CovarianceMatrix(svd(V), k, V)
end

## ------------------------------------------------------
## HC
## Note: adjfactor is always n/(n-k) for H1-H5 types
## ------------------------------------------------------
adjfactor(k::HC0, m, x) = one(first(x))
adjfactor(k::HC1, m, x) = numobs(m)./dof_resid(m)

function adjfactor(k::T, m, x) where T<:Union{HC2, HC3, HC4, HC4m, HC5}
    mm = typeof(m)
    @warn "$k variance is not defined for type $mm. Using `HC1`."
    numobs(m)./dof_resid(m)
end

adjust!(m, adj::AbstractFloat) = m
adjust!(m, adj::AbstractMatrix) = m.*sqrt.(adj)

StatsBase.vcov(k::HC, m::T; scale::Real=1) where T = _vcov(k, m, scale)

vcovmatrix(k::HC, m::T, factorization = Cholesky; scale::Real = 1)  where T = _vcovmatrix(k, m, scale, factorization)

function _vcovmatrix(k::HC, m::T, scale::Real, ::Type{Cholesky}) where T
    V = _vcov(k, m, scale)
    CovarianceMatrix(cholesky(V, check=true), k, V)
end

function _vcovmatrix(k::HC, m::T, scale::Real, ::Type{SVD}) where T
    V = _vcov(k, m, scale)
    CovarianceMatrix(svd(V), k, V)
end

## This needs to be here - because for GLM
## we used modmatrix and some other indirection
## to get care of weights and dispersion
## TODO: It could be fixed by more cleverly dealing
##       with the weights situation
function _vcov(k::HC, m::T, scale) where T
    B   = invpseudohessian(m)
    mm  = momentmatrix(m)
    adj = adjfactor(k, m, modelmatrix(m))  ## <- here `modelmatrix` in GLM.jl `modmatrix`
    mm  = adjust!(mm, adj)
    scale *= length(adj) > 1 ? one(first(adj)) : adj
    A = mm'*mm
    return Symmetric((B*A*B).*scale)
end

"""
Calculate sandwich type variance-covariance of an estimator

"""
function sandwich(k::RobustVariance, B::Matrix, momentmatrix::Matrix; demean = false, prewhite = false, dof = 0.0)
    K, m = size(B)
    n, m_m = size(momentmatrix)
    @assert m_m == m "number of rows of `momentmatrix` must be equal to the number of column of `B`"
    scale = n^2/(n-dof)
    A = lrvar(k, momentmatrix, demean = demean, scale = scale, prewhite = prewhite)   # df adjustment is built into vcov
    Symmetric(B*A*B)
end

## Standard errors
StatsBase.stderror(k::RobustVariance, m; kwargs...) = sqrt.(diag(vcov(k, m; kwargs...)))
StatsBase.stderror(k::Smoothed, m; kwargs...) = sqrt.(diag(vcov(k, m; kwargs...)))
StatsBase.stderror(v::CovarianceMatrix) = sqrt.(diag(v.V))


text:
function get_max_growth(obj::MultistepIntegrator)
    return jcall(obj, "getMaxGrowth", jdouble, ())
end

function get_max_step(obj::AdaptiveStepsizeIntegrator)
    return jcall(obj, "getMaxStep", jdouble, ())
end

function get_min_reduction(obj::MultistepIntegrator)
    return jcall(obj, "getMinReduction", jdouble, ())
end

function get_min_step(obj::AdaptiveStepsizeIntegrator)
    return jcall(obj, "getMinStep", jdouble, ())
end

function get_n_steps(obj::MultistepIntegrator)
    return jcall(obj, "getNSteps", jint, ())
end

function get_safety(obj::MultistepIntegrator)
    return jcall(obj, "getSafety", jdouble, ())
end

function get_starter_integrator(obj::MultistepIntegrator)
    return jcall(obj, "getStarterIntegrator", ODEIntegrator, ())
end

function initialize_step(obj::AdaptiveStepsizeIntegrator, arg0::jboolean, arg1::jint, arg2::Vector{jdouble}, arg3::ODEStateAndDerivative, arg4::EquationsMapper)
    return jcall(obj, "initializeStep", jdouble, (jboolean, jint, Vector{jdouble}, ODEStateAndDerivative, EquationsMapper), arg0, arg1, arg2, arg3, arg4)
end

function set_initial_step_size(obj::AdaptiveStepsizeIntegrator, arg0::jdouble)
    return jcall(obj, "setInitialStepSize", void, (jdouble,), arg0)
end

function set_max_growth(obj::MultistepIntegrator, arg0::jdouble)
    return jcall(obj, "setMaxGrowth", void, (jdouble,), arg0)
end

function set_min_reduction(obj::MultistepIntegrator, arg0::jdouble)
    return jcall(obj, "setMinReduction", void, (jdouble,), arg0)
end

function set_safety(obj::MultistepIntegrator, arg0::jdouble)
    return jcall(obj, "setSafety", void, (jdouble,), arg0)
end

function set_starter_integrator(obj::MultistepIntegrator, arg0::ODEIntegrator)
    return jcall(obj, "setStarterIntegrator", void, (ODEIntegrator,), arg0)
end

function set_step_size_control(obj::AdaptiveStepsizeIntegrator, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::Vector{jdouble})
    return jcall(obj, "setStepSizeControl", void, (jdouble, jdouble, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2, arg3)
end

function set_step_size_control(obj::AdaptiveStepsizeIntegrator, arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return jcall(obj, "setStepSizeControl", void, (jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end


text:
export Abstract_BC_Solver_Conf

@doc raw"""
```julia
abstract type Abstract_BC_Solver_Conf end
```

Abstract solver configuration. These are the solvers to be used to
solve bound constrained nonlinear least squares:

```math
\min\limits_{\theta_l \le \theta \le \theta_u } \frac{1}{2}\|r(\theta)\|^2
```

Implementations:
- [`LevenbergMarquardt_BC_Conf`](@ref)
"""
abstract type Abstract_BC_Solver_Conf end


text:
# run batch_molar_mass.jl first
# re-run this script if there are network issues.
#   The GISSMO server might be blocking you for downloading so many entries.
#   This script will skip downloading entries that were already downloaded.

#include("../src/GISSMOReader.jl")
import GISSMOReader
#import DelimitedFiles

#include("../src/coupling/parse2.jl")

entries = GISSMOReader.getGISSMOentriesall()
#entries = getGISSMOentries(["L-Lysine";])

# storage for the experiment 1D 1H files.
#base_dir = "/home/roy/MEGAsync/data/molecules"
base_dir = "/home/roy/Documents/repo/NMRData/src/input/GISSMO_data/"

# storage for the extract J-coupling and chemical shift values, stored in a JLD file per GISSMO entry.
save_dir = "/home/roy/Documents/repo/NMRData/src/input/molecules/"

δ_lb = 0.1
δ_ub = 0.1
unique_cs_tol = 1e-6

GISSMOReader.downloadGISSMOentries(entries,
base_dir,
save_dir;
δ_lb = δ_lb,
δ_ub = δ_ub,
unique_cs_tol = unique_cs_tol)

# julia batch_graphs.jl > /home/roy/MEGAsync/inputs/NMR/molecules/output.txt


text:
module jubart

include("jubart_functions.jl")

export juBart, predict_juBart, sim_friedman_simple

end # module


text:
sum1 = 0
sum2 = 0
for i = 1:100
    global sum1+=i
    global sum2+=i^2
end
sum1=sum1^2
println(sum1-sum2)


text:
module ContComplianceDemo1

using TopOpt, LinearAlgebra, StatsFuns
# using Makie
# using TopOpt.TopOptProblems.Visualization: visualize

E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0 # downward force
rmin = 4.0 # filter radius
xmin = 0.0001 # minimum density
problem_size = (60, 20)
V = 0.5 # maximum volume fraction
p = 4.0 # penalty
x0 = fill(V, prod(problem_size)) # initial design

problem = HalfMBB(Val{:Linear}, problem_size, (1.0, 1.0), E, v, f)

solver = FEASolver(Direct, problem; xmin=xmin)
cheqfilter = DensityFilter(solver; rmin=rmin)
comp = TopOpt.Compliance(problem, solver)

function obj(x)
    # minimize compliance
    return comp(cheqfilter(x))
end
function constr(x)
    # volume fraction constraint
    return sum(cheqfilter(x)) / length(x) - V
end

m = Model(obj)
addvar!(m, zeros(length(x0)), ones(length(x0)))
Nonconvex.add_ineq_constraint!(m, constr)

options = MMAOptions(; maxiter=1000, tol=Tolerance(; kkt=1e-4, x=1e-4, f=1e-4))
TopOpt.setpenalty!(solver, p)
# Method of Moving Asymptotes
@time r = Nonconvex.optimize(m, MMA87(), x0; options=options)

@show obj(r.minimizer)
@show constr(r.minimizer)
topology = cheqfilter(r.minimizer)
# fig = visualize(problem, solver.u;
#     topology = topology, default_exagg_scale=0.0, scale_range=10.0)
# Makie.display(fig)

end


text:
using MAT, Reachability, MathematicalSystems

SUITE["ISS"] = BenchmarkGroup()

# ==============================
# Load model
# ==============================
file = matopen(@relpath "iss.mat")
A, B, C = sparse(read(file, "A")), read(file, "B"), Matrix(read(file, "C")[3, :]')
n = size(A, 1)
Cvec = C[:]
time_horizon = 20.0
X0 = BallInf(zeros(n), 0.0001)

# ==============================
# Time-varying input
# ==============================
U = Hyperrectangle(low=[0.0, 0.8, 0.9], high=[0.1, 1., 1.])
S = ConstrainedLinearControlContinuousSystem(A, eye(A), nothing, ConstantInput(B * U))
problem_TV = InitialValueProblem(S, X0)

# ==============================
# ISU01 and ISS01
# ==============================
ISU01 = LinearConstraintProperty([Clause(LinearConstraint(Cvec, 0.0005)), Clause(LinearConstraint(-Cvec, 0.0005))])
ISS01 = LinearConstraintProperty([Clause(LinearConstraint(Cvec, 0.0007)), Clause(LinearConstraint(-Cvec, 0.0007))])

sol = solve(problem_TV, :T=>time_horizon, :δ=>5e-3, :vars=>136:270, :mode=>"check", :property=>ISU01, :projection_matrix=>C, :assume_sparse=>true)
@assert !sol.satisfied
SUITE["ISS"]["ISU01", "dense"] = @benchmarkable solve($problem_TV, :T=>$time_horizon, :δ=>5e-3, :vars=>136:270, :mode=>"check", :property=>$ISU01, :projection_matrix=>$C, :assume_sparse=>true)

sol = solve(problem_TV, :T=>time_horizon, :δ=>6e-4, :vars=>136:270, :mode=>"check", :property=>ISS01, :projection_matrix=>C, :assume_sparse=>true, :lazy_inputs_interval=>-1, :partition=>[1:135, 136:270])
@assert sol.satisfied
SUITE["ISS"]["ISS01", "dense"] = @benchmarkable solve($problem_TV, :T=>$time_horizon, :δ=>6e-4, :vars=>136:270, :mode=>"check", :property=>$ISS01, :projection_matrix=>$C, :assume_sparse=>true, :lazy_inputs_interval=>-1, :partition=>[1:135, 136:270])

sol = solve(problem_TV, :approx_model=>"nobloating", :T=>time_horizon, :δ=>5e-3, :vars=> 136:270, :mode=>"check", :property=>ISU01, :projection_matrix=>C, :assume_sparse=>true)
@assert !sol.satisfied
SUITE["ISS"]["ISU01", "discrete"] = @benchmarkable solve($problem_TV, :approx_model=>"nobloating", :T=>$time_horizon, :δ=>5e-3, :vars=> 136:270, :mode=>"check", :property=>$ISU01, :projection_matrix=>$C, :assume_sparse=>true)

sol = solve(problem_TV, :approx_model=>"nobloating", :T=>time_horizon, :δ=>5e-3, :vars=>136:270, :mode=>"check", :property=>ISS01, :projection_matrix=>C, :assume_sparse=>true, :lazy_inputs_interval=>-1, :partition=>[1:135, 136:270])
@assert sol.satisfied
SUITE["ISS"]["ISS01", "discrete"] = @benchmarkable solve($problem_TV, :approx_model=>"nobloating", :T=>$time_horizon, :δ=>5e-3, :vars=>136:270, :mode=>"check", :property=>$ISS01, :projection_matrix=>$C, :assume_sparse=>true, :lazy_inputs_interval=>-1, :partition=>[1:135, 136:270])

# ==============================
# Constant input
# ==============================
import Reachability.add_dimension
A = sparse(read(file, "A"))
Aext = add_dimension(A, 3)
Aext[1:270, 271:273] = B
S = LinearContinuousSystem(Aext)
X0 = X0 * U
problem_CONST = InitialValueProblem(S, X0)
C = hcat(C, [0.0 0.0 0.0])
Cvec = C[1, :]

# ==============================
# ISU02 and ISS02
# ==============================
ISU02 = LinearConstraintProperty([Clause(LinearConstraint(Cvec, 0.00017)), Clause(LinearConstraint(-Cvec, 0.00017))])
ISS02 = LinearConstraintProperty([Clause(LinearConstraint(Cvec, 0.0005)), Clause(LinearConstraint(-Cvec, 0.0005))])

sol = solve(problem_CONST, :T=>time_horizon, :δ=>5e-3, :vars=>136:270, :mode=>"check", :property=>ISU02, :projection_matrix=>C, :assume_sparse=>true, :assume_homogeneous=>true)
@assert !sol.satisfied
SUITE["ISS"]["ISU02", "dense"] = @benchmarkable solve($problem_CONST, :T=>$time_horizon, :δ=>5e-3, :vars=>136:270, :mode=>"check", :property=>$ISU02, :projection_matrix=>$C, :assume_sparse=>true, :assume_homogeneous=>true)

sol = solve(problem_CONST, :T=>time_horizon, :δ=>5e-3, :vars=> 136:270, :mode=>"check", :property=>ISS02, :projection_matrix=>C, :assume_sparse=>true, :assume_homogeneous=>true, :lazy_inputs_interval=>-1, :partition=>[1:135, 136:270])
@assert sol.satisfied
SUITE["ISS"]["ISS02", "dense"] = @benchmarkable solve($problem_CONST, :T=>$time_horizon, :δ=>5e-3, :vars=> 136:270, :mode=>"check", :property=>$ISS02, :projection_matrix=>$C, :assume_sparse=>true, :assume_homogeneous=>true, :lazy_inputs_interval=>-1, :partition=>[1:135, 136:270])

sol = solve(problem_CONST, :approx_model=>"nobloating", :T=>time_horizon, :δ=>5e-3, :vars=>136:270, :mode=>"check", :property=>ISU02, :projection_matrix=>C, :assume_sparse=>true, :assume_homogeneous=>true)
@assert !sol.satisfied
SUITE["ISS"]["ISU02", "discrete"] = @benchmarkable solve($problem_CONST, :approx_model=>"nobloating", :T=>$time_horizon, :δ=>5e-3, :vars=>136:270, :mode=>"check", :property=>$ISU02, :projection_matrix=>$C, :assume_sparse=>true, :assume_homogeneous=>true)

sol = solve(problem_CONST, :approx_model=>"nobloating", :T=>time_horizon, :δ=>5e-3, :vars=> 136:270, :mode=>"check", :property=>ISS02, :projection_matrix=>C, :assume_sparse=>true, :assume_homogeneous=>true, :lazy_inputs_interval=>-1, :partition=>[1:135, 136:270])
@assert sol.satisfied
SUITE["ISS"]["ISS02", "discrete"] = @benchmarkable solve($problem_CONST, :approx_model=>"nobloating", :T=>$time_horizon, :δ=>5e-3, :vars=> 136:270, :mode=>"check", :property=>$ISS02, :projection_matrix=>$C, :assume_sparse=>true, :assume_homogeneous=>true, :lazy_inputs_interval=>-1, :partition=>[1:135, 136:270])


text:
using Test
using MIMEs

sub(s) = SubString(s, 1)

@test mime_from_path("a/foo.txt") === MIME"text/plain"()
@test mime_from_path(sub("a/foo.txt")) === MIME"text/plain"()
@test mime_from_path("a/fo🌟oê.txt") === MIME"text/plain"()
@test mime_from_path("a/foo.txt") === MIME"text/plain"()
@test mime_from_path("a/foo.json") === MIME"application/json"()
@test mime_from_extension(".json") === MIME"application/json"()
@test mime_from_extension("json") === MIME"application/json"()
@test mime_from_extension("JSON") === MIME"application/json"()
@test mime_from_extension(sub("JSON")) === MIME"application/json"()
@test extension_from_mime(MIME"application/json"()) == ".json"
@test extension_from_mime(MIME"application/x-asfdafd"()) == ""
@test extension_from_mime(MIME"application/x-asfdafd"(), ".a") == ".a"

@test charset_from_mime(MIME"application/x-asfdafd"()) == nothing
@test charset_from_mime(MIME"text/x-asfdafd"()) == "UTF-8"
@test charset_from_mime(MIME"text/plain"()) == "UTF-8"
@test charset_from_mime(MIME"text/html"()) == "UTF-8"
@test contenttype_from_mime(MIME"application/x-asfdafd"()) == "application/x-asfdafd"
@test contenttype_from_mime(MIME"text/html"()) == "text/html; charset=utf-8"
@test contenttype_from_mime(MIME"text/asdfasdfasdf"()) == "text/asdfasdfasdf; charset=utf-8"

@test mime_from_extension("js") == MIME"text/javascript"()
@test charset_from_mime(mime_from_extension("js")) == "UTF-8"
@test mime_from_extension("jl") == MIME"text/julia"()
@test charset_from_mime(mime_from_extension("jl")) == "UTF-8"

@test compressible_from_mime(MIME"text/html"())
@test compressible_from_mime(MIME"text/css"())
@test !compressible_from_mime(MIME"text/x-sadfa"())

@test contenttype_from_mime(MIME"application/json"()) == "application/json; charset=utf-8"
@test contenttype_from_mime(MIME"application/x-bogus"()) == "application/x-bogus"
@test contenttype_from_mime(mime_from_extension(".png", MIME"application/octet-stream"())) == "image/png"

# from https://developer.mozilla.org/en-US/docs/Web/HTTP/Basics_of_HTTP/MIME_types/Common_types
const mdn = Dict(
    ".bin" => "application/octet-stream",
    ".bmp" => "image/bmp",
    ".css" => "text/css",
    ".csv" => "text/csv",
    ".eot" => "application/vnd.ms-fontobject",
    ".gz" => "application/gzip",
    ".gif" => "image/gif",
    ".htm" => "text/html",
    ".html" => "text/html",
    ".ico" => "image/vnd.microsoft.icon",
    ".jpeg" => "image/jpeg",
    ".jpg" => "image/jpeg",
    ".js" => "text/javascript",
    ".json" => "application/json",
    ".jsonld" => "application/ld+json",
    ".mjs" => "text/javascript",
    ".mp3" => "audio/mpeg",
    ".mp4" => "video/mp4",
    ".mpeg" => "video/mpeg",
    ".ogx" => "application/ogg",  ".otf" => "font/otf",
    ".png" => "image/png",
    ".pdf" => "application/pdf",
    ".rtf" => "application/rtf",
    ".sh" => "application/x-sh",
    ".svg" => "image/svg+xml",
    ".tar" => "application/x-tar",
    ".tif" => "image/tiff",
    ".tiff" => "image/tiff",
    ".ttf" => "font/ttf",
    ".txt" => "text/plain",
    ".weba" => "audio/webm",
    ".webm" => "video/webm",
    ".webp" => "image/webp",
    ".woff" => "font/woff",
    ".woff2" => "font/woff2",
    ".xhtml" => "application/xhtml+xml",
    ".xml" => "application/xml",
    ".xul" => "application/vnd.mozilla.xul+xml",
    ".zip" => "application/zip",
    ".wasm" => "application/wasm",
)

for (ex, ms) in mdn
    @test mime_from_extension(ex) === MIME(ms)
end

# Mismatched with MDN, but that's fine:
# ".aac" => "audio/aac",
# ".wav" => "audio/wav",
# ".oga" => "audio/ogg", ".ogv" => "video/ogg",
# ".opus" => "audio/opus",

max_stars_repo_name: fonsp/MIMEs

max_issues_repo_name: fonsp/MIMEs

max_forks_repo_name: fonsp/MIMEs


text:
export Handler,
    Context

const EVENT_TYPES = Dict{String, Symbol}(
    "READY"                       => :Ready,
    "RESUMED"                     => :Resumed,
    "CHANNEL_CREATE"              => :ChannelCreate,
    "CHANNEL_UPDATE"              => :ChannelUpdate,
    "CHANNEL_DELETE"              => :ChannelDelete,
    "CHANNEL_PINS_UPDATE"         => :ChannelPinsUpdate,
    "GUILD_CREATE"                => :GuildCreate,
    "GUILD_UPDATE"                => :GuildUpdate,
    "GUILD_DELETE"                => :GuildDelete,
    "GUILD_BAN_ADD"               => :GuildBanAdd,
    "GUILD_BAN_REMOVE"            => :GuildBanRemove,
    "GUILD_EMOJIS_UPDATE"         => :GuildEmojisUpdate,
    "GUILD_INTEGRATIONS_UPDATE"   => :GuildIntegrationsUpdate,
    "GUILD_MEMBER_ADD"            => :GuildMemberAdd,
    "GUILD_MEMBER_REMOVE"         => :GuildMemberRemove,
    "GUILD_MEMBER_UPDATE"         => :GuildMemberUpdate,
    "GUILD_MEMBERS_CHUNK"         => :GuildMembersChunk,
    "GUILD_ROLE_CREATE"           => :GuildRoleCreate,
    "GUILD_ROLE_UPDATE"           => :GuildRoleUpdate,
    "GUILD_ROLE_DELETE"           => :GuildRoleDelete,
    "MESSAGE_CREATE"              => :MessageCreate,
    "MESSAGE_UPDATE"              => :MessageUpdate,
    "MESSAGE_DELETE"              => :MessageDelete,
    "MESSAGE_DELETE_BULK"         => :MessageDeleteBulk,
    "MESSAGE_REACTION_ADD"        => :MessageReactionAdd,
    "MESSAGE_REACTION_REMOVE"     => :MessageReactionRemove,
    "MESSAGE_REACTION_REMOVE_ALL" => :MessageReactionRemoveAll,
    "INTERACTION_CREATE"          => :InteractionCreate,
    "PRESENCE_UPDATE"             => :PresenceUpdate,
    "TYPING_START"                => :TypingStart,
    "USER_UPDATE"                 => :UserUpdate,
    "VOICE_STATE_UPDATE"          => :VoiceStateUpdate,
    "VOICE_SERVER_UPDATE"         => :VoiceServerUpdate,
    "WEBHOOKS_UPDATE"             => :WebhooksUpdate,
)

"""
    Handler(
        f::Function
        d::Dict{Symbol, Any}
    )

Handler is a wrapper for a `Dict{Symbol, Any}` that also contains a function.
"""
struct Handler
    f::Function
    d::Dict{Symbol, Any}
end

"""
    Handler(; kwargs...) -> Handler

Generates a handler based on kwargs and a function.
"""
Handler(f; kwargs...) = Handler(f, Dict(kwargs))

"""
    Context(
        data::Dict{Symbol, Any}
    )

Context is a wrapper for a `Dict{Symbol, Any}` with some special functionality.
"""
struct Context
    data::Dict{Symbol, Any}
end
"""
    Context(; kwargs...) -> Context

Generates a context based on kwargs.
"""
Context(; kwargs...) = Context(Dict(kwargs))

quickdoc(name::String, ar::String) = "
    $(name)(
        f::Function
        c::Client
    )

Adds a handler for the $(replace(ar, "_"=>"\\_")) gateway event.
The `f` parameter's signature should be:
```
    (ctx::Context) -> Any
```
"

context(data::Dict{Symbol, Any}) = Context(data)
"""
    context(t::Symbol, data::Dict{Symbol, Any}) -> Context

Checks if the Context needs to be created in a special way based on the event provided by `t`.\n
Then, returns the generated context.
"""
function context(t::Symbol, data::Dict{Symbol, Any})
    t == :OnMessageCreate && return Context(; message=Message(data))
    t ∈ [:OnGuildCreate, :OnGuildUpdate] && return Context(; guild=Guild(data))
    t ∈ [:OnMessageReactionAdd, :OnMessageReactionRemove] && return Context(;
        emoji=make(Emoji, data, :emoji),
        channel=DiscordChannel(; id=data[:channel_id]),
        message=Message(; id=data[:message_id], channel_id=data[:channel_id]),
        data...)
    t  == :OnReady && return Context(; user=make(User, data, :user), data...)
    t == :OnInteractionCreate && return Context(; interaction=Interaction(data))
    Context(data)
end

make(::Type{T}, data::Dict{Symbol, Any}, k::Symbol) where T <: DiscordObject = T(pop!(data, k, missing))

for (k, v) in EVENT_TYPES
    nm = Symbol("On"*String(v))
    hm = Symbol("on_"*lowercase(k)*"!")
    @eval ($nm)(f::Function; kwargs...) = Handler(f; type=Symbol($nm), kwargs...)
    k = quote
        @doc quickdoc(string($hm), $k)
        ($hm)(f::Function, c; kwargs...) = add_handler!(c, ($nm)(f, kwargs...))
        export $hm
    end
    eval(k)
end

# Deprecated
# ----------------------------------------
# |    Removed: 1.0+                     |
# |    Added 0.3                         |
# |    Replaced by on_message_create!    |
# ----------------------------------------
const on_message! = on_message_create!
const on_reaction_add! = on_message_reaction_add!
const on_reaction_remove! = on_message_reaction_remove!

Base.getproperty(ctx::Context, sym::Symbol) = getfield(ctx, :data)[sym]
Base.hasproperty(ctx::Context, sym::Symbol) = haskey(getfield(ctx, :data), sym)

Base.getproperty(h::Handler, sym::Symbol) = sym != :f ? getfield(h, :d)[sym] : getfield(h, sym)
Base.hasproperty(h::Handler, sym::Symbol) = haskey(getfield(h, :d), sym)

handlerkind(h::Handler) = h.type
handlerargs(f::Function) = method_args(first(methods(f)))


text:
"""
    rstar(
        rng=Random.GLOBAL_RNG,
        classifier,
        samples::AbstractMatrix,
        chain_indices::AbstractVector{Int};
        subset::Real=0.8,
        verbosity::Int=0,
    )

Compute the ``R^*`` convergence statistic of the `samples` with shape (draws, parameters)
and corresponding chains `chain_indices` with the `classifier`.

This implementation is an adaption of algorithms 1 and 2 described by Lambert and Vehtari.

The `classifier` has to be a supervised classifier of the MLJ framework (see the
[MLJ documentation](https://alan-turing-institute.github.io/MLJ.jl/dev/list_of_supported_models/#model_list)
for a list of supported models). It is trained with a `subset` of the samples. The training
of the classifier can be inspected by adjusting the `verbosity` level.

If the classifier is deterministic, i.e., if it predicts a class, the value of the ``R^*``
statistic is returned (algorithm 1). If the classifier is probabilistic, i.e., if it outputs
probabilities of classes, the scaled Poisson-binomial distribution of the ``R^*`` statistic
is returned (algorithm 2).

!!! note
    The correctness of the statistic depends on the convergence of the `classifier` used
    internally in the statistic.

# Examples

```jldoctest rstar; setup = :(using Random; Random.seed!(100))
julia> using MLJBase, MLJXGBoostInterface, Statistics

julia> samples = fill(4.0, 300, 2);

julia> chain_indices = repeat(1:3; outer=100);
```

One can compute the distribution of the ``R^*`` statistic (algorithm 2) with the
probabilistic classifier.

```jldoctest rstar
julia> distribution = rstar(XGBoostClassifier(), samples, chain_indices);

julia> isapprox(mean(distribution), 1; atol=0.1)
true
```

For deterministic classifiers, a single ``R^*`` statistic (algorithm 1) is returned.
Deterministic classifiers can also be derived from probabilistic classifiers by e.g.
predicting the mode. In MLJ this corresponds to a pipeline of models.

```jldoctest rstar
julia> xgboost_deterministic = Pipeline(XGBoostClassifier(); operation=predict_mode);

julia> value = rstar(xgboost_deterministic, samples, chain_indices);

julia> isapprox(value, 1; atol=0.2)
true
```

# References

Lambert, B., & Vehtari, A. (2020). ``R^*``: A robust MCMC convergence diagnostic with uncertainty using decision tree classifiers.
"""
function rstar(
    rng::Random.AbstractRNG,
    classifier::MLJModelInterface.Supervised,
    x::AbstractMatrix,
    y::AbstractVector{Int};
    subset::Real=0.8,
    verbosity::Int=0,
)
    # checks
    size(x, 1) != length(y) && throw(DimensionMismatch())
    0 < subset < 1 || throw(ArgumentError("`subset` must be a number in (0, 1)"))

    # randomly sub-select training and testing set
    N = length(y)
    Ntrain = round(Int, N * subset)
    0 < Ntrain < N ||
        throw(ArgumentError("training and test data subsets must not be empty"))
    ids = Random.randperm(rng, N)
    train_ids = view(ids, 1:Ntrain)
    test_ids = view(ids, (Ntrain + 1):N)

    # train classifier on training data
    ycategorical = MLJModelInterface.categorical(y)
    fitresult, _ = MLJModelInterface.fit(
        classifier, verbosity, Tables.table(x[train_ids, :]), ycategorical[train_ids]
    )

    # compute predictions on test data
    xtest = Tables.table(x[test_ids, :])
    predictions = MLJModelInterface.predict(classifier, fitresult, xtest)

    # compute statistic
    ytest = ycategorical[test_ids]
    result = _rstar(predictions, ytest)

    return result
end

function rstar(
    classif::MLJModelInterface.Supervised,
    x::AbstractMatrix,
    y::AbstractVector{Int};
    kwargs...,
)
    return rstar(Random.GLOBAL_RNG, classif, x, y; kwargs...)
end

# R⋆ for deterministic predictions (algorithm 1)
function _rstar(predictions::AbstractVector{T}, ytest::AbstractVector{T}) where {T}
    length(predictions) == length(ytest) ||
        error("numbers of predictions and targets must be equal")
    mean_accuracy = Statistics.mean(p == y for (p, y) in zip(predictions, ytest))
    nclasses = length(MLJModelInterface.classes(ytest))
    return nclasses * mean_accuracy
end

# R⋆ for probabilistic predictions (algorithm 2)
function _rstar(predictions::AbstractVector, ytest::AbstractVector)
    length(predictions) == length(ytest) ||
        error("numbers of predictions and targets must be equal")

    # create Poisson binomial distribution with support `0:length(predictions)`
    distribution = Distributions.PoissonBinomial(map(Distributions.pdf, predictions, ytest))

    # scale distribution to support in `[0, nclasses]`
    nclasses = length(MLJModelInterface.classes(ytest))
    scaled_distribution = (nclasses//length(predictions)) * distribution

    return scaled_distribution
end


text:
# generated from SQSTypeDict via AWS/src/codegen.jl
# do not edit, edit specification file and regenerate instead

type BatchResultErrorEntry <: AbstractAWSType
    code::Union{Void,String}
    id::Union{Void,String}
    message::Union{Void,String}
    senderFault::Union{Void,Bool}
    function BatchResultErrorEntry(;code=nothing, id=nothing, message=nothing, senderFault=nothing)
        new(code, id, message, senderFault)
    end
    BatchResultErrorEntry(pd) = parse_from_xml(BatchResultErrorEntry, SQSTypeDict, pd)
end # BatchResultErrorEntry
export BatchResultErrorEntry

type ChangeMessageVisibilityBatchRequestEntry <: AbstractAWSType
    id::Union{Void,String}
    receiptHandle::Union{Void,String}
    visibilityTimeout::Union{Void,Int64}
    function ChangeMessageVisibilityBatchRequestEntry(;id=nothing, receiptHandle=nothing, visibilityTimeout=nothing)
        new(id, receiptHandle, visibilityTimeout)
    end
    ChangeMessageVisibilityBatchRequestEntry(pd) = parse_from_xml(ChangeMessageVisibilityBatchRequestEntry, SQSTypeDict, pd)
end # ChangeMessageVisibilityBatchRequestEntry
export ChangeMessageVisibilityBatchRequestEntry

type ChangeMessageVisibilityBatchResultEntry <: AbstractAWSType
    id::Union{Void,String}
    function ChangeMessageVisibilityBatchResultEntry(;id=nothing)
        new(id)
    end
    ChangeMessageVisibilityBatchResultEntry(pd) = parse_from_xml(ChangeMessageVisibilityBatchResultEntry, SQSTypeDict, pd)
end # ChangeMessageVisibilityBatchResultEntry
export ChangeMessageVisibilityBatchResultEntry

type DeleteMessageBatchRequestEntry <: AbstractAWSType
    id::Union{Void,String}
    receiptHandle::Union{Void,String}
    function DeleteMessageBatchRequestEntry(;id=nothing, receiptHandle=nothing)
        new(id, receiptHandle)
    end
    DeleteMessageBatchRequestEntry(pd) = parse_from_xml(DeleteMessageBatchRequestEntry, SQSTypeDict, pd)
end # DeleteMessageBatchRequestEntry
export DeleteMessageBatchRequestEntry

type DeleteMessageBatchResultEntry <: AbstractAWSType
    id::Union{Void,String}
    function DeleteMessageBatchResultEntry(;id=nothing)
        new(id)
    end
    DeleteMessageBatchResultEntry(pd) = parse_from_xml(DeleteMessageBatchResultEntry, SQSTypeDict, pd)
end # DeleteMessageBatchResultEntry
export DeleteMessageBatchResultEntry

type Message <: AbstractAWSType
    attribute::Union{Void,Dict{String,String}}
    body::Union{Void,String}
    mD5OfBody::Union{Void,String}
    mD5OfMessageAttributes::Union{Void,String}
    messageAttribute::Union{Void,Dict{String,String}}
    messageId::Union{Void,String}
    receiptHandle::Union{Void,String}
    function Message(;attribute=nothing, body=nothing, mD5OfBody=nothing, mD5OfMessageAttributes=nothing, messageAttribute=nothing, messageId=nothing, receiptHandle=nothing)
        new(attribute, body, mD5OfBody, mD5OfMessageAttributes, messageAttribute, messageId, receiptHandle)
    end
    Message(pd) = parse_from_xml(Message, SQSTypeDict, pd)
end # Message
export Message

type MessageAttributeValue <: AbstractAWSType
    binaryListValue::Union{Void,Array{Array{UInt8,1},1}}
    binaryValue::Union{Void,Array{UInt8,1}}
    dataType::Union{Void,String}
    stringListValue::Union{Void,Array{String,1}}
    stringValue::Union{Void,String}
    function MessageAttributeValue(;binaryListValue=nothing, binaryValue=nothing, dataType=nothing, stringListValue=nothing, stringValue=nothing)
        new(binaryListValue, binaryValue, dataType, stringListValue, stringValue)
    end
    MessageAttributeValue(pd) = parse_from_xml(MessageAttributeValue, SQSTypeDict, pd)
end # MessageAttributeValue
export MessageAttributeValue

type SendMessageBatchRequestEntry <: AbstractAWSType
    delaySeconds::Union{Void,Int64}
    id::Union{Void,String}
    messageAttribute::Union{Void,Dict{String,MessageAttributeValue}}
    messageBody::Union{Void,String}
    messageDeduplicationId::Union{Void,String}
    messageGroupId::Union{Void,String}
    function SendMessageBatchRequestEntry(;delaySeconds=nothing, id=nothing, messageAttribute=nothing, messageBody=nothing, messageDeduplicationId=nothing, messageGroupId=nothing)
        new(delaySeconds, id, messageAttribute, messageBody, messageDeduplicationId, messageGroupId)
    end
    SendMessageBatchRequestEntry(pd) = parse_from_xml(SendMessageBatchRequestEntry, SQSTypeDict, pd)
end # SendMessageBatchRequestEntry
export SendMessageBatchRequestEntry

type SendMessageBatchResultEntry <: AbstractAWSType
    id::Union{Void,String}
    mD5OfMessageAttributes::Union{Void,String}
    mD5OfMessageBody::Union{Void,String}
    messageId::Union{Void,String}
    sequenceNumber::Union{Void,String}
    function SendMessageBatchResultEntry(;id=nothing, mD5OfMessageAttributes=nothing, mD5OfMessageBody=nothing, messageId=nothing, sequenceNumber=nothing)
        new(id, mD5OfMessageAttributes, mD5OfMessageBody, messageId, sequenceNumber)
    end
    SendMessageBatchResultEntry(pd) = parse_from_xml(SendMessageBatchResultEntry, SQSTypeDict, pd)
end # SendMessageBatchResultEntry
export SendMessageBatchResultEntry


text:
module TestUtil
using Test
using OptimizationAlgorithms: CircularBuffer, _circshift!

@testset "CircularBuffer" begin
    cb = CircularBuffer{Float64}(2)
    push!(cb, 1.)
    push!(cb, 2.)
    @test _circshift!(cb) ≈ [2., 1.]
end

end # TestUtil


text:
# issue #11551 (#11004,#10959)
function tstcvt(strUTF8::UTF8Str, strUTF16::UTF16Str, strUTF32::UTF32Str)
    @test utf16(strUTF8) == strUTF16
    @test utf32(strUTF8) == strUTF32
    @test utf8(strUTF16) == strUTF8
    @test utf32(strUTF16) == strUTF32
    @test utf8(strUTF32)  == strUTF8
    @test utf16(strUTF32) == strUTF16
end

# Create some ASCII, UTF8, UTF16, and UTF32 strings
let strAscii = to_ascii("abcdefgh"),
    strA_UTF8 = utf8(("abcdefgh\uff")[1:8]),
    strL_UTF8 = utf8("abcdef\uff\uff"),
    str2_UTF8 = utf8("abcd\uff\uff\u7ff\u7ff"),
    str3_UTF8 = utf8("abcd\uff\uff\u7fff\u7fff"),
    str4_UTF8 = utf8("abcd\uff\u7ff\u7fff\U7ffff"),
    #strS_UTF8 = UTF8Str([0xc3, 0xbf, 0xdf, 0xbf, 0xe7, 0xbf, 0xbf, 0xed, 0xa0, 0x80, 0xed, 0xb0, 0x80])
    strC_UTF8 = UTF8Str(vcat([0xc3, 0xbf, 0xdf, 0xbf, 0xe7, 0xbf, 0xbf], codeunits("\U10000")))
    #strz_UTF8 = UTF8Str([0xc3, 0xbf, 0xdf, 0xbf, 0xe7, 0xbf, 0xbf, 0x00])
    #strZ      = ['a'%UInt8, 'b'%UInt8, 'c'%UInt8, 'd'%UInt8, 0xc3, 0xbf, 0xdf, 0xbf, 0xe7, 0xbf, 0xbf, 0xc0, 0x80]

    strA_UTF16 = utf16(strA_UTF8)
    strL_UTF16 = utf16(strL_UTF8)
    str2_UTF16 = utf16(str2_UTF8)
    str3_UTF16 = utf16(str3_UTF8)
    str4_UTF16 = utf16(str4_UTF8)
    #strS_UTF16 = utf16(strS_UTF8)

    strA_UTF32 = utf32(strA_UTF8)
    strL_UTF32 = utf32(strL_UTF8)
    str2_UTF32 = utf32(str2_UTF8)
    str3_UTF32 = utf32(str3_UTF8)
    str4_UTF32 = utf32(str4_UTF8)
    #strS_UTF32 = utf32(strS_UTF8)

    @test utf8(strAscii) == strAscii
    @test utf16(strAscii) == strAscii
    @test utf32(strAscii) == strAscii

    tstcvt(strA_UTF8,strA_UTF16,strA_UTF32)
    tstcvt(strL_UTF8,strL_UTF16,strL_UTF32)
    tstcvt(str2_UTF8,str2_UTF16,str2_UTF32)
    tstcvt(str3_UTF8,str3_UTF16,str3_UTF32)
    tstcvt(str4_UTF8,str4_UTF16,str4_UTF32)

    #=
    # Test converting surrogate pairs
    @test utf16(strS_UTF8) == strC_UTF8
    @test utf32(strS_UTF8) == strC_UTF8
    @test utf8(strS_UTF16) == strC_UTF8
    @test utf32(strS_UTF16) == strC_UTF8
    @test utf8(strS_UTF32)  == strC_UTF8
    @test utf16(strS_UTF32) == strC_UTF8

    # Test converting overlong \0
    @test utf8(strZ)  == strz_UTF8
    @test utf16(UTF8Str(strZ)) == strz_UTF8
    @test utf32(UTF8Str(strZ)) == strz_UTF8
    =#
end