text:
"""
    gcrange(sys; zeros = "none", coinner = false, atol = 0, atol1 = atol, atol2 = atol, rtol,
            fast = true, offset = sqrt(ϵ)) -> (sysr, sysx, info)

Compute for the descriptor system `sys = (A-λE,B,C,D)` with the transfer function matrix `G(λ)`,
the proper descriptor system `sysr = (Ar-λEr,Br,Cr,Dr)` with a full row rank
transfer function matrix `R(λ)` such that `Coimage(G(λ)) = Coimage(R(λ))` and the
descriptor system `sysx = (A-λE,B,Cx,Dx)` with the full column rank transfer function matrix `X(λ)`,
which satisfies

     G(λ) = X(λ)*R(λ) ,

representing a full rank factorization of `G(λ)`.
The number of rows of `R(λ)` is the normal rank `r` of `G(λ)`.
The rows of `R(λ)` form a rational basis of the coimage space of the rational matrix `G(λ)`.
A selected set of zeros of `G(λ)` are included as zeros of `R(λ)`.

The resulting named triple `info` contains `(nrank, nfuz, niuz) `, where `info.nrank = r`,
the normal rank of `G(λ)`, `info.nfuz` is the number of finite zeros of `sys` on
the boundary of the stability domain `Cs`, and `info.niuz` is the number of infinite zeros of `sys` in
the continuous-time case and is set to `0` in the discrete-time case.

The following options can be selected via the keyword parameter `zeros` for which zeros of `G(λ)`
to be included in `R(λ)`:

     "none"       - include no zeros (default)
     "all"        - include all zeros of `sys`
     "unstable"   - include all unstable zeros of `sys`
     "s-unstable" - include all strictly unstable zeros of `sys`, both finite and infinite
     "stable"     - include all stable zeros of `sys`
     "finite"     - include all finite zeros of `sys`
     "infinite"   - include all infinite zeros of `sys`

If `coinner = true`, the resulting basis `R(λ)` is _coinner_, i.e., `R(λ)*R(λ)' = I`, where `R(s)' = transpose(R(-s))` for a
continuous-time system with `λ = s` and `R(z)' = transpose(R(1/z))` for a discrete-time system with `λ = z`.
This option can be used only in conjunction with `zeros = "none"` or `zeros = "unstable"`.

For a continuous-time system `sys`, the stability domain `Cs` is defined as the set of
complex numbers with real parts at most `-β`,
while for a discrete-time system `sys`, `Cs` is the set of complex numbers with
moduli at most `1-β` (i.e., the interior of a disc of radius `1-β` centered in the origin).
The boundary offset `β` to be used to assess the stability of zeros and their number
on the boundary of `Cs` can be specified via the keyword parameter `offset = β`.
Accordingly, for a continuous-time system,
the boundary of `Cs` contains the complex numbers with real parts within the interval `[-β,β]`,
while for a discrete-time system, the boundary of `Cs` contains
the complex numbers with moduli within the interval `[1-β,1+β]`.
The default value used for `β` is `sqrt(ϵ)`, where `ϵ` is the working machine precision.

The keyword arguments `atol1`, `atol2` and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of `A`, `B`, `C` and `D`,
the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `E`, `B`, `C` and `D`.
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`.

For the assessment of zeros, the dual system pencil `transpose([A-λE B; C D])` is reduced to a
special Kronecker-like form (see [2]). In this reduction, the
performed rank decisions are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

_Method:_  The range computation method described in [1], is applied to the
dual descriptor system realization corresponding to the transpose of the
rational matrix `G(λ)`. The underlying pencil reduction algorithm of [2],
has been adapted to deal with several zero selection options. The computation of the involved
Kronecker-like form is based on the algorithm of [3].

_References:_

[1] Varga, A.
    A note on computing the range of rational matrices.
    arXiv:1707.0048, [https://arxiv.org/abs/1707.0048](https://arxiv.org/abs/1707.004), 2017.

[2] C. Oara.
    Constructive solutions to spectral and inner–outer factorizations
    respect to the disk. Automatica,  41, pp. 1855–1866, 2005.

[3] C. Oara and P. Van Dooren.
    An improved algorithm for the computation of structural invariants of a system pencil and related geometric aspects.
    Syst. Control Lett., 30:39–48, 1997.
"""
function gcrange(sys::DescriptorStateSpace{T}; zeros = "none", coinner::Bool = false,
                 atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
                 rtol::Real = sys.nx*eps(real(float(one(T))))*iszero(max(atol1,atol2)),
                 offset::Real = sqrt(eps(float(real(T)))), fast::Bool = true) where T
   sysr, sysx, info = grange(gdual(sys); zeros = zeros, inner = coinner, offset = offset,
                      atol = atol, atol1 = atol1, atol2 = atol2, rtol = rtol, fast = fast)
   return gdual(sysr), gdual(sysx), info
end
"""
    grange(sys; zeros = "none", atol = 0, atol1 = atol, atol2 = atol, rtol,
           fast = true, offset = sqrt(ϵ)) -> (sysr, sysx, info)

Compute for the descriptor system `sys = (A-λE,B,C,D)` with the transfer function matrix `G(λ)`,
the proper descriptor system `sysr = (Ar-λEr,Br,Cr,Dr)` with a full column rank
transfer function matrix `R(λ)` such that `Range(G(λ)) = Range(R(λ))` and the
descriptor system `sysx = (A-λE,B,Cx,Dx)` with the full row rank transfer function matrix `X(λ)`,
which satisfies

     G(λ) = R(λ)*X(λ) ,

representing a full rank factorization of `G(λ)`.
The number of columns of `R(λ)` is the normal rank `r` of `G(λ)`.
The columns of `R(λ)` form a rational basis of the range (or image) space of the rational matrix `G(λ)`.
A selected set of zeros of `G(λ)` are included as zeros of `R(λ)`.

The resulting named triple `info` contains `(nrank, nfuz, niuz) `, where `info.nrank = r`,
the normal rank of `G(λ)`, `info.nfuz` is the number of finite zeros of `sys` on
the boundary of the stability domain `Cs`, and `info.niuz` is the number of infinite zeros of `sys` in
the continuous-time case and is set to `0` in the discrete-time case.

Depending on the value of the keyword parameter `zeros`, the following options can be selected
for the zeros of `G(λ)` to be included in `R(λ)`:

     "none"       - include no zeros (default)
     "all"        - include all zeros of `sys`
     "unstable"   - include all unstable zeros of `sys`
     "s-unstable" - include all strictly unstable zeros of `sys`, both finite and infinite
     "stable"     - include all stable zeros of `sys`
     "finite"     - include all finite zeros of `sys`
     "infinite"   - include all infinite zeros of `sys`

If `inner = true`, the resulting basis `R(λ)` is _inner_, i.e., `R(λ)'*R(λ) = I`, where `R(s)' = transpose(R(-s))` for a
continuous-time system with `λ = s` and `R(z)' = transpose(R(1/z))` for a discrete-time system with `λ = z`.
This option can be used only in conjunction with `zeros = "none"` or `zeros = "unstable"`.

For a continuous-time system `sys`, the stability domain `Cs` is defined as the set of
complex numbers with real parts at most `-β`,
while for a discrete-time system `sys`, `Cs` is the set of complex numbers with
moduli at most `1-β` (i.e., the interior of a disc of radius `1-β` centered in the origin).
The boundary offset  `β` to be used to assess the stability of zeros and their number
on the boundary of `Cs` can be specified via the keyword parameter `offset = β`.
Accordingly, for a continuous-time system,
the boundary of `Cs` contains the complex numbers with real parts within the interval `[-β,β]`,
while for a discrete-time system, the boundary of `Cs` contains
the complex numbers with moduli within the interval `[1-β,1+β]`.
The default value used for `β` is `sqrt(ϵ)`, where `ϵ` is the working machine precision.

The keyword arguments `atol1`, `atol2` and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of `A`, `B`, `C` and `D`,
the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `E`, `B`, `C` and `D`.
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`, `atol2 = atol`.

For the assessment of zeros, the system pencil `[A-λE B; C D]` is reduced to a
special Kronecker-like form (see [2]). In this reduction, the
performed rank decisions are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

_Method:_  The range computation method is described in [1] and is based on
the reduction algorithm of [2], which has been adapted to deal with
several zero selection options. The computation of the involved
Kronecker-like form is based on the algorithm of [3].

_References:_

[1] Varga, A.
    A note on computing the range of rational matrices.
    arXiv:1707.0048, [https://arxiv.org/abs/1707.0048](https://arxiv.org/abs/1707.004), 2017.

[2] C. Oara.
    Constructive solutions to spectral and inner–outer factorizations
    respect to the disk. Automatica,  41, pp. 1855–1866, 2005.

[3] C. Oara and P. Van Dooren.
    An improved algorithm for the computation of structural invariants of a system pencil and related geometric aspects.
    Syst. Control Lett., 30:39–48, 1997.
"""
function grange(sys::DescriptorStateSpace{T}; zeros = "none", inner::Bool = false,
                atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
                rtol::Real = sys.nx*eps(real(float(one(T))))*iszero(max(atol1,atol2)),
                offset::Real = sqrt(eps(float(real(T)))), fast::Bool = true) where T

   disc = !iszero(sys.Ts)

   # enforce controllability of sys
   sys = gir(sys, atol1 = atol1, atol2 = atol2, obs = false)

   #  Reduce the system matrix pencil to the special Kronecker-like form
   #
   #                 ( Arg-λ*Erg     *        *   *   )
   #                 (   0        Abl-λ*Ebl  Bbl  *   )
   #       At-λ*Et = (   0           0        0   Bn  )
   #                 (--------------------------------)
   #                 (   0          Cbl      Dbl  *   )
   #
   # where the subpencil
   #                           ( Abl-λ*Ebl  Bbl )
   #                           (   Cbl      Ddl )
   #
   # is full column rank, the pair (Abl-λ*Ebl,Bbl) is stabilizable, and
   # Abl-λ*Ebl contains either no zeros of sys if
   # zeros = "none", or contains zeros according to the specified option in zeros.

   At, Et, _, Z, dimsc, nmszer, nizer = gsklf(dssdata(sys)...; disc = disc,
                jobopt = zeros,
                offset = offset, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol,
                withQ = false, withZ = true)

   n, m = size(sys.B); p = size(sys.D,1);

   # form the reduced system (Abl-λ*Ebl,Bbl,Cbl,Dbl)
   nr = dimsc[1]; nric = dimsc[2]; mric = dimsc[3]; nsinf = dimsc[4];
   il = n-nsinf-nric+1:n-nsinf; jal = nr+1:nr+nric; jbl = nr+nric+1:nr+nric+mric; icl = n+1:n+p
   # A = At[il,jal]; E = Et[il,jal]; B = At[il,jbl]; C = At[icl,jal]; D = At[icl,jbl];
   A = view(At,il,jal); E = view(Et,il,jal); B = view(At,il,jbl); C = view(At,icl,jal); D = view(At,icl,jbl)

   if inner && (zeros == "none" || zeros == "unstable")
      # move unstable zeros to stable positions and compress the reduced system to one
      # with full row rank transfer function matrix
      if disc
         # X, _,F, = gared(A,E,B,D'*D,C'*C,C'*D);
         X, _,F, = try
            gared(A,E,B,D'*D,C'*C,C'*D)
         catch err
            findfirst("dichotomic",string(err)) === nothing ? error("$err") :
               error("Solution of the DARE failed: Symplectic pencil has eigenvalues on the unit circle")
         end
         H = cholesky(Hermitian(D'*D+B'*X*B)).U
      else
         # X, _, F, = garec(A,E,B,D'*D,C'*C,C'*D);
         X, _, F, = try
            garec(A,E,B,D'*D,C'*C,C'*D)
         catch err
            findfirst("dichotomic",string(err)) === nothing ? error("$err") :
               error("Solution of the CARE failed: Hamiltonian pencil has jw-axis eigenvalues")
         end
         H = qr(D).R[1:mric,:]
      end
      # construct the inner basis
      sysr = dss(A-B*F, E, B/H, C-D*F, D/H, Ts = sys.Ts)
      # construct the full row rank factor
      CDt = [H*F H]*Z[:,nr+1:n+m-nsinf]'
      sysx = dss(sys.A, sys.E, sys.B, CDt[:,1:n], CDt[:,n+1:n+m], Ts = sys.Ts);
   else
      sysr = dss(A, E, B, C, D, Ts = sys.Ts)
      sysx = dss(sys.A, sys.E, sys.B, Z[1:n,jbl]', Z[n+1:n+m,jbl]', Ts = sys.Ts);
   end
   info = (nrank = mric, nfuz = nmszer, niuz = nizer)

   return sysr, sysx, info

end
"""
    glnull(sys; polynomial = false, simple = false, coinner = false, fast = true, poles = missing, sdeg = missing,
           atol = 0, atol1 = atol, atol2 = atol, rtol, offset = sqrt(ϵ) ) -> (syslnull, info)

Determine for the descriptor systems `sys = (A-λE,B,C,D)` with the `p x m` transfer function matrix `G(λ)`,
the descriptor system `syslnull = (Al-λEl,Bl,Cl,Dl)` with the transfer function matrix `Nl(λ)`
such that `Nl(λ)` is a minimal rational left nullspace basis of `G(λ)` and satisfies `Nl(λ)*G(λ) = 0`.

For the call with

    glnull(sys, m2; polynomial = false, simple = false, coinner = false, fast = true, poles = missing, sdeg = missing,
           atol = 0, atol1 = atol, atol2 = atol, rtol, offset = sqrt(ϵ) ) -> (syslnull, info)

`sys` contains the compound system `sys = [sys1 sys2]`, with `G(λ)`, the transfer function matrix of `sys1`, and
`G2(λ)`, the transfer function matrix of `sys2`, and has the descriptor realization `sys = (A-λE,[B B2],C,[D D2])`,
where `sys2` has `m2` inputs. The resulting `syslnull` contains the compound system
`[syslnull1 syslnull1*sys2] = (Al-λEl,[Bl Bl2],Cr,[Dl Dl2])`, where
`syslnull1 = (Al-λEl,Bl,Cl,Dl)` has the transfer function matrix `Nl(λ)`, which is a rational left nullspace basis of `G(λ)`
satisfying `Nl(λ)*G(λ) = 0` and `syslnull1*sys2 = (Al-λEl,Bl2,Cl,Dl2)` has the transfer function matrix `Nl(λ)*G2(λ)`.

The returned named tuple `info` has the components `info.nrank`, `info.stdim`, `info.degs`, `info.fnorm` and `info.tcond`.

If `polynomial = false`, the resulting `syslnull` has a proper transfer function matrix,
while for `polynomial = true` the resulting `syslnull` has a polynomial transfer function matrix.
The resulting basis `Nl(λ)` contains `p-r` basis vectors, where `r = rank G(λ)`. The rank `r` is returned in `info.nrank`.
If `simple = true`, the resulting basis is _simple_ and satisfies the condition that the sum of the
number of poles of the `p-r` basis vectors is equal to the number of poles of `Nl(λ)` (i.e., its McMillan degree) .

For a non-simple proper basis, the realization `(Al-λEl,Bl,Cl,Dl)` is observable and the pencil `[Al-λEl; Cl]`
is in an observable staircase form. The row dimensions of the full column rank diagonal blocks are returned in `info.stdim`
and the corresponding left Kronecker indices are returned in `info.degs`.
For a simple basis, the regular pencil `Al-λEl` is block diagonal, with the `i`-th block of size `info.stdim[i]`.
The increasing numbers of poles of the basis vectors are returned in `info.degs`.
For the `i`-th basis vector `vi(λ)` (i.e., the `i`-th row of `Nl(λ)`)
a minimal realization can be explicitly constructed as `(Ali-λEli,Bl,Cli,Dl[i,:])`, where `Ali`, `Eli` and `Cli` are
the `i`-th diagonal blocks of `Al`, `El`, and `Cl`, respectively, and `Dl[i,:]` is the `i`-th row of `Dl`.
The corresponding realization of `vi(λ)*G2(λ)` can be constructed as `(Ali-λEli,Bl2,Cl2,Dl2[i,:])`, where
`Dl2[i,:]` is the `i`-th row of `Dl2`.

For a proper basis, the poles of `Nl(λ)` can be freely assigned, by assigning the eigenvalues of the pencil `Al-λEl`.
The vector `poles`, specified as a keyword argument, can be used to specify the desired eigenvalues,
alternatively to or jointly with enforcing a desired stability degree `sdeg` of the real parts of the eigenvalues,
for a continuous-time system, or the moduli of eigenvalues, for a discrete-time system.
If `coinner = true`, the resulting basis `Nl(λ)` is _coinner_, i.e., `Nl(λ)*Nl(λ)' = I`, where `Nl(s)' = transpose(Nl(-s))` for a
continuous-time system with `λ = s` and `Nl(z)' = transpose(Nl(1/z))` for a discrete-time system with `λ = z`.
If the proper basis is simple, each of the resulting individual basis vector is inner.
If `sys2` has poles on the boundary of the appropriate stability domain `Cs`, which are not poles of `sys1` too,
then there exists no inner `Nl(λ)` such that `Nl(λ)*G2(λ)` is stable. An offset can be specified via the keyword parameter `offset = β`
to be used to assess the existence of zeros on the stability domain boundary. Accordingly, for a continuous-time system,
the boundary of `Cs` contains the complex numbers with real parts within the interval `[-β,β]`,
while for a discrete-time system, the boundary of `Cs` contains
the complex numbers with moduli within the interval `[1-β,1+β]`.
The default value used for `β` is `sqrt(ϵ)`, where `ϵ` is the working machine precision.

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

The rank determinations in the performed reductions
 are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.
The computation of simple bases involves the solution of several Type 1 minimum dynamic cover
problems. This computation involves using
non-orthogonal transformations whose worst condition number is returned in `info.tcond`, in conjunction with
using feedback gains, whose norms are returned in `info.fnorm`. High values of these quantities
indicate a potential loss of numerical stability of computations.

_Note:_ The resulting realization of `syslnull` is minimal provided the realization of `sys` is minimal.
However, `syslnull1` is a minimal basis only if the realization (A-lambda E,B,C,D) of `sys1` is
minimal. In this case, `info.degs` are the degrees of the vectors of a minimal polynomial basis or,
if `simple = true`, of the resulting minimal simple proper basis.

_Method:_ The computation method for the computation of a right nullspace basis is applied
to the dual of descriptor system `sys`.
The computation of a minimal proper right nullspace basis is based
on [1]; see also [2]. For the computation of a minimal simple proper
right nullspace basis the method of [3] is emloyed to compute a simple basis from a
minimal proper basis. For the computation of an inner proper right nullspace basis,
the inner factor of an inner-outer factorization of `Nl(λ)` is explicitly
constructed using formulas given in [4].

_References:_

[1] T.G.J. Beelen.
    New algorithms for computing the Kronecker structure of a pencil
    with applications to systems and control theory.
    Ph. D. Thesis, Technical University Eindhoven, 1987.

[2] A. Varga.
    On computing least order fault detectors using rational nullspace bases.
    IFAC SAFEPROCESS'03 Symposium, Washington DC, USA, 2003.

[3] A. Varga.
    On computing nullspace bases – a fault detection perspective.
    Proc. IFAC 2008 World Congress, Seoul, Korea, pages 6295–6300, 2008.

[4] K. Zhou, J. C. Doyle, and K. Glover.
    Robust and Optimal Control. Prentice Hall, 1996.
"""
function glnull(sys::DescriptorStateSpace{T}, m2::Int = 0; polynomial::Bool = false,
                poles::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
                offset::Real = sqrt(eps(float(real(T)))),
                atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
                rtol::Real = sys.nx*eps(real(float(one(T))))*iszero(max(atol1,atol2)),
                fast::Bool = true, coinner::Bool = false, simple::Bool = false) where T
   p, m = size(sys)
   (m2 <= m && m2 >= 0) || throw(DimensionMismatch("m2 must be at most $m, got $m2"))
   polynomial && coinner && (coinner = false; @warn "Coinner option not allowed for polynomial basis")
   sysnull, info1 = grnull(gdual(sys), m2; polynomial = polynomial, poles = poles, sdeg = sdeg, offset = offset,
                         atol = atol, atol1 = atol1, atol2 = atol2, rtol = rtol, fast = fast, inner = coinner, simple = simple)
   info = (nrank = info1.nrank, stdim = reverse(info1.stdim), degs = info1.degs, tcond = info1.tcond, fnorm = info1.fnorm)
   return gdual(sysnull[:,p-info.nrank:-1:1],rev = true), info
end
"""
    grnull(sys; polynomial = false, simple = false, inner = false, fast = true, poles = missing, sdeg = missing,
           atol = 0, atol1 = atol, atol2 = atol, rtol, offset = sqrt(ϵ) ) -> (sysrnull, info)

Determine for the descriptor systems `sys = (A-λE,B,C,D)` with the `p x m` transfer function matrix `G(λ)`,
the descriptor system `sysrnull = (Ar-λEr,Br,Cr,Dr)` with the transfer function matrix `Nr(λ)`
such that `Nr(λ)` is a minimal rational right nullspace basis of `G(λ)` and satisfies `G(λ)*Nr(λ) = 0`.

For the call with

    grnull(sys, p2; polynomial = false, simple = false, inner = false, fast = true, poles = missing, sdeg = missing,
           atol = 0, atol1 = atol, atol2 = atol, rtol, offset = sqrt(ϵ) ) -> (sysrnull, info)

`sys` contains the compound system `sys = [sys1; sys2]`, with `G(λ)`, the transfer function matrix of `sys1`, and
`G2(λ)`, the transfer function matrix of `sys2`, and has the descriptor realization `sys = (A-λE,B,[C;C2],[D;D2])`,
where `sys2` has `p2` outputs. The resulting `sysrnull` contains the compound system
`[sysrnull1; sys2*sysrnull1] = (Ar-λEr,Br,[Cr;Cr2],[Dr;Dr2])`, where
`sysrnull1 = (Ar-λEr,Br,Cr,Dr)` has the transfer function matrix `Nr(λ)`, which is a rational right nullspace basis of `G(λ)`
satisfying `G(λ)*Nr(λ) = 0` and `sys2*sysrnull1 = (Ar-λEr,Br,Cr2,Dr2)` has the transfer function matrix `G2(λ)*Nr(λ)`.

The returned named tuple `info` has the components `info.nrank`, `info.stdim`, `info.degs`, `info.fnorm` and `info.tcond`.

If `polynomial = false`, the resulting `sysrnull` has a proper transfer function matrix,
while for `polynomial = true` the resulting `sysrnull` has a polynomial transfer function matrix.
The resulting basis `Nr(λ)` contains `m-r` basis vectors, where `r = rank G(λ)`. The rank `r` is returned in `info.nrank`.
If `simple = true`, the resulting basis is _simple_ and satisfies the condition that the sum of the
number of poles of the `m-r` basis vectors is equal to the number of poles of `Nr(λ)` (i.e., its McMillan degree) .

For a non-simple proper basis, the realization `(Ar-λEr,Br,Cr,Dr)` is controllable and the pencil `[Br Ar-λEr]`
is in a controllable staircase form. The column dimensions of the full row rank diagonal blocks are returned in `info.stdim`
and the corresponding right Kronecker indices are returned in `info.degs`.
For a simple basis, the regular pencil `Ar-λEr` is block diagonal, with the `i`-th block of size `info.deg[i]`
(the `i`-th right Kronecker index) for a proper basis and  `info.deg[i]+1` for a polynomial basis.
The dimensions of the diagonal blocks are returned in this case in `info.stdim`, while the increasing numbers of poles of
the basis vectors are returned in `info.degs`. For the `i`-th basis vector `vi(λ)` (i.e., the `i`-th column of `Nr(λ)`)
a minimal realization can be explicitly constructed as `(Ari-λEri,Bri,Cr,Dr[:,i])`, where `Ari`, `Eri` and `Bri` are
the `i`-th diagonal blocks of `Ar`, `Er`, and `Br`, respectively, and `Dr[:,i]` is the `i`-th column of `Dr`.
The corresponding realization of `G2(λ)*vi(λ)` can be constructed as `(Ari-λEri,Bri,Cr2,Dr2[:,i])`, where
`Dr2[:,i]` is the `i`-th column of `Dr2`.

For a proper basis, the poles of `Nr(λ)` can be freely assigned, by assigning the  eigenvalues of the pencil `Ar-λEr`.
The vector `poles`, specified as a keyword argument, can be used to specify the desired eigenvalues,
alternatively to or jointly with enforcing a desired stability degree `sdeg` of the real parts of the eigenvalues,
for a continuous-time system, or the moduli of eigenvalues, for a discrete-time system.
If `inner = true`, the resulting basis `Nr(λ)` is _inner_, i.e., `Nr(λ)'*Nr(λ) = I`, where `Nr(s)' = transpose(Nr(-s))` for a
continuous-time system with `λ = s` and `Nr(z)' = transpose(Nr(1/z))` for a discrete-time system with `λ = z`.
If the proper basis is simple, each of the resulting individual basis vector is inner.
If `sys2` has poles on the boundary of the appropriate stability domain `Cs`, which are not poles of `sys1` too,
then there exists no inner `Nr(λ)` such that `G2(λ)*Nr(λ)` is stable. An offset can be specified via the keyword parameter `offset = β`
to be used to assess the existence of zeros on the stability domain boundary. Accordingly, for a continuous-time system,
the boundary of `Cs` contains the complex numbers with real parts within the interval `[-β,β]`,
while for a discrete-time system, the boundary of `Cs` contains
the complex numbers with moduli within the interval `[1-β,1+β]`.
The default value used for `β` is `sqrt(ϵ)`, where `ϵ` is the working machine precision.

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

The rank determinations in the performed reductions
 are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.
The computation of simple bases involves the solution of several Type 1 minimum dynamic cover
problems. This computation involves using
non-orthogonal transformations whose worst condition number is returned in `info.tcond`, in conjunction with
using feedback gains, whose norms are returned in `info.fnorm`. High values of these quantities
indicate a potential loss of numerical stability of computations.

_Note:_ The resulting realization of `sysrnull` is minimal provided the realization of `sys` is minimal.
However, `sysrnull1` is a minimal basis only if the realization (A-lambda E,B,C,D) of `sys1` is
minimal. In this case, `info.degs` are the degrees of the vectors of a minimal polynomial basis or,
if `simple = true`, of the resulting minimal simple proper basis.

_Method:_ The computation of a minimal proper right nullspace basis is based
on [1]; see also [2]. For the computation of a minimal simple proper
right nullspace basis the method of [3] is emloyed to compute a simple basis from a
minimal proper basis. For the computation of an inner proper right nullspace basis,
the inner factor of an inner-outer factorization of `Nr(λ)` is explicitly
constructed using formulas given in [4].

_References:_

[1] T.G.J. Beelen.
    New algorithms for computing the Kronecker structure of a pencil
    with applications to systems and control theory.
    Ph. D. Thesis, Technical University Eindhoven, 1987.

[2] A. Varga.
    On computing least order fault detectors using rational nullspace bases.
    IFAC SAFEPROCESS'03 Symposium, Washington DC, USA, 2003.

[3] A. Varga.
    On computing nullspace bases – a fault detection perspective.
    Proc. IFAC 2008 World Congress, Seoul, Korea, pages 6295–6300, 2008.

[4] K. Zhou, J. C. Doyle, and K. Glover.
    Robust and Optimal Control. Prentice Hall, 1996.
"""
function grnull(sys::DescriptorStateSpace{T}, p2::Int = 0; polynomial::Bool = false,
                poles::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
                offset::Real = sqrt(eps(float(real(T)))),
                atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
                rtol::Real = ((size(sys.A,1)+1)*eps(real(float(one(T)))))*iszero(max(atol1,atol2)),
                fast::Bool = true, inner::Bool = false, simple::Bool = false) where T

   p, m = size(sys)
   (p2 <= p && p2 >= 0) || throw(DimensionMismatch("p2 must be at most $p, got $p2"))
   Ts = sys.Ts;
   polynomial && inner && (inner = false; @warn "Inner option not allowed for polynomial basis")

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

   disc = !iszero(Ts)
   sdeg_nomissing = !ismissing(sdeg)
   poles_nomissing = !ismissing(poles)

   stabilize = sdeg_nomissing || poles_nomissing

   # set default values of sdeg if poles = missing
   # sdeg_nomissing || (sdeg = disc ? real(T1)(0.95) : real(T1)(-0.05))

   # sort desired eigenvalues
   if poles_nomissing
      tempc = poles[imag.(poles) .> 0]
      if !isempty(tempc)
         tempc1 = conj(poles[imag.(poles) .< 0])
         isequal(tempc[sortperm(real(tempc))],tempc1[sortperm(real(tempc1))]) ||
                 error("poles must be a self-conjugated complex vector")
      end
      # check that all eigenvalues are inside of the stability region
      sdeg_nomissing && ( ((disc && any(abs.(poles) .> sdeg) )  || (!disc && any(real.(poles) .> sdeg)))  &&
            error("The elements of poles must lie in the stability region of interest") )
   end

   # maximum allowed condition number
   maxtcond = 1.e4;

   At, Et, Bt, Ct, Dt = dssdata(T1,sys);

   n = size(At,1)
   p1 = p-p2
   ip1 = 1:p1; ip2 = p1+1:p;

   M, N, Q, Z, νr, μr, = sklf_right(At, Et, Bt, Ct[ip1,:], Dt[ip1,:]; fast = true, atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = true, withZ = true)

   # separate the left Kronecker structure using the Kronecker-like form of
   # the system pencil S(λ) = [ At-λEt Bt; Ct1 Dt1]:
   #
   #        Q'*S(λ)*Z = [ Br Ar - λ Er      *     ]
   #                    [ 0     0       Al - λ El ]
   #
   # A = [a b ; c(ip1,:) d(ip1,:)]; E = [e zeros(n,m);zeros(p1,n+m)];
   # mode = 2; # perform separation of infinite part
   # qz = 2;   # accumulate only Z
   # [A,E,~,Z,nur,mur] = sl_klf(A,E,tolc,mode,qz);

   # determine main dimensional parameters
   nr = sum(νr);    # order of left structure (also order of Ar)
   mr = sum(μr)-nr; # number of basis elements
   nrank = m - mr;   # normal rank of the transfer matrix SYS1

   tcond = 1;     # set condition number of used transformations
   fnorm = 0;

   degs = zeros(Int,mr)
   if nrank == m
      # full row rank: the null space is empty
      sysrnull = dss(zeros(T,m+p2,0), Ts = Ts)
      info = (nrank = nrank, stdim = zeros(Int,0), degs = degs, tcond = 1, fnorm = 0)
      return sysrnull, info
   end
   # in the case nrank = 0, a null space basis is the identity matrix
   # and therefore sysrnull = [eye(p); sys(ip2,:)] and degs = zeros(1,m);
   # however, to handle properly this case, the right Kronecker indices
   # must be determined in degs; the basis is generally not minimal

   # compute the right Kronecker indices, which are also the degrees of
   # the vectors of a polynomial right nullspace basis
   nvi = μr-νr;    # there are nvi(i) vectors of degree i-1
   nb = length(νr); # number of blocks in [ Br Ar ]
   k = 1;
   for i = 1:nb
       for j = 1:nvi[i]
           degs[k] = i-1
           k += 1
       end
   end

   # normalize right part
   klf_right_refineut!(νr, μr, M, N, Q, Z, withQ = false, withZ = true)

   # form  CR1 = [0 I]*Z = [ Dr1 Cr1 * ] and
   #       CR2 = [ C2 D2]*Z = [ Dr2 Cr2 * ]
   CR1 =  Z[n+1:n+m,:]
   CR2 = [ Ct[ip2,:] Dt[ip2,:]]*Z;

   iar = 1:nr; jar = mr+1:mr+nr; jbr = 1:mr;

   ar = view(M,iar,jar); cr = [ CR1[:,jar]; CR2[:,jar] ];  er = view(N,iar,jar);
   br = view(M,iar,jbr); dr = [ CR1[:,jbr]; CR2[:,jbr] ];
   #sklf_rightnorm!(νr,ar,er,br,cr,dr)
   sysrnull = dss(ar,er,br,cr,dr)
   if maximum(degs) == 0
      info = (nrank = nrank, stdim = simple ? zeros(Int,mr) : νr[νr .> 0], degs = degs, tcond = 1, fnorm = 0)
      return dss(ar,er,br,cr,dr,Ts = Ts), info
   end

   if inner && p2 > 0
      # check stabilizability of SYSF*NR1
      zer = gzero(dss(ar,er,br,cr[1:m,:],dr[1:m,:],Ts=Ts), atol1 = atol1,atol2 = atol2, rtol = rtol)
      if disc
         zf = abs.(zer[isfinite.(zer)])
         nostab = any((zf .>= 1-offset) .& (zf .<= 1+offset))
      else
         zf = real.(zer[isfinite.(zer)])
         nostab = any(isinf.(zer)) || any((zf .>= -offset) .& (zf .<= offset))
      end
      if nostab
         mes = "No inner nullspace Nr exists such that SYSF*Nr is stable - \n  standard nullspace computed instead"
         @warn mes
         inner = false;  # perform standard nullspace computation
      end
   end

   if !simple || mr == 1
      if polynomial
         # pr = size(dr,1)
         f, g, SF, = salocinf(ar, er, br; fast = fast, atol1 = atol1, atol2 = atol2, atol3 = atol1)
         # build (ar+br*f-λ(er+br*g), br, cr+dr*f+λdr*g, dr)
         # ar1 = [I zeros(T1,pr,nr); zeros(T1,nr,pr) SF.S]
         # er1 = [zeros(T1,pr,pr) dr*g*SF.Z; zeros(T1,nr,pr) SF.T]
         # br1 = [zeros(T1,pr,mr); SF.Q'*br]
         # cr1 = [-I (cr+dr*f)*SF.Z]
         # sysrnull = dss(ar1,er1,br1,cr1,dr,Ts = Ts)
         sysrnull = dss(SF.S, SF.T, SF.Q'*br, missing, (cr+dr*f)*SF.Z, dr*g*SF.Z, dr, missing, Ts = Ts,
                        compacted = true, atol1 = atol1, atol2 = atol2, rtol = rtol)
         fnorm = max(fnorm, norm(g))
         simple &&  (blkdims = [order(sysrnull)])
      else
         if inner
            # compute the outer factor Nro of the inner-outer factorization
            # NR1 = Nri*Nro and form explicitly [ NR1; SYS2*NR1]*inv(Nro)
            C = cr[1:m,:]; D = dr[1:m,:];
            if disc
               X, _,f, = try
                  gared(ar,er,br,D'*D,C'*C,C'*D)
               catch err
                  findfirst("dichotomic",string(err)) === nothing ? error("$err") :
                     error("Solution of the DARE failed: Symplectic pencil has eigenvalues on the unit circle")
               end
               H = cholesky(Hermitian(D'*D+br'*X*br)).U
            else
               X, _, f, = try
                  garec(ar,er,br,D'*D,C'*C,C'*D)
               catch err
                  findfirst("dichotomic",string(err)) === nothing ? error("$err") :
                     error("Solution of the CARE failed: Hamiltonian pencil has jw-axis eigenvalues")
               end
               H = qr(D).R[1:mr,:];
            end
            fnorm = norm(f)
            sysrnull = dss(ar-br*f,er,br/H,cr-dr*f,dr/H,Ts = Ts)
         else
            fnorm = 0
            if stabilize
               # make poles stable
               if nr > 0
                  f, SF, = saloc(ar, er, br; evals = poles, sdeg = sdeg, disc = disc, fast = fast, atol1 = atol1, atol2 = atol2, atol3 = atol1)
                  sysrnull = dss(SF.S, SF.T, SF.Q'*br, (cr+dr*f)*SF.Z, dr, Ts = Ts)
                  # mul!(ar, br, f, ONE, ONE)
                  # mul!(cr, dr, f, ONE, ONE)
                  fnorm = norm(f)
               else
                  sysrnull = dss(ar,er,br,cr,dr,Ts = Ts)
               end
            else
               sysrnull = dss(ar,er,br,cr,dr,Ts = Ts)
            end
         end
         simple &&  (blkdims = [size(ar,1)])
      end
   else
      # compute a simple basis
      sysrnull = dss(ar,er,br,cr,dr,Ts = Ts)
      sysn = dss(zeros(T1,m+p2,0),Ts = Ts)
      tcond = 1;
      blkdims = zeros(Int,mr)
      for i = 1:mr
          imt = [[i]; Vector(1:i-1); Vector(i+1:mr)];  # row permutation indices
          sysni, _, info1 = grmcover1(sysrnull[:,imt],1,atol1 = atol1, atol2 = atol2, rtol = rtol, fast = fast)
          ndegi = order(sysni)
          ndegi == degs[i] ||
               (@warn "Resulting McMillan degree of $i-th vector is $ndegi, expected $(degs[i])")
          tcond = max(tcond, info1.fnorm, info1.tcond)
          if polynomial
             if ndegi > 0
               a1, e1, b1, c1, d1 = dssdata(sysni)
               f, g, SF, = salocinf(a1, e1, b1; fast = fast, atol1 = atol1, atol2 = atol2, atol3 = atol1)
               sysni = dss(SF.S, SF.T, SF.Q'*b1, missing, (c1+d1*f)*SF.Z, d1*g*SF.Z, d1, missing, Ts = Ts, compacted = true,
                           atol1 = atol1, atol2 = atol2, rtol = rtol)
               # f, g, SF, = salocinf(sysni.A, sysni.E, sysni.B; fast = fast, atol1 = atol1, atol2 = atol2, atol3 = atol1)
               # sysni = dss(SF.S, SF.T, SF.Q'*sysni.B, missing, (sysni.C+sysni.D*f)*SF.Z, sysni.D*g*SF.Z, sysni.D, missing, Ts = Ts)
               #  den = fromroots(eigvals(sysni.A,sysni.E))
               #  T <: Complex || (den = real(den))
               #  sysni = gir(sysni*dss(den,Ts=Ts),atol1=atol1, atol2=atol2, rtol = rtol, fast = fast, infinite = false)
               ndegi = order(sysni)
             end
          else
             ar,er,br,cr,dr = dssdata(sysni)
             if inner
                C = cr[1:m,:]; D = dr[1:m,:];
                if disc
                   X, _,f, = try
                     gared(ar,er,br,D'*D,C'*C,C'*D)
                   catch err
                     findfirst("dichotomic",string(err)) === nothing ? error("$err") :
                        error("Solution of the DARE failed: Symplectic pencil has eigenvalues on the unit circle")
                   end
                   H = cholesky(Hermitian(D'*D+br'*X*br)).U
                else
                   X, _, f, = try
                     garec(ar,er,br,D'*D,C'*C,C'*D)
                   catch err
                     findfirst("dichotomic",string(err)) === nothing ? error("$err") :
                        error("Solution of the CARE failed: Hamiltonian pencil has jw-axis eigenvalues")
                   end
                   H = qr(D).R[1:1,1:1]
                 end
                 fnorm = max(fnorm, norm(f))
                 sysni = dss(ar-br*f,er,br/H,cr-dr*f,dr/H,Ts = Ts)
              else
                 if stabilize
                    f, SF, = saloc(ar, er, br; evals = poles, sdeg = sdeg, disc = disc, fast = fast, atol1 = atol1, atol2 = atol2, atol3 = atol1)
                    sysni = dss(SF.S, SF.T, SF.Q'*br, (cr+dr*f)*SF.Z, dr, Ts = Ts)
                    #sysni.A[:,:] = ar + br*f; sysni.C[:,:] = sysni.C + sysni.D*f
                    fnorm = max(fnorm, norm(f))
                 end
              end
          end
          blkdims[i] = ndegi
          sysn = [sysn sysni]
      end
      tcond > maxtcond && (@warn "Possible loss of numerical stability due to ill-conditioned transformations")
      sysrnull = sysn
   end

   info = (nrank = nrank, stdim = simple ? blkdims : νr[νr .> 0], degs = degs, tcond = tcond, fnorm = fnorm)
   return sysrnull, info

   # end GRNULL

end


text:
### Fetch required packages and reaction networks ###
using Catalyst, OrdinaryDiffEq, Random, Test
using ModelingToolkit: get_states, get_ps
include("test_networks.jl")

using StableRNGs
rng = StableRNG(12345)

### Tests various ways to input u0 and p for various functions ###

# Tests for reaction_networks_standard[7]
test_network = reaction_networks_standard[7]
@parameters t p1 p2 p3 k1 k2 k3 v1 K1 d1 d2 d3 d4 d5
@variables X1(t) X2(t) X3(t) X4(t) X5(t) X6(t) X(t)

for factor = [1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3]
    u0_1 = factor*rand(rng,length(get_states(test_network)))
    u0_2 = [X1=>u0_1[1], X2=>u0_1[2], X3=>u0_1[3], X4=>u0_1[4], X5=>u0_1[5]]
    u0_3 = [X2=>u0_1[2], X5=>u0_1[5], X4=>u0_1[4], X3=>u0_1[3], X1=>u0_1[1]]
    p_1 = 0.01 .+ factor*rand(rng,length(get_ps(test_network)))
    p_2 = [p1=>p_1[1], p2=>p_1[2], p3=>p_1[3], k1=>p_1[4], k2=>p_1[5], k3=>p_1[6], v1=>p_1[7], K1=>p_1[8], d1=>p_1[9], d2=>p_1[10], d3=>p_1[11], d4=>p_1[12], d5=>p_1[13]]
    p_3 = [k2=>p_1[5], k3=>p_1[6], v1=>p_1[7], d5=>p_1[13], p2=>p_1[2], p1=>p_1[1], d2=>p_1[10], K1=>p_1[8], d1=>p_1[9], d4=>p_1[12], d3=>p_1[11], p3=>p_1[3], k1=>p_1[4]]

    sols = []
    push!(sols,solve(ODEProblem(test_network,u0_1,(0.,10.),p_1),Rosenbrock23()))
    push!(sols,solve(ODEProblem(test_network,u0_1,(0.,10.),p_2),Rosenbrock23()))
    push!(sols,solve(ODEProblem(test_network,u0_1,(0.,10.),p_3),Rosenbrock23()))
    push!(sols,solve(ODEProblem(test_network,u0_2,(0.,10.),p_1),Rosenbrock23()))
    push!(sols,solve(ODEProblem(test_network,u0_2,(0.,10.),p_2),Rosenbrock23()))
    push!(sols,solve(ODEProblem(test_network,u0_2,(0.,10.),p_3),Rosenbrock23()))
    push!(sols,solve(ODEProblem(test_network,u0_3,(0.,10.),p_1),Rosenbrock23()))
    push!(sols,solve(ODEProblem(test_network,u0_3,(0.,10.),p_2),Rosenbrock23()))
    push!(sols,solve(ODEProblem(test_network,u0_3,(0.,10.),p_3),Rosenbrock23()))

    ends = map(sol -> sol.u[end],sols)
    for i in 1:length(u0_1)
        @test abs(maximum(getindex.(ends,1))-minimum(first.(getindex.(ends,1)))) < 1e-5
    end

    u0_1 = rand(rng,1:Int64(factor*100),length(get_states(test_network)))
    u0_2 = [X1=>u0_1[1], X2=>u0_1[2], X3=>u0_1[3], X4=>u0_1[4], X5=>u0_1[5]]
    u0_3 = [X2=>u0_1[2], X5=>u0_1[5], X4=>u0_1[4], X3=>u0_1[3], X1=>u0_1[1]]

    discrete_probs = []
    push!(discrete_probs,DiscreteProblem(test_network,u0_1,(0.,1.),p_1))
    push!(discrete_probs,DiscreteProblem(test_network,u0_1,(0.,1.),p_2))
    push!(discrete_probs,DiscreteProblem(test_network,u0_1,(0.,1.),p_3))
    push!(discrete_probs,DiscreteProblem(test_network,u0_2,(0.,1.),p_1))
    push!(discrete_probs,DiscreteProblem(test_network,u0_2,(0.,1.),p_2))
    push!(discrete_probs,DiscreteProblem(test_network,u0_2,(0.,1.),p_3))
    push!(discrete_probs,DiscreteProblem(test_network,u0_3,(0.,1.),p_1))
    push!(discrete_probs,DiscreteProblem(test_network,u0_3,(0.,1.),p_2))
    push!(discrete_probs,DiscreteProblem(test_network,u0_3,(0.,1.),p_3))

    for i in 2:9
        @test discrete_probs[1].p == discrete_probs[i].p
        @test discrete_probs[1].u0 == discrete_probs[i].u0
    end
end

ext: jl


text:
module FileSession

using Genie
import Serialization, Logging

SESSIONS_PATH = "sessions"

function sessions_path(path::String)
  SESSIONS_PATH = normpath(path) |> abspath
end
function sessions_path()
  SESSIONS_PATH
end

"""
    write(session::Genie.Sessions.Session) :: Genie.Sessions.Session

Persists the `Session` object to the file system, using the configured sessions folder and returns it.
"""
function write(session::Genie.Sessions.Session) :: Genie.Sessions.Session
  if ! isdir(sessions_path())
    @warn "Sessions folder $(sessions_path()) does not exist"
    @info "Creating sessions folder at $(sessions_path())"

    try
      mkpath(sessions_path())
    catch ex
      @error "Can't create session storage path $(sessions_path())"
      @error ex
    end
  end

  try
    write_session(session)

    return session
  catch ex
    @error "Failed to store session data"
    @error ex
  end

  try
    @warn "Resetting session"

    session = Genie.Sessions.Session(Genie.Sessions.id())
    Genie.Cookies.set!(Genie.Router.params(Genie.PARAMS_RESPONSE_KEY), Genie.config.session_key_name, session.id, Genie.config.session_options)
    write_session(session)
    Genie.Router.params(Genie.PARAMS_SESSION_KEY, session)

    return session
  catch ex
    @error "Failed to regenerate and store session data. Giving up."
    @error ex
  end

  session
end

function write_session(session::Genie.Sessions.Session)
  open(joinpath(SESSIONS_PATH, session.id), "w") do io
    Serialization.serialize(io, session)
  end
end

"""
    read(session_id::Union{String,Symbol}) :: Union{Nothing,Genie.Sessions.Session}
    read(session::Genie.Sessions.Session) :: Union{Nothing,Genie.Sessions.Session}

Attempts to read from file the session object serialized as `session_id`.
"""
function read(session_id::Union{String,Symbol}) :: Union{Nothing,Genie.Sessions.Session}
  isfile(joinpath(SESSIONS_PATH, session_id)) || return nothing

  try
    open(joinpath(SESSIONS_PATH, session_id), "r") do (io)
      Serialization.deserialize(io)
    end
  catch ex
    @error "Can't read session"
    @error ex
  end
end

function read(session::Genie.Sessions.Session) :: Union{Nothing,Genie.Sessions.Session}
  read(session.id)
end

#===#
# IMPLEMENTATION

"""
    persist(s::Session) :: Session

Generic method for persisting session data - delegates to the underlying `SessionAdapter`.
"""
function Genie.Sessions.persist(req::Genie.Sessions.HTTP.Request, res::Genie.Sessions.HTTP.Response, params::Dict{Symbol,Any}) :: Tuple{Genie.Sessions.HTTP.Request,Genie.Sessions.HTTP.Response,Dict{Symbol,Any}}
  write(params[Genie.PARAMS_SESSION_KEY])

  req, res, params
end
function Genie.Sessions.persist(s::Genie.Sessions.Session) :: Genie.Sessions.Session
  write(s)
end

"""
    load(session_id::String) :: Session

Loads session data from persistent storage.
"""
function Genie.Sessions.load(session_id::String) :: Genie.Sessions.Session
  session = read(session_id)

  session === nothing ? Genie.Sessions.Session(session_id) : (session)
end

end


text:
# cfg.all_joints = [      [0, 5], [1, 4], [2, 3], [6, 11], [7, 10],  [8, 9],     [12],      [13]]
# cfg.all_joints_names = ['ankle', 'knee', 'hip', 'wrist', 'elbow', 'shoulder', 'chin', 'forehead']
# cfg.num_joints = 14

# TODO: add this to config
path_to_processed_mat = "/userfiles/gsoykan20/mpii_human_pose/cropped/dataset.mat"
path_to_single_person_mat = "/userfiles/gsoykan20/mpii_human_pose/cropped/annolist-singlePerson-h400.mat"
path_to_full_mat = "/userfiles/gsoykan20/mpii_human_pose/cropped/annolist-full-h400.mat"
path_to_multi_person_mat = "/userfiles/gsoykan20/mpii_human_pose/cropped/annolist-multPerson-h400.mat"
global_num_joints = 14
# all_joints = [[1, 6], [2, 5], [3, 4], [7, 12], [8, 11], [9, 10], [13], [14]]
symmetric_joint_dict = Dict(
    1 => 6,
    6 => 1,
    2 => 5,
    5 => 2,
    3 => 4,
    4 => 3,
    7 => 12,
    12 => 7,
    8 => 11,
    11 => 8,
    9 => 10,
    10 => 9,
    13 => 13,
    14 => 14 )

# Total image count => 28883
# multi-person image count => 9698
# single-person image count => 19185 ( 17440 - 1744 )
test_image_count = 1000
validation_image_count = 64 #28883 - 25600 - 1000
train_image_count =  1024 # 25600
read_image_w = 512
read_image_h = 512
max_image_number_to_read = validation_image_count + train_image_count + test_image_count

global_scale = 0.8452830189
scale_jitter_interval = [0.85, 1.15]
preprocess_stride = 8
pos_dist_thresh = 17

weight_decay = 0.0001
momentum_gamma = 0.9

# TODO: output consider threshold might not be needed from now on
output_consider_threshold = 0.0
# TODO: learn how this was computed
global_locref_stdev = 7.2801
locref_loss_weight = 0.05
PCKh_range=0.5
whole_dataset_count = 28883
single_dataset_count = 19185
mean_pixel = [123.68 / 255, 116.779 / 255, 103.939 / 255]
reshaped_mean_pixel = reshape(mean_pixel, (1, 1, 3, 1));

mean_pixel_255 = [123.68, 116.779, 103.939]
reshaped_mean_pixel_255 = reshape(mean_pixel_255, (1, 1, 3, 1));

use_locref_mask_weights = false
global_add_random_mirroring = false

pre_full_path = "/kuacc/users/gsoykan20/comp541_term_project/deeper-cut/results/"


text:
# export container object for GPU temporary arrays
export PlinkGPUVariables

type PlinkGPUVariables{T <: Float, V <: cl.Buffer}
    df_buff     :: V
    x_buff      :: cl.Buffer{Int8}
    y_buff      :: V
    mask_buff   :: cl.Buffer{Int}
    device      :: cl.Device
    ctx         :: cl.Context
    queue       :: cl.CmdQueue
    m_buff      :: V
    p_buff      :: V
    red_buff    :: V
    xtyk        :: cl.Kernel
    rxtyk       :: cl.Kernel
    reset_x     :: cl.Kernel
    wg_size     :: Int
    y_chunks    :: Int
    r_chunks    :: Int
    n           :: Int
    p           :: Int
    p2          :: Int
    n32         :: Int32
    p32         :: Int32
    y_chunks32  :: Int32
    blocksize32 :: Int32
    wg_size32   :: Int32
    y_blocks32  :: Int32
    r_length32  :: Int32
    genofloat   :: cl.LocalMem{T}
end

function PlinkGPUVariables{T <: Float}(
    df_buff     :: cl.Buffer{T},
    x_buff      :: cl.Buffer{Int8},
    y_buff      :: cl.Buffer{T},
    mask_buff   :: cl.Buffer{Int},
    device      :: cl.Device,
    ctx         :: cl.Context,
    queue       :: cl.CmdQueue,
    m_buff      :: cl.Buffer{T},
    p_buff      :: cl.Buffer{T},
    red_buff    :: cl.Buffer{T},
    xtyk        :: cl.Kernel,
    rxtyk       :: cl.Kernel,
    reset_x     :: cl.Kernel,
    wg_size     :: Int,
    y_chunks    :: Int,
    r_chunks    :: Int,
    n           :: Int,
    p           :: Int,
    p2          :: Int,
    n32         :: Int32,
    p32         :: Int32,
    y_chunks32  :: Int32,
    blocksize32 :: Int32,
    wg_size32   :: Int32,
    y_blocks32  :: Int32,
    r_length32  :: Int32,
    genofloat   :: cl.LocalMem{T}
)
    PlinkGPUVariables{T, eltype(y_buff)}(df_buff, x_buff, y_buff, mask_buff, device, ctx, queue, m_buff, p_buff, red_buff, xtyk, rxtyk, reset_x, wg_size, y_chunks, r_chunks, n, p, p2, n32, p32, y_chunks32, blocksize32, wg_size32, y_blocks32, r_length32, genofloat)
end

function PlinkGPUVariables{T <: Float}(
    z      :: DenseVector{T},
    x      :: BEDFile{T},
    y      :: DenseVector{T},
    kern   :: String = PLINK.gpucode64,
    mask_n :: Vector{Int} = ones(Int, size(y))
)
    n,p         = size(x.geno)
    wg_size     = 512
    y_chunks    = div(n, wg_size) + (n % wg_size != 0 ? 1 : 0)
    y_blocks    = div(y_chunks, wg_size) + (y_chunks % wg_size != 0 ? 1 : 0)
    r_chunks    = div(p*y_chunks, wg_size) + ((p*y_chunks) % wg_size != 0 ? 1 : 0)
    wg_size32   = convert(Int32, wg_size)
    n32         = convert(Int32, n)
    p32         = convert(Int32, p)
    y_chunks32  = convert(Int32, y_chunks)
    y_blocks32  = convert(Int32, y_blocks)
    blocksize32 = convert(Int32, x.geno.blocksize)
    r_length32  = convert(Int32, p*y_chunks)
    device      = last(cl.devices(:gpu))
    ctx         = cl.Context(device) :: cl.Context
    queue       = cl.CmdQueue(ctx) :: cl.CmdQueue
    program     = cl.Program(ctx, source=kern) :: cl.Program
    cl.build!(program)
    xtyk        = cl.Kernel(program, "compute_xt_times_vector")
    rxtyk       = cl.Kernel(program, "reduce_xt_vec_chunks")
    reset_x     = cl.Kernel(program, "reset_x")
    x_buff      = cl.Buffer(Int8, ctx, (:r,  :copy), hostbuf = sdata(x.geno.x)) :: cl.Buffer{Int8}
    y_buff      = cl.Buffer(T,    ctx, (:r,  :copy), hostbuf = sdata(y)) :: cl.Buffer{T}
    m_buff      = cl.Buffer(T,    ctx, (:r,  :copy), hostbuf = sdata(x.means)) :: cl.Buffer{T}
    p_buff      = cl.Buffer(T,    ctx, (:r,  :copy), hostbuf = sdata(x.precs)) :: cl.Buffer{T}
    df_buff     = cl.Buffer(T,    ctx, (:rw, :copy), hostbuf = sdata(z)) :: cl.Buffer{T}
    red_buff    = cl.Buffer(T,    ctx, (:rw), p*y_chunks) :: cl.Buffer{T}
    mask_buff   = cl.Buffer(Int,  ctx, (:r,  :copy), hostbuf = sdata(mask_n)) :: cl.Buffer{Int}
    genofloat   = cl.LocalMem(T, wg_size)

    PlinkGPUVariables{T, cl.Buffer{T}}(df_buff, x_buff, y_buff, mask_buff, device, ctx, queue, m_buff, p_buff, red_buff, xtyk, rxtyk, reset_x, wg_size, y_chunks, r_chunks, x.geno.n, x.geno.p, x.covar.p, n32, p32, y_chunks32, blocksize32, wg_size32, y_blocks32, r_length32, genofloat)
end

"""
    df_x2!(snp, df, x, r, mask_n)

This subroutine calculates the gradient for one nongenetic covariate.

Arguments:

- `snp` is the SNP to use in calculations.
- `df` is the gradient vector.
- `x` is the `BEDFile` object with nongenetic covariates
- `r` is the vector of residuals.
- `mask_n` is the bitmask on the data.
"""
function df_x2!{T <: Float}(
    snp     :: Int,
    df      :: DenseVector{T},
    x       :: BEDFile{T},
    y       :: DenseVector{T},
    mask_n  :: Vector{Int}
)
    m = x.means[x.geno.p+snp]
    s = x.precs[x.geno.p+snp]
    df[x.geno.p+snp] = zero(T)
    @inbounds for case = 1:x.geno.n
        if mask_n[case] == 1
            @fastmath df[x.geno.p+snp] += y[case] * (x.covar.x[case,snp] - m) * s
        end
    end
    return nothing
end

function df_x2!{T <: Float}(
    df      :: DenseVector{T},
    x       :: BEDFile{T},
    y       :: DenseVector{T},
    mask_n  :: Vector{Int}
)
    for i = 1:x.covar.p
        df_x2!(i, df, x, y, mask_n)
    end
    return nothing
end

function Base.At_mul_B!{T <: Float}(
    df          :: SharedVector{T},
    df_buff     :: cl.Buffer,
    x           :: BEDFile{T},
    x_buff      :: cl.Buffer,
    y           :: SharedVector{T},
    y_buff      :: cl.Buffer,
    mask_n      :: DenseVector{Int},
    mask_buff   :: cl.Buffer,
    q           :: cl.CmdQueue,
    m_buff      :: cl.Buffer,
    p_buff      :: cl.Buffer,
    red_buff    :: cl.Buffer,
    xtyk        :: cl.Kernel,
    rxtyk       :: cl.Kernel,
    reset_x     :: cl.Kernel,
    wg_size     :: Int,
    y_chunks    :: Int,
    r_chunks    :: Int,
    n           :: Int,
    p           :: Int,
    p2          :: Int,
    n32         :: Int32,
    p32         :: Int32,
    y_chunks32  :: Int32,
    blocksize32 :: Int32,
    wg_size32   :: Int32,
    y_blocks32  :: Int32,
    r_length32  :: Int32,
    genofloat   :: cl.LocalMem
)
    cl.wait(cl.copy!(q, y_buff, sdata(y)))
    cl.wait(queue(reset_x, (wg_size*r_chunks,1,1), nothing, red_buff, r_length32, zero(T)))
    cl.wait(queue(xtyk, (wg_size*y_chunks,p,1), (wg_size,1,1), n32, p32, y_chunks32, blocksize32, wg_size32, x_buff, red_buff, y_buff, m_buff, p_buff, mask_buff, genofloat))
    cl.wait(queue(rxtyk, (wg_size,p,1), (wg_size,1,1), n32, y_chunks32, y_blocks32, wg_size32, red_buff, df_buff, genofloat))
    cl.wait(cl.copy!(q, sdata(df), df_buff))
    df_x2!(df, x, y, mask_n)
    return nothing
end

# wrapper functions for At_mul_B!
function copy_y!{T <: Float}(v::PlinkGPUVariables{T}, y::SharedVector{T})
    cl.copy!(v.queue, v.y_buff, sdata(y)) :: cl.NannyEvent
    #cl.queue(cl.copy!, v.y_buff, sdata(y)) :: cl.NannyEvent
end

function reset_x!{T <: Float}(v::PlinkGPUVariables{T})
    #cl.call(v.queue, v.reset_x, (v.wg_size*v.r_chunks,1,1), nothing, v.red_buff, v.r_length32, zero(T)) :: cl.Event
    v.queue(v.reset_x, (v.wg_size*v.r_chunks,1,1), nothing, v.red_buff, v.r_length32, zero(T)) :: cl.Event
end

function xty!{T <: Float}(v::PlinkGPUVariables{T})
    #cl.call(v.queue, v.xtyk, (v.wg_size*v.y_chunks, v.p, 1), (v.wg_size, 1, 1), v.n32, v.p32, v.y_chunks32, v.blocksize32, v.wg_size32, v.x_buff, v.red_buff, v.y_buff, v.m_buff, v.p_buff, v.mask_buff, v.genofloat) :: cl.Event
    v.queue(v.xtyk, (v.wg_size*v.y_chunks, v.p, 1), (v.wg_size, 1, 1), v.n32, v.p32, v.y_chunks32, v.blocksize32, v.wg_size32, v.x_buff, v.red_buff, v.y_buff, v.m_buff, v.p_buff, v.mask_buff, v.genofloat) :: cl.Event
end

function xty_reduce!{T <: Float}(v::PlinkGPUVariables{T})
    #cl.call(v.queue, v.rxtyk, (v.wg_size,v.p,1), (v.wg_size,1,1), v.n32, v.y_chunks32, v.y_blocks32, v.wg_size32, v.red_buff, v.df_buff, v.genofloat) :: cl.Event
    v.queue(v.rxtyk, (v.wg_size,v.p,1), (v.wg_size,1,1), v.n32, v.y_chunks32, v.y_blocks32, v.wg_size32, v.red_buff, v.df_buff, v.genofloat) :: cl.Event
end

function copy_xty!{T <: Float}(xty::DenseVector{T}, v::PlinkGPUVariables{T})
    cl.copy!(v.queue, sdata(xty), v.df_buff) :: cl.NannyEvent
    #cl.queue(cl.copy!, sdata(xty), v.df_buff) :: cl.NannyEvent
end

"""
    At_mul_B!(xty, x::BEDFile, y, mask_n, v)

Compute `x' * y` for a `BEDFile` object `x` using a GPU command queue configured in the `PlinkGPUVariables` object `v`.
"""
function Base.At_mul_B!{T <: Float}(
    xty    :: SharedVector{T},
    x      :: BEDFile{T},
    y      :: SharedVector{T},
    mask_n :: Vector{Int},
    v      :: PlinkGPUVariables{T}
)
    copy_y!(v, y)
#   cl.fill!(queue, red_buff, zero(T))    # only works with OpenCL 1.2
    reset_x!(v)
    xty!(v)
    xty_reduce!(v)
    copy_xty!(xty, v)
    df_x2!(xty, x, y, mask_n)
    return nothing
end

"""
    At_mul_B(x::BEDFile, y, kernfile, mask_n)

If called with a kernel file, then `At_mul_B()` will attempt to accelerate computations with a GPU.
"""
function Base.At_mul_B{T <: Float}(
    x      :: BEDFile{T},
    y      :: SharedVector{T},
    v      :: PlinkGPUVariables{T};
    mask_n :: Vector{Int} = ones(Int, size(y)),
    pids   :: Vector{Int} = procs(x),
)
    xty = SharedArray{T}((size(x,2),), pids=pids) :: SharedVector{T}
    At_mul_B!(xty, x, y, v)
    return xty
end

max_stars_repo_name: klkeys/PLINK

max_issues_repo_name: klkeys/PLINK


text:
using SchwarzSolver
using Test

@testset "SchwarzSolver.jl" begin
    # Write your own tests here.
end


text:
"""Julia program to implement Merge sort
Merge Sort is a Divide and Conquer algorithm. It divides the input array into two halves, calls itself for the two halves, and then merges the two sorted halves.
First divide the list into the smallest unit (1 element), then compare each element with the adjacent list to sort and merge the two adjacent lists.
Finally all the elements are sorted and merged.
"""
function merge_sort(arr)
    if length(arr) <= 1
        return arr
    end
    middle = trunc(Int, length(arr) / 2) # Finding the mid of the array
    L = arr[1:middle]
    R = arr[middle+1:end]

    merge_sort(L) # Sorting the first half
    merge_sort(R) # Sorting the second half

    # Copy data to temp arrays L[] and R[]
    i = j = k = 1
    while i <= length(L) && j <= length(R)
        if L[i] < R[j]
            arr[k] = L[i]
            i+=1
        else
            arr[k] = R[j]
            j+=1
         end
         k+=1
     end
    # Checking if any element was left
     while i <= length(L)
          arr[k] = L[i]
          i+=1
          k+=1
      end
      while j <= length(R)
          arr[k] = R[j]
          j+=1
          k+=1
      end
      arr
end

print("How many numbers do you want to sort? ")
n = readline()
n = parse(Int64, n)
if (n <= 0)
    println("No numbers to sort!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
merge_sort(arr)
print("The numbers in sorted order is ")
print(arr)

"""
Time Complexity - O(n*Log n), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT
How many numbers do you want to sort? 4
Enter the numbers: 8 2 7 6
The numbers in sorted order is [2, 6, 7, 8]
"""

max_issues_repo_name: adarshnjena/NeoAlgo

max_forks_repo_name: adarshnjena/NeoAlgo


text:
using WinRPM
WinRPM.update()

max_stars_repo_name: wookay/TestWinRPM

max_issues_repo_name: wookay/TestWinRPM

max_forks_repo_name: wookay/TestWinRPM


text:
module DirectedRoundings

export RoundNearest, RoundUp, RoundDown,
       RoundToZero, RoundFromZero,
       RoundLoHi, RoundValue

import Base: RoundNearest, RoundUp, RoundDown, RoundToZero, RoundFromZero
             +, -, *, /, \, hypot, sqrt, cbrt,
             div, fld, cld, mod, rem, divrem, fldmod,
             abs, flipsign, copysign,

export  +, -, *, /, \, hypot, sqrt, cbrt,
        div, fld, cld, mod, rem, divrem, fldmod,
        abs, flipsign, copysign

using Base.Rounding

const RoundLoHi  = RoundingMode{:HiLo}()
const RoundValue = RoundingMode{:Value}()

#=
    •    round the significand to use fewer bits

        •    RoundNearest (default)
        •    RoundUp
        •    RoundDown

    •    round the floating point value to an integer

        •    RoundToZero
        •    RoundFromZero
=#

# allow multiparameter operators taking 1, 2, 3, or 4 args

@inline RoundNearest(fn::Function, a::T) where {T<:AbstractFloat} =
    rounded(fn, a, RoundNearest)
@inline RoundNearest(fn::Function, a::T, b::T) where {T<:AbstractFloat} =
    rounded(fn, a, b, RoundNearest)
@inline RoundNearest(fn::Function, a::T, b::T, c::T) where {T<:AbstractFloat} =
    rounded(fn, a, b, c, RoundNearest)
@inline RoundNearest(fn::Function, a::T, b::T, c::T, d::T) where {T<:AbstractFloat} =
    rounded(fn, a, b, c, d, RoundNearest)

@inline RoundUp(fn::Function, a::T) where {T<:AbstractFloat} =
    rounded(fn, a, RoundUp)
@inline RoundUp(fn::Function, a::T, b::T) where {T<:AbstractFloat} =
    rounded(fn, a, b, RoundUp)
@inline RoundUp(fn::Function, a::T, b::T, c::T) where {T<:AbstractFloat} =
    rounded(fn, a, b, c, RoundUp)
@inline RoundUp(fn::Function, a::T, b::T, c::T, d::T) where {T<:AbstractFloat} =
    rounded(fn, a, b, c, d, RoundUp)

@inline RoundDown(fn::Function, a::T) where {T<:AbstractFloat} =
    rounded(fn, a, RoundDown)
@inline RoundDown(fn::Function, a::T, b::T) where {T<:AbstractFloat} =
    rounded(fn, a, b, RoundDown)
@inline RoundDown(fn::Function, a::T, b::T, c::T) where {T<:AbstractFloat} =
    rounded(fn, a, b, c, RoundDown)
@inline RoundDown(fn::Function, a::T, b::T, c::T, d::T) where {T<:AbstractFloat} =
    rounded(fn, a, b, c, d, RoundDown)

@inline RoundToZero(fn::Function, a::T) where {T<:AbstractFloat} =
    rounded(fn, a, RoundToZero)
@inline RoundToZero(fn::Function, a::T, b::T) where {T<:AbstractFloat} =
    rounded(fn, a, b, RoundToZero)
@inline RoundToZero(fn::Function, a::T, b::T, c::T) where {T<:AbstractFloat} =
    rounded(fn, a, b, c, RoundToZero)
@inline RoundToZero(fn::Function, a::T, b::T, c::T, d::T) where {T<:AbstractFloat} =
    rounded(fn, a, b, c, d, RoundToZero)

# support RoundFromZero

@inline RoundFromZero(fn::Function, a::T) where {T<:BigFloat} =
    rounded(fn, a, RoundFromZero)
@inline RoundFromZero(fn::Function, a::T, b::T) where {T<:BigFloat} =
    rounded(fn, a, b, RoundFromZero)
@inline RoundFromZero(fn::Function, a::T, b::T, c::T) where {T<:BigFloat} =
    rounded(fn, a, b, c, RoundFromZero)
@inline RoundFromZero(fn::Function, a::T, b::T, c::T, d::T) where {T<:BigFloat} =
    rounded(fn, a, b, c, d, RoundFromZero)

@inline function RoundFromZero(fn::Function, a::T) where {T<:AbstractFloat}
    result = RoundUp(a, b)
    if result < 0
       result = RoundDown(a, b)
    end
    return result
end
@inline function RoundFromZero(fn::Function, a::T, b::T) where {T<:AbstractFloat}
    result = RoundUp(a, b)
    if result < 0
       result = RoundDown(a, b)
    end
    return result
end
@inline function RoundFromZero(fn::Function, a::T, b::T, c::T) where {T<:AbstractFloat}
    result = RoundUp(a, b, c)
    if result < 0
       result = RoundDown(a, b, c)
    end
    return result
end
@inline function RoundFromZero(fn::Function, a::T, b::T, c::T, d::T) where {T<:AbstractFloat}
    result = RoundUp(a, b, c, d)
    if result < 0
       result = RoundDown(a, b, c, d)
    end
    return result
end

# directed rounding in a functional context

@inline function rounded(fn::Function, a::T,
                         mode::RoundingMode) where {T<:AbstractFloat}
     setrounding(T, mode) do
         fn(a)
     end
end

@inline function rounded(fn::Function, a::T, b::T,
                         mode::RoundingMode) where {T<:AbstractFloat}
     setrounding(T, mode) do
         fn(a, b)
     end
end

@inline function rounded(fn::Function, a::T, b::T, c::T,
                         mode::RoundingMode) where {T<:AbstractFloat}
     setrounding(T, mode) do
         fn(a, b, c)
     end
end

@inline function rounded(fn::Function, a::T, b::T, c::T, d::T,
                         mode::RoundingMode) where {T<:AbstractFloat}
     setrounding(T, mode) do
         fn(a, b, c, d)
     end
end

# bidirectional rounding in a functional context

@inline function RoundLoHi(fn::Function, a::T) where {T<:AbstractFloat}
     hi = rounded(fn, a, RoundUp)
     lo = rounded(fn, a, RoundDown)
     return minmax(hi, lo)
end

@inline function RoundLoHi(fn::Function, a::T, b::T) where {T<:AbstractFloat}
     hi = rounded(fn, a, b, RoundUp)
     lo = rounded(fn, a, b, RoundDown)
     return minmax(hi, lo)
end

@inline function RoundLoHi(fn::Function, a::T, b::T, c::T) where {T<:AbstractFloat}
     hi = rounded(fn, a, b, c, RoundUp)
     lo = rounded(fn, a, b, c, RoundDown)
     return minmax(hi, lo)
end

@inline function RoundLoHi(fn::Function, a::T, b::T, c::T, d::T) where {T<:AbstractFloat}
     hi = rounded(fn, a, b, c, d, RoundUp)
     lo = rounded(fn, a, b, c, d, RoundDown)
     return minmax(hi, lo)
end

# obtain the most informing, least misleading representation of value

const twopowex = [2.0, 4.0, 8.0, 16.0, 32.0, 64.0]

function RoundValue(fn::Function, a::T) where {T<:AbstractFloat}
   hi, lo = RoundLoHi(fn, a)
   frhi, exhi = frexp(hi)
   frlo, exlo = frexp(lo)
   if exhi !== exlo
      dex = exhi - exlo
      fex = dex < 7 ? twopowex[dex] : 2.0^dex
      frlo = frlo * fex
      frex = frex - dex
   end
   signifhi = frhi
   signiflo = frlo
   sigbits  = 53
   while sigbits > 0 && (signifhi !== signiflo)
      sigbits -= 1
      signifhi = signif(hi, sigbits, 2)
      signiflo = signif(lo, sigbits, 2)
   end
   if sigbits == 0 && (signifhi !== signiflo)
       span = hi - lo
       throw(ErrorException("interval [$lo, $hi] spans $span, which does not resolve"))
   end
   return signiflo
end

end # DirectedRoundings module


text:
using Oceananigans.Architectures: device_event, arch_array

"""
    BatchedTridiagonalSolver

A batched solver for large numbers of triadiagonal systems.
"""
struct BatchedTridiagonalSolver{A, B, C, T, G, R, P}
               a :: A
               b :: B
               c :: C
               t :: T
            grid :: G
    architecture :: R
      parameters :: P
end

"""
    BatchedTridiagonalSolver(grid; lower_diagonal, diagonal, upper_diagonal, parameters=nothing)

Construct a solver for batched tridiagonal systems on `grid` of the form

                           bⁱʲ¹ ϕⁱʲ¹ + cⁱʲ¹ ϕⁱʲ²   = fⁱʲ¹,  k = 1
           aⁱʲᵏ⁻¹ ϕⁱʲᵏ⁻¹ + bⁱʲᵏ ϕⁱʲᵏ + cⁱʲᵏ ϕⁱʲᵏ⁺¹ = fⁱʲᵏ,  k = 2, ..., N-1
           aⁱʲᴺ⁻¹ ϕⁱʲᴺ⁻¹ + bⁱʲᴺ ϕⁱʲᴺ               = fⁱʲᴺ,  k = N

where `a` is the `lower_diagonal`, `b` is the `diagonal`, and `c` is the `upper_diagonal`.
`ϕ` is the solution and `f` is the right hand side source term passed to `solve!(ϕ, tridiagonal_solver, f)`

`a`, `b`, `c`, and `f` can be specified in three ways:

1. A 1D array means that `aⁱʲᵏ = a[k]`.

2. A 3D array means that `aⁱʲᵏ = a[i, j, k]`.

3. Otherwise, `a` is assumed to be callable:
    * If `isnothing(parameters)` then `aⁱʲᵏ = a(i, j, k, grid, args...)`.
    * If `!isnothing(parameters)` then `aⁱʲᵏ = a(i, j, k, grid, parameters, args...)`.
    where `args...` are `Varargs` passed to `solve_batched_tridiagonal_system!(ϕ, solver, args...)`.
"""
function BatchedTridiagonalSolver(arch, grid;
                                  lower_diagonal,
                                  diagonal,
                                  upper_diagonal,
                                  scratch = arch_array(arch, zeros(eltype(grid), grid.Nx, grid.Ny, grid.Nz)),
                                  parameters = nothing)

    return BatchedTridiagonalSolver(lower_diagonal, diagonal, upper_diagonal,
                                    scratch, grid, arch, parameters)
end

@inline get_coefficient(a::AbstractArray{T, 1}, i, j, k, grid, p, args...) where {T} = @inbounds a[k]
@inline get_coefficient(a::AbstractArray{T, 3}, i, j, k, grid, p, args...) where {T} = @inbounds a[i, j, k]
@inline get_coefficient(a::Base.Callable, i, j, k, grid, p, args...) = a(i, j, k, grid, p, args...)
@inline get_coefficient(a::Base.Callable, i, j, k, grid, ::Nothing, args...) = a(i, j, k, grid, args...)

"""
    solve!(ϕ, solver, rhs, args...; dependencies = device_event(solver.architecture))

Solve the batched tridiagonal system of linear equations with right hand side
`rhs` and lower diagonal, diagonal, and upper diagonal coefficients described by the
`BatchedTridiagonalSolver` `solver`. `BatchedTridiagonalSolver` uses a modified
TriDiagonal Matrix Algorithm (TDMA).

The result is stored in `ϕ` which must have size `(grid.Nx, grid.Ny, grid.Nz)`.

Reference implementation per Numerical Recipes, Press et. al 1992 (§ 2.4).
"""
function solve!(ϕ, solver::BatchedTridiagonalSolver, rhs, args...; dependencies = device_event(solver.architecture))

    a, b, c, t, parameters = solver.a, solver.b, solver.c, solver.t, solver.parameters
    grid = solver.grid

    event = launch!(solver.architecture, grid, :xy,
                    solve_batched_tridiagonal_system_kernel!, ϕ, a, b, c, rhs, t, grid, parameters, args...,
                    dependencies = dependencies)

    wait(device(solver.architecture), event)

    return nothing
end

@inline float_eltype(ϕ::AbstractArray{T}) where T <: AbstractFloat = T
@inline float_eltype(ϕ::AbstractArray{<:Complex{T}}) where T <: AbstractFloat = T

@kernel function solve_batched_tridiagonal_system_kernel!(ϕ, a, b, c, f, t, grid, p, args...)
    Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz

    i, j = @index(Global, NTuple)

    @inbounds begin
        β  = get_coefficient(b, i, j, 1, grid, p, args...)
        f₁ = get_coefficient(f, i, j, 1, grid, p, args...)
        ϕ[i, j, 1] = f₁ / β

        @unroll for k = 2:Nz
            cᵏ⁻¹ = get_coefficient(c, i, j, k-1, grid, p, args...)
            bᵏ   = get_coefficient(b, i, j, k,   grid, p, args...)
            aᵏ⁻¹ = get_coefficient(a, i, j, k-1, grid, p, args...)

            t[i, j, k] = cᵏ⁻¹ / β
            β = bᵏ - aᵏ⁻¹ * t[i, j, k]

            fᵏ = get_coefficient(f, i, j, k, grid, p, args...)

            # If the problem is not diagonally-dominant such that `β ≈ 0`,
            # the algorithm is unstable and we elide the forward pass update of ϕ.
            definitely_diagonally_dominant = abs(β) > 10 * eps(float_eltype(ϕ))
            !definitely_diagonally_dominant && break
            ϕ[i, j, k] = (fᵏ - aᵏ⁻¹ * ϕ[i, j, k-1]) / β
        end

        @unroll for k = Nz-1:-1:1
            ϕ[i, j, k] -= t[i, j, k+1] * ϕ[i, j, k+1]
        end
    end
end


text:
immutable NoncentralChisq <: ContinuousUnivariateDistribution
    df::Float64
    ncp::Float64
    function NonCentralChisq(d::Real, nc::Real)
    	if d >= 0.0 && nc >= 0.0
    		new(float64(d), float64(nc))
    	else
    		error("df and ncp must be non-negative")
    	end
    end
end

@_jl_dist_2p NoncentralChisq nchisq

insupport(d::NoncentralChisq, x::Number) = isreal(x) && isfinite(x) && 0.0 < x


text:
module TestExamplesDGMulti1D

using Test
using Trixi

include("test_trixi.jl")

EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_1d")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "DGMulti 1D" begin
  @trixi_testset "elixir_euler_flux_diff.jl " begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_flux_diff.jl"),
      cells_per_dimension = (16,),
      l2 = [7.853842541289665e-7, 9.609905503440606e-7, 2.832322219966481e-6],
      linf = [1.5003758788711963e-6, 1.802998748523521e-6, 4.83599270806323e-6]
    )
  end

  @trixi_testset "elixir_euler_flux_diff.jl (convergence)" begin
    mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "elixir_euler_flux_diff.jl"), 3)
    @test isapprox(mean_convergence[:l2], [4.1558759698638434, 3.977911306037128, 4.041421206468769], rtol=0.05)
  end

end

# Clean up afterwards: delete Trixi output directory
@test_nowarn isdir(outdir) && rm(outdir, recursive=true)

end # module


text:
using DataStructures
using Base.Meta

struct EClass
    id::Int64
end

"""
Check if an expr is an enode ⟺ all args are e-classes
"""
function isenode(e::Expr)
    return all(x -> x isa EClass, getfunargs(e))
end
# literals are enodes
isenode(x::EClass) = false
isenode(x) = true

### Definition 2.3: canonicalization
iscanonical(U::IntDisjointSets, n::Expr) = n == canonicalize(U, n)
iscanonical(U::IntDisjointSets, e::EClass) = find_root!(U, e.id) == e.id

# canonicalize an e-term n
# throws a KeyError from find_root! if any of the child classes
# was not found as the representative element in a set in U
function canonicalize(U::IntDisjointSets, n::Expr)
    @assert isenode(n)
    ne = copy(n)
    setfunargs!(ne, [EClass(find_root!(U, x.id)) for x ∈ getfunargs(ne)])
    @debug("canonicalized ", n, " to ", ne)
    return ne
end

# canonicalize in place
function canonicalize!(U::IntDisjointSets, n::Expr)
    @assert isenode(n)
    setfunargs!(n, [EClass(find_root!(U, x.id)) for x ∈ getfunargs(n)])
    @debug("canonicalized ", n)
    return n
end

# literals are already canonical
canonicalize(U::IntDisjointSets, n) = n
canonicalize!(U::IntDisjointSets, n) = n


text:
Base.copy!(dest::AxisArray{TT, N, AA, As},
           source::AxisArray{T, N, A, As}) where {TT, N, AA, T, A, As} = begin
    copy!(dest.data, source.data)
    dest
end

"""
    copy!(dest::AxisArray, source::AxisArray, ignore_axes::Bool=false) -> AxisArray

Copies data from one AxisArray to another. If `ignore_axes` is `false` (default), then all
axis names should match, though they may be permutated. If `ignore_axes` is `true`, then the
function copies the underlying arrays, using the behaviour standard for the types of the
underlying arrays.
"""
Base.copy!(dest::AxisArray{TT, N, AA, AAs}, source::AxisArray{T, N, A, As},
           ignore_axes::Bool=false) where {TT, N, AA, AAs, T, A, As} = begin
    if ignore_axes
        copy!(dest.data, source.data)
    else
        perm = indexin(collect(axisnames(dest)), collect(axisnames(source)))
        any(iszero(x) for x in perm) &&
                throw(ArgumentError("""
                    Axes of source and destination do not match.
                    Use copy!(A, B, false) if you want to ignore axis information.
                    """))

        permutedims!(dest.data, source.data, perm)
    end
    dest
end

"""
    copy!(dest::AxisArray, daxes::Tuple, src::AxisArray, saxes::Tuple, ignore_axes=false)

Copies a block from one array to a block of another. The tuples identifying the blocks can
be anything that `view` accepts. Hence it may include axis information.
"""
Base.copy!(dest::AxisArray, daxes::Tuple{Axis, Vararg{Axis}},
           src::AxisArray, saxes::Tuple{Axis, Vararg{Axis}}) =
    copy!(view(dest, daxes...), view(src, saxes...))

max_stars_repo_name: mdavezac/DFTShims

max_issues_repo_name: mdavezac/DFTShims

max_forks_repo_name: mdavezac/DFTShims


text:
using Infinities, Base64, Base.Checked, Test

@testset "InfiniteCardinal" begin
    @testset "basics" begin
        @test !isone(ℵ₀)
        @test !iszero(ℵ₀)
        @test sign(ℵ₀) ≡ 1 && !signbit(ℵ₀)
        @test angle(ℵ₀) ≡ 0
        @test Integer(∞) ≡ convert(Integer,∞) ≡ Integer(ℵ₀) ≡ convert(Integer, ℵ₀) ≡ ℵ₀
        @test abs(ℵ₀) ≡ ℵ₀
        @test zero(ℵ₀) ≡ 0
        @test one(ℵ₀) ≡ one(InfiniteCardinal{0}) ≡ oneunit(ℵ₀) ≡ oneunit(InfiniteCardinal{0}) ≡ 1
        @test isinf(ℵ₀) && !isfinite(ℵ₀)
        @test Integer(RealInfinity()) ≡ Integer(ComplexInfinity()) ≡ ℵ₀
        @test_throws InexactError Integer(-∞)
        @test_throws InexactError Integer(exp(0.1im)*∞)
    end

    @testset "equality" begin
        @test ℵ₀ == ℵ₀
        @test ℵ₀ ≠ ℵ₁
        @test ℵ₀ == ∞ && ∞ == ℵ₀
        @test ℵ₁ ≠ ∞ && ∞ ≠ ℵ₁
        @test ℵ₀ ≠ -∞ && -∞ ≠ ℵ₀
        @test ℵ₀ == Inf && Inf == ℵ₀
        @test ℵ₀ ≠ 5 && 5 ≠ ℵ₀
    end

    @testset "inequalities" begin
        @test !isless(ℵ₀, ℵ₀)
        @test isless(ℵ₀, ℵ₁)
        @test isless(Inf, ℵ₁) && !isless(ℵ₁, Inf)
        @test !isless(Inf, ℵ₀) && !isless(ℵ₀, Inf)
        @test isless(5, ℵ₀) && !isless(ℵ₀, 5)
        @test isless(5, ℵ₁) && !isless(ℵ₁, 5)

        @test !(ℵ₀ < ℵ₀) && !(ℵ₀ > ℵ₀)
        @test ℵ₀ ≤ ℵ₀ && ℵ₀ ≥ ℵ₀
        @test ℵ₀ < ℵ₁ && ℵ₀ ≤ ℵ₁
        @test !(ℵ₀ > ℵ₁) && !(ℵ₀ ≥ ℵ₁)
        @test ℵ₁ ≤ ℵ₁ && ℵ₁ ≥ ℵ₁ && !(ℵ₁ > ℵ₁) && !(ℵ₁ < ℵ₁)

        @test !(∞ < ℵ₀) && ∞ ≤ ℵ₀
        @test !(∞ > ℵ₀) && ∞ ≥ ℵ₀
        @test ∞ < ℵ₁ && ∞ ≤ ℵ₁
        @test !(∞ > ℵ₁) && !(∞ ≥ ℵ₁)
        @test ℵ₀ ≤ ∞
        @test !(ℵ₁ < ∞) && !(ℵ₁ ≤ ∞)
        @test !(ℵ₀ > ∞)
        @test ℵ₁ > ∞ && ℵ₁ ≥ ∞

        @test -∞ < ℵ₀ && -∞ ≤ ℵ₀
        @test -∞ < ℵ₁ && -∞ ≤ ℵ₁
        @test !(-∞ > ℵ₀) && RealInfinity() ≥ ℵ₀
        @test !(-∞ > ℵ₁) && !(RealInfinity() ≥ ℵ₁)
        @test !(-∞ > ℵ₀) && !(-∞ ≥ ℵ₀)
        @test ℵ₀ > -∞
        @test ℵ₁ > -∞ && ℵ₁ ≥ -∞
        @test ℵ₀ ≤ RealInfinity()
        @test !(ℵ₀ ≤ -∞)
        @test !(ℵ₁ ≤ RealInfinity())

        @test Inf < ℵ₁ && !(ℵ₁ < Inf)
        @test !(Inf < ℵ₀) && !(ℵ₀ < Inf)
        @test (5 < ℵ₀) && !(ℵ₀ < 5)
        @test 5 ≤ ℵ₁ && !(5 ≥ ℵ₁) && !(ℵ₁ ≤ 5) && ℵ₁ > 5
        @test 5 < ℵ₀ && 5 ≤ ℵ₀
        @test !(ℵ₀ < 5) && !(ℵ₀ ≤ 5)
        @test ℵ₀ > 5 && ℵ₀ ≥ 5
        @test !(5 > ℵ₀) && !(5 ≥ ℵ₀)
    end

    @testset "min/max" begin
        @test @inferred(min(ℵ₀,ℵ₁)) ≡ ℵ₀
        @test @inferred(min(ℵ₀,ℵ₁)) ≡ min(ℵ₁,ℵ₀) ≡ ℵ₀
        @test @inferred(max(ℵ₀,ℵ₁)) ≡ max(ℵ₁,ℵ₀) ≡ ℵ₁
        @test min(∞,ℵ₀) ≡ min(ℵ₀,∞) ≡ ∞
        @test max(∞,ℵ₀) ≡ max(ℵ₀,∞) ≡ ℵ₀
        @test min(-∞,ℵ₀) ≡ min(ℵ₀,-∞) ≡ -∞
        @test max(-∞,ℵ₀) ≡ max(ℵ₀,-∞) ≡ ℵ₀
        @test min(5,ℵ₀) ≡ min(ℵ₀,5) ≡ 5
        @test max(5,ℵ₀) ≡ max(ℵ₀,5) ≡ ℵ₀
    end

    @testset "algebra" begin
        @test +(ℵ₀) ≡ ℵ₀
        @test ℵ₀ + ℵ₁ ≡ ℵ₁ + ℵ₀ ≡ ℵ₁
        @test 5 + ℵ₀ ≡ ℵ₀ + 5 ≡ ℵ₀
        @test ℵ₀ - 5 ≡ ℵ₀
        @test 5 - ℵ₀ ≡ -∞
        @test ℵ₀ - ℵ₀ ≡ NotANumber()
        @test -ℵ₀ ≡ -∞

        @test *(ℵ₀) ≡ ℵ₀
        @test ℵ₀ * ℵ₀ ≡ ℵ₀
        @test ℵ₀ * ∞ ≡ ∞ * ℵ₀ ≡ ℵ₀
        @test 5 * ℵ₀ ≡ ℵ₀ * 5 ≡ ℵ₀
        @test 5.0ℵ₀ ≡ ℵ₀*5.0 ≡ RealInfinity()
        @test_throws ArgumentError (-5) * ℵ₀

        @test ℵ₀ - 5.1 ≡ ∞
        @test ℵ₀ + 5.1 ≡ ∞
        @test 5.1 + ℵ₀ ≡ ∞
        @test 5.1 - ℵ₀ ≡ -∞

        @test ℵ₀ + ∞ == ∞
        @test ∞ + ℵ₀ == ∞
        @test ℵ₀ + RealInfinity() == ∞
        @test ℵ₀ - (-∞) == ∞
        @test  RealInfinity() + ℵ₀ == ∞
    end

    @testset "fld/cld/div/mod" begin
        @test ℵ₀ ÷ 5 ≡ ℵ₀
        @test ℵ₀ ÷ ℵ₀ ≡ NotANumber()
        @test 5 ÷ ℵ₀ ≡ 0
        @test fld(-5, ℵ₀) ≡ -1
        @test cld(5, ℵ₀) ≡ 1
        @test mod(ℵ₀,ℵ₀) ≡ NotANumber()
        @test mod(ℵ₀,6) ≡ NotANumber()
        @test mod(5,ℵ₀) ≡ 5
        @test_throws ArgumentError mod(-1,ℵ₀)
    end

    @test string(ℵ₀) == stringmime("text/plain", ℵ₀) == "ℵ₀"
    @test string(ℵ₁) == stringmime("text/plain", ℵ₁) == "ℵ₁"
    @test Base.to_index(ℵ₀) ≡ Base.to_shape(ℵ₀) ≡ ℵ₀
    @test Base.to_shape((∞,)) ≡ (ℵ₀,)

    @testset "Set" begin
        s = Set([ℵ₀,ℵ₁,∞,1])
        @test 1 ∈ s
        @test ∞ ∈ s
        @test ℵ₀ ∈ s
        @test ℵ₁ ∈ s
        @test 2 ∉ s
    end

    @testset "checked" begin
        @test checked_add(5, ℵ₀) ≡ checked_add(ℵ₀, 5) ≡ ℵ₀
        @test checked_sub(ℵ₀, 5) ≡ ℵ₀
        @test checked_sub(5, ℵ₀) ≡ -∞
        @test checked_mul(-5, ℵ₀) ≡ checked_mul(ℵ₀, -5) ≡ -∞
    end

    @testset "indexing" begin
        @test_throws BoundsError randn(3)[ℵ₀]
        @test_throws ErrorException Base._unsafe_getindex(IndexCartesian(),permutedims(1:3)',ℵ₀)
        @test_throws BoundsError view(randn(3),1:2)[ℵ₀]
        # check ambiguity is not introduced
        @test view(collect(1:5),2)[] == 2
    end
end


text:
@testset "649.dota2-senate.jl" begin
    @test predict_party_victory("RDD") == "Dire"
    @test predict_party_victory("RD") == "Radiant"
end


text:
using JLD
using PyPlot
using StatsBase
using Distributions
using VBMatrixFactorization

function toy_matrix(L, M, H, std)
    # B is a random matrix whose columns we select from
    B = randn(L, H)
    # A is the sparse matrix that selects the columns from A
    A = zeros(M, H)
    for m in 1:M
        A[m, rand(1:H)] = 1
    end

    Y = B*A' + rand(Normal(0, std), L, M)
    return Y, A, B
end

L = 10
M = 20
H = 2
Y_toy, A_toy, B_toy = toy_matrix(L, M, H, 0.05);

# first, initialize a variable of the vbmf_parameters type
# then use it to run the basic vbmf algorithm
# here, the whole run is saved into the data_path folder
# if desc is not specified, then the subfolder with the actual data
# is named using date and time

println(" ----------- Basic VB Matrix factorization ---------------- ")
println("")
params_init = VBMatrixFactorization.vbmf_init(Y_toy, H, ca = 0.1, cb = 0.1, sigma2 = 0.1);
data_path = "./data"
res_vbmf = VBMatrixFactorization.vbmf(Y_toy, params_init, 100, est_covs = true, est_var = true, verb = true,
    logdir = data_path, desc = "vbmf_test");
err = norm(Y_toy - res_vbmf.YHat)
println("||Y - Yhat|| = $err")

# this is how the log is loaded and used
vbmf_log, data, priors = VBMatrixFactorization.load_log(string(data_path,"/vbmf_test"));
params_vbmf = VBMatrixFactorization.vbmf_parameters(); # dummy variable to store individual step data

iter = 3
time_slice = VBMatrixFactorization.extract_params!(vbmf_log, iter, params_vbmf);
err = norm(Y_toy - time_slice.YHat)
println("in iteration number $iter, the error was ||Y - Yhat|| = $err")
println("")

# vbmf_sparse usage
println(" ----------- VB Matrix factorization with sparse A ---------------- ")
println("")
params_sparse_init = VBMatrixFactorization.vbmf_sparse_init(Y_toy, H, ca = 0.1, cb = 0.1, sigma = 0.1);
params_sparse_init.AHat = A_toy
sparse_vbmf = VBMatrixFactorization.vbmf_sparse(Y_toy, params_sparse_init, 100, diag_var = false, verb = true,
    logdir = data_path, desc = "sparse_test", full_cov = true);
err = norm(Y_toy - sparse_vbmf.YHat)
println("||Y - Yhat|| = $err")
println("")
# compare the results

println(" original A          basic reconstruction          sparse reconstruction with true start")
for m in 1:M
    println(A_toy[m,:],"   ", res_vbmf.AHat[m,:], "   ", sparse_vbmf.AHat[m,:])
end


text:
using Dates
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(ESDL.DAT.getchunkoffsets), ESDL.DAT.DATConfig{1, 1}})
    precompile(Tuple{typeof(ESDL.DAT.getLoopCacheSize), Int64, Array{Int64, 1}, Float64, Array{NamedTuple{(:iloopax, :cs, :iscompressed, :innerleap), Tuple{Int64, Int64, Bool, Int64}}, 1}})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##s74#101")), Int, Int, Int, Int, Int, Int, Int, Int, Int, Int, Int, Int, Int, Int, Int, Int, Int, Int, Int, Int, Int, Int})
    precompile(Tuple{typeof(ESDL.CubeAPI.getLonLatsToRead), ESDL.CubeAPI.CubeConfig, Tuple{Int64, Int64}, Tuple{Int64, Int64}})
    precompile(Tuple{typeof(ESDL.CubeAPI.getCubeData), ESDL.CubeAPI.Cube, Array{String, 1}, Tuple{Dates.Date, Dates.Date}, Tuple{Int64, Int64}, Tuple{Int64, Int64}})
    precompile(Tuple{typeof(ESDL.DAT.allocatecachebuf), ESDL.DAT.InputCube{4}, Array{Int64, 1}})
    precompile(Tuple{typeof(ESDL.DAT.runLoop), ESDL.DAT.DATConfig{1, 1}})
    precompile(Tuple{typeof(ESDL.DAT.analyzeAxes), ESDL.DAT.DATConfig{1, 1}})
    precompile(Tuple{typeof(ESDL.DAT.allocatecachebuf), ESDL.DAT.OutputCube, Array{Int64, 1}})
    precompile(Tuple{typeof(ESDL.CubeAPI.readFromDataYear), ESDL.CubeAPI.Cube, Base.SubArray{Float32, 3, Array{Float32, 4}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64}, true}, Base.SubArray{UInt8, 3, Array{UInt8, 4}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64}, true}, String, Int64, Int64, Int64, Int64, Int64, Int64, Int64, Int64, Int64})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##90#95")), ESDL.DAT.OutputCube})
    precompile(Tuple{typeof(ESDL.DAT.updateinars), ESDL.DAT.DATConfig{1, 1}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}}})
    precompile(Tuple{getfield(ESDL.CubeAPI, Symbol("##getCubeData#7")), Array{Int64, 1}, Array{Any, 1}, Tuple{Int64, Int64}, Tuple{Int64, Int64}, Array{Any, 1}, typeof(identity), (ESDL.CubeAPI).Cube})
    precompile(Tuple{typeof(ESDL.Cubes._read), ESDL.CubeAPI.SubCubeVPerm{Float32}, Tuple{Base.SubArray{Float32, 4, Array{Float32, 4}, Tuple{Base.IteratorsMD.CartesianIndices{4, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}}}}, false}, Base.SubArray{UInt8, 4, Array{UInt8, 4}, Tuple{Base.IteratorsMD.CartesianIndices{4, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}}}}, false}}, Base.IteratorsMD.CartesianIndices{4, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}}}})
    precompile(Tuple{typeof(ESDL.Cubes.caxes), ESDL.CubeAPI.SubCubeV{Float32, ESDL.CubeAPI.Cube}})
    precompile(Tuple{typeof(ESDL.Proc.MSC.fillmsc), Int64, Array{Float64, 1}, Array{Int64, 1}, Array{Float32, 1}, Int64})
    precompile(Tuple{typeof(ESDL.DAT.__init__)})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##104#117")), Int64})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##105#118")), Int64})
    precompile(Tuple{typeof(ESDL.Cubes.chunkoffset), ESDL.CubeAPI.SubCubeV{Float32, ESDL.CubeAPI.Cube}})
    precompile(Tuple{typeof(ESDL.Cubes.cubechunks), ESDL.CubeAPI.SubCubeVPerm{Float32}})
    precompile(Tuple{typeof(ESDL.DAT.getRetCubeType), ESDL.DAT.OutputCube, Bool, Float64})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:ref_time}})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##114#127")), Int64})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##128#129")), Int64})
    precompile(Tuple{typeof(ESDL.Cubes.chunkoffset), ESDL.CubeAPI.SubCubeVPerm{Float32}})
    precompile(Tuple{typeof(ESDL.ESDLTools.freshworkermodule)})
    precompile(Tuple{typeof(ESDL.CubeAPI.expandknownvars), Array{String, 1}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseConfig), Array{Base.SubString{String}, 1}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:end_time}})
    precompile(Tuple{typeof(((ESDL.Cubes).Axes).getOutAxis), Tuple{((ESDL.Cubes).Axes).ByName}, Array{(((ESDL.Cubes).Axes).CubeAxis{T, S} where S) where T, 1}, Tuple{(ESDL.CubeAPI).SubCubeV{Float32, (ESDL.CubeAPI).Cube}}, Tuple{Int64, Array{Float64, 1}, Array{Int64, 1}}, typeof(identity)})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:calendar}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:start_time}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:file_format}})
    precompile(Tuple{typeof((ESDL.DAT).generateOutCube), Type{((ESDL.Cubes).CubeMem{T, N} where N) where T}, Type{Int}, (ESDL.DAT).OutputCube, Array{Int64, 1}})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##113#126")), Int64})
    precompile(Tuple{typeof(ESDL.DAT.getCacheSizes), ESDL.DAT.DATConfig{1, 1}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:model_version}})
    precompile(Tuple{typeof(ESDL.CubeAPI.getTimesToRead), Dates.Date, Dates.Date, ESDL.CubeAPI.CubeConfig})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##110#123")), Int64})
    precompile(Tuple{typeof(ESDL.Cubes.cubechunks), ESDL.CubeAPI.SubCubeV{Float32, ESDL.CubeAPI.Cube}})
    precompile(Tuple{typeof(ESDL.Cubes.Axes.findAxis), ESDL.Cubes.Axes.ByName, Array{ESDL.Cubes.Axes.CubeAxis{T, S} where S where T, 1}})
    precompile(Tuple{typeof(ESDL.Cubes.Axes.findAxis), ESDL.Cubes.Axes._RangeAxis{Dates.Date, :Time, ESDL.Cubes.Axes.YearStepRange}, Array{ESDL.Cubes.Axes.CubeAxis{T, S} where S where T, 1}})
    precompile(Tuple{typeof(ESDL.Cubes.Axes.findAxis), ESDL.Cubes.Axes._RangeAxis{Float64, :Lat, Base.StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, ESDL.CubeAPI.SubCubeVPerm{Float32}})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##mapCube#54")), Float64, (ESDL.DAT).InDims, (ESDL.DAT).OutDims, Bool, String, Bool, Bool, Bool, (Base.Iterators).Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, typeof((ESDL.DAT).mapCube), typeof(((ESDL.Proc).MSC).removeMSC), Tuple{(ESDL.CubeAPI).SubCubeV{Float32, (ESDL.CubeAPI).Cube}}, Int64, Int})
    precompile(Tuple{typeof(ESDL.Cubes.Axes.findAxis), ESDL.Cubes.Axes._RangeAxis{Float64, :Lat, Base.StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, Array{ESDL.Cubes.Axes.CubeAxis{T, S} where S where T, 1}})
    precompile(Tuple{typeof(ESDL.Cubes.Axes.findAxis), ESDL.Cubes.Axes._RangeAxis{Float64, :Lon, Base.StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, Array{ESDL.Cubes.Axes.CubeAxis{T, S} where S where T, 1}})
    precompile(Tuple{typeof(ESDL.Cubes.Axes.findAxis), ESDL.Cubes.Axes.CategoricalAxis{String, :Variable, Array{String, 1}}, ESDL.CubeAPI.SubCubeVPerm{Float32}})
    precompile(Tuple{typeof(ESDL.Cubes.Axes.findAxis), ESDL.Cubes.Axes._RangeAxis{Float64, :Lon, Base.StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, ESDL.CubeAPI.SubCubeVPerm{Float32}})
    precompile(Tuple{typeof(ESDL.Cubes.Axes.findAxis), ESDL.Cubes.Axes.CategoricalAxis{String, :Variable, Array{String, 1}}, Array{ESDL.Cubes.Axes.CubeAxis{T, S} where S where T, 1}})
    precompile(Tuple{getfield(ESDL.CubeAPI, Symbol("#kw##getCubeData")), NamedTuple{(:longitude, :latitude), Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, typeof(ESDL.CubeAPI.getCubeData), ESDL.CubeAPI.Cube})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##106#119")), Int64})
    precompile(Tuple{typeof(ESDL.CubeAPI.getLandSeaMask!), Array{UInt8, 4}, ESDL.CubeAPI.Cube, Int64, Int64, Int64, Int64})
    precompile(Tuple{typeof(ESDL.Proc.MSC.removeMSC), ESDL.Cubes.Mask.MaskArray{Float32, 1, Array{Float32, 1}}, ESDL.Cubes.Mask.MaskArray{Float32, 1, Array{Float32, 1}}, Int64, Array{Float64, 1}, Array{Int64, 1}})
    precompile(Tuple{getfield(ESDL.ESDLTools, Symbol("##s12#5")), Int, Int, Int, Int, Int, Int, Int, Int})
    precompile(Tuple{typeof(ESDL.DAT.reOrderInCubes), ESDL.DAT.DATConfig{1, 1}})
    precompile(Tuple{typeof(ESDL.DAT.getFrontPerm), ESDL.CubeAPI.SubCubeV{Float32, ESDL.CubeAPI.Cube}, Array{ESDL.Cubes.Axes.CubeAxis{T, S} where S where T, 1}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:chunk_sizes}})
    precompile(Tuple{getfield((ESDL.Proc).MSC, Symbol("##removeMSC#3")), (Base.Iterators).Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, typeof(identity), (ESDL.CubeAPI).SubCubeV{Float32, (ESDL.CubeAPI).Cube}})
    precompile(Tuple{typeof((ESDL.DAT).createworkarray), Type{Int}, Tuple{Int64}})
    precompile(Tuple{typeof(ESDL.ESDLTools.mypermutedims!), Base.SubArray{Float32, 4, Array{Float32, 4}, Tuple{Base.IteratorsMD.CartesianIndices{4, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}}}}, false}, Array{Float32, 4}, Type{Base.Val{(3, 1, 2, 4)}}})
    precompile(Tuple{typeof(ESDL.DAT.runLoop), ESDL.DAT.DATConfig{1, 1}, Base.Iterators.ProductIterator{Tuple{Array{Base.UnitRange{Int64}, 1}, Array{Base.UnitRange{Int64}, 1}, Array{Base.UnitRange{Int64}, 1}}}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:static_data}})
    precompile(Tuple{typeof(ESDL.ESDLTools.mypermutedims!), Base.SubArray{UInt8, 4, Array{UInt8, 4}, Tuple{Base.IteratorsMD.CartesianIndices{4, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}}}}, false}, Array{UInt8, 4}, Type{Base.Val{(3, 1, 2, 4)}}})
    precompile(Tuple{typeof(ESDL.Cubes.caxes), ESDL.CubeAPI.SubCubeVPerm{Float32}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:compression}})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##88#93")), ESDL.DAT.InputCube{4}})
    precompile(Tuple{typeof(ESDL.CubeAPI.getvartype), ESDL.CubeAPI.Cube, String})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("#kw##mapCube")), NamedTuple{(:indims, :outdims), Tuple{(ESDL.DAT).InDims, (ESDL.DAT).OutDims}}, typeof((ESDL.DAT).mapCube), typeof(identity), (ESDL.CubeAPI).SubCubeV{Float32, (ESDL.CubeAPI).Cube}, Int64, Int})
    precompile(Tuple{typeof((ESDL.DAT).setSubRange2), ((ESDL.Cubes).Mask).MaskArray{Float32, 1, Array{Float32, 1}}, Tuple{Array{Float32, 4}, Array{UInt8, 4}}, typeof(identity), Int})
    precompile(Tuple{typeof(ESDL.DAT.innerLoop), typeof(ESDL.Proc.MSC.removeMSC), Tuple{Tuple{Array{Float32, 4}, Array{UInt8, 4}}}, Tuple{Tuple{Array{Float32, 4}, Array{UInt8, 4}}}, Tuple{Tuple{ESDL.DAT.AllMissing}}, ESDL.DAT.InnerObj{(1,), (1,), ((1, 2, 3), (1, 2, 3)), true, (), false}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}}, Tuple{ESDL.Cubes.Mask.MaskArray{Float32, 1, Array{Float32, 1}}}, Tuple{ESDL.Cubes.Mask.MaskArray{Float32, 1, Array{Float32, 1}}}, Tuple{ESDL.Cubes.Axes._RangeAxis{Float64, :Lon, Base.StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, ESDL.Cubes.Axes._RangeAxis{Float64, :Lat, Base.StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, ESDL.Cubes.Axes.CategoricalAxis{String, :Variable, Array{String, 1}}}, Tuple{Int64, Array{Float64, 1}, Array{Int64, 1}}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
    precompile(Tuple{typeof((ESDL.DAT).getSubRange), Tuple{Array{Float32, 4}, Array{UInt8, 4}}, typeof(identity), Int})
    precompile(Tuple{typeof((ESDL.DAT).getSubRange2), ((ESDL.Cubes).Mask).MaskArray{Float32, 1, Array{Float32, 1}}, Tuple{Array{Float32, 4}, Array{UInt8, 4}}, typeof(identity), Int})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:temporal_res}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:comp_level}})
    precompile(Tuple{typeof(ESDL.CubeAPI.defaulttime), ESDL.CubeAPI.Cube, Array{Int64, 1}})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##mapCube#52")), (Base.Iterators).Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:indims, :outdims), Tuple{(ESDL.DAT).InDims, (ESDL.DAT).OutDims}}}, typeof(identity), typeof(identity), (ESDL.CubeAPI).SubCubeV{Float32, (ESDL.CubeAPI).Cube}, Int64, Int})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:spatial_res}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:grid_x0}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:grid_width}})
    precompile(Tuple{typeof(ESDL.Proc.getNpY), ESDL.CubeAPI.SubCubeV{Float32, ESDL.CubeAPI.Cube}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:grid_y0}})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:variables}})
    precompile(Tuple{typeof(ESDL.Cubes._read), ESDL.CubeAPI.SubCubeV{Float32, ESDL.CubeAPI.Cube}, Tuple{Array{Float32, 4}, Array{UInt8, 4}}, Base.IteratorsMD.CartesianIndices{4, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}}}})
    precompile(Tuple{getfield(ESDL.ESDLTools, Symbol("##s13#18")), Int, Int, Int, Int})
    precompile(Tuple{typeof(ESDL.CubeAPI.parseEntry), ESDL.CubeAPI.CubeConfig, ESDL.CubeAPI.ConfigEntry{:grid_height}})
    precompile(Tuple{typeof(ESDL.Cubes.iscompressed), ESDL.CubeAPI.Cube})
    precompile(Tuple{typeof((ESDL.DAT).mysizeof), Type{Int}})
    precompile(Tuple{typeof(((ESDL.Cubes).Axes).getOutAxis), Tuple{}, Array{(((ESDL.Cubes).Axes).CubeAxis{T, S} where S) where T, 1}, Tuple{(ESDL.CubeAPI).SubCubeV{Float32, (ESDL.CubeAPI).Cube}}, Tuple{Int64, Array{Float64, 1}, Array{Int64, 1}}, typeof(identity)})
    precompile(Tuple{typeof(ESDL.Proc.MSC.removeMSC), ESDL.CubeAPI.SubCubeV{Float32, ESDL.CubeAPI.Cube}})
    precompile(Tuple{typeof(ESDL.Cubes.iscompressed), ESDL.CubeAPI.SubCubeVPerm{Float32}})
    precompile(Tuple{typeof(ESDL.DAT.getCubeHandles), ESDL.DAT.DATConfig{1, 1}})
    precompile(Tuple{getfield(ESDL.CubeAPI, Symbol("##Cube#1")), String, Type{Int}})
    precompile(Tuple{typeof(ESDL.DAT.generateworkarrays), ESDL.DAT.DATConfig{1, 1}})
    precompile(Tuple{typeof(ESDL.CubeAPI.readAllyears), ESDL.CubeAPI.SubCubeV{Float32, ESDL.CubeAPI.Cube}, Array{Float32, 4}, Array{UInt8, 4}, Int64, Int64, Int64, Int64, Int64, Int64, Int64, Int64, Int64, Int64})
    precompile(Tuple{typeof(ESDL.DAT.updatears), ESDL.DAT.DATConfig{1, 1}, Tuple{ESDL.DAT.InputCube{4}}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}}, typeof(ESDL.Cubes._read)})
    precompile(Tuple{typeof(ESDL.DAT.generateOutCubes), ESDL.DAT.DATConfig{1, 1}})
    precompile(Tuple{typeof(ESDL.DAT.mustReorder), ESDL.CubeAPI.SubCubeV{Float32, ESDL.CubeAPI.Cube}, Array{ESDL.Cubes.Axes.CubeAxis{T, S} where S where T, 1}})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##97#99")), Int64, Int64, Int64})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("##89#94")), Array{Int64, 1}, ESDL.DAT.InputCube{4}})
    precompile(Tuple{typeof(ESDL.DAT.setworkarray), ESDL.DAT.OutputCube})
    precompile(Tuple{typeof(ESDL.Cubes.Axes.axname), Type{ESDL.Cubes.Axes._RangeAxis{Float64, :Lat, Base.StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}})
    precompile(Tuple{typeof(ESDL.Cubes._write), ESDL.Cubes.CubeMem{Float32, 4}, Tuple{Base.SubArray{Float32, 4, Array{Float32, 4}, Tuple{Base.IteratorsMD.CartesianIndices{4, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}}}}, false}, Base.SubArray{UInt8, 4, Array{UInt8, 4}, Tuple{Base.IteratorsMD.CartesianIndices{4, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}}}}, false}}, Base.IteratorsMD.CartesianIndices{4, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}, Base.UnitRange{Int64}}}})
    precompile(Tuple{typeof(ESDL.Cubes.Axes.axname), Type{ESDL.Cubes.Axes._RangeAxis{Float64, :Lon, Base.StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}})
    precompile(Tuple{typeof(ESDL.DAT.wrapWorkArray), ESDL.DAT.AsArray, ESDL.Cubes.Mask.MaskArray{Float32, 1, Array{Float32, 1}}, Array{ESDL.Cubes.Axes.CubeAxis{T, S} where S where T, 1}})
    precompile(Tuple{getfield(ESDL.DAT, Symbol("#kw##mapCube")), NamedTuple{(:indims, :outdims), Tuple{(ESDL.DAT).InDims, (ESDL.DAT).OutDims}}, typeof((ESDL.DAT).mapCube), typeof(identity), Tuple{(ESDL.CubeAPI).SubCubeV{Float32, (ESDL.CubeAPI).Cube}}, Int64, Int})
    precompile(Tuple{typeof(ESDL.DAT.distributeLoopRanges), Tuple{Int64, Int64, Int64}, Tuple{Int64, Int64, Int64}, Tuple{Int64, Int64, Int64}})
    precompile(Tuple{typeof(ESDL.Cubes.Axes.axname), Type{ESDL.Cubes.Axes.CategoricalAxis{String, :Variable, Array{String, 1}}}})
    precompile(Tuple{typeof(ESDL.DAT.setworkarray), ESDL.DAT.InputCube{4}})
    precompile(Tuple{typeof(ESDL.DAT.generateOutCube), ESDL.DAT.OutputCube, Bool, Float64, Array{Int64, 1}})
    precompile(Tuple{typeof(ESDL.DAT.mask2miss), Tuple{Base.SubArray{Float32, 1, Array{Float32, 4}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64, Int64, Int64}, true}, Base.SubArray{UInt8, 1, Array{UInt8, 4}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64, Int64, Int64}, true}}, ESDL.Cubes.Mask.MaskArray{Float32, 1, Array{Float32, 1}}})
    precompile(Tuple{typeof(ESDL.DAT.miss2mask!), Tuple{Base.SubArray{Float32, 1, Array{Float32, 4}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64, Int64, Int64}, true}, Base.SubArray{UInt8, 1, Array{UInt8, 4}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64, Int64, Int64}, true}}, ESDL.Cubes.Mask.MaskArray{Float32, 1, Array{Float32, 1}}})
end