text:
"""

     H = freqresp(sys, ω)

Compute the frequency response `H` of the descriptor system `sys = (A-λE,B,C,D)` at the real frequencies `ω`.

For a system with `p` outputs and `m` inputs, and for `N` frequency values in `ω`,
`H` is a `p × m × N` array, where `H[:,:,i]` contains the `i`-th value of the frequency response computed as
`C*((w[i]*E - A)^-1)*B + D`, where `w[i] = im*ω[i]` for a continuous-time system and `w[i] = exp(im*ω[i]*|Ts|)`
for a discrete-time system with sampling time `Ts`.

_Method:_ For an efficient evaluation of `C*((w[i]*E - A)^-1)*B + D` for many values of `w[i]`, a preliminary
orthogonal coordinate transformation is performed such that for the input-output equivalent transformed
system `sys = (At-λEt,Bt,Ct,Dt)`, the matrix `w[i]*Et - At` is upper Hessenberg.
This allows an efficient computation of the frequency response using the Hessenberg-form based
approach proposed in [1].

_Reference:_

[1] Laub, A.J., "Efficient Multivariable Frequency Response Computations",
    IEEE Transactions on Automatic Control, AC-26 (1981), pp. 407-408.
"""
function freqresp(sys::DescriptorStateSpace{T}, ω::Union{AbstractVector{<:Real},Real}) where T
    T1 = T <: BlasFloat ? T : promote_type(Float64,T)
    ONE = one(T1)
    typeof(ω) <: Real && (ω = [ω])
    N = length(ω)
    a, e, b, c, d = dssdata(T1,sys)
    if isempty(a)
        # overcome bug in orghr! for null dimension #43680
        H = similar(d, eltype(d), sys.ny, sys.nu, N)
        for i = 1:N
           copyto!(view(H,:,:,i),d)
       end
       return complex(H)
    end

    Ts = abs(sys.Ts)
    disc = !iszero(Ts)
    # create the complex frequency scaling
    sw = disc ? im*abs(sys.Ts) : im
    desc = !(e == I)
    # reduce to complex Hessenberg form
    ac, ec, bc, cc, dc = chess(a, e, b, c, d)
    H = similar(dc, eltype(dc), sys.ny, sys.nu, N)
    bct = similar(bc)
    for i = 1:N
        if isinf(ω[i])
           # exceptional call to evalfr
           H[:,:,i] = evalfr(sys, fval=Inf)
        else
           Hi = view(H,:,:,i)
           copyto!(Hi,dc)
           copyto!(bct,bc)
           w = disc ? -exp(sw*ω[i]) : -sw*ω[i]
           desc ? ldiv!(UpperHessenberg(ac+w*ec),bct) : ldiv!(ac,bct,shift = w)
           mul!(Hi,cc,bct,-ONE,ONE)
        end
    end
    return H
end
function chess(a::MT, e::Union{MT,UniformScaling}, b::MT, c::MT, d::MT) where {T <: BlasFloat, MT<: AbstractMatrix{T}}
    # reduce the descriptor system (A-λE,B,C,D) to an equivalent complex descriptor system
    # (Ac-λEc,Bc,Cc,Dc) such that Ac-λEc is in an upper Hessenberg form
    # Note: This function is only used in the context of efficient frequency response computations.
    isempty(a) && (return complex(a), e == I ? e : complex(e), complex(b), complex(c), complex(d))
    desc = !(e == I)
    if desc
        # Reduce (A,E) to (generalized) upper-Hessenberg form for
        # fast frequency response computation
        at = copy(a)
        et = copy(e)
        bt = copy(b)
        # first reduce E to upper triangular form
        _, tau = LinearAlgebra.LAPACK.geqrf!(et)
        T <: Complex ? tran = 'C' : tran = 'T'
        LinearAlgebra.LAPACK.ormqr!('L', tran, et, tau, at)
        LinearAlgebra.LAPACK.ormqr!('L', tran, et, tau, bt)
        # reduce A to Hessenberg form and keep E upper triangular
        _, _, Q, Z = MatrixPencils.gghrd!('I', 'I',  1, size(at,1), at, et, similar(at),similar(et))
        if T <: Complex
           bc = Q'*bt; cc = c*Z; ac = at; ec = et; dc = d;
        else
           bc = complex(Q'*bt); cc = complex(c*Z); ac = complex(at); ec = complex(et); dc = complex(d);
        end
    else
        # Reduce A to upper Hessenberg form for
        # fast frequency response computation
        Ha = hessenberg(a)
        ac = Ha.H
        if T <: Complex
           bc = Ha.Q'*b; cc = c*Ha.Q; dc = d;
        else
           bc = complex(Ha.Q'*b); cc = complex(c*Ha.Q); dc = complex(d);
        end
        ec = e
    end
    return ac, ec, bc, cc, dc
end
"""
    nx = order(sys)

Return the order `nx` of the descriptor system `sys` as the dimension of the state variable vector.
For improper or non-minimal systems, `nx` is less than the McMillan degree of the system.
"""
function order(sys::DescriptorStateSpace)
    return sys.nx
end
function Base.ndims(::DescriptorStateSpace)
# """
#     ndims(sys)

# Return the number of dimensions of the descriptor system array `sys` (1 for a single model).
# """
   return 1
end
"""
    size(sys) -> (p,m)
    size(sys,1) -> p
    size(sys,2) -> m

Return the number of outputs `p` and the number of inputs `m` of a descriptor system `sys`.
"""
function Base.size(sys::DescriptorStateSpace)
    size(sys.D)
end
function Base.size(sys::DescriptorStateSpace, d::Integer)
    d <= 2 ? size(sys)[d] : 1
end
function Base.length(sys::DescriptorStateSpace)
# """
#     length(sys)

# Return the number of elements of the descriptor system array `sys` (1 for a single model).
# """
   return 1
end
"""
     iszero(sys; atol = 0, atol1 = atol, atol2 = atol, rtol, fastrank = true)

Return `true` if the transfer function matrix of the descriptor system `sys` is zero.
For a descriptor system `sys = (A-λE,B,C,D)` with the transfer function matrix `G(λ)` it is checked
that the normal rank of `G(λ)` is zero, or equivalently (see [1]), that the normal rank of
the system matrix pencil

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

is equal to `n`, the order of the system `sys`.

If `fastrank = true`, the normal rank of `S(λ)` is evaluated by counting how many singular values of `S(γ)` have magnitudes
greater than `max(max(atol1,atol2), rtol*σ₁)`, where `σ₁` is the largest singular value of `S(γ)` and `γ` is a randomly generated value.
If `fastrank = false`, the rank is evaluated as `nr + ni + nf + nl`, where `nr` and `nl` are the sums of right and left Kronecker indices,
respectively, while `ni` and `nf` are the number of infinite and finite eigenvalues, respectively. The sums `nr+ni` and
`nf+nl`, are determined from an appropriate Kronecker-like form of the pencil `S(λ)`, exhibiting the spliting of the right and left structures.

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

[1] A. Varga,  On checking null rank conditions of rational matrices, 2018.
[arXiv:2006.06825](https://arxiv.org/pdf/1812.11396).
"""
function iszero(sys::DescriptorStateSpace{T}; atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
                rtol::Real =  ((max(size(sys.A)...))+1)*eps(real(float(one(real(T)))))*iszero(max(atol1,atol2)),
                fastrank::Bool = true) where T
    return (sprank(dssdata(sys)..., atol1 = atol1, atol2 = atol2, rtol = rtol, fastrank = fastrank) == sys.nx)
end
"""
    Gval = evalfr(sys, val; atol1 = 0, atol2 = 0, rtol, fast = true)

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

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of matrices `A`, `B`, `C`, `D`, the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`.

The computation of minimal realizations of individual input-output channels relies on pencil manipulation algorithms,
which employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.
"""
function evalfr(SYS::DescriptorStateSpace, val::Number; kwargs...)
    return lseval(dssdata(SYS)..., val; kwargs...)
end
"""
    Gval = evalfr(sys; fval = 0, atol = 0, atol1 = atol, atol2 = atol, rtol, fast = true)

Evaluate for a descriptor system `sys = (A-λE,B,C,D)` with the transfer function matrix `G(λ)`,
`Gval`, the value of the rational matrix `G(λ) = C*inv(λE-A)*B+D` for `λ = val`, where `val = im*fval`
for a continuous-time system or `val = exp(im*fval*sys.Ts)` for a discrete-time system.
The computed `Gval` has infinite entries if `val` is a pole (finite or infinite) of `G(λ)`.
If `val` is finite and `val*E-A` is singular or if `val = Inf` and `E` is singular,
then the entries of `Gval` are evaluated separately for minimal realizations of each input-output channel.

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

The computation of minimal realizations of individual input-output channels relies on pencil manipulation algorithms,
which employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.
"""
function evalfr(SYS::DescriptorStateSpace{T}; fval::Real = 0, atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
    rtol::Real =  ((max(size(SYS.A)...))+1)*eps(real(float(one(real(T)))))*iszero(max(atol1,atol2)), fast::Bool = true) where T
    val = SYS.Ts == 0 ? im*abs(fval) : exp(im*abs(fval*SYS.Ts))
    return lseval(dssdata(SYS)..., val; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
end
"""
    Gval = dcgain(sys; atol1, atol2, rtol, fast = true)

Evaluate for a descriptor system `sys = (A-λE,B,C,D)` with the transfer function matrix `G(λ)`,
`Gval`, the DC (or steady-state) gain. `Gval` is the value of the
rational matrix `G(λ)` for `λ = val`, where for a continuous-time system `val = 0` and
for a discrete-time system `val = 1`. The computed `Gval` has infinite entries if `val` is a pole of `G(λ)`.
In this case (i.e., `val*E-A` is singular), the entries of `Gval` are evaluated separately for minimal realizations
of each input-output channel.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of matrices `A`, `B`, `C`, `D`, the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`.

The computation of minimal realizations of individual input-output channels relies on pencil manipulation
algorithms, which employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.
"""
function dcgain(SYS::DescriptorStateSpace; kwargs...)
    SYS.Ts == 0 ? (return lseval(dssdata(SYS)..., 0; kwargs...) ) :
                  (return lseval(dssdata(SYS)..., 1; kwargs...) )
end
"""
     opnorm(sys[, p = Inf]; kwargs...)
     opnorm(sys, 2; kwargs...) -> sysnorm
     opnorm(sys, Inf; kwargs...) -> (sysnorm, fpeak)
     opnorm(sys; kwargs...) -> (sysnorm, fpeak)

Compute for a descriptor system `sys = (A-λE,B,C,D)` with the transfer function matrix `G(λ)`,
the `L2` or `L∞` system norm `sysnorm` induced by the vector `p`-norm, where valid values of `p` are `2` or `Inf`.
For the `L∞` norm, the frequency `fpeak` is also returned, where `G(λ)` achieves its peak gain.
See [`gh2norm`](@ref) and [`ghinfnorm`](@ref) for a description of the allowed keyword arguments.
"""
function opnorm(SYS::DescriptorStateSpace{T}, p::Real=Inf; fast::Bool = true, offset::Real = sqrt(eps(float(real(T)))),
    atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol, atolinf::Real = atol, rtolinf::Real = real(T)(0.001),
    rtol::Real = SYS.nx*eps(real(float(one(T))))*iszero(min(atol1,atol2))) where T
    if p == 2
        return gl2norm(SYS, fast = fast, offset = offset, atol1 = atol1, atol2 = atol2, atolinf = atolinf, rtol = rtol)
    elseif p == Inf
        return glinfnorm(SYS, rtolinf = rtolinf, fast = fast, offset = offset, atol1 = atol1, atol2 = atol2, rtol = rtol)
    else
        throw(ArgumentError("invalid p-norm p=$p. Valid: p = 2 and p = Inf"))
    end
end
"""
    sys = rss(n, p, m; disc = false, T = Float64, stable = false, nuc = 0, nuo = 0, randt = true)

Generate a randomized `n+nuc+nuo`-th order standard state-space system `sys = (A,B,C,D)` with `p` outputs and `m` inputs, with all matrices
randomly generated of type `T`.
The resulting `sys` is a continuous-time system if `disc = false` and a discrete-time system if `disc = true`.
If `stable = true`, the resulting system is stable, with `A` having all eigenvalues with negative real parts for a continuous-time system,
or with moduli less than one for a discrete-time system.
If `nuc+nuo > 0`, the system `sys` is non-minimal, with `A` having `nuc` uncontrollable and `nuo` unobservable eigenvalues.
If  `randt = true`, a randomly generated orthogonal or unitary similarity transformation is additionally applied.
If `randt = false`, the system matrices `A`, `B`, and `C` result in block stuctured forms exhibitting the
uncontrollable and unobservable eigenvalues of `A`:

    A = diag(A1, A2, A3),  B = [B1; 0; B3], C = [C1 C2 0]

with the diagonal blocks `A1`, `A2`, `A3` of orders `n`, `nuc`, and `nuo`, respectively.
"""
function rss(n::Int, p::Int, m::Int; disc::Bool = false, T::Type = Float64,
             nuc::Int = 0, nuo::Int = 0, randt = true, stable::Bool = false)
    nf = n; Af = randn(T,nf,nf); Bf = randn(T,nf,m); Cf = randn(T,p,nf);
    n > 0 && stable && (disc ? rmul!(Af,1/(rand()+opnorm(Af,1))) : Af -= I*(rand()+maximum(real(eigvals(Af)))) )
    Afuc = randn(T,nuc,nuc); Bfuc = zeros(T,nuc,m); Cfuc = randn(T,p,nuc);
    stable && nuc > 0 && (disc ? rmul!(Afuc,1/(rand()+opnorm(Afuc,1))) : Afuc -= I*(rand()+maximum(real(eigvals(Afuc)))) )
    Afuo = randn(T,nuo,nuo); Bfuo = randn(T,nuo,m); Cfuo = zeros(T,p,nuo);
    stable && nuo > 0 && (disc ? rmul!(Afuo,1/(rand()+opnorm(Afuo,1))) : Afuo -= I*(rand()+maximum(real(eigvals(Afuo)))) )
    nx = nf+nuc+nuo
    Q = randt ? qr!(rand(nx,nx)).Q : I
    return DescriptorStateSpace{T}(Q'*blockdiag(Af,Afuc,Afuo)*Q, I,
                                   Q'*[Bf;Bfuc;Bfuo], [ Cf Cfuc Cfuo ]*Q, rand(T,p,m),
                                   disc ? -one(real(T)) : zero(real(T)))
end
"""
    sys = rdss(n, p, m; id = [ ], disc = false, T = Float64, stable = false, nfuc = 0, iduc = [ ],
               nfuo = 0, iduo = [ ], randlt = true, randrt = true)

Generate a randomized descriptor state-space system `sys = (A-λE,B,C,D)` with `p` outputs and `m` inputs, with all matrices
randomly generated of type `T`.
The resulting `sys` is a continuous-time system if `disc = false` and a discrete-time system if `disc = true`.

If the vector `id` is nonempty, then `id[i]` specifies the order of the `i`-th infinite elementary divisor of the
resulting pencil `A-λE`, which thus has `n` finite eigenvalues and `ni = sum(id)` infinite eigenvalues which are
controllable and observable.
If `nfuc+nfuo > 0`, the system `sys` is non-minimal, with `A` having `nfuc` uncontrollable and `nfuo` unobservable finite eigenvalues.
If the vector `iduc` is a nonempty, then `iduc[i]` specifies the order of the `i`-th infinite elementary divisor
with uncontrollable infinite eigenvalues of the resulting pencil `A-λE`, which thus has `niuc = sum(iduc)` uncontrollable infinite eigenvalues.
If the vector `iduo` is a nonempty, then `iduo[i]` specifies the order of the `i`-th infinite elementary divisor with
unobservable infinite eigenvalues of the resulting pencil `A-λE`, which thus has `niuo = sum(iduo)` unobservable infinite eigenvalues.
If `niuc+niuo > 0`, the system `sys` is non-minimal, with `A` having `niuc` uncontrollable and `niuo` unobservable infinite eigenvalues.

It follows, that the resulting pencil `A-λE` has  `n+nfuc+nfuo` finite eigenvalues and `ni+niuc+niuo` infinite eigenvalues.
If `stable = true`, the proper part of the system `sys` is stable, with `A` having all finite eigenvalues with negative real parts
for a continuous-time system, or with moduli less than one for a discrete-time system.

If  `randlt = true`, a randomly generated orthogonal or unitary transformation is additionally applied to `A`, `E`, and `B` from the left.
If  `randrt = true`, a randomly generated orthogonal or unitary transformation is additionally applied to `A`, `E`,  and `C` from the right.
If `randlt = false` and `randrt = false`, the system matrices `A`, `E`, `B`, and `C` result in block stuctured forms exhibitting the
uncontrollable and unobservable finite and infinite eigenvalues of `A-λE`:

    A-λE = diag(A1-λE1, A2-λE2, A3-λE3, A4-λE4, A5-λE5, A6-λE6),
    B = [B1; B2; 0; 0; B5; B6 ],
    C = [C1 C2 C3 C4 0 0]

with the diagonal blocks `A1`, `A2`, `A3`, `A4`, `A5`, `A6` of orders `n`, `ni`, `nfuc`, `niuc`, `nfuo` and `niuo`, respectively.
"""
function rdss(n::Int, p::Int, m::Int; disc::Bool = false, stable::Bool = false, T::Type = Float64, id::Vector{Int} = Int[],
                                      nfuc::Int = 0, nfuo::Int = 0, iduc::Vector{Int} = Int[], iduo::Vector{Int} = Int[],
                                      randlt = true, randrt = true)
    nf = n; Af = randn(T,nf,nf); Ef = randn(T,nf,nf); Bf = randn(T,nf,m); Cf = randn(T,p,nf);
    n > 0 && stable && (disc ? Af = Ef*rmul!(Af,1/(rand()+opnorm(Af,1))) : Af = Ef*(Af-I*(rand()+maximum(real(eigvals(Af))))) )
    ni = sum(id); Ai = Matrix{T}(I,ni,ni); Ei = jordanblockdiag(zero(T),id); Bi = randn(T,ni,m); Ci = randn(T,p,ni);
    Afuc = randn(T,nfuc,nfuc); Efuc = randn(T,nfuc,nfuc); Bfuc = zeros(T,nfuc,m); Cfuc = randn(T,p,nfuc);
    stable && nfuc > 0 && (disc ? Afuc = Efuc*rmul!(Afuc,1/(rand()+opnorm(Afuc,1))) : Afuc = Efuc*(Afuc-I*(rand()+maximum(real(eigvals(Afuc))))) )
    niuc = sum(iduc); Aiuc = Matrix{T}(I,niuc,niuc); Eiuc = jordanblockdiag(zero(T),iduc); Biuc = zeros(T,niuc,m); Ciuc = randn(T,p,niuc);
    Afuo = randn(T,nfuo,nfuo); Efuo = randn(T,nfuo,nfuo); Bfuo = randn(T,nfuo,m); Cfuo = zeros(T,p,nfuo);
    stable && nfuo > 0 && (disc ? Afuo = Efuo*rmul!(Afuo,1/(rand()+opnorm(Afuo,1))) : Afuo = Efuo*(Afuo-I*(rand()+maximum(real(eigvals(Afuo))))) )
    niuo = sum(iduo); Aiuo = Matrix{T}(I,niuo,niuo); Eiuo = jordanblockdiag(zero(T),iduo); Biuo = randn(T,niuo,m); Ciuo = zeros(T,p,niuo);
    nx = nf+ni+nfuc+nfuo+niuc+niuo
    Q = randlt ? qr!(rand(nx,nx)).Q : I
    Z = randrt ? qr!(rand(nx,nx)).Q : I
    return DescriptorStateSpace{T}(Q*blockdiag(Af,Ai,Afuc,Aiuc,Afuo,Aiuo)*Z, Q*blockdiag(Ef,Ei,Efuc,Eiuc,Efuo,Eiuo)*Z,
                                   Q*[Bf;Bi;Bfuc;Biuc;Bfuo;Biuo], [ Cf Ci Cfuc Ciuc Cfuo Ciuo]*Z, rand(T,p,m),
                                   disc ? -one(real(T)) : zero(real(T)))
end
"""
    sysr = dsxvarsel(sys,ind)

Construct for the descriptor system `sys = (A-λE,B,C,D)` of order `n` the descriptor system
`sysr = (A[ind,ind]-λE[ind,ind],B[ind,:],C[:,ind],D)` of order `nr = length(ind)`,
by selecting the state variables of `sys` with indices specified by `ind`.
If `ind` is a permutation vector of length `n`, then `sysr` has the same transfer function matrix as `sys`
and permuted state variables.
"""
function dsxvarsel(SYS::DescriptorStateSpace{T},ind::Union{UnitRange{Int64},StepRange{Int64, Int64}, Array{Int64,1}}) where T
    isempty(ind) &&  (return DescriptorStateSpace{T}(zeros(T,0,0),I,zeros(T,0,SYS.nu),zeros(T,SYS.ny,0),SYS.D,SYS.Ts))
    (minimum(ind) < 1 || maximum(ind) > SYS.nx) && error("BoundsError: selected indices $ind out of range $(1:SYS.nx)")
    return DescriptorStateSpace{T}(SYS.A[ind,ind], SYS.E == I ? I : SYS.E[ind,ind], SYS.B[ind,:], SYS.C[:,ind], SYS.D, SYS.Ts)
end

"""
     rtfbilin(type = "c2d"; Ts = T, Tis = Ti, a = val1, b = val2, c = val3, d = val4) -> (g, ginv)

Build the rational transfer functions of several commonly used bilinear transformations and their inverses.
The resulting `g` describes the rational transfer function `g(δ)` in the bilinear transformation `λ = g(δ)` and
`ginv` describes its inverse transformation `ginv(λ)` in the bilinear transformation `δ = ginv(λ)`.
In accordance with the values of `type` and the keyword argument values `Ts`, `Tis`, `a`, `b`, `c`, and `d`,
the resulting `g` and `ginv` contain first order rational transfer functions of the form `g(δ) = (a*δ+b)/(c*δ+d)`
and `ginv(λ) = (d*λ-b)/(-c*λ+a)`, respectively, which satisfy `g(ginv(λ)) = λ` and `ginv(g(δ)) = δ`.

Depending on the value of `type`, the following types of transformations can be generated
in conjunction with parameters specified in `Ts`, `Tis`, `a`, `b`, `c`, and `d`:

    "Cayley" - Cayley transformation: `s = g(z) = (z-1)/(z+1)` and `z = ginv(s) = (s+1)/(-s+1)`;
               the sampling time `g.Ts` is set to the value `T ≠ 0` (default `T = -1`), while `ginv.Ts = 0`;
               g(z) and ginv(s) are also known as the continuous-to-discrete and discrete-to-continuous transformations, respectively;

    "c2d"    - is alias to "Cayley"

    "Tustin" - Tustin transformation (also known as trapezoidal integration): `s = g(z) = (2*z-2)/(T*z+T)` and `z = ginv(s) = (T*s+2)/(-T*s+2)`;
               the sampling time `g.Ts` is set to the value `T ≠ 0` (default `T = -1`), while `ginv.Ts = 0`;
               a nonzero prewarping frequency `freq` can be specified using the keyword parameter `prewarp_freq = freq`;

    "Euler"  - Euler integration (or forward Euler integration): `s = g(z) = s = (z-1)/T` and `z = ginv(s) = T*s+1`;
               the sampling time `g.Ts` is set to the value `T ≠ 0` (default `T = -1`), while `ginv.Ts = 0`;

    "BEuler" - Backward Euler integration (or backward Euler integration): `s = g(z) = (z-1)/(T*z)` and `z = ginv(s) = 1/(-T*s+1)`;
               the sampling time `g.Ts` is set to the value `T ≠ 0` (default `T = -1`), while `ginv.Ts = 0`;

    "Moebius" - general (Moebius ) bilinear transformation: `λ = g(δ) = (a*δ+b)/(c*δ+d)` and `δ = ginv(λ) = (d*λ-b)/(-c*λ+a)`;
                the sampling times `g.Ts` and `ginv.Ts` are  set to the values `T` and `Ti`, respectively;
                the default values of the parameters `a`, `b`, `c`, and `d` are `a = 1`, `b = 0`, `c = 0`, and `d = 1`.
                Some useful particular Moebius transformations correspond to the following choices of parameters:
                - _translation_: `a = 1`, `b ≠ 0`, `c = 0`, `d = 1` (`λ = δ+b`)
                - _scaling_: `a ≠ 0`, `b = 0`, `c = 0`, `d = 1` (`λ = a*δ`)
                - _rotation_: `|a| = 1`, `b = 0`, `c = 0`, `d = 1` (`λ = a*δ`)
                - _inversion_: `a = 0`, `b = 1`, `c = 1`, `d = 0`  (`λ = 1/δ`)

    "lft"     - alias to "Moebius" (linear fractional transformation).
"""
function rtfbilin(type::String = "c2d"; Ts::Union{Real,Missing} = missing, Tsi::Union{Real,Missing} = missing,
                  prewarp_freq::Real = 0,  a::Number = 1, b::Number = 0, c::Number = 0, d::Number = 1 )
    s = rtf('s');
    type = lowercase(type)
    if type == "c2d" || type == "cayley"
        ismissing(Ts) && (Ts = -1)
        z = rtf('z',Ts = Ts)
        g = (z-1)/(z+1); ginv = (s+1)/(-s+1)
    elseif type == "d2c"
        ismissing(Tsi) && (Tsi = -1)
        z = rtf('z',Ts = Tsi);
        g = (s+1)/(-s+1); ginv = (z-1)/(z+1);
    elseif type ==  "tustin"
        ismissing(Ts) && (Ts = 2)
        Ts > 0 || error("sampling time must be positive")
        z = rtf('z',Ts = Ts);
        prewarp_freq == 0 ? t = Ts : t = 2*tan(prewarp_freq*Ts/2)/prewarp_freq
        g = (2*z-2)/(t*z+t); ginv = (t*s+2)/(-t*s+2);
    elseif type ==  "euler"
        ismissing(Ts) && (Ts = 1)
        Ts > 0 || error("sampling time must be positive")
        z = rtf('z',Ts = Ts);
        g = (z-1)/Ts; ginv = Ts*s+1;
    elseif type ==  "beuler"
        ismissing(Ts) && (Ts = 1)
        Ts > 0 || error("sampling time must be positive")
        z = rtf('z',Ts = Ts);
        g = (z-1)/(Ts*z); ginv = 1/(-Ts*s+1);
    elseif type ==  "lft" || type ==  "moebius"
        c == 0 && d == 0 &&  error("the denominator of g must be nonzero")
        a*d ≈ b*c &&  error("the McMillan degree of g must be one")
        ismissing(Ts) && (Ts = 0)
        if Ts == 0
            g = (a*s+b)/(c*s+d)
        else
            z = rtf('z',Ts = Ts)
            g = (a*z+b)/(c*z+d)
        end
        ismissing(Tsi) && (Tsi = 0)
        if Tsi == 0
            ginv = (d*s-b)/(-c*s+a)
       else
            z = rtf('z',Ts = Tsi)
            ginv = (d*z-b)/(-c*z+a)
        end
    else
        error("Unknown type")
    end
    return g, ginv
    # end RTFBILIN
end
"""
     res = dssubset(sys, subsys, rows, cols; minimal = true, fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol = nϵ)

Assign for a given descriptor system `sys = (A-λE,B,C,D)` its subsystem `sys[row,cols]` to `subsys` and
return the modified system in `res`. `rows` and `cols` are indices, vectors of indices, index ranges, `:` or
any combinations of them.

If `minimal = true` (default), an irreducible realization of the resulting system `res` is computed, otherwise a
possibly non-minimal realization is returned if `minimal = false`.

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

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

"""
function dssubset(sys::DescriptorStateSpace{T1}, subsys::DescriptorStateSpace{T2}, inds...;
                        atol::Real = zero(real(T1)), atol1::Real = atol, atol2::Real = atol,
                        rtol::Real =  sys.nx*eps(real(float(one(real(T1)))))*iszero(max(atol1,atol2)),
                        fast::Bool = true, minimal::Bool = true  ) where {T1, T2}
    size(inds, 1) != 2 &&
        error("Must specify 2 indices to index descriptor state-space models")
    rows, cols = index2range(inds...)
    p, m = size(sys)
    rows == Colon() && (rows = 1:p)
    cols == Colon() && (cols = 1:m)
    pred, mred = size(subsys)
    length(rows) == pred  || error("number of row indices must be equal to the number of outputs of subsys")
    length(cols) == mred  || error("number of column indices must be equal to the number of inputs of subsys")
    pred > 0 && maximum(rows) > p && error("output indices must not exceed $p")
    mred > 0 && maximum(cols) > m && error("input indices must not exceed $m")
    irows = setdiff(Vector(1:p),Vector(rows))
    jcols = setdiff(Vector(1:m),Vector(cols))
    if minimal
        res = gir([sys[irows,[jcols;cols]]; sys[rows,jcols] subsys]; atol, atol1, atol2, rtol, fast )
    else
        res = [sys[irows,[jcols;cols]]; sys[rows,jcols] subsys]
    end
    ip = sortperm([irows;rows])
    jp = sortperm([jcols;cols])
    return DescriptorStateSpace{eltype(res)}(res.A, res.E, res.B[:, jp], res.C[ip, :], res.D[ip, jp], sys.Ts)
end
dssubset(sys::DST, subsys::AbstractNumOrArray, inds...; kwargs...) where {T, DST <: DescriptorStateSpace{T}} = dssubset(sys, dss(T.(subsys); Ts = sys.Ts), inds...; kwargs...)
"""
     res = dszeros(sys, rows, cols; minimal = true, fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol = nϵ)

Set for a given descriptor system `sys = (A-λE,B,C,D)` its subsystem `sys[row,cols]` to zero and
return the modified system in `res`. `rows` and `cols` are indices, vectors of indices, index ranges, `:` or
any combinations of them.

If `minimal = true` (default), an irreducible realization of the resulting system `res` is computed, otherwise a
possibly non-minimal realization is returned if `minimal = false`.

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

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

"""
function dszeros(sys::DescriptorStateSpace{T}, inds...;
                        atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
                        rtol::Real =  sys.nx*eps(real(float(one(real(T)))))*iszero(max(atol1,atol2)),
                        fast::Bool = true, minimal::Bool = true  ) where T
    size(inds, 1) != 2 &&
        error("Must specify 2 indices to index descriptor state-space models")
    rows, cols = index2range(inds...)
    p, m = size(sys)
    rows == Colon() && (rows = 1:p)
    cols == Colon() && (cols = 1:m)
    irows = setdiff(Vector(1:p),Vector(rows))
    jcols = setdiff(Vector(1:m),Vector(cols))
    subsys = dss(zeros(T,length(rows),length(cols)); Ts = sys.Ts)
    if minimal
        res = gir([sys[irows,[jcols;cols]]; sys[rows,jcols] subsys]; atol, atol1, atol2, rtol, fast )
    else
        res = [sys[irows,[jcols;cols]]; sys[rows,jcols] subsys]
    end
    ip = sortperm([irows;rows])
    jp = sortperm([jcols;cols])
    return DescriptorStateSpace{eltype(res)}(res.A, res.E, res.B[:, jp], res.C[ip, :], res.D[ip, jp], sys.Ts)
end
"""
     res = dssubsel(sys, S; minimal = true, fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol = nϵ)

Select for a given descriptor system `sys = (A-λE,B,C,D)` its subsystem `res` corresponding to the zero-one pattern specified by
the binary structure matrix `S`. If `G(λ)` is the transfer function matrix of `sys`, then the transfer function matrix of `res`
is the element-wise product `S .* G(λ)`.

If `minimal = true` (default), an irreducible realization of the resulting system `res` is computed, otherwise a
possibly non-minimal realization is returned if `minimal = false`.

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

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of matrices `A`, `B`, `C`, `D`, the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`.
The default relative tolerance is `nϵ`, where `ϵ` is the working _machine epsilon_
and `n` is the order of the system `sys`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.
"""
function dssubsel(sys::DescriptorStateSpace{T}, S::Union{BitMatrix,Matrix{Bool}};
                        atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
                        rtol::Real =  sys.nx*eps(real(float(one(real(T)))))*iszero(max(atol1,atol2)),
                        fast::Bool = true, minimal::Bool = true ) where T
    p, m = size(sys)
    (p, m ) == size(S) || error("missmatch between the sizes of sys and S")
    syst = dss(zeros(T,p,m); Ts = sys.Ts)
    cols = collect(1:m)
    for i = 1:p
        colsi = cols[S[i,:]]
        syst = dssubset(syst,sys[i,colsi],i,colsi; atol1, atol2, rtol, fast, minimal = false)
    end
    return minimal ? gir(syst; atol, atol1, atol2, rtol, fast) : syst
end
"""
     sysdiag = dsdiag(sys, k)

Build for a given descriptor system `sys` and non-negative integer `k` a descriptor system `sysdiag`
obtained as `sysdiag = diag( sys, ..., sys )` (i.e., `k`-times repeated application of `append`).
"""
function dsdiag(sys::DescriptorStateSpace{T}, k::Int) where T
    p, m = size(sys)
    k >=0 || error("k must be a non-negative integer")
    n = sys.nx
    N = k*n
    M = k*m
    P = k*p
    a, e, b, c, d = dssdata(sys)
    desc = !(e == I)
    A = zeros(T, N, N)
    E = desc ? zeros(T, N, N) : I
    B = zeros(T, N, M)
    C = zeros(T, P, N)
    D = zeros(T, P, M)
    ia = 1
    jb = 1
    ic = 1
    for i = 1:k
        ia1 = ia:ia+n-1
        jb1 = jb:jb+m-1
        ic1 = ic:ic+p-1
        A[ia1,ia1] = a
        desc && (E[ia1,ia1] = e)
        B[ia1,jb1] = b
        C[ic1,ia1] = c
        D[ic1,jb1] = d
        ia += n
        jb += m
        ic += p
    end
    return DescriptorStateSpace{T}(A, E, B, C, D, sys.Ts)
end
"""
     syscl = feedback(sys, K, inp1, out1; negative = true)

Build for a given descriptor system `sys` with input vector `u` and output vector `y` and
a static output feedback gain `K` the closed-loop descriptor system `syscl` is constructed
corresponding to the static output feedback `u[inp] = -K*y[out] + v`, where `inp` and `out` are
are indices, vectors of indices, index ranges, `:` or any combinations of them. Only distinct indices
can be specified. If `negative = false`, a positive feedback `u[inp] = K*y[out] + v` is used.
"""
function feedback(sys::DescriptorStateSpace{T1}, K::Matrix{T2}, inds...; negative::Bool = true) where {T1, T2}
    size(inds, 1) != 2 &&
        error("Must specify 2 indices to index descriptor state-space models")
    u1, y1 = index2range(inds...)
    p, m = size(sys)
    u1 == Colon() && (u1 = 1:m)
    y1 == Colon() && (y1 = 1:p)
    maximum(u1) > m && error("input indices must not exceed $m")
    maximum(y1) > m && error("output indices must not exceed $p")
    allunique(u1) || error("all input indices must be distinct")
    allunique(y1) || error("all output indices must be distinct")
    m1, p1 = size(K)
    length(u1) == m1  || error("number of row indices must be equal to the number of rows of K")
    length(y1) == p1  || error("number of column indices must be equal to the number of columns of K")
    y2 = setdiff(Vector(1:p),Vector(y1))
    u2 = setdiff(Vector(1:m),Vector(u1))

    A, E, B, C, D = dssdata(sys)
    B1 = view(B,:,u1)
    B2 = view(B,:,u2)
    C1 = view(C,y1,:)
    C2 = view(C,y2,:)
    D11 = view(D,y1,u1)
    D12 = view(D,y1,u2)
    D21 = view(D,y2,u1)
    D22 = view(D,y2,u2)
    nullD11 = iszero(D11)
    Li = nullD11 ? I : (negative ? inv(I+K*D11) : inv(I-K*D11))
    KT = nullD11 ? (negative ? -K : K) : (negative ? -Li*K : Li*K)
    ip = sortperm([y1;y2])
    jp = sortperm([u1;u2])

    return dss(A+B1*KT*C1, E, [B1*Li (B2+B1*KT*D12)][:,jp], [C1 + D11*KT*C1 ; C2 + D21*KT*C1][ip,:],
                        [ D11*Li D12+D11*KT*D12; D21*Li D22+D21*KT*D12 ][ip,jp]; Ts = sys.Ts)

end
feedback(sys::DescriptorStateSpace{T1}, K::Matrix{T2}; negative::Bool = true) where {T1, T2} = feedback(sys,K,1:size(K,1),1:size(K,2); negative)


text:
# started conversion from Python to julia
using Revise
using MccDaqHats
using Infiltrator
includet(joinpath(@__DIR__, "scan_utils.jl"))

# Constants
CURSOR_BACK_2 = "\x1b[2D"
ERASE_TO_END_OF_LINE = "\x1b[0K"

"""
MCC 172 Functions Demonstrated:
    mcc172_trigger_config
    mcc172_iepe_config_write
    mcc172_a_in_clock_config_write
    mcc172_a_in_clock_config_read
    mcc172_a_in_scan_start
    mcc172_a_in_scan_read
    mcc172_a_in_scan_stop

Purpose:
    Perform a triggered finite acquisition on 1 or more channels.

Description:
    Waits for an external trigger to occur and then acquires blocks
    of analog input data for a user-specified group of channels.  The
    RMS voltage for each channel is displayed for each block
    of data received from the device.  The acquisition is stopped when
    the specified number of samples is acquired for each channel.
"""
function finite_scan_with_trigger()
    # Store the channels in a list and convert the list to a channel mask that
    # can be passed as a parameter to the MCC 172 functions.
    channels = [0, 1]
    channel_mask = chan_list_to_mask(channels)
    num_channels = length(channels)

    samples_per_channel = UInt32(10000)
    scan_rate = 10240.0
    options = Set([:OPTS_EXTTRIGGER])
    trigger_mode = :TRIG_RISING_EDGE

    try
        # Select an MCC 172 HAT device to use.
        hats = select_hat_devices(:MCC_172, 1)
        address = hats[1].address

        print("\nSelected MCC 172 HAT device at address $address")

        # Turn on IEPE supply?
        iepe_enable = get_iepe()

        mcc172_open(address)
        for channel in channels
            mcc172_iepe_config_write(address, channel, iepe_enable)
        end

        # Configure the clock and wait for sync to complete.
        mcc172_a_in_clock_config_write(address, :SOURCE_LOCAL, scan_rate)

        synced = false
        actual_scan_rate = 0.0  # initialize to Float64
        while !synced
            (_source_type, actual_scan_rate, synced) = mcc172_a_in_clock_config_read(address)
            if !synced
                sleep(0.005)
            end
        end

        println("\nMCC 172 continuous scan example")
        println("    Functions demonstrated:")
        println("         mcc172_trigger_mode")
        println("         mcc172_a_in_clock_config_write")
        println("         mcc172_a_in_clock_config_read")
        println("         mcc172_a_in_scan_start")
        println("         mcc172_a_in_scan_read")
        println("         mcc172_a_in_scan_stop")
        println("         mcc172_a_in_scan_cleanup")
        println("    IEPE power: $(iepe_enable ? "on" : "off")")
        println("    Channels: $(join(channels, ", "))")
        println("    Requested scan rate: $scan_rate")
        println("    Actual scan rate: $(round(actual_scan_rate, digits=3))")
        println("    Samples per channel $samples_per_channel")
        println("    Options: $(options)")
        println("    Trigger Mode: $trigger_mode")

        try
            println("\nPress 'Enter' to continue")
            readline()
        catch
            error("^C to end program")
        end

        mcc172_trigger_config(address, :SOURCE_LOCAL, trigger_mode)

        # Configure and start the scan.
        mcc172_a_in_scan_start(address, channel_mask, samples_per_channel, options)

        try
            # wait for the external trigger to occur
            println("\nWaiting for trigger ... hit Ctrl-C to cancel the trigger")
            wait_for_trigger(address)

            println("\nStarting scan ... Press Ctrl-C to stop\n")

            # Display the header row for the data table.
            print("Samples Read    Scan Count")
            for chan in channels
                print("       Channel  $chan")
            end
            println()

            read_and_display_data(address, samples_per_channel, num_channels)

        catch e # KeyboardInterrupt
            if isa(e, InterruptException)
                # Clear the "^C" from the display.
                println("$CURSOR_BACK_2, $ERASE_TO_END_OF_LINE, Aborted\n")
            else
                println("\n $e")
            end
            mcc172_a_in_scan_stop(address)
        end

        mcc172_a_in_scan_cleanup(address)

    catch err  # (HatError, ValueError) as err
        println("\n $err")
    end
end

"""
Reads data from the specified channels on the specified DAQ HAT devices
and updates the data on the terminal display.  The reads are executed in a
loop that continues until either the scan completes or an overrun error
is detected.

Args:
hat (mcc172): The mcc172 HAT device object.
samples_per_channel: The number of samples to read for each channel.
num_channels (int): The number of channels to display.

Returns:
Nothing
"""
function read_and_display_data(address::Integer, samples_per_channel::Integer, num_channels::Integer)
    total_samples_read = 0
    read_request_size = 1000
    timeout = 5.0

    # Since the read_request_size is set to a specific value, a_in_scan_read()
    # will block until that many samples are available or the timeout is
    # exceeded.

    # Continuously update the display value until Ctrl-C is pressed
    # or the number of samples requested has been read.
    while total_samples_read < samples_per_channel
        resultcode, statuscode, result, samples_read =
        mcc172_a_in_scan_read(address, Int32(read_request_size), num_channels, timeout)
        # @show(length(result), samples_read)

        # Check for an overrun error
        status = mcc172_status_decode(statuscode)
        if status.hardwareoverrun
            println("\n\nHardware overrun\n")
            break
        elseif status.bufferoverrun
            println("\n\nBuffer overrun\n")
            break
        end

        samples_read_per_channel = length(result) ÷ num_channels
        total_samples_read += samples_read_per_channel

        print("\r$(samples_read_per_channel)         $(total_samples_read)")

        # Display the RMS voltage for each channel.
        if samples_read_per_channel > 0
            for i in 1:num_channels
                value = calc_rms(result, i, num_channels, samples_read_per_channel)
                print("      $(round(value, digits=5)) Vrms   ")
            end
            println()
        end
    end
    print("\n")
end

ext: jl

max_stars_repo_name: Spectrum-Tec/MccDaqHats

max_issues_repo_name: Spectrum-Tec/MccDaqHats

max_forks_repo_name: Spectrum-Tec/MccDaqHats


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCell
include("../../test/useful_routines.jl")

function step_interface(levelset, levelsetcoeffs, mesh, levelsetspeed, dt)
    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    paddedmesh = CutCell.BoundaryPaddedMesh(cutmesh, 1)
    refseedpoints, spatialseedpoints, seedcellids =
        CutCell.seed_zero_levelset(2, levelset, levelsetcoeffs, cutmesh)
    paddedlevelset = CutCell.BoundaryPaddedLevelset(
        paddedmesh,
        refseedpoints,
        spatialseedpoints,
        seedcellids,
        levelset,
        levelsetcoeffs,
        cutmesh,
        1e-10,
    )
    return CutCell.step_first_order_levelset(paddedlevelset,levelsetspeed,dt)
end

function grid_size(mesh)
    w = CutCell.widths(mesh)
    nn = CutCell.nodes_per_mesh_side(mesh)

    return w ./ (nn .- 1)
end

function time_step_size(levelsetspeed,mesh;CFL=0.9)
    dx = minimum(grid_size(mesh))
    s = maximum(abs.(levelsetspeed))
    return CFL*dx/s
end

function run_time_steps(levelset,initialcondition,mesh,levelsetspeed,dt,nsteps)
    levelsetcoeffs = [copy(initialcondition) for i = 1:nsteps+1]
    for i = 1:nsteps
        levelsetcoeffs[i+1] = step_interface(levelset,levelsetcoeffs[i],mesh,levelsetspeed,dt)
    end
    return levelsetcoeffs
end

x0 = [0.0, 0.0]
L, W = 1.0, 1.0
nelmts = 10
numghostlayers = 1
polyorder = 2

xI = [0.05, 0.0]
normal = [1.0, 0.0]
tol = 1e-8
speed = 1.0
stoptime = 0.8

basis = TensorProductBasis(2, polyorder)
mesh = CutCell.Mesh(x0, [L, W], [nelmts, nelmts], basis)

levelset = InterpolatingPolynomial(1, basis)
initiallevelset =
    CutCell.levelset_coefficients(x -> plane_distance_function(x, normal, xI), mesh)
levelsetspeed = speed*ones(length(initiallevelset))

dt = time_step_size(levelsetspeed,mesh)
nsteps = ceil(Int,stoptime/dt)
levelsetcoeffs = run_time_steps(levelset,initiallevelset,mesh,levelsetspeed,stoptime)

actualstoptime = dt*nsteps
xT = xI + [actualstoptime*speed,0.0]

quad = tensor_product_quadrature(2,3)
err = uniform_mesh_L2_error(levelsetcoeffs[end]',x->plane_distance_function(x,normal,xT),basis,quad,mesh)
@test allapprox(err,[0.0],1e2eps())

function grid_range(mesh)
    x0 = CutCell.reference_corner(mesh)
    w = CutCell.widths(mesh)
    nfmside = CutCell.nodes_per_mesh_side(mesh)

    x = range(x0[1],stop=x0[1]+w[1],length=nfmside[1])
    y = range(x0[2],stop=x0[2]+w[2],length=nfmside[2])

    return x,y
end

using Plots
x,y = grid_range(mesh)
Z = reshape(levelsetcoeffs[1],length(y),:)
Plots.contour(x,y,Z,levels=-1:0.25:1,color="black",linewidth=2,legend=false)

anim = @animate for i = 1:length(levelsetcoeffs)
    Z = reshape(levelsetcoeffs[i],length(y),:)
    contour(x, y, Z, levels = -1:0.25:1,color="black",linewidth=2,legend=false)
end

gif(anim,fps=5)


text:
# Use this simple script to record an interaction on an opened scene

framerate = 10
total_time = 30 # in seconds
framen = framerate*total_time

record(scene, "name.mp4"; framerate = framerate) do io
    for i = 1:framen
        sleep(1/framerate)
        recordframe!(io)
    end
end


text:
transform_ldiv(A::AbstractQuasiArray{T}, f::AbstractQuasiArray{V}, ::Tuple{<:Any,InfiniteCardinal{0}}) where {T,V}  =
    adaptivetransform_ldiv(convert(AbstractQuasiArray{promote_type(T,V)}, A), f)

pad(c::AbstractVector{T}, ax::Union{OneTo,OneToInf}) where T = [c; Zeros(length(ax)-length(c))]
pad(c, ax...) = PaddedArray(c, ax)

padchop!(cfs, tol, ax...) = pad(chop!(cfs, tol), ax...)
padchop(cfs, tol, ax...) = pad(chop(cfs, tol), ax...)

# ax will impose block structure for us
padchop!(cfs::PseudoBlockVector, tol, ax...) = padchop!(cfs.blocks, tol, ax...)

increasingtruncations(::OneToInf) = oneto.(2 .^ (4:∞))
increasingtruncations(::BlockedUnitRange) = broadcast(n -> Block.(oneto(n)), (2 .^ (4:∞)))
function adaptivetransform_ldiv(A::AbstractQuasiArray{U}, f::AbstractQuasiVector{V}) where {U,V}
    T = promote_type(eltype(U),eltype(V))

    r = checkpoints(A)
    fr = f[r]
    maxabsfr = norm(fr,Inf)

    tol = 20eps(real(T))
    ax = axes(A,2)

    for jr in increasingtruncations(ax)
        An = A[:,jr]
        cfs = An \ f
        maxabsc = maximum(abs, cfs)
        if maxabsc == 0 && maxabsfr == 0
            return pad(similar(cfs,0), ax)
        end

        m = length(cfs)
        if maximum(abs,@views(cfs[max(1,m-2):m])) < 10tol*maxabsc
            n = length(cfs)
            c = padchop!(cfs, tol, ax)
            un = A * c
            if all(norm.(un[r] - fr, 1) .< tol * n * maxabsfr*1000)
                return c
            end
        end
    end
    error("Have not converged")
end

function adaptivetransform_ldiv(A::AbstractQuasiArray{U}, f::AbstractQuasiMatrix{V}) where {U,V}
    T = promote_type(eltype(U),eltype(V))

    bx = axes(f,2)
    r = checkpoints(A)
    fr = f[r,:]
    maxabsfr = norm(fr,Inf)

    tol = 20eps(real(T))
    ax = axes(A,2)

    for jr = increasingtruncations(ax)
        An = A[:,jr]
        cfs = An \ f
        maxabsc = maximum(abs, cfs)
        if maxabsc == 0 && maxabsfr == 0
            return pad(similar(cfs,0,size(cfs,2)), ax, bx)
        end

        if maximum(abs,@views(cfs[end-2:end,:])) < 10tol*maxabsc
            n = size(cfs,1)
            c = padchop(cfs, tol, ax, bx)
            un = A * c # expansion
        # we allow for transformed coefficients being a different size
        ##TODO: how to do scaling for unnormalized bases like Jacobi?
            if all(norm.(un[r,:] - fr, 1) .< tol * n * maxabsfr*1000)
                return c
            end
        end
    end
    error("Have not converged")
end


text:
const _rtol = 1e-10
const _atol = 1e-10

_to_psd(A::Matrix{<:Real}) = A * A' + I
_to_psd(a::Vector{<:Real}) = exp.(a) .+ 1
_to_psd(σ::Real) = exp(σ) + 1

Base.length(::Nothing) = 0

function print_adjoints(adjoint_ad, adjoint_fd, rtol, atol)
    @show typeof(adjoint_ad), typeof(adjoint_fd)
    adjoint_ad, adjoint_fd = to_vec(adjoint_ad)[1], to_vec(adjoint_fd)[1]
    println("atol is $atol, rtol is $rtol")
    println("ad, fd, abs, rel")
    abs_err = abs.(adjoint_ad .- adjoint_fd)
    rel_err = abs_err ./ adjoint_ad
    display([adjoint_ad adjoint_fd abs_err rel_err])
    return println()
end

# # AbstractArrays.
# function to_vec(x::ColVecs{<:Real})
#     x_vec, back = to_vec(x.X)
#     return x_vec, x_vec -> ColVecs(back(x_vec))
# end

Base.zero(d::Dict) = Dict([(key, zero(val)) for (key, val) in d])
Base.zero(x::Array) = zero.(x)

# My version of isapprox
function fd_isapprox(x_ad::Nothing, x_fd, rtol, atol)
    return fd_isapprox(x_fd, zero(x_fd), rtol, atol)
end
function fd_isapprox(x_ad::AbstractArray, x_fd::AbstractArray, rtol, atol)
    return all(fd_isapprox.(x_ad, x_fd, rtol, atol))
end
function fd_isapprox(x_ad::Real, x_fd::Real, rtol, atol)
    return isapprox(x_ad, x_fd; rtol=rtol, atol=atol)
end
function fd_isapprox(x_ad::NamedTuple, x_fd, rtol, atol)
    f = (x_ad, x_fd) -> fd_isapprox(x_ad, x_fd, rtol, atol)
    return all([f(getfield(x_ad, key), getfield(x_fd, key)) for key in keys(x_ad)])
end
function fd_isapprox(x_ad::Tuple, x_fd::Tuple, rtol, atol)
    return all(map((x, x′) -> fd_isapprox(x, x′, rtol, atol), x_ad, x_fd))
end
function fd_isapprox(x_ad::Dict, x_fd::Dict, rtol, atol)
    return all([
        fd_isapprox(get(() -> nothing, x_ad, key), x_fd[key], rtol, atol) for
        key in keys(x_fd)
    ])
end

function adjoint_test(
    f, ȳ, x...; rtol=_rtol, atol=_atol, fdm=central_fdm(5, 1), print_results=false
)
    # Compute forwards-pass and j′vp.
    y, back = Zygote.pullback(f, x...)
    adj_ad = back(ȳ)
    adj_fd = j′vp(fdm, f, ȳ, x...)

    # Check that forwards-pass agrees with plain forwards-pass.
    @test y ≈ f(x...)

    # Check that ad and fd adjoints (approximately) agree.
    print_results && print_adjoints(adj_ad, adj_fd, rtol, atol)
    @test fd_isapprox(adj_ad, adj_fd, rtol, atol)
end

"""
    mean_function_tests(m::MeanFunction, X::AbstractVector)

Test _very_ basic consistency properties of the mean function `m`.
"""
function mean_function_tests(m::MeanFunction, x::AbstractVector)
    @test AbstractGPs._map(m, x) isa AbstractVector
    @test length(ew(m, x)) == length(x)
end

"""
    differentiable_mean_function_tests(
        m::MeanFunction,
        ȳ::AbstractVector,
        x::AbstractVector,
    )

Ensure that the gradient w.r.t. the inputs of `MeanFunction` `m` are approximately correct.
"""
function differentiable_mean_function_tests(
    m::MeanFunction,
    ȳ::AbstractVector{<:Real},
    x::AbstractVector{<:Real};
    rtol=_rtol,
    atol=_atol,
)
    # Run forward tests.
    mean_function_tests(m, x)

    # Check adjoint.
    @assert length(ȳ) == length(x)
    return adjoint_test(x -> ew(m, x), ȳ, x; rtol=rtol, atol=atol)
end

# function differentiable_mean_function_tests(
#     m::MeanFunction,
#     ȳ::AbstractVector{<:Real},
#     x::ColVecs{<:Real};
#     rtol=_rtol,
#     atol=_atol,
# )
#     # Run forward tests.
#     mean_function_tests(m, x)

#     @assert length(ȳ) == length(x)
#     adjoint_test(X->ew(m, ColVecs(X)), ȳ, x.X; rtol=rtol, atol=atol)
# end

function differentiable_mean_function_tests(
    rng::AbstractRNG, m::MeanFunction, x::AbstractVector; rtol=_rtol, atol=_atol
)
    ȳ = randn(rng, length(x))
    return differentiable_mean_function_tests(m, ȳ, x; rtol=rtol, atol=atol)
end

"""
    abstractgp_interface_tests(
        f::AbstractGP,
        x::AbstractVector,
        z::AbstractVector;
        atol=1e-12,
    )

Check that the `AbstractGP` interface is at least implemented for `f` and is
self-consistent. `x` and `z` must be valid inputs for `f`. For tests to pass, the minimum
eigenvalue of `cov(f, x)` must be greater than `-eig_tol`.
"""
function abstractgp_interface_tests(
    f::AbstractGP, x::AbstractVector, z::AbstractVector; eig_tol::Real=1e-12, σ²::Real=1e-9
)
    @assert length(x) ≠ length(z)

    # Verify that `mean` works and is the correct length and type.
    m = mean(f, x)
    @test m isa AbstractVector{<:Real}
    @test length(m) == length(x)

    # Verify that cov(f, x, z) works, is the correct size and type.
    C_xy = cov(f, x, z)
    @test C_xy isa AbstractMatrix{<:Real}
    @test size(C_xy) == (length(x), length(z))

    # Reversing arguments transposes the return.
    @test C_xy ≈ cov(f, z, x)'

    # Verify cov(f, x) works, is the correct size and type.
    C_xx = cov(f, x)
    @test size(C_xx) == (length(x), length(x))

    # Check that C_xx is positive definite.
    @test eigmin(Symmetric(C_xx)) > -eig_tol

    # Check that C_xx is consistent with cov(f, x, x).
    @test C_xx ≈ cov(f, x, x)

    # Check that var(f, x) works, is the correct size and type.
    C_xx_diag = var(f, x)
    @test C_xx_diag isa AbstractVector{<:Real}
    @test length(C_xx_diag) == length(x)

    # Check C_xx_diag is consistent with cov(f, x).
    @test C_xx_diag ≈ diag(C_xx)

    # Check that mean_and_cov is consistent.
    let
        m, C = mean_and_cov(f, x)
        @test m ≈ mean(f, x)
        @test C ≈ cov(f, x)
    end

    # Check that mean_and_var is consistent.
    let
        m, c = mean_and_var(f, x)
        @test m ≈ mean(f, x)
        @test c ≈ var(f, x)
    end

    # Construct a FiniteGP, and check that all standard methods defined on it at least run.
    fx = f(x, σ²)
    fz = f(z, σ²)
    @test mean(fx) ≈ mean(f, x)
    @test cov(fx) ≈ cov(f, x) + fx.Σy
    @test cov(fx, fz) ≈ cov(f, x, z)
    @test first(mean_and_cov(fx)) ≈ mean(f, x)
    @test last(mean_and_cov(fx)) ≈ cov(f, x)
    @test mean.(marginals(fx)) ≈ mean(f, x)
    @test var.(marginals(fx)) ≈ var(f, x) .+ diag(fx.Σy)

    # Generate, compute logpdf, compare against VFE and DTC.
    y = rand(fx)
    @test length(y) == length(x)
    @test logpdf(fx, y) isa Real
    @test elbo(fx, y, f(x)) ≈ logpdf(fx, y) rtol = 1e-5 atol = 1e-5
    @test dtc(fx, y, f(x)) ≈ logpdf(fx, y) rtol = 1e-5 atol = 1e-5
end


text:
using AdventOfCodeSolutions
using Test

function input(puzzle::Puzzle{2020, 1, n}) where n
    io = openInput(puzzle)
    inputs = split(read(io, String), '\n', keepempty=false)
    return parse.(Int, inputs)
end

function solveForProduct(product)
    solution = Iterators.filter(
        ns -> sum(ns) == 2020,
        product
    ) |> first

    prod(solution)
end

function solve(::Puzzle{2020, 1, 1}, numbers)
    product = Iterators.product(numbers, numbers)
    solveForProduct(product)
end

function solve(::Puzzle{2020, 1, 2}, numbers)
    product = Iterators.product(numbers, numbers, numbers)
    solveForProduct(product)
end

@testset "Day 1" begin
    @test(solve(Puzzle(1, 1), input(Puzzle(1, 1))) == 73371)
    @test(solve(Puzzle(1, 2), input(Puzzle(1, 2))) == 127642310)
end


text:
using GraphBLAS
using Documenter

DocMeta.setdocmeta!(GraphBLAS, :DocTestSetup, :(using GraphBLAS); recursive=true)

makedocs(;
    modules=[GraphBLAS],
    authors="Will Kimmerer <wrkimmerer@outlook.com> and contributors",
    repo="https://github.com/Wimmerer/GraphBLAS.jl/blob/{commit}{path}#{line}",
    sitename="GraphBLAS.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Wimmerer.github.io/GraphBLAS.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Wimmerer/GraphBLAS.jl",
)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# NOTES on range reduction
# [1] compute numbers near pi: http://www.cs.berkeley.edu/~wkahan/testpi/nearpi.c
# [2] range reduction: http://hal-ujm.ccsd.cnrs.fr/docs/00/08/69/04/PDF/RangeReductionIEEETC0305.pdf
# [3] precise addition, see Add22: http://ftp.nluug.nl/pub/os/BSD/FreeBSD/distfiles/crlibm/crlibm-1.0beta3.pdf

# Examples
# ΓΓ = 6411027962775774 / 2^45  # see [2] above, section 1.2
# julia> mod(ΓΓ, 2pi)    # "naive" way - easily wrong
# 7.105427357601002e-15
# julia> mod2pi(ΓΓ)      # using function provided here
# 2.475922546353431e-18
# Wolfram Alpha: mod(6411027962775774 / 2^45, 2pi)
# 2.475922546353430800060268586243862383453213646146648435... × 10^-18

# Test Cases. Each row contains: x and x mod 2pi (as from Wolfram Alpha)
# The values x are:
# -pi/2, pi/2, -pi, pi, 2pi, -2pi
#   (or rather, the Float64 approx to those numbers.
#   Thus, x mod pi will result in a small, but positive number)
# ΓΓ = 6411027962775774 / 2^47
#   from [2], section 1.2:
#   the Float64 greater than 8, and less than 2**63 − 1 closest to a multiple of π/4 is
#   Γ = 6411027962775774 / 2^48. We take ΓΓ = 2*Γ to get cancellation with pi/2 already
# 3.14159265359, -3.14159265359
# pi/16*k +/- 0.00001 for k in [-20:20] # to cover all quadrants
# numerators of continuous fraction approximations to pi
#   see http://oeis.org/A002485
#   (reason: for max cancellation, we want x = k*pi + eps for small eps, so x/k ≈ pi)

testCases = [
      -1.5707963267948966         4.71238898038469
       1.5707963267948966       1.5707963267948966
       -3.141592653589793       3.1415926535897936
        3.141592653589793        3.141592653589793
        6.283185307179586        6.283185307179586
       -6.283185307179586   2.4492935982947064e-16
          45.553093477052       1.5707963267948966
            3.14159265359            3.14159265359
           -3.14159265359       3.1415926535895866
      -3.9269808169872418        2.356204490192345
        -3.73063127613788       2.5525540310417068
      -3.5342817352885176        2.748903571891069
       -3.337932194439156        2.945253112740431
      -3.1415826535897935        3.141602653589793
      -2.9452331127404316        3.337952194439155
      -2.7488835718910694       3.5343017352885173
      -2.5525340310417075        3.730651276137879
       -2.356184490192345       3.9270008169872415
      -2.1598349493429834        4.123350357836603
      -1.9634854084936209        4.319699898685966
      -1.7671358676442588        4.516049439535328
      -1.5707863267948967         4.71239898038469
      -1.3744367859455346        4.908748521234052
      -1.1780872450961726        5.105098062083414
      -0.9817377042468104        5.301447602932776
      -0.7853881633974483       5.4977971437821385
      -0.5890386225480863          5.6941466846315
      -0.3926890816987242        5.890496225480862
      -0.1963395408493621       6.0868457663302244
                   1.0e-5                   1.0e-5
      0.19635954084936205      0.19635954084936205
       0.3927090816987241       0.3927090816987241
       0.5890586225480862       0.5890586225480862
       0.7854081633974482       0.7854081633974482
       0.9817577042468103       0.9817577042468103
       1.1781072450961723       1.1781072450961723
       1.3744567859455343       1.3744567859455343
       1.5708063267948964       1.5708063267948964
       1.7671558676442585       1.7671558676442585
       1.9635054084936205       1.9635054084936205
        2.159854949342982        2.159854949342982
       2.3562044901923445       2.3562044901923445
       2.5525540310417063       2.5525540310417063
       2.7489035718910686       2.7489035718910686
       2.9452531127404304       2.9452531127404304
       3.1416026535897927       3.1416026535897927
       3.3379521944391546       3.3379521944391546
        3.534301735288517        3.534301735288517
       3.7306512761378787       3.7306512761378787
        3.927000816987241        3.927000816987241
      -3.9270008169872415        2.356184490192345
      -3.7306512761378796        2.552534031041707
      -3.5343017352885173       2.7488835718910694
      -3.3379521944391555        2.945233112740431
       -3.141602653589793       3.1415826535897935
      -2.9452531127404313       3.3379321944391553
       -2.748903571891069       3.5342817352885176
       -2.552554031041707       3.7306312761378795
       -2.356204490192345       3.9269808169872418
       -2.159854949342983        4.123330357836603
      -1.9635054084936208       4.3196798986859655
      -1.7671558676442587        4.516029439535328
      -1.5708063267948966         4.71237898038469
      -1.3744567859455346        4.908728521234052
      -1.1781072450961725        5.105078062083414
      -0.9817577042468104        5.301427602932776
      -0.7854081633974483        5.497777143782138
      -0.5890586225480863        5.694126684631501
     -0.39270908169872415        5.890476225480862
     -0.19635954084936208        6.086825766330224
                  -1.0e-5        6.283175307179587
      0.19633954084936206      0.19633954084936206
      0.39268908169872413      0.39268908169872413
       0.5890386225480861       0.5890386225480861
       0.7853881633974482       0.7853881633974482
       0.9817377042468103       0.9817377042468103
       1.1780872450961724       1.1780872450961724
       1.3744367859455344       1.3744367859455344
       1.5707863267948965       1.5707863267948965
       1.7671358676442586       1.7671358676442586
       1.9634854084936206       1.9634854084936206
       2.1598349493429825       2.1598349493429825
       2.3561844901923448       2.3561844901923448
       2.5525340310417066       2.5525340310417066
        2.748883571891069        2.748883571891069
       2.9452331127404308       2.9452331127404308
        3.141582653589793        3.141582653589793
        3.337932194439155        3.337932194439155
        3.534281735288517        3.534281735288517
        3.730631276137879        3.730631276137879
       3.9269808169872413       3.9269808169872413
                     22.0       3.1504440784612404
                    333.0       6.2743640266615035
                    355.0       3.1416227979431572
                 103993.0        6.283166177843807
                 104348.0        3.141603668607378
                 208341.0        3.141584539271598
                 312689.0    2.9006993893361787e-6
                 833719.0       3.1415903406703767
               1.146408e6       3.1415932413697663
               4.272943e6        6.283184757600089
               5.419351e6       3.1415926917902683
              8.0143857e7        6.283185292406739
             1.65707065e8       3.1415926622445745
             2.45850922e8        3.141592647471728
             4.11557987e8    2.5367160519636766e-9
            1.068966896e9         3.14159265254516
            2.549491779e9    4.474494938161497e-10
            6.167950454e9        3.141592653440059
          1.4885392687e10   1.4798091093322177e-10
          2.1053343141e10         3.14159265358804
        1.783366216531e12    6.969482408757582e-13
        3.587785776203e12        3.141592653589434
        5.371151992734e12       3.1415926535901306
        8.958937768937e12        6.283185307179564
      1.39755218526789e14          3.1415926535898
      4.28224593349304e14       3.1415926535897927
     5.706674932067741e15    4.237546464512562e-16
     6.134899525417045e15        3.141592653589793
]

function testModPi()
    numTestCases = size(testCases,1)
    modFns = [mod2pi]
    xDivisors = [2pi]
    errsNew, errsOld = Vector{Float64}(), Vector{Float64}()
    for rowIdx in 1:numTestCases
        xExact = testCases[rowIdx,1]
        for colIdx in 1:1
            xSoln = testCases[rowIdx,colIdx+1]
            xDivisor = xDivisors[colIdx]
            modFn = modFns[colIdx]
            # 2. want: xNew := modFn(xExact)  ≈  xSoln  <--- this is the crucial bit, xNew close to xSoln
            # 3. know: xOld := mod(xExact,xDivisor) might be quite a bit off from xSoln - that's expected
            xNew = modFn(xExact)
            xOld = mod(xExact,xDivisor)

            newDiff  = abs(xNew - xSoln)  # should be zero, ideally (our new function)
            oldDiff  = abs(xOld - xSoln)  # should be zero in a perfect world, but often bigger due to cancellation
            oldDiff  = min(oldDiff, abs(xDivisor - oldDiff)) # we are being generous here:
            # if xOld happens to end up "on the wrong side of 0", eg
            # if xSoln = 3.14 (correct), but xOld reports 0.01,
            # we don't take the long way around the circle of 3.14 - 0.01, but the short way of 3.1415.. - (3.14 - 0.1)
            push!(errsNew,abs(newDiff))
            push!(errsOld,abs(oldDiff))
        end
    end
    sort!(errsNew)
    sort!(errsOld)
    totalErrNew = sum(errsNew)
    totalErrOld = sum(errsOld)
    @test totalErrNew ≈ 0.0
end
testModPi()

# 2pi
@test mod2pi(10) ≈ mod(10,2pi)
@test mod2pi(-10) ≈ mod(-10,2pi)
@test mod2pi(355) ≈ 3.1416227979431572
@test mod2pi(Int32(355)) ≈ 3.1416227979431572
@test mod2pi(355.0) ≈ 3.1416227979431572
@test mod2pi(355.0f0) ≈ 3.1416228f0
@test mod2pi(Int64(2)^60) == mod2pi(2.0^60)
@test_throws ArgumentError mod2pi(Int64(2)^60-1)

@testset "rem_pio2_kernel" begin
    # test worst case
    x = 6381956970095103.0 * 2.0^797
    a = setprecision(BigFloat, 4096) do
        rem(big(x), big(pi)/2, RoundNearest)
    end

    n, yrem = Base.Math.rem_pio2_kernel(x)
    y=yrem.hi+yrem.lo
    @test a-y<nextfloat(y)/2
    # The following has easy and hard cases in each interval. A hard case is one
    # where x ≈ k*pi/2 for some integer k.
    cases = [0.0, pi/6, # -π/4 <= x <= π/4
             2*pi/4-0.1, 2*pi/4, # -2π/4 <= x <= 2π/4
             3*pi/4-0.1, # -3π/4 <= x <= 3π/4
             pi-0.1, Float64(pi), # -4π/4 <= x <= 4π/4
             5*pi/4-0.1, # -5π/4 <= x <= 5π/4
             6*pi/4-0.1, 6*pi/4, # -6π/4 <= x <= 6π/4
             7*pi/4-0.1, # -7π/4 <= x <= 7π/4
             2*pi, 2*pi-0.1, # -8π/4 <= x <= 8π/4
             9*pi/4-0.1, # -9π/4 <= x <= 9π/4
             2.0^10*pi/4, # -2.0^20π/2 <= x <= 2.0^20π/2
             2.0^30*pi/4, # |x| >= 2.0^20π/2, idx < 0
             2.0^80*pi/4] # |x| >= 2.0^20π/2, idx > 0-0.22370138542135648

     # ieee754_rem_pio2_return contains the returned value from the ieee754_rem_pio2
     # function in openlibm: https://github.com/JuliaLang/openlibm/blob/0598080ca09468490a13ae393ba17d8620c1b201/src/e_rem_pio2.c
     ieee754_rem_pio2_return = adjoint([1.5707963267948966 1.5707963267948966;
                                1.0471975511965979 -1.0471975511965979;
                                0.10000000000000014 -0.10000000000000014;
                                6.123233995736766e-17 -6.123233995736766e-17;
                                -0.6853981633974481 0.6853981633974481;
                                0.10000000000000021 -0.10000000000000021;
                                1.2246467991473532e-16 -1.2246467991473532e-16;
                                -0.6853981633974481 0.6853981633974481;
                                0.09999999999999983 -0.09999999999999983;
                                1.8369701987210297e-16 -1.8369701987210297e-16;
                                -0.6853981633974484 0.6853981633974484;
                                2.4492935982947064e-16 -2.4492935982947064e-16;
                                0.0999999999999999 -0.0999999999999999;
                                -0.6853981633974484 0.6853981633974484;
                                3.135095805817224e-14 -3.135095805817224e-14;
                                3.287386219680602e-8 -3.287386219680602e-8;
                                -0.1757159771004682 0.1757159771004682])

    for (i, case) in enumerate(cases)
        # negative argument
        n, ret = Base.Math.rem_pio2_kernel(-case)
        ret_sum = ret.hi+ret.lo
        ulp_error = (ret_sum-ieee754_rem_pio2_return[1, i])/eps(ieee754_rem_pio2_return[1, i])
        @test ulp_error <= 0.5
        diff = Float64(mod(big(-case), big(pi)/2))-(ret.hi+ret.lo)
        @test abs(diff) in (0.0, 1.5707963267948966, 1.5707963267948968)
        # positive argument
        n, ret = Base.Math.rem_pio2_kernel(case)
        ret_sum = ret.hi+ret.lo
        ulp_error = (ret_sum-ieee754_rem_pio2_return[2, i])/eps(ieee754_rem_pio2_return[2, i])
        @test ulp_error <= 0.5
        diff = Float64(mod(big(case), big(pi)/2))-(ret.hi+ret.lo)
        @test abs(diff) in (0.0, 1.5707963267948966, 1.5707963267948968)
    end
end
@testset "rem_pio2_kernel and mod2pi" begin
    for int in (3632982096228748, 1135326194816)
        bignum = int*big(pi)/2+0.00001
        bigrem = rem(bignum, big(pi)/2, RoundDown)
        fnum = Float64(bignum)
        n, ret = Base.Math.rem_pio2_kernel(fnum)
        @test mod2pi(fnum) == (ret.hi+ret.lo)
    end
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
import Base: iterate
import Base.Iterators.take

"""
Structure to hold the iterates of the dual conditional gradient method.

DualCGIterable(M, b, A; optTol=√eps)

Construct a dual conditional gradient iterator for the problem

    minimize  ½‖Mx-b‖²  subject to gauge(A|x) ≤ 1.

The absolute tolerance `optTol` is used to determine the stopping condition:

    gap ≤ optTol(1+‖b‖₂)

"""
mutable struct DualCGIterable{T1<:Float64, T2<:Vector{T1}, T3<:AbstractArray{T1}, T4<:AbstractLinearOp, T5<:AbstractAtomicSet, T6<:AbstractAtom}
    M::T4           # linear operator
    Mt::T4          # adjoint of the linear operator
    b::T2           # observation
    A::T5           # atomic set
    r::T2           # residual
    z::T3           # negative gradient
    a::T6           # exposed atom
    Ma::T3          # image of the atom
    Mx::T3          # image of primal variable
    Δr::T3          # residual step
    τ::T1           # current scale of atomic set
    optTol::T1      # tolerance
    function DualCGIterable(M::AbstractLinearOp, b::Vector{Float64}, A::AbstractAtomicSet, optTol = sqrt(eps(1.0)))
        r = copy(b)
        Mt = M'
        z = Mt*r
        a = expose(A, z)
        Ma = M*a
        Mx = copy(Ma); fill!(Mx, 0.0)
        Δr = Ma - Mx
        τ = 0.0
        new{Float64, Vector{Float64}, AbstractArray{Float64}, AbstractLinearOp, AbstractAtomicSet, AbstractAtom}(M, Mt, b, A, r, z, a, Ma, Mx, Δr, τ, optTol)
    end
end

function oracleExit(p::DualCGIterable, gap::Float64)
    return gap < p.optTol
end

"""
Defines one complete iteration of the dual conditional gradient method.
"""
function iterate(p::DualCGIterable, k::Int=0)

    # Compute  M*a  where  a = expose(τA, z)
    expose!(p.A, p.z, p.a)
    # p.Ma .= p.M*p.a
    mul!(p.Ma, p.M, p.a)
    rmul!(p.Ma, p.τ)

    # Compute search direction Δr = M(a - x)
    @. p.Δr = p.Ma - p.Mx

    # Compute gap and possible exit
    gap = sumdot(p.Δr, p.r)
    oracleExit(p, gap) && return (gap, p.r), k+1

    if iszero(p.Mx)
        # Mx ← Ma
        p.Mx .= p.Ma
        # r ← b - Mx
        p.r .= p.b - sum(p.Mx, dims=2)[:]
    else
        # Linesearch
        α = linesearch(p.Δr, p.r)
        # r ← r - Δr*α
        mul!(p.r, p.Δr, α, -1, 1)
        # Mx ← Mx + Δr*Diagonal(α)
        broadcast!(*, p.Δr, p.Δr, α')
        BLAS.axpy!(1, p.Δr, p.Mx)
    end

    # Update gradient: z ← M'r
    mul!(p.z, p.Mt, p.r)

    return (gap, p.r), k+1
end

"""
    sumdot(X::Matrix, y::Vector)
Compute the sum of the dotproducts of each column of `X`
with the vector `y`.
"""
function sumdot(X::AbstractMatrix{Float64}, y::Vector{Float64})
    fdot = x -> dot(x, y)
    return sum(fdot, eachcol(X))
end

sumdot(x::AbstractVector{Float64}, y::Vector{Float64}) = dot(x, y)

"""
    setRadius!(p::DualCGIterable, τ::Float64)

Increase to `τ` the radius of the current Lasso problem. This routine
reflects the increase in the radius on various quantities of the dual
CG method.

If the current radius 0, then this routine only sets the new parameter,
but doesn't rescale any of the dual iterates.
"""
function setRadius!(p::DualCGIterable, τ::Float64)
    if p.τ > 0
        rescaleIterates!(p, τ/p.τ)
    end
    p.τ = τ
end

function rescaleIterates!(p::DualCGIterable, s::Float64)
    p.Mx .*= s
    Mxs = vec(sum(p.Mx, dims=2))
    p.r .= p.b .- Mxs
    p.z .= p.M'*p.r
end

function replaceIterates!(p::DualCGIterable, M::AbstractMatrix{Float64}, F::AbstractFace, c::Vector{Float64})
    MF = M*F
    p.Mx .= MF*c
    p.Mx .*= p.τ/sum(c)
    Mxs = vec(sum(p.Mx, dims=2))
    p.r .= p.b .- Mxs
    p.z .= p.M'*p.r
end

getNegGradient(p::DualCGIterable) = p.z
getResidual(p::DualCGIterable) = p.r

"""
    linesearch(Δr::Vector, r::Vector)

Solve univariate least-squares problem

    minimize_{θ∈[0,1]} ||θ⋅Δr-r||₂,

which has the solution

    θ = min(1, dot(Δr,Δr)/dot(Δr,Δr))
"""
function linesearch(Δr::AbstractVector{Float64}, r::Vector{Float64})
    return min(1.0, dot(Δr,r)/dot(Δr,Δr))
end

"""
    linesearch(Δr::Matrix, r::Vector)

Solve the box-constrained least squares problem

    minimize_{θ∈[0,1]ⁿ} ||Δr⋅θ - r||₂
"""
function linesearch(Δr::AbstractMatrix{Float64}, r::Vector{Float64})
    n = size(Δr, 2)
    bl = zeros(n)
    bu = ones(n)
    Q = Quadratic(Δr'*Δr, Δr'r)
    lbfgsq = LBFGSBQuad(Q, bl, bu, m=2)
    θ = boxquad(lbfgsq)
    return θ
end

"""
    primal recover
"""
function primalrecover!(p::DualCGIterable, sol::Solution, α::Float64, feaTol::Float64, exitFlag)
    flag = exitFlag
    s = support(p.A, p.z)
    face!(p.A, p.z/s, sol.Fnew)
    ϵ = copy(getResidual(p)); ϵ .*= sqrt(2*α)/norm(ϵ)
    face_project!(p.M, sol.Fnew, p.b-ϵ, sol.cnew, sol.r)
    sol.feasnew = norm(sol.r + ϵ)^2/2
    if sol.feasnew  ≤ sol.feas
        update!(sol)
    end
    if sol.feas - α ≤ feaTol
        flag = :feasible
    end
    return flag
end


text:
"""
    theme(name; kwrgs...)

Offer themes support. NAME is the theme name. So far the three options are:

- `modern`: - This is the default theme (same as GMT modern theme but with thinner FRAME_PEN [0.75p])
- `classic`: - The GMT classic theme
- `dark`: - A modern theme variation with dark background.
- `A0|2[XY|XX|YY][atg][ag][g][H][V][NT|nt][ITit][Graph][Dark]` Make a composition of these to select a theme.
   The main condition is that it starts with an A (Annotate). Hence `A2` means annotate two axis
   and `A0` means no axes at all. `XY` means to plot only left and bottom axes, `YY` only left and right
   and `XX` bottom and top. `atg` (or `afg`) means annotate, tick and grid lines. `ag` does not tick.
   `H` and `V` means grid lines will only be horizontal or vertical. Note, these require `atg` or `ag`.
   `NT` stands for not ticks at all and `IT` plots the ticks inside the axes. `Graph` adds a vector and
   to the end of each axis (sets `XY`), and `Dark` put the background in dark mode.
   - Example: `A2YYg` -> plot left and right axes (only) and add grid lines.
   - Example: `A2Graph` -> plot left and right axes (only) and adds arrows at the end of them

On top of the modern mode variations (so far `dark` only) one can set the following `kwargs` options:

- `noticks` or `no_ticks`: Axes will have annotations but no tick marks
- `inner_ticks` or `innerticks`: - Ticks will be drawn inside the axes instead of outside.
- `gray_grid` or `graygrid`: - When drawing grid line use `gray` instead of `black`
- `save`: - Save the name in the directory printed in shell by ``gmt --show-userdir`` and make it permanent.
- `reset`: - Remove the saved theme name and return to the default `modern` theme.

Note: Except `save` and `reset`, the changes operated by the `kwargs` are temporary and operate only until
an image is `show`(n) or saved.

This function can be called alone, e.g. ``theme("dark")`` or as an option in the ``plot()`` module.
"""
function theme(name="modern"; kwargs...)
	# Provide the support for themes
	global API

	(GMTver < v"6.2.0") && return nothing
	(!isa(API, Ptr{Nothing}) || API == C_NULL) && (API = GMT_Create_Session("GMT", 2, GMT_SESSION_NOEXIT + GMT_SESSION_EXTERNAL + GMT_SESSION_COLMAJOR))

	d = KW(kwargs)
	font = ((val = find_in_dict(d, [:font])[1]) !== nothing) ? string(val) : ""
	bg_color = ((val = find_in_dict(d, [:bg_color])[1]) !== nothing) ? string(val) : ""
	color = ((val = find_in_dict(d, [:fg_color])[1]) !== nothing) ? string(val) : ""

	# Some previous calls may have changed these and a new theme option may be caught with the pens down
	def_fig_axes[1]  = def_fig_axes_bak		# So that we always start with the defaults
	def_fig_axes3[1] = def_fig_axes3_bak

	_name = string(name)
	if (_name == "classic")
		theme_classic()
	else
		theme_modern()			# All other themes are variations over the GMT_MODERN theme
		if (_name == "dark" || contains(_name, "Dark"))
			helper_theme_fonts_colors(font, color, bg_color, true)
		elseif ((_name != "modern") || (font != "" || color != ""))	# bg_color alone wont trigger next call
			helper_theme_fonts_colors(font, color, bg_color, false)
		end
		parse_theme_names(_name)
	end

	(find_in_dict(d, [:noticks :no_ticks])[1] !== nothing) && helper_theme_noticks()		# No ticks
	(find_in_dict(d, [:inner_ticks :innerticks])[1] !== nothing) &&  helper_theme_inticks()	# Inner ticks
	if (find_in_dict(d, [:gray_grid :graygrid])[1] !== nothing)
		gmtlib_setparameter(API, "MAP_GRID_PEN_PRIMARY", "auto,gray")
		gmtlib_setparameter(API, "MAP_GRID_PEN_SECONDARY", "auto,gray")
	end

	isOn = true				# Means this theme will be reset to the default (modern) in showfig()
	if (haskey(d, :save))
		f = joinpath(readlines(`$(joinpath("$(GMT_bindir)", "gmt")) --show-userdir`)[1], "theme_jl.txt")
		(isfile(f)) && rm(f)
		(_name == "reset") ? theme_modern() : write(f, _name)
		isOn = false		# So we wont reset defaults in showfig()
	end
	ThemeIsOn[1] = isOn
	return nothing
end

# ---------------------------------------------------------------------------------------------------
function parse_theme_names(name::String)
	# Basically only short forms for setting -B option: A2atg, A2afg, A2ag, A2agHnt, A2agV, A2XY..., A0...
	# A2 means annotate two axes (default to WSrt)
	# If the name contains XY plot only WS (no EN). Containing "YY" or "XX" plots only WE or SN, respectively.
	# "afg" and "ag" have the normal meaning of 'a'nnotate, 't'icks & 'g'rid lines. 'f' can be used instead of 't'
	# "H" in name plots only the horizontal grid lines and 'V' only the vertical (requires "afg" or "ag")
	# "nt" or "NT" anywhere after "A2" means no ticks
	# "it" or "IT" means inside ticks
	# "A0" means no axes, but it can have annotations, ticks and gridlines
	if (name[1] == 'A')			# Name should either start with A2 or A0
		t1 = "a"
		if (contains(name, "atg") || contains(name, "afg"))
			t1 = "afg"
			contains(name, "H") && (t1 = "xaf -Byafg")		# Only horizontal grid lines
			contains(name, "V") && (t1 = "xafg -Byaf")		# Only vertical grid lines
		elseif (contains(name, "ag") || contains(name, "g"))
			t1 = "ag"
			contains(name, "H") && (t1 = "xa -Byag")		# Only horizontal grid lines
			contains(name, "V") && (t1 = "xag -Bya")		# Only vertical grid lines
		end
		if contains(name, "Graph")							# Add a vector to the end of each axis
			gmtlib_setparameter(API, "MAP_FRAME_TYPE", "graph")
			name *= "XY"									# Must be if already there, no problems
		end
		t2 = contains(name, "XY") ? "WS" : (contains(name, "YY")) ? "WE" : (contains(name, "XX")) ? "SN" : "WSrt"
		t1, t2 = " -B" * t1, " -B" * t2

		if (name[2] == '0')
			if (length(name) == 2)  t1, t2 = " ", ""		# Really no axes no annotations/ticks
			else                    gmtlib_setparameter(API, "MAP_FRAME_PEN", "0.001,white@100") # No axes, but need this sad trick
			end
		end

		def_fig_axes[1] = t1 * t2
		if     (contains(name, "nt") || contains(name, "NT"))  helper_theme_noticks()	# No ticks
		elseif (contains(name, "it") || contains(name, "IT"))  helper_theme_inticks()	# Inside ticks
		end
		gmtlib_setparameter(API, "MAP_GRID_PEN_PRIMARY", "auto,gray")
		gmtlib_setparameter(API, "MAP_GRID_PEN_SECONDARY", "auto,gray")
		return true
	end
	return false
end

# ---------------------------------------------------------------------------------------------------
function helper_theme_inticks()
	# Inside ticks. Fixed length because GMT has no way to set it auto
	gmtlib_setparameter(API, "MAP_TICK_LENGTH_PRIMARY", "-4p")
	gmtlib_setparameter(API, "MAP_TICK_LENGTH_SECONDARY", "-12p")
end

# ---------------------------------------------------------------------------------------------------
function helper_theme_noticks()
	gmtlib_setparameter(API, "MAP_TICK_LENGTH_PRIMARY", "0/0")
	gmtlib_setparameter(API, "MAP_TICK_LENGTH_SECONDARY", "0/0")
end

# ---------------------------------------------------------------------------------------------------
function helper_theme_fonts_colors(font, color, bg_color, dark::Bool)
	fonts  = (font == "")  ? ["Helvetica", "Helvetica", "Helvetica"] : [font, font, font]
	colors = (color == "") ? (dark ? ["gray92", "gray86"] : ["black", "black"]) : [color, color]
	(bg_color == "") && (bg_color = (dark) ? "5/5/35" : "white")
	fonts_colors_settings(fonts, colors, bg_color)
end

# ---------------------------------------------------------------------------------------------------
function fonts_colors_settings(fonts, colors, bg_color)
	# This function serves mainly the dark theme but can be used by any other theme that
	# wants to change font and/or colors.
	gmtlib_setparameter(API, "FONT_ANNOT_PRIMARY", "auto,$(fonts[1]),$(colors[1])")
	gmtlib_setparameter(API, "FONT_ANNOT_SECONDARY", "auto,$(fonts[1]),$(colors[1])")
	gmtlib_setparameter(API, "FONT_HEADING", "auto,$(fonts[2]),$(colors[1])")
	gmtlib_setparameter(API, "FONT_LABEL", "auto,$(fonts[1]),$(colors[1])")
	gmtlib_setparameter(API, "FONT_LOGO", "auto,$(fonts[3]),$(colors[1])")
	gmtlib_setparameter(API, "FONT_TAG", "auto,$(fonts[1]),$(colors[1])")
	gmtlib_setparameter(API, "FONT_TITLE", "auto,$(fonts[2]),$(colors[1])")
	gmtlib_setparameter(API, "MAP_DEFAULT_PEN", "0.25p,$(colors[1])")
	gmtlib_setparameter(API, "MAP_FRAME_PEN", "0.75,$(colors[1])")
	gmtlib_setparameter(API, "MAP_GRID_PEN_PRIMARY", "auto,$(colors[1])")
	gmtlib_setparameter(API, "MAP_GRID_PEN_SECONDARY", "auto,$(colors[2])")
	gmtlib_setparameter(API, "MAP_TICK_PEN_PRIMARY", "auto,$(colors[1])")
	gmtlib_setparameter(API, "MAP_TICK_PEN_SECONDARY", "auto,$(colors[2])")
	gmtlib_setparameter(API, "PS_PAGE_COLOR", "$bg_color")
end

# ---------------------------------------------------------------------------------------------------
function theme_modern()
	# Set the MODERN mode settings
	(GMTver < v"6.2.0") && return nothing
	swapmode(API, classic=false)		# Set GMT->current.setting.run_mode = GMT_MODERN
	reset_defaults(API)					# Set the modern mode settings
	gmtlib_setparameter(API, "MAP_FRAME_PEN", "0.75")
	!IamModern[1] && swapmode(API, classic=true)	# Reset GMT->current.setting.run_mode = GMT_CLASSIC
	return nothing
end

# ---------------------------------------------------------------------------------------------------
function theme_classic()
	(GMTver < v"6.2.0") && return nothing
	swapmode(API, classic=true)			# Set GMT->current.setting.run_mode = GMT_CLASSIC
	reset_defaults(API)					# Set the classic mode settings
end

# ---------------------------------------------------------------------------------------------------
function swapmode(API; classic::Bool=true)
	# GMT6.2 shot out the access the gmtinit_conf_classic() and gmtinit_conf_modern() functions
	# (declared them GMT_LOCAL) and direct access to members of the API mega-structure is,
	# for the time being, impractical so we must resort to tricks. This function changes the state
	# of GMT->current.setting.run_mode to MODERN or CLASSIC, which in turn calls the wished conf function.
	if (classic == false)
		GMT_Create_Options(API, 0, "-pdf lix")			# Cheat to make GMT go into MODERN mode
	else
		gmtlib_setparameter(API, "GMT_VERBOSE", "q")	# Shut up messages about session dir not found
		gmt_manage_workflow(API, GMT_END_WORKFLOW, NULL)		# Force reset to CLASSIC
		gmtlib_setparameter(API, "GMT_VERBOSE", "w")
	end
	return nothing
end


text:
function compare_images(testfn::AbstractString, reffn::AbstractString; sigma = (1,1), tol = 0.02)

    result = VisualTestResult(testfn, nothing, reffn, nothing, PROCESSING_ERROR, 1.0, nothing)

    # load test and ref images... if error, return immediately
    try
        result.testImage = load(testfn)
        result.refImage  = load(reffn)
    catch err
        result.err = err
        return result
    end

    # we loaded both images, now do the comparison
    result.diff = blurdiff(result.testImage, result.refImage, sigma)

    if result.diff > tol
        result.status = DOES_NOT_MATCH
        result.err = "Images differ.  Difference: $(result.diff)  tolerance: $tol"
    else
        @info "Reference image $reffn matches.  Difference: $(result.diff)"
        result.status = (result.diff == 0) ? EXACT_MATCH : CLOSE_MATCH
        result.err = nothing
    end

    result
end

function test_images(testfn::AbstractString, reffn::AbstractString; popup=isinteractive(), newfn = reffn, kw...)
    result = compare_images(testfn, reffn; kw...)

    if !success(result)
        @warn "Image did not match reference image $reffn. err: $(result.err)"

        if popup
            # open a popup and give us a chance to examine the images,
            # then ask to replace the reference
            @warn "Should we make this the new reference image?"
            replace_refimg_dialog(testfn, reffn; reffn_new = newfn)
        end
    end

    result
end

function test_images(visualtest::VisualTest; popup=isinteractive(), kw...)
    visualtest.testFunction(visualtest.testFilename, visualtest.args...; visualtest.kw...)
    test_images(visualtest.testFilename, visualtest.refFilename; popup=popup, kw...)
end

function test_images(visualtests::AbstractVector{VisualTest}; popup=isinteractive(), kw...)
    VisualTestResult[test_images(vtest, popup=popup, kw...) for vtest in visualtests]
end


text:
struct Bad {
  int x;
};

int main() {
  Bad b = new Bad;
  b->x = true;
  return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
@info "Running show tests"

const TEMPLATES_DIR = contractuser(PT.default_file())
const LICENSES_DIR = joinpath(TEMPLATES_DIR, "licenses")

function test_show(expected::AbstractString, observed::AbstractString)
    if expected == observed
        @test true
    else
        print_diff(expected, observed)
        @test :expected == :observed
    end
end

@testset "Show methods" begin
    @testset "Plugins" begin
        expected = """
            Readme:
              file: "$(joinpath(TEMPLATES_DIR, "README.md"))"
              destination: "README.md"
              inline_badges: false
              badge_order: DataType[Documenter{GitHubActions}, Documenter{GitLabCI}, Documenter{TravisCI}, GitHubActions, GitLabCI, TravisCI, AppVeyor, DroneCI, CirrusCI, Codecov, Coveralls, BlueStyleBadge, ColPracBadge, PkgEvalBadge]
              badge_off: DataType[]
            """
        test_show(rstrip(expected), sprint(show, MIME("text/plain"), Readme()))
    end

    @testset "Template" begin
        expected = """
            Template:
              authors: ["$USER"]
              dir: "$(contractuser(Pkg.devdir()))"
              host: "github.com"
              julia: v"1.0.0"
              user: "$USER"
              plugins:
                CompatHelper:
                  file: "$(joinpath(TEMPLATES_DIR, "github", "workflows", "CompatHelper.yml"))"
                  destination: "CompatHelper.yml"
                  cron: "0 0 * * *"
                Git:
                  ignore: String[]
                  name: nothing
                  email: nothing
                  branch: "main"
                  ssh: false
                  jl: true
                  manifest: false
                  gpgsign: false
                License:
                  path: "$(joinpath(LICENSES_DIR, "MIT"))"
                  destination: "LICENSE"
                ProjectFile:
                  version: v"0.1.0"
                Readme:
                  file: "$(joinpath(TEMPLATES_DIR, "README.md"))"
                  destination: "README.md"
                  inline_badges: false
                  badge_order: DataType[Documenter{GitHubActions}, Documenter{GitLabCI}, Documenter{TravisCI}, GitHubActions, GitLabCI, TravisCI, AppVeyor, DroneCI, CirrusCI, Codecov, Coveralls, BlueStyleBadge, ColPracBadge, PkgEvalBadge]
                  badge_off: DataType[]
                SrcDir:
                  file: "$(joinpath(TEMPLATES_DIR, "src", "module.jl"))"
                TagBot:
                  file: "$(joinpath(TEMPLATES_DIR, "github", "workflows", "TagBot.yml"))"
                  destination: "TagBot.yml"
                  trigger: "JuliaTagBot"
                  token: Secret("GITHUB_TOKEN")
                  ssh: Secret("DOCUMENTER_KEY")
                  ssh_password: nothing
                  changelog: nothing
                  changelog_ignore: nothing
                  gpg: nothing
                  gpg_password: nothing
                  registry: nothing
                  branches: nothing
                  dispatch: nothing
                  dispatch_delay: nothing
                Tests:
                  file: "$(joinpath(TEMPLATES_DIR, "test", "runtests.jl"))"
                  project: false
            """
        # `with_clean_gitconfig` requires Git to be installed, but if Git is not installed,
        # then we probably don't need to worry about any conflicting Git config files.
        f = () -> test_show(
            rstrip(expected),
            sprint(show, MIME("text/plain"), tpl(; authors=USER)),
        )
        if PT.git_is_installed()
            with_clean_gitconfig() do
                run(`git config user.name Tester`)
                run(`git config user.email te@st.er`)
                f()
            end
        else
            f()
        end
    end

    @testset "show as serialization" begin
        t1 = tpl()
        t2 = eval(Meta.parse(sprint(show, t1)))
        @test t1 == t2

        foreach((NoDeploy, GitHubActions)) do T
            p1 = Documenter{T}()
            p2 = eval(Meta.parse(sprint(show, p1)))
            @test p1 == p2
        end
    end
end


text:
struct Exp{X} <: AbstractFloat
    log::X
end

Base.iszero(::Exp) = false

Base.log(x::Exp) = x.log

function Base.show(io::IO, x::Exp)
    io = IOContext(io, :compact => true)
    print(io, "Exp(", log(x), ")")
end

Base.:*(a::Exp, b::Exp) = Exp(log(a) + log(b))

Base.:/(a::Number, b::Exp) = a * Exp(-log(b))

Base.inv(a::Exp) = Exp(-log(a))

Base.:^(a::Exp, b::Number) = Exp(b * log(a))

Base.:^(a::Exp, b::Integer) = Exp(b * log(a))

Base.:^(a::Exp, b) = Exp(b * log(a))
Base.:^(a::Exp, b::Rational) = Exp(b * log(a))

Base.promote_rule(T::Type{<:Real}, ::Type{Exp}) = T
Base.promote_rule(::Type{Exp}, T::Type{<:Real}) = T

Base.convert(::Type{<:Real}, x::Exp) = exp(log(x))

# Exp(3) / Exp(2)

# Exp(3) * Exp(2)

# 3 + Exp(2)

# 3 / Exp(4)

# Exp(2) ^ 3

# 1 / Exp(2)

# 2 + Exp(3) ^ 0


text:
# develop cascade tree init

using Caesar, RoME
using Pkg

## everything in csm2
# canCliqMargRecycle_StateMachine
# canCliqIncrRecycle_StateMachine
# isCliqUpSolved_StateMachine
# buildCliqSubgraph_StateMachine
# canCliqMargSkipUpSolve_StateMachine
# blockUntilChildrenHaveStatus_StateMachine
# trafficRedirectConsolidate459_StateMachine
# maybeNeedDwnMsg_StateMachine
# doesParentNeedDwn_StateMachine
# determineCliqNeedDownMsg_StateMachine
# slowIfChildrenNotUpSolved_StateMachine

## first steps in csm5
# canCliqMargRecycle_StateMachine
# canCliqIncrRecycle_StateMachine
# isCliqUpSolved_StateMachine
# buildCliqSubgraph_StateMachine
# canCliqMargSkipUpSolve_StateMachine
# blockUntilChildrenHaveStatus_StateMachine
# trafficRedirectConsolidate459_StateMachine
# maybeNeedDwnMsg_StateMachine
# doesParentNeedDwn_StateMachine
# determineCliqNeedDownMsg_StateMachine
# towardUpOrDwnSolve_StateMachine

allCSMFunctions = [
  doCliqDownSolve_StateMachine,
  cleanupAfterDownSolve_StateMachine,
  specialCaseRootDownSolve_StateMachine,
  canCliqDownSolve_StateMachine,
  checkUpsolveFinished_StateMachine,
  prepInitUp_StateMachine,
  doCliqUpSolveInitialized_StateMachine,
  rmUpLikeliSaveSubFg_StateMachine,
  wipRedirect459Dwn_StateMachine,
  waitChangeOnParentCondition_StateMachine,
  slowOnPrntAsChildrNeedDwn_StateMachine,
  towardUpOrDwnSolve_StateMachine,
  canCliqMargSkipUpSolve_StateMachine,
  tryDwnInitCliq_StateMachine,
  rmMsgLikelihoodsAfterDwn_StateMachine,
  blockSiblingStatus_StateMachine,
  slowIfChildrenNotUpSolved_StateMachine,
  blockUntilChildrenHaveStatus_StateMachine,
  dwnInitSiblingWaitOrder_StateMachine,
  collectDwnInitMsgFromParent_StateMachine,
  trafficRedirectConsolidate459_StateMachine,
  doAllSiblingsNeedDwn_StateMachine,
  maybeNeedDwnMsg_StateMachine,
  determineCliqNeedDownMsg_StateMachine,
  decideUpMsgOrInit_StateMachine,
  doesParentNeedDwn_StateMachine,
  attemptCliqInitUp_StateMachine,
  sendCurrentUpMsg_StateMachine,
  buildCliqSubgraph_StateMachine,
  buildCliqSubgraphForDown_StateMachine,
  isCliqUpSolved_StateMachine,
  checkChildrenAllUpRecycled_StateMachine,
  canCliqIncrRecycle_StateMachine,
  canCliqMargRecycle_StateMachine,
  exitStateMachine,
]

runfile = joinpath(pathof(Caesar) |> dirname |> dirname, "examples", "dev", "devCSMHexDelaysTest.jl")
include(runfile)

for csmDlyFnc in allCSMFunctions, cliqId in 1:6

  injectDelayBefore=[cliqId=>(csmDlyFnc=>5); ]
  # injectDelayBefore=[1=>(exitStateMachine=>0.1); ]
  # injectDelayBefore = nothing

  fg, didFail = csmHexDelayTest(injectDelayBefore)

  fid = open("/tmp/caesar/dbgCSMHexDelayMatrix2.log","a")
  println(fid,getLogPath(fg), " -- solve_passed=$(!didFail), ", injectDelayBefore)
  close(fid)

end

#


text:
# EXCLUDE FROM TESTING
# this example requires BenchmarkTools, which can be hard to install

using BenchmarkTools
include("reduce.jl")

dev = CuDevice(0)
@assert(capability(dev) >= v"3.0", "this example requires a newer GPU")

len = 10^7
input = ones(Int32, len)

## CUDAnative

# PTX generation
open(joinpath(@__DIR__, "reduce.jl.ptx"), "w") do f
    code_ptx(f, reduce_grid,
                        Tuple{typeof(+), CuDeviceArray{Int32,1},
                              CuDeviceArray{Int32,1}, Int32};
                        cap=v"6.1.0")
end

ctx = CuContext(dev)
benchmark_gpu = @benchmarkable begin
        gpu_reduce(+, gpu_input, gpu_output)
        val = Array(gpu_output)[1]
    end setup=(
        val = nothing;
        gpu_input = CuArray($input);
        gpu_output = similar(gpu_input)
    ) teardown=(
        gpu_input = nothing;
        gpu_output = nothing;
        gc()
    )
println(run(benchmark_gpu))
destroy(ctx)

## CUDA

# Entry-point wrappers
lib = Libdl.dlopen(joinpath(@__DIR__, "reduce.so"))
setup_cuda(input)    = ccall(Libdl.dlsym(lib, "setup"), Ptr{Void},
                             (Ptr{Cint}, Csize_t), input, length(input))
run_cuda(state)      = ccall(Libdl.dlsym(lib, "run"), Cint,
                             (Ptr{Void},), state)
teardown_cuda(state) = ccall(Libdl.dlsym(lib, "teardown"), Void,
                             (Ptr{Void},), state)

# Correctness check (not part of verify.jl which is meant to run during testing)
let
    using Base.Test
    cuda_state = setup_cuda(input)
    cuda_val = run_cuda(cuda_state)
    teardown_cuda(cuda_state)
    @assert cuda_val == reduce(+, input)
end

benchmark_cuda = @benchmarkable begin
        val = run_cuda(state)
    end setup=(
        val = nothing;
        state = setup_cuda($input);
    ) teardown=(
        teardown_cuda(state)
    )
println(run(benchmark_cuda))