text:
"""
    grasol(sysg, sysf[, γ]; L2sol = false, nehari = false, reltol = 0.0001, mindeg = false, poles, sdeg,
           fast = true, offset = β, atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (sysx, info)

Determine for the descriptor systems `sysg = (Ag-λEg,Bg,Cg,Dg)` and `sysf = (Af-λEf,Bf,Cf,Df)`
with the transfer function matrices `G(λ)` and `F(λ)`,
respectively, the descriptor system `sysx` with the transfer function matrix `X(λ)`
such that `X(λ)` is the approximate solution of the linear rational equation `G(λ)X(λ) = F(λ)`,
which achieves the minimum error norm ``{\\small \\text{mindist} := \\min \\|G(λ)X(λ) - F(λ)\\|}``.
The resulting `X(λ)` has all poles stable or lying on the boundary of the stability domain `Cs`.
If `L2sol = false` (default) then the `L∞`-norm optimal solution is computed, while if `L2sol = true` the
`L2`-norm optimal solution is computed.
`sysg` and `sysf` must not have poles on the boundary of the stability domain `Cs`.

If  `γ > 0` is a desired sub-optimality degree, then the `γ`-suboptimal model-matching problem

```math
     \\text{mindist} := \\|G(λ)X(λ) - F(λ) \\| < γ
```
is solved and `mindist` is the achieved suboptimal distance.

The call with

    grasol(sysgf[, mf[, γ]]; L2sol = false, nehari = false, reltol = 0.0001, mindeg = false, poles, sdeg,
           fast = true, offset = β, atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (sysx, info)

uses the compound descriptor system `sysgf = (A-λE,[Bg Bf],C,[Dg Df])`,
where `Bf` and `Df` have `mf` columns, to define
the descriptor systems `sysg = (A-λE,Bg,C,Dg)` and `sysf = (A-λE,Bf,C,Df)`
(i.e., `Ag-λEg = Af-λEf = A-λE` and `Cg = Cf = C`).
`sysgf` must not have poles on the boundary of the stability domain `Cs`.

If `nehari = true`, the optimal or suboptimal Nehari approximation is used to
compute a `L∞`-suboptimal solution of the underlying _least-distance problem_ (`LDP`).
If `nehari = false` (default), the `L∞`-optimal solution is computed using the `γ`-iteration
in the underlying `LDP` [2].

If `mindeg = true`, a minimum order solution is determined (if possible),
while if `mindeg = false` (default) a particular solution of non-minimal order is determined.

The resulting named tuple `info` contains additional information:
`info.nrank` is the normal rank of `G(λ)`,
`info.nr` is the number of freely assignable poles of the solution `X(λ)`,
`info.mindist` is the achieved approximation error norm and
`info.nonstandard` is `true` for a non-standard problem, with `G(λ)`
having zeros on the boundary of the stability domain, and `false` for
a standard problem, when `G(λ)` has no zeros on the boundary
of the stability domain.

The keyword argument `reltol` specifies the relative tolerance for the desired accuracy of
the `γ`-iteration employed to solve the underlying least-distance problem.
The iterations are performed until the current estimations of maximum ``γ_u`` and minimum ``γ_l`` of
the optimal distance satisfies
``{\\small γ_u-γ_l \\leq \\text{reltol} * \\text{gap}}``, where `gap` is the initial estimation of the error gap.

To assess the presence of poles on the boundary of the stability domain `Cs`, a boundary offset  `β`
can be specified via the keyword parameter `offset = β`.
Accordingly, for a continuous-time setting,
the boundary of `Cs` contains the complex numbers with real parts within the interval `[-β,β]`,
while for a discrete-time setting, 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 vector `poles` specified as a keyword argument, can be used to specify the desired poles
of `sysx` alternatively to or jointly with enforcing a desired stability degree `sdeg` of poles.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `Ag`, `Af`, `A`, `Bg`, `Bf`, `Cg`, `Cf`, `Dg`, `Df`,
the absolute tolerance for the nonzero elements of `Eg`, `Ef`, and the relative tolerance
for the nonzero elements of all above matrices. The default relative tolerance is `n*ϵ`, where `ϵ` is the working machine epsilon
and `n` is the order of the system `sysg`. The keyword argument `atol` can be used
to simultaneously set `atol1 = atol` and `atol2 = atol`.

The rank decisions in the underlying pencil manipulation algorithms are
based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

_Method:_  An extension of the approach of [1] to descriptor systems is used.

_References:_

[1]  B. A. Francis. A Course in H-infinity Theory,
       Vol. 88 of Lecture Notes in Control and Information Sciences,
       Springer-Verlag, New York, 1987.

[2] C.-C. Chu, J. C. Doyle, and E. B. Lee.
    The general distance problem in H∞  optimal control theory,
    Int. J. Control, vol 44, pp. 565-596, 1986.

"""
function grasol(sysgf::DescriptorStateSpace{T}, mf::Int, γ::Union{Real,Missing} = missing;
                L2sol::Bool = false, nehari::Bool = false, offset::Real = sqrt(eps(float(real(T)))),
                poles::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
                atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
                rtol::Real = (sysgf.nx*eps(real(float(one(T)))))*iszero(max(atol1,atol2)),
                reltol::Real = 0.0001, fast::Bool = true, mindeg::Bool = false) where T
   p, mgf = size(sysgf);
   (mf <= mgf && mf >= 0) || throw(DimensionMismatch("mf must be at most $mgf, got $mf"))
   m = mgf-mf
   return grasol(sysgf[:,1:m], sysgf[:,m+1:mgf], γ; L2sol = L2sol, nehari = nehari, offset = offset,
                 reltol = reltol, poles = poles, sdeg = sdeg, fast = fast,
                 mindeg = mindeg, atol1 = atol1, atol2 = atol2, rtol = rtol)
end
function grasol(sysg::DescriptorStateSpace{T1}, sysf::DescriptorStateSpace{T2},
                γ::Union{Real,Missing} = missing; L2sol::Bool = false, nehari::Bool = false,
                offset::Real = sqrt(eps(float(real(T1)))), reltol::Real = 0.0001,
                poles::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
                atol::Real = zero(real(T1)), atol1::Real = atol, atol2::Real = atol,
                rtol::Real = (sysg.nx*eps(real(float(one(T1)))))*iszero(max(atol1,atol2)),
                fast::Bool = true, mindeg::Bool = false) where {T1,T2}
   T = promote_type(T1,T2)
   T <: BlasFloat || ( T = promote_type(Float64,T) )

   p, m = size(sysg)
   pf, mf = size(sysf)
   p == pf|| throw(DimensionMismatch("sysg and sysf must have the same number of outputs"))
   Ts = promote_Ts(sysg.Ts,sysf.Ts)
   disc = !iszero(Ts)

   #compute the extended inner-quasi-outer factorization
   Gi, Go, info1 = giofac(sysg, atol1 = atol1, atol2 = atol2, atol3 = atol1, rtol = rtol,
                          fast = fast, minphase = true, offset = offset)
   ro = info1.nrank

   #detect nonstandard problem
   nonstandard = (info1.nfuz + info1.niuz > 0)

   #define and solve the LDP
   # ensure stabilizability of F
   F = gir(Gi'*sysf, atol1 = atol1, atol2 = atol2, rtol = rtol)
   if L2sol
      #set tolerance for feedthrough matrix
      atol1 > 0 ? told = atol1 : told = 1.e4*eps()
      disc || norm(F.D[ro+1:end,:],Inf) > told || (F.D[ro+1:end,:] = zeros(T1,p-ro,mf))
      #solve the H2-LDP min||[ F1-Xt; F2 ] ||_2
      Xt, Xtu = gsdec(F[1:ro,:], job = "stable")
      gopt = gl2norm(grcfid(gir([Xtu; F[ro+1:end,:] ], fast = fast, atol1 = atol1, atol2 = atol2))[1],
                     atolinf = told)
   else
      #solve the H_inf-LDP min ||[ F1-Xt; F2 ] ||_inf
      nehari || (nehari = iszero(F[ro+1:end,:]; atol1, atol2, rtol))
      Yt, gopt = glinfldp(gdual(F), p-ro, γ; nehari = nehari, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, offset = offset, reltol = reltol);
      Xt = gdual(Yt)
   end

   if ro == m
      nr = 0
      sysx = gir(Go\Xt, noseig = true, atol1 = atol1, atol2 = atol2)
   else
      sysx, info1, = grsol(Go, Xt, mindeg = mindeg, poles = poles, sdeg = sdeg,
                           atol1 = atol1, atol2 = atol2, rtol = rtol)
      if min(sysx.ny,sysx.nu) > 0
         ev = gpole(sysx, atol1 = atol1, atol2 = atol2);
         disc || (ev = ev[isfinite.(ev)])
         if !isempty(ev) && !nonstandard
            if (!disc && maximum(real.(ev)) > -offset) || (disc && maximum(abs.(ev)) > 1-offset)
               ismissing(poles) && ismissing(sdeg) && (sdeg = disc ? 0.95 : -0.05)
               sysx = grsol(Go, Xt, mindeg = false, poles = poles, sdeg = sdeg,
                            atol1 = atol1, atol2 = atol2, rtol = rtol)[1]
            end
         end
      end
      nr = info1.nr
   end

   info = (nrank = ro, nr = nr, mindist = gopt, nonstandard = nonstandard)
   return sysx, info

   #end GRASOL
end

"""
    glasol(sysg, sysf[, γ]; L2sol = false, nehari = false, reltol = 0.0001, mindeg = false, poles, sdeg,
           fast = true, offset = β, atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (sysx, info)

Determine for the descriptor systems `sysg = (Ag-λEg,Bg,Cg,Dg)` and
`sysf = (Af-λEf,Bf,Cf,Df)` with the transfer function matrices `G(λ)` and `F(λ)`,
respectively, the descriptor system `sysx` with the transfer function matrix `X(λ)`
such that `X(λ)` is the approximate solution of the linear rational equation `X(λ)G(λ) = F(λ)`,
which achieves the minimum error norm ``{\\small \\text{mindist} := \\min \\|X(λ)G(λ) - F(λ)\\|}``.
The resulting `X(λ)` has all poles stable or lying on the boundary of the stability domain `Cs`.
If `L2sol = false` (default) then the `L∞`-norm optimal solution is computed, while if `L2sol = true` the
`L2`-norm optimal solution is computed.
`sysg` and `sysf` must not have poles on the boundary of the stability domain `Cs`.

If  `γ > 0` is a desired sub-optimality degree, then the `γ`-suboptimal model-matching problem

```math
     \\text{mindist} := \\|X(λ)G(λ) - F(λ) \\| < γ
```

is solved and `mindist` is the achieved suboptimal distance.

The call with

    glasol(sysgf[, pf[, γ]]; L2sol = false, nehari = false, reltol = 0.0001, mindeg = false, poles, sdeg,
           fast = true, offset = β, atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (sysx, info)

uses the compound descriptor system `sysgf = (A-λE,B,[Cg; Cf],[Dg; Df])`,
where `Cf` and `Df` have `pf` rows, to define
the descriptor systems `sysg = (A-λE,B,Cg,Dg)` and `sysf = (A-λE,B,Cf,Df)`
(i.e., `Ag-λEg = Af-λEf = A-λE` and `Bg = Bf = B`).
`sysgf` must not have poles on the boundary of the stability domain `Cs`.

If `nehari = true`, the optimal or suboptimal Nehari approximation is used to solve the
 underlying _least-distance problem_ (`LDP`).
If `nehari = false` (default), the optimal solution is computed using the `γ`-iteration
in the underlying `LDP` [2].

If `mindeg = true`, a minimum order solution is determined (if possible),
while if `mindeg = false` (default) a particular solution of non-minimal order is determined.

The resulting named tuple `info` contains additional information:
`info.nrank` is the normal rank of `G(λ)`,
`info.nl` is the number of freely assignable poles of the solution `X(λ)`,
`info.mindist` is the achieved approximation error norm and
`info.nonstandard` is `true` for a non-standard problem, with `G(λ)`
having zeros on the boundary of the stability domain, and `false` for
a standard problem, when `G(λ)` has no zeros on the boundary
of the stability domain.

The keyword argument `reltol` specifies the relative tolerance for the desired accuracy of
the `γ`-iteration employed to solve the underlying least-distance problem.
The iterations are performed until the current estimations of maximum ``γ_u`` and minimum ``γ_l`` of
the optimal distance satisfies
``{\\small γ_u-γ_l < \\text{reltol}* \\text{gap}}``, where `gap` is the initial estimation of the error gap.

To assess the presence of poles on the boundary of the stability domain `Cs`, a boundary offset  `β`
can be specified via the keyword parameter `offset = β`.
Accordingly, for a continuous-time setting,
the boundary of `Cs` contains the complex numbers with real parts within the interval `[-β,β]`,
while for a discrete-time setting, 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 vector `poles` specified as a keyword argument, can be used to specify the desired poles
of `sysx` alternatively to or jointly with enforcing a desired stability degree `sdeg` of poles.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `Ag`, `Af`, `A`, `Bg`, `Bf`, `Cg`, `Cf`, `Dg`, `Df`,
the absolute tolerance for the nonzero elements of `Eg`, `Ef`, and the relative tolerance
for the nonzero elements of all above matrices. The default relative tolerance is `n*ϵ`, where `ϵ` is the working machine epsilon
and `n` is the order of the system `sysg`. The keyword argument `atol` can be used
to simultaneously set `atol1 = atol` and `atol2 = atol`.

The rank decisions in the underlying pencil manipulation algorithms are
based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

_Method:_  An extension of the approach of [1] to descriptor systems is used.

_References:_

[1]  B. A. Francis. A Course in H-infinity Theory,
       Vol. 88 of Lecture Notes in Control and Information Sciences,
       Springer-Verlag, New York, 1987.

[2] C.-C. Chu, J. C. Doyle, and E. B. Lee.
    The general distance problem in H∞  optimal control theory,
    Int. J. Control, vol 44, pp. 565-596, 1986.

"""
function glasol(sysgf::DescriptorStateSpace{T}, pf::Int, γ::Union{Real,Missing} = missing;
                L2sol::Bool = false, nehari::Bool = false, offset::Real = sqrt(eps(float(real(T)))),
                poles::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
                atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
                rtol::Real = (sysgf.nx*eps(real(float(one(T)))))*iszero(max(atol1,atol2)),
                reltol::Real = 0.0001, fast::Bool = true, mindeg::Bool = false) where T
   pgf, m = size(sysgf);
   (pf <= pgf && pf >= 0) || throw(DimensionMismatch("pf must be at most $pgf, got $pf"))
   p = pgf-pf
   return glasol(sysgf[1:p,:], sysgf[p+1:pgf,:], γ; L2sol = L2sol, nehari = nehari, offset = offset,
                 reltol = reltol, poles = poles, sdeg = sdeg, fast = fast,
                 mindeg = mindeg, atol1 = atol1, atol2 = atol2, rtol = rtol)
end
function glasol(sysg::DescriptorStateSpace{T1}, sysf::DescriptorStateSpace{T2},
                γ::Union{Real,Missing} = missing; L2sol::Bool = false, nehari::Bool = false,
                offset::Real = sqrt(eps(float(real(T1)))), reltol::Real = 0.0001,
                poles::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
                atol::Real = zero(real(T1)), atol1::Real = atol, atol2::Real = atol,
                rtol::Real = (sysg.nx*eps(real(float(one(T1)))))*iszero(max(atol1,atol2)),
                fast::Bool = true, mindeg::Bool = false) where {T1,T2}
   T = promote_type(T1,T2)
   T <: BlasFloat || ( T = promote_type(Float64,T) )

   p, m = size(sysg)
   pf, mf = size(sysf)
   m == mf|| throw(DimensionMismatch("sysg and sysf must have the same number of inputs"))
   # isstable(sysf, atol1 = atol1, atol2 = atol2, rtol = rtol, offset = offset) ||
   #          error("sysf must be a stable system")
   Ts = promote_Ts(sysg.Ts,sysf.Ts)
   disc = !iszero(Ts)

   #compute the extended quasi-co-outer-co-inner factorization
   Gi, Go, info1 = goifac(sysg; atol1, atol2, atol3 = atol1, rtol, fast, minphase = true, offset)
   ro = info1.nrank

   #detect nonstandard problem
   #nonstandard = order(gcrange(Go,struct('tol',tol,'zeros','s-unstable'))) > 0;
   nonstandard = (info1.nfuz + info1.niuz > 0)

   #define and solve the LDP
   # ensure detectability of F
   F = gir(sysf*Gi', atol1 = atol1, atol2 = atol2, rtol = rtol)
   if L2sol
      #set tolerance for feedthrough matrix
      atol1 > 0 ? told = atol1 : told = 1.e4*eps()
      disc || norm(F.D[:,ro+1:end],Inf) > told || (F.D[:,ro+1:end] = zeros(T1,pf,m-ro))
      #solve the H2-LDP min||[ F1-Xt F2 ] ||_2
      Xt, Xtu = gsdec(F[:,1:ro],job = "stable")
      Y = glcfid(gir([Xtu F[:,ro+1:end] ], atol1 = atol1, atol2 = atol2))[1]
      gopt = gl2norm(glcfid(gir([Xtu F[:,ro+1:end] ], fast = fast, atol1 = atol1, atol2 = atol2))[1],
                     atolinf = told)
   else
      #solve the H_inf-LDP min ||[ F1-Xt F2 ] ||_inf
      nehari || (nehari = iszero(F[:,ro+1:end]; atol1, atol2, rtol))
      Xt, gopt = glinfldp(F, m-ro, γ; nehari = nehari, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol, offset = offset, reltol = reltol);
   end

   if ro == p
      nl = 0
      sysx = gir(Xt/Go, noseig = true, atol1 = atol1, atol2 = atol2)
   else
      sysx, info1, = glsol(Go, Xt, mindeg = mindeg, poles = poles, sdeg = sdeg,
                    atol1 = atol1, atol2 = atol2, rtol = rtol)
      if min(sysx.ny,sysx.nu) > 0
         ev = gpole(sysx, atol1 = atol1, atol2 = atol2);
         disc || (ev = ev[isfinite.(ev)])
         if !isempty(ev) && !nonstandard
            if (!disc && maximum(real.(ev)) > -offset) || (disc && maximum(abs.(ev)) > 1-offset)
               ismissing(poles) && ismissing(sdeg) && (sdeg = disc ? 0.95 : -0.05)
               sysx = glsol(Go, Xt, mindeg = false, poles = poles, sdeg = sdeg,
                            atol1 = atol1, atol2 = atol2, rtol = rtol)[1]
            end
         end
      end
      nl = info1.nl
   end

   info = (nrank = ro, nl = nl, mindist = gopt, nonstandard = nonstandard)
   return sysx, info

   #end GLASOL
end
"""
    glinfldp(sys1, sys2, [, γ]; nehari = false, reltol = 0.0001, fast = true, offset = β,
             atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (sysx, mindist)

Determine for the descriptor systems `sys1 = (A1-λE1,B1,C1,D1)` and
`sys2 = (A2-λE2,B2,C2,D2)` with the transfer function matrices ``G_1(λ)`` and ``G_2(λ)``,
respectively, the descriptor system `sysx` with the transfer function matrix ``X(λ)``
such that ``X(λ)`` is the solution of the 2-block `L∞` _least distance problem_ (`LDP`)

```math
      \\text{mindist} := \\min \\|G_1(λ)-X(λ) \\mid   G_2(λ) \\|_\\infty
```

`mindist` is the achieved minimum distance corresponding to the optimal solution.
If `sys2 = []`, an 1-block `LDP` is solved.
`sys1` and `sys2` must not have poles on the boundary of the stability domain `Cs`.
It is assumed that `sys2` has no poles on the boundary of the stability domain (see below).
The resulting solution is stable, provided `sys1` has no poles on the boundary
of the stability domain as well.

If  ``{\\small γ > \\|G_2(λ)\\|_\\infty}`` is a desired sub-optimality degree, then the
`γ`-suboptimal `LDP`

```math
     \\text{mindist} := \\|G_1(λ)-X(λ) \\mid G_2(λ) \\|_\\infty < γ
```

is solved and `mindist` is the achieved suboptimal distance.

The call with

    glinfldp(sys[, m2[, γ]]; nehari = false, fast = true, offset = β,
             atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (sysx, mindist)

uses the compound descriptor system `sys = (A-λE,[B1 B2],C,[D1 D2])`,
where `B2` has `m2` columns, to define
the descriptor systems `sys1 = (A-λE,B1,C,D1)` and `sys2 = (A-λE,B2,C,D2)`
(i.e., `A1-λE1 = A2-λE2 = A-λE` and `C1 = C2 = C`).
If `m2 = 0`, an 1-block `LDP` is solved.
`sys2` must not have poles on the boundary of the stability domain `Cs`.

If `nehari = true`, the optimal or suboptimal Nehari approximation is used to solve the `LDP`.
If `nehari = false` (default), the optimal solution is computed using the `γ`-iteration [1].

The keyword argument `reltol` specifies the relative tolerance for the desired accuracy of `γ`-iteration.
The iterations are performed until the current estimations of maximum ``γ_u`` and minimum ``γ_l`` of
the optimal distance ``γ_o``, ``γ_l \\leq γ_o \\leq γ_u``, satisfies

```math
     γ_u-γ_l \\leq \\text{reltol} * \\text{gap} ,
```
where `gap` is the original gap (internally determined).

To assess the presence of poles on the boundary of the stability domain `Cs`, a boundary offset  `β`
can be specified via the keyword parameter `offset = β`.
Accordingly, for a continuous-time setting,
the boundary of `Cs` contains the complex numbers with real parts within the interval `[-β,β]`,
while for a discete-time setting, 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`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` 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 decisions in the underlying pencil manipulation algorithms are
based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

_Method:_ The approach of [1] is used for the solution of the 2-block least distance problem.
Necessary conditions for solvability of the LDP with a stable solution is that `sys1` and `sys2` have no poles
on the boundary of the stability domain `Cs`.

_References:_
[1] C.-C. Chu, J. C. Doyle, and E. B. Lee
    The general distance problem in H∞  optimal control theory,
    Int. J. Control, vol 44, pp. 565-596, 1986.
"""
function glinfldp(sys1::DescriptorStateSpace{T1}, sys2::DescriptorStateSpace{T2},
                  γ::Union{Real,Missing} = missing; nehari::Bool = false, offset::Real = sqrt(eps(float(real(T1)))),
                  atol::Real = zero(real(T1)), atol1::Real = atol, atol2::Real = atol,
                  rtol::Real = (sys1.nx*eps(real(float(one(T1)))))*iszero(max(atol1,atol2)),
                  reltol::Real = 0.0001, fast::Bool = true) where {T1,T2}
   size(sys1,1) == size(sys2,1) || throw(DimensionMismatch("sys1 and sys2 must have the same number of outputs"))
   return glinfldp(gir([sys1 sys2]; fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol), size(sys2,2),γ;
                   nehari = nehari, offset = offset, reltol = reltol, fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
end
function glinfldp(sys::DescriptorStateSpace{T}, m2::Int, γ::Union{Real,Missing} = missing;  nehari::Bool = false,
                  offset::Real = sqrt(eps(float(real(T)))), reltol::Real = 0.0001,
                  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) where T

   T1 = T <: BlasFloat ? T : promote_type(Float64,T)
   p, m = size(sys);
   (m2 <= m && m2 >= 0) || throw(DimensionMismatch("m2 must be at most $m, got $m2"))
   m1 = m-m2
   tol = sqrt(eps(real(float(one(T1)))))
   disc = (sys.Ts != 0)

   # set stability margin to separate stable and unstable parts
   smarg = disc ? 1+offset : offset

   if nehari
      # compute a suboptimal solution as a Nehari approximation of G1
      sysx, s1 = gnehari(sys[:,1:m1], γ; offset = -offset, atol1, atol2, rtol, fast)
      mindist = glinfnorm(gir(sys-sysx*eye(m1,m); atol1, atol2, rtol, fast); offset, atol1, atol2, rtol)[1]
      isinf(mindist) && error("sys2 has poles on the boundary of stability domain")
      return sysx, mindist
   end

   subopt = !ismissing(γ)
   gl1 = m2 > 0 ? (glinfnorm(sys[:,m1+1:end]; offset, atol1, atol2, rtol)[1]) : 0 # prevent failure of opnorm
   isinf(gl1) && error("sys2 has poles on the boundary of stability domain")

   # address constant case
   if order(sys) == 0 || m1 == 0
      # solve a 2-block Parrott problem
      # handle opnorm bug for 1 x 0 matrices
      return sys[:,1:m1], m1 == m ? 0 : opnorm(sys.D[:,m1+1:end])
   end

   subopt && gl1 > γ && error("γ must be chosen greater than $gl1")
   if disc
      hlu = ghanorm(gsdec(sys[:,1:m1]; job = "unstable", smarg, atol1, atol2, rtol, fast)[1]')[1]
   else
      hlu = ghanorm(gsdec(gsdec(sys[:,1:m1]; job = "finite", atol1, atol2, rtol, fast)[1]; job = "unstable", smarg, atol1, atol2, rtol, fast)[1]')[1]
   end
   if m2 == 0 || hlu <= atol || gl1 <= tol*hlu
      # solve one block LDP
      if subopt
         # solve sub-optimal Nehari problem
         sysx,  = gnehari(sys[:,1:m1], γ; offset, atol1, atol2, rtol, fast)
         mindist = glinfnorm(sys[:,1:m1]-sysx; offset, atol1, atol2, rtol)[1]
      else
        # solve optimal Nehari problem
        sysx, mindist = gnehari(sys[:,1:m1]; offset = -offset, atol1, atol2, rtol, fast)
      end
   else
      # solve 2-block LDP
      if subopt
         # solve sub-optimal LDP
         gam = max(γ,gl1*1.01);
         V = glsfg(sys[:,m1+1:end], gam; stabilize = true, offset, atol1, atol2, rtol, fast)
         syse = V\sys[:,1:m1]
         hlu = ghanorm(gsdec(syse; job = "unstable", smarg, atol1, atol2, rtol, fast)[1]')[1];
         # solve the sub-optimal Nehari problem
         if hlu < 1
            sysx =  gbalmr(V*gnehari(syse, hlu; offset = -offset, atol1, atol2, rtol, fast)[1];
                           offset, atol1, atol2, rtol, fast)[1];
         else
            subopt = false;
         end
         mindist = glinfnorm(sys-sysx*eye(m1,m); offset, atol1, atol2, rtol)[1]
      end
      if !subopt
         # solve optimal LDP using gamma-iteration
         # initialize the gamma-iteration
         # compute upper bound and inital gap
         # hlu = ghanorm(gsdec(sys[:,1:m1], job = "unstable", atol1 = atol1, atol2 = atol2, rtol = rtol, fast = fast)[1]')[1];
         gu = hypot(hlu,gl1); gl = max(gl1,hlu)
         gap = max(1,gu-gl);
         gam = (gl+gu)/2;
         # stabilize G2 only once
         G2s = grcfid(sys[:,m1+1:end], mininf = true, fast = fast, offset = offset,
                      atol1 = atol1, atol2 = atol2, atol3 = atol1, rtol = rtol)[1]
         # perform the gamma-iteration
         iter = 0; iterl = 0; g0 = gam; hlu = 1;
         while gu-gl > reltol*gap && iter <= 50
            iter += 1
            V = glsfg(G2s, gam, stabilize = false, fast = fast, offset = offset,
                      atol1 = atol1, atol2 = atol1, rtol = rtol)
            # fail = false
            # V = try
            #    V = glsfg(G2s, gam, stabilize = false, fast = fast, offset = offset,
            #           atol1 = atol1, atol2 = atol1, rtol = rtol)
            #    g0 = gam
            #    fail = false
            #    V
            # catch
            #    if iter > 1
            #       gam = g0
            #       V = glsfg(G2s, gam, stabilize = false, fast = fast, offset = offset,
            #                 atol1 = atol1, atol2 = atol1, rtol = rtol)
            #       fail = true
            #       @warn "Desired accuracy of γ-iteration not achievable"
            #    else
            #       error("Failure of the minimum-phase spectral factorization")
            #    end
            # end

            # compute inv(V)*G1
            syse = gminreal(V\sys[:,1:m1],atol=1.e-7)
            hlu = ghanorm(gsdec(syse; job = "unstable", smarg, atol1, atol2, rtol, fast)[1]')[1];
            #fail && break
            if hlu < 1
               gu = gam;           # gamma > gamma_opt
            else
               gl = gam;           # gamma <= gamma_opt
               iterl += 1
            end
            gam = (gl+gu)/2;
         end
         if iter > 0
            if hlu >= 1
              # recompute spectral factorization if gamma < gamma_opt
              V = glsfg(G2s, gu, stabilize = false)
              # compute inv(V)*G1
              syse = V\sys[:,1:m1]
              gam = gu
            end
         else
            # compute spectral factorization for increased gamma
            gam = min(gu, gam*(1+reltol))
            V = glsfg(G2s, gam, stabilize = false, offset = offset,
                     atol1 = atol1, atol2 = atol1, rtol = rtol)
            # compute inv(V)*G1
            syse = V\sys[:,1:m1]
         end
         if iterl == 0 && gam - gl < reltol*gl1
            # handle the case gam = norm(G2)
            @warn "Optimal LDP solution has poles near to the boundary of stability domain"
            # set offset to machine precision
            offset = eps(float(real(T1)))
         end

         # solve the Nehari problem for the unstable part and compute solution
         # t1, t2 = gsdec(syse, job = "unstable", fast = fast, atol1 = atol1, atol2 = atol2, rtol = rtol)
         # tt1 = gnehari(t1; offset, atol1, atol2, rtol, fast)[1]
         # sysx = gminreal(V*(tt1+t2); atol1, atol2, rtol, fast)
         # sysx = gminreal(V*gnehari(gminreal(syse,atol=1.e-7), offset = offset, atol1 = atol1, atol2 = atol2, rtol = rtol, fast = fast)[1],
         #                atol1 = atol1, atol2 = atol2, rtol = rtol, fast = fast);
         sysx = gminreal(V*gnehari(syse; offset = -offset, atol1, atol2, rtol, fast)[1]; atol1, atol2, rtol, fast)
         mindist = glinfnorm(sys-sysx*eye(m1,m), offset = offset, atol1 = atol1, atol2 = atol2, rtol = rtol)[1]

      end
   end
   return sysx, mindist
   # end of GLINFLDP
end
glinfldp(sys::DescriptorStateSpace{T}, γ::Union{Real,Missing} = missing;  kwargs...) where T =
         glinfldp(sys, 0, γ; kwargs...)
"""
    gnehari(sys[, γ]; fast = true, offset = β,
                    atol = 0, atol1 = atol, atol2 = atol, rtol = n*ϵ) -> (sysx, σ1)

Compute for the descriptor system `sys = (A-λE,B,C,D)` with the transfer function matrix `G(λ)`
the optimal or suboptimal Nehari approximation `sysx = (Ax-λEx,Bx,Cx,Dx)`
with the transfer function matrix `X(λ)`. The optimal Nehari approximation  `X(λ)` satisfies

```math
     \\| G(\\lambda) - X(\\lambda) \\|_\\infty = \\| G^{*}_u(\\lambda) \\|_H := \\sigma_1,
```

where ``{\\small G_u(\\lambda)}`` is the proper antistable part of `G(λ)`.
The resulting ``σ_1`` is the Hankel-norm of ``{\\small G^{*}_u(\\lambda)}``
(also the L∞-norm of the optimal approximation error).
For a given ``γ > σ_1``, the suboptimal approximation satisfies

```math
     \\| G(\\lambda) - X(\\lambda) \\|_\\infty \\leq \\gamma .
```
The Nehari approximation is stable if  `G(λ)` has no poles on the boundary of the stability domain
`Cs`. In the continuous-time, `Cs` is the set of complex numbers with negative real parts and its boundary
is the extended imaginary axis (containig also the poit at infinity),
while, in the discrete-time case, `Cs` is the set complex number with moduli less than one and its boundary
is the unit circle centered in the origin.

To assess the presence of poles in `Cs`, a boundary offset  `β`
can be specified via the keyword parameter `offset = β`.
Accordingly, for a continuous-time system, the stable poles in
`Cs` have real parts less than or equal to `β`,
while for a discrete-time system, they have moduli less than or equal to 1+β`.
The default value used for `β` is `sqrt(ϵ)`, where `ϵ` is the working machine precision.
For a negative values of `β` (e.g., `β = -sqrt(ϵ)`), an extended stability domain corresponding to the closure of `Cs`
is used instead of `Cs`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` 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 separation of the finite and infinite eigenvalues is performed using
rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

_Method:_ Let `G(λ) = Gs(λ) + Gu(λ)` be an additive decomposition of `G(λ)`
such that `Gs(λ)` has only poles in the closure of `Cs` (or its closure if `β < 0`)
and `Gu(λ)` is the antistable part having only
unstable poles.   The Hankel-norm approximation methods of [1] and [2],
with extensions for descriptor systems, are used to compute a stable Nehari approximation `Gn(λ)`
of the unstable part `Gu(λ)` and the resulting solution is computed as `X(λ) = Gs(λ) + Gn(λ)`.

_References:_

[1] K. Glover. All optimal Hankel-norm approximations of linear
       multivariable systems and their L∞ error bounds,
       Int. J. Control, vol. 39, pp. 1115-1193, 1984.

[2] M. G. Safonov, R. Y. Chiang, and D. J. N. Limebeer.
       Optimal Hankel model reduction for nonminimal systems.
       IEEE Trans. Automat. Control, vol. 35, pp. 496–502, 1990.
"""
function gnehari(sys::DescriptorStateSpace{T}, γ::Union{Real,Missing} = missing;
                 offset::Real = sqrt(eps(float(real(T)))), atol::Real = zero(float(real(T))),
                 atol1::Real = atol, atol2::Real = atol, atol3::Real = atol,
                 rtol::Real = sys.nx*eps(real(float(one(T))))*iszero(max(atol1,atol2,atol3)),
                 fast::Bool = true) where T
   sys.nx == 0 && (return sys, 0)
   subopt_flg = !ismissing(γ)

   ONE = one(T)
   disc = (sys.Ts != 0)

   # set stability margin to separate stable and unstable parts
   smarg = disc ? 1-offset : -offset

   # Stable/antistable separation; possible eigenvalues on the boundary of
   # the stability domain are included in the unstable part
   #sysx, sysu = gsdec(sys; job = "stable", smarg = smarg, fast = fast,  atol1 = atol1, atol2 = atol2, rtol = rtol)
   if disc
      sysu, sysx = gsdec(sys; job = "unstable", smarg, fast, atol1, atol2, rtol)
   else
      tfin, t1 = gsdec(sys; job = "finite", fast, atol1, atol2, rtol)
      sysu, t2 = gsdec(tfin; job = "unstable", smarg, fast, atol1, atol2, rtol)
      sysx = t1+t2
      sysx.D[:,:] = sysx.D + sysu.D
      sysu.D[:,:] = zeros(T, sys.ny, sys.nu)
   end
   s2eps = sqrt(eps(real(float(one(T)))))
   # only proper systems are handled in the continuous-time case
   # if !disc && sysu.E != I && (norm(sysu.E) < atol2 || rcond(sysu.E) < s2eps)
   #    sysu = gss2ss(sysu, atol2 = s2eps)[1]
   #    if sysu.E != I && rcond(sysu.E) < s2eps
   #        error("gnehari: improper continuous-time system sys")
   #    else
   #        sysx.D[:,:] = sysx.D + sysu.D
   #        sysu.D[:,:] = zeros(T, sys.ny, sys.nu)
   #    end
   # end

   if offset > 0
      #e == I ? (isqtriu(a) ? ev = ordeigvals(a) : ev = eigvals(a)) : ev = ordeigvals(a,e)[1]
      sysx.E == I ? (isqtriu(sysx.A) ? ev = ordeigvals(sysx.A) : ev = eigvals(sysx.A)) : ev = ordeigvals(sysx.A,sysx.E)[1]
      if disc
         any(abs.(ev) .>= 1-offset) &&
             error("gnehari: the system has possibly poles on the unit circle")
      else
         any(isinf.(ev)) && error("gnehari: the system is possibly improper")
         any(real.(ev) .>= -offset) &&
             error("gnehari: the system has possibly poles on the imaginary axis")
      end
   end

   # exit if the system has only stable part
   order(sysu) == 0 && (return sysx, zero(float(real(T))))

   a, e, b, c, d = dssdata(sysu)

   if disc
      # use the bilinear transformation to compute an equivalent unstable,
      # possibly non-minimal, continuous-time system
      sqrt2 = sqrt(2)
      esave = copy(e)
      e += a; a -= esave
      c = c/e
      mul!(d,c,b,-ONE,ONE)
      lmul!(sqrt2,b); c = (sqrt2*c)*esave;
   end
   # compute a balanced minimal realization of antistable part and the
   # Hankel singular values
   sysr, hs = gbalmr(dss(-a,e,-b,c,d), balance = true, atolhsv = s2eps, atol1 = atol1, atol2 = atol2, rtol = rtol, fast = fast)
   a, e, b, c, d = dssdata(sysr)

   na = size(a,1)
   p, m = size(d)

   s1 = hs[1];  # expected L-inf norm of the optimal approximation error

   # Determine the type of required approximation
   if subopt_flg
      if abs(γ-s1) <= s1*s2eps
         subopt_flg = false; # enforce optimal approximation
      elseif γ < s1
        error("Suboptimal antistable Hankel approximation cannot be computed")
      end
   else
      γ = s1
   end

   if subopt_flg
      sk = γ; r = 0
   else
      sk = s1; ep = s1*s2eps;
      if na > 0
         r = 1;
         for i = 2:na
             abs(sk-hs[i]) > ep && break
             r += 1;
         end
      else
         r = 0;
      end
   end
   #          ^
   #  Compute F(s)
   ns = na-r; sk2 = sk*sk; i1 = r+1:na;
   # make system square by padding with zeros rows or column
   Tw = eltype(a)
   pm = max(p,m)
   if m < p
       b = [b zeros(Tw,na,p-m)]
   elseif m > p
       c = [c; zeros(Tw,m-p,na)]
   end
   if subopt_flg
      u = zeros(Tw,pm,pm);
      b1 = b; c1 = c;
   else
      i2 = 1:r;
      SV1 = svd(view(b,i2,:),full=true); # b[i2,:] = U1*Σ1*V1'
      SV2 = svd(view(c,:,i2)*SV1.U,full=true)
      #sig = [diag(SV2.Vt) .< 0; falses(pm-r)]
      sig = [real.(diag(SV2.Vt)) .< 0; falses(pm-r)]
      SV2.U[:,sig] = -SV2.U[:,sig]
      u = -SV2.U*SV1.Vt
      b1 = b[i1,:]; c1 = c[:,i1]
   end
   jm = 1:m; jp = 1:p;
   if ns > 0
      # optimal solution
      hsd = Diagonal(view(hs,i1))
      e = Matrix(hsd*hsd - sk2*I)
      a = -(sk2*a[i1,i1]' + hsd*a[i1,i1]*hsd - sk*c1'*u*b1')
      b = -(hsd*b1[:,jm] + sk*c1'*u[:,jm])
      c = c1[jp,:]*hsd + sk*u[jp,:]*b1'
   else
      a = zeros(Tw,0,0); e = I; b = zeros(Tw,0,m); c = zeros(Tw,p,0);
   end
   if na > 0
      d -= sk*u[jp,jm]
   end

   if disc
      # compute the equivalent discrete-time solution
      esave = copy(e)
      e -= a; a += esave
      c = c/e
      mul!(d,c,b,ONE,ONE)
      lmul!(sqrt2,b); c = sqrt2*c*esave;
   end

   return sysx + dss(a, e, b, c, d, Ts = sys.Ts), s1

   # end GNEHARI
end


text:
"""
    MultiScaleDeepEquilibriumNetwork(main_layers::Tuple, mapping_layers::Matrix, solver;
                                     post_fuse_layers::Union{Tuple,Nothing}=nothing, p=nothing,
                                     sensealg=SteadyStateAdjoint(0.1f0, 0.1f0, 10), kwargs...)

Multiscale Deep Equilibrium Network as proposed in [baimultiscale2020](@cite)

## Arguments

* `main_layers`: Tuple of Explicit Neural Networks. The first network needs to take 2 inputs, the other ones only take 1 input
* `mapping_layers`: Matrix of Explicit Neural Networks. The ``(i, j)^{th}`` network takes the output of ``i^{th}`` `main_layer`
                    and passes it to the ``j^{th}`` `main_layer`
* `solver`: Solver for the optimization problem (See: [`ContinuousDEQSolver`](@ref) & [`DiscreteDEQSolver`](@ref))
* `post_fuse_layers`: Tuple of Explicit Neural Networks. Applied after the `mapping_layers` (Default: `nothing`)
* `p`: Optional parameters for the `model`
* `sensealg`: See [`SteadyStateAdjoint`](@ref)
* `kwargs`: Additional Parameters that are directly passed to `solve`

## Example

```julia
model = MultiScaleDeepEquilibriumNetwork(
    (
        Parallel(+, Dense(4, 4, tanh_fast), Dense(4, 4, tanh_fast)),
        Dense(3, 3, tanh_fast), Dense(2, 2, tanh_fast),
        Dense(1, 1, tanh_fast)
    ),
    [
        NoOpLayer() Dense(4, 3, tanh_fast) Dense(4, 2, tanh_fast) Dense(4, 1, tanh_fast);
        Dense(3, 4, tanh_fast) NoOpLayer() Dense(3, 2, tanh_fast) Dense(3, 1, tanh_fast);
        Dense(2, 4, tanh_fast) Dense(2, 3, tanh_fast) NoOpLayer() Dense(2, 1, tanh_fast);
        Dense(1, 4, tanh_fast) Dense(1, 3, tanh_fast) Dense(1, 2, tanh_fast) NoOpLayer()
    ],
    ContinuousDEQSolver(VCABM3(); abstol=0.01f0, reltol=0.01f0),
)

model(rand(Float32, 4, 1))
```

See also: [`DeepEquilibriumNetwork`](@ref), [`SkipDeepEquilibriumNetwork`](@ref), [`MultiScaleSkipDeepEquilibriumNetwork`](@ref)
"""
struct MultiScaleDeepEquilibriumNetwork{N,M1<:Parallel,M2<:Union{Chain,FChain},RE1,RE2,P,A,K,S} <:
       AbstractDeepEquilibriumNetwork
    main_layers::M1
    mapping_layers::M2
    main_layers_re::RE1
    mapping_layers_re::RE2
    p::P
    ordered_split_idxs::NTuple{N,Int}
    solver::A
    kwargs::K
    sensealg::S
    stats::DEQTrainingStats

    function MultiScaleDeepEquilibriumNetwork(main_layers::Parallel, mapping_layers::Union{Chain,FChain}, re1, re2,
                                              p, ordered_split_idxs, solver::A, kwargs::K, sensealg::S, stats) where {A,K,S}
        @assert length(mapping_layers) == 2
        @assert mapping_layers[1] isa MultiParallelNet

        p_main_layers, re_main_layers = destructure_parameters(main_layers)
        p_mapping_layers, re_mapping_layers = destructure_parameters(mapping_layers)

        ordered_split_idxs = tuple(cumsum([0, length(p_main_layers), length(p_mapping_layers)])...)

        p = p === nothing ? vcat(p_main_layers, p_mapping_layers) : convert(typeof(p_main_layers), p)

        return new{length(ordered_split_idxs),
                   typeof.((main_layers, mapping_layers, re_main_layers, re_mapping_layers, p))...,
                   A,K,S}(main_layers, mapping_layers, re_main_layers, re_mapping_layers, p, ordered_split_idxs,
                          solver, kwargs, sensealg, stats)
    end
end

Flux.@functor MultiScaleDeepEquilibriumNetwork

function MultiScaleDeepEquilibriumNetwork(main_layers::Tuple, mapping_layers::Matrix, solver;
                                          post_fuse_layers::Union{Tuple,Nothing}=nothing, p=nothing,
                                          sensealg=SteadyStateAdjoint(0.1f0, 0.1f0, 10), kwargs...)
    mapping_layers = if post_fuse_layers === nothing
        @assert size(mapping_layers, 1) == size(mapping_layers, 2) == length(main_layers)
        FChain(MultiParallelNet(Parallel.(+, map(x -> tuple(x...), eachcol(mapping_layers)))...), NoOpLayer())
    else
        @assert size(mapping_layers, 1) == size(mapping_layers, 2) == length(main_layers) == length(post_fuse_layers)
        FChain(MultiParallelNet(Parallel.(+, map(x -> tuple(x...), eachcol(mapping_layers)))...),
               Parallel(flatten_merge, post_fuse_layers...))
    end

    main_layers = Parallel(flatten_merge, main_layers...)

    return MultiScaleDeepEquilibriumNetwork(main_layers, mapping_layers, nothing, nothing, p,
                                            nothing, solver, kwargs, sensealg, DEQTrainingStats(0))
end

function (mdeq::MultiScaleDeepEquilibriumNetwork)(x::AbstractArray{T}) where {T}
    current_nfe = mdeq.stats.nfe

    z = zero(x)
    initial_conditions = Zygote.@ignore map(l -> l(z), map(l -> l.layers[1], mdeq.mapping_layers[1].layers))
    u_sizes = Zygote.@ignore size.(initial_conditions)
    u_split_idxs = Zygote.@ignore vcat(0, cumsum(length.(initial_conditions) .÷ size(x, ndims(x)))...)
    u0 = Zygote.@ignore vcat(Flux.flatten.(initial_conditions)...)

    N = length(u_sizes)
    update_is_variational_hidden_dropout_mask_reset_allowed(false)

    function dudt_(u, _p)
        mdeq.stats.nfe += 1

        uₛ = split_array_by_indices(u, u_split_idxs)
        p1, p2 = split_array_by_indices(_p, mdeq.ordered_split_idxs)

        u_reshaped = ntuple(i -> reshape(uₛ[i], u_sizes[i]), N)

        main_layers_output = mdeq.main_layers_re(p1)((u_reshaped[1], x), u_reshaped[2:end]...)

        return mdeq.mapping_layers_re(p2)(main_layers_output)
    end

    dudt(u, _p, t) = vcat(Flux.flatten.(dudt_(u, _p))...) .- u

    ssprob = SteadyStateProblem(dudt, u0, mdeq.p)
    res = solve(ssprob, mdeq.solver; u0=u0, sensealg=mdeq.sensealg, mdeq.kwargs...).u

    x_ = dudt_(res, mdeq.p)

    residual = Zygote.@ignore Tuple(map((iu) -> reshape(iu[2], u_sizes[iu[1]]),
                                    enumerate(split_array_by_indices(dudt(res, mdeq.p, nothing), u_split_idxs))))
    update_is_variational_hidden_dropout_mask_reset_allowed(true)

    return x_, DeepEquilibriumSolution(x_, initial_conditions, residual, T(0), mdeq.stats.nfe - current_nfe)
end

ext: jl


text:
struct Mesh
    vertices::Array{Float64, 2}
    topology::Dict{Tuple{Int8, Int8}, Array{Int64, 2}}
    reference_cell
end

"""

    unit_square_mesh(nx, ny[,cell=ufl.triangle])

Build unit square mesh with lexicographic vertex ordering.

"""
function unit_square_mesh(nx::Int64, ny::Int64; cell=ufl.triangle)::Mesh
    vertices = [[x, y] for x=LinRange(0, 1, nx + 1), y=LinRange(0, 1, ny + 1)]
    vertices = collect(transpose(hcat(vertices...)))

    if cell == ufl.triangle
        cell_vert_conn = zeros(Int64, 2 * nx * ny, 3)

        for iy in 1:ny
            for ix in 1:nx
                v0 = (iy - 1) * (nx + 1) + ix
                v1 = v0 + 1
                v2 = v0 + (nx + 1)
                v3 = v1 + (nx + 1)

                c0 = 2 * ((iy - 1) * nx + ix - 1)
                cell_vert_conn[c0 + 1, :] = [v0, v1, v2]
                cell_vert_conn[c0 + 2, :] = [v1, v2, v3]
            end
        end

        ref_cell = FIAT.reference_element.ufc_cell("triangle")

    elseif cell == ufl.quadrilateral
        cell_vert_conn = zeros(Int64, nx * ny, 4)

        for iy in 1:ny
            for ix in 1:nx
                v0 = (iy - 1) * (nx + 1) + ix
                v1 = v0 + 1
                v2 = v0 + (nx + 1)
                v3 = v1 + (nx + 1)

                cell_vert_conn[(iy - 1) * nx + ix, :] = [v0, v1, v2, v3]
            end
        end

        ref_cell = FIAT.reference_element.ufc_cell("quadrilateral")
    end

    return Mesh(vertices, Dict((2, 0) => cell_vert_conn), ref_cell)
end

"""Number of mesh entities of given dimension"""
function num_entities(mesh::Mesh, dim::Int64)::Int64
    if dim == 0
        return length(mesh.vertices[:, 1])
    else
        return length(mesh.topology[(dim, 0)][:, 1])
    end
end

"""Compute connectivity between (tdim, d) and (d, 0) mesh entities"""
function compute_connectivity_tdim_d_0!(mesh::Mesh, d::Int64)
    tdim = mesh.reference_cell.get_dimension()::Int64
    cell_vertex_conn = mesh.topology[(tdim, 0)]

    # Get local (FIAT) connectivity
    fiat_connectivity = mesh.reference_cell.get_connectivity()
    ent_vert_conn_local = fiat_connectivity[(d, 0)]
    cell_ent_conn_local  = fiat_connectivity[(tdim, d)]

    ent_per_cell = length(ent_vert_conn_local)
    vertices_per_ent = length(ent_vert_conn_local[1])
    num_cells = num_entities(mesh, tdim)

    # Allocate space for global connectivities
    ent_vert_conn = zeros(Int64, num_cells * ent_per_cell, vertices_per_ent)
    cell_ent_conn = zeros(Int64, num_cells, ent_per_cell)

    for cell_idx in 1:num_cells
        for ent_idx in 1:(ent_per_cell)
            for vert_idx in 1:(vertices_per_ent)
                # Get local index for vertex on this entity
                # +1 for python --> Julia indexing shift
                vert_for_entity = ent_vert_conn_local[ent_idx][vert_idx] + 1
                # Fill global connectivities
                ent_vert_conn[(cell_idx - 1) * ent_per_cell + ent_idx, vert_idx] = cell_vertex_conn[cell_idx, vert_for_entity]
            end
        end
    end

    ent_vert_conn, cell_ent_conn = unique_inverse(ent_vert_conn)
    cell_ent_conn = reshape(cell_ent_conn, num_cells, ent_per_cell)

    mesh.topology[(tdim, d)] = cell_ent_conn
    mesh.topology[(d, 0)] = ent_vert_conn
end

function unique_inverse(arr::Array{Int64, 2})

    inverse = Dict()
    for (i, row) in enumerate(eachrow(arr))
        if haskey(inverse, row)
            # If the row is already in keys then append to the
            # indices under which it is found
            push!(inverse[row], i)
        else
            inverse[row] = [i]
        end
    end

    # Populate arrays
    unique = zeros(Int64, length(inverse), size(arr)[2])
    unique_inverse = zeros(Int64, size(arr)[1])

    for (k, (row, I)) in enumerate(inverse)
        unique[k, :] = row
        for iI in I
            unique_inverse[iI] = k
        end
    end

    return (unique, unique_inverse)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -29.23432407510253 <= q <= 25.005662503891426)

                     add_NL_constraint(m, :(sigmoid(-0.9584800931585375 + 0.8291526833328144*sigmoid(0.45833851557375827 + -0.07893784992449904*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.05269389376863742*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.9157122134504139*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.14234137821989812*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + -0.47308336785041494*sigmoid(0.8182392255084516 + 0.5838043616725095*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.4648893036822681*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.819336300807536*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.1444146117738878*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + 0.6108443940953072*sigmoid(-0.7015763069538381 + -0.6574215519391942*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.13767395548269734*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.21605537563661414*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.3991258949611094*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + -0.4993930565410927*sigmoid(-0.21931284624450997 + 0.9567749529797225*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.8630566689068395*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.8170122079846647*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.11471605851894884*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))) + sigmoid(0.475745519197857 + -0.20728788527253217*sigmoid(0.45833851557375827 + -0.07893784992449904*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.05269389376863742*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.9157122134504139*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.14234137821989812*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + 0.5536591664216424*sigmoid(0.8182392255084516 + 0.5838043616725095*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.4648893036822681*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.819336300807536*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.1444146117738878*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + 0.5772469178059145*sigmoid(-0.7015763069538381 + -0.6574215519391942*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.13767395548269734*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.21605537563661414*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.3991258949611094*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + 0.14258534816037427*sigmoid(-0.21931284624450997 + 0.9567749529797225*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.8630566689068395*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.8170122079846647*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.11471605851894884*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))) + sigmoid(-0.19885854025662475 + 0.9363781836864602*sigmoid(0.45833851557375827 + -0.07893784992449904*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.05269389376863742*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.9157122134504139*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.14234137821989812*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + -0.2868856241387201*sigmoid(0.8182392255084516 + 0.5838043616725095*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.4648893036822681*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.819336300807536*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.1444146117738878*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + 0.6989952318088708*sigmoid(-0.7015763069538381 + -0.6574215519391942*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.13767395548269734*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.21605537563661414*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.3991258949611094*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + 0.2430276118759056*sigmoid(-0.21931284624450997 + 0.9567749529797225*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.8630566689068395*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.8170122079846647*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.11471605851894884*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))) + sigmoid(-0.4748268724301301 + -0.943475605090133*sigmoid(0.45833851557375827 + -0.07893784992449904*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.05269389376863742*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.9157122134504139*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.14234137821989812*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + -0.504486324487766*sigmoid(0.8182392255084516 + 0.5838043616725095*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.4648893036822681*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.819336300807536*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.1444146117738878*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + -0.8621266685661562*sigmoid(-0.7015763069538381 + -0.6574215519391942*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.13767395548269734*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.21605537563661414*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.3991258949611094*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2])))) + -0.8772206856506037*sigmoid(-0.21931284624450997 + 0.9567749529797225*sigmoid(-0.22994750449881307 + -0.3590609593879259*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.10337862983418544*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.321850031280575*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.612120130564894*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.8630566689068395*sigmoid(0.19384674006188307 + -0.9239566552941776*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + -0.4571782009479768*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + -0.7401514915717242*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.3444562251116623*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + 0.8170122079846647*sigmoid(0.4065836733363821 + 0.7200826994240335*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.11568748742914314*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.5958961905666036*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + -0.0402986153245215*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))) + -0.11471605851894884*sigmoid(-0.18406387526761447 + 0.28290146144347617*sigmoid(-0.9795230057859095 + -0.9823842606294351*$(x[1]) + 0.3621989809961321*$(x[2])) + 0.5991916639069861*sigmoid(0.2857383741568795 + 0.17154081707188507*$(x[1]) + -0.2093148258397095*$(x[2])) + 0.7799932719867426*sigmoid(-0.5451083144913764 + -0.518481685587207*$(x[1]) + 0.024249687825614874*$(x[2])) + 0.858692422726886*sigmoid(-0.02703791931361943 + -0.7859299073841055*$(x[1]) + 0.7832525238663837*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# each_word_iteration.jl
# ======================
#
# An iterator over arrays of binary fetching 31 bits at a time.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/WAHVectors.jl/blob/master/LICENSE.md

immutable Every31Bits{T}
    data::T
end

# General methods.

Base.eltype(iter::Every31Bits)= UInt32

@inline function Base.start(iter::Every31Bits)
    return 1
end

# Specific methods.

# BitVectors.

@inline function Base.next(iter::Every31Bits{BitVector}, state::Int)
    chunk = 0x00000000
    lv = length(iter.data)
    niter = ifelse(lv >= (state + 30), 31, (lv - state) + 1)
    @inbounds for i in 1:niter
        bit_i = UInt32(iter.data[state])
        chunk |= (bit_i << (31 - i))
        state += 1
    end
    return chunk, state
end

@inline function Base.done(iter::Every31Bits{BitVector}, state::Int)
    return state > length(iter.data)
end

@inline function Base.length(iter::Every31Bits{BitVector})
    return UInt64(ceil(length(iter.data) / 31))
end

# WAHVectors

@inline function Base.start(iter::Every31Bits{WAHVector})
    return (1,nwords(iter.data.data[1]))
end

@inline function Base.next(iter::Every31Bits{WAHVector}, state)
    elem = iter.data.data[state[1]]
end

@inline function Base.length(iter::Every31Bits{WAHVector})
    return nwords(iter.data)
end


text:
struct SplineSparse{T,I,n_chunks,chunk_len} <: AbstractSparseMatrix{T,I}
    ncol::Int
    vals::Vector{T}
    cols::Vector{I}

    function SplineSparse{T,I,n_chunks,chunk_len}(
            col, vals, cols
        ) where {T,I,n_chunks,chunk_len}
        if length(cols)*chunk_len != length(vals)
            error("vals and cols not conformable")
        end
        new{T,I,n_chunks,chunk_len}(col, vals, cols)
    end
end

function SplineSparse(N::Int, L::Int, ncol::Int, vals::AbstractVector{T},
                      cols::AbstractVector{I}) where {T,I}
    SplineSparse{T,I,N,L}(ncol, vals, cols)
end

n_chunks(::Type{SplineSparse{T,I,N,L}}) where {T,I,N,L} = N
chunk_len(::Type{SplineSparse{T,I,N,L}}) where {T,I,N,L} = L
Base.eltype(::Type{SplineSparse{T,I,N,L}}) where {T,I,N,L} = T
ind_type(::Type{SplineSparse{T,I,N,L}}) where {T,I,N,L} = I
Base.copy(x::SplineSparse{T,I,N,L}) where {T,I,N,L} = SplineSparse{T,I,N,L}(x.ncol, x.vals, x.cols)

for f in [:n_chunks, :chunk_len, :(Base.eltype), :ind_type]
    @eval $(f)(s::SplineSparse) = $(f)(typeof(s))
end

@inline val_ix(s::SplineSparse{T,I,N,L}, row, chunk, n) where {T,I,N,L} =
    N*L*(row-1) + L*(chunk-1) + n

@inline col_ix(s::SplineSparse{T,I,N}, row, chunk) where {T,I,N} =
    N*(row-1) + chunk

function Base.Array(s::SplineSparse{T,I,N,L}) where {T,I,N,L}
    nrow = _nrows(s)
    out = zeros(T, nrow, s.ncol)

    for row in 1:nrow
        for chunk in 1:N
            first_col = s.cols[col_ix(s, row, chunk)]
            @simd for n in 1:L
                out[row, first_col+n-1] = s.vals[val_ix(s, row, chunk, n)]
            end
        end
    end

    out
end

function SparseArrays.findnz(s::SplineSparse{T,I,N,L}) where {T,I,N,L}
    nrow = _nrows(s)
    rows = repeat(collect(I, 1:nrow), inner=[N*L])
    cols = Array{I}(undef, length(s.vals))

    for row in 1:nrow
        for chunk in 1:N
            first_col = s.cols[col_ix(s, row, chunk)]
            for n in 1:L
                cols[val_ix(s, row, chunk, n)] = first_col+n-1
            end
        end
    end

    rows, cols, s.vals
end

function Base.convert(::Type{SparseMatrixCSC}, s::SplineSparse)
    I, J, V = findnz(s)
    sparse(I, J, V, _nrows(s), s.ncol)
end

SparseArrays.sparse(s::SplineSparse) = convert(SparseMatrixCSC, s)

function Base.getindex(s::SplineSparse{T,I,N,L}, row::Integer, cols::Integer) where {T,I,N,L}
    for chunk in 1:N
        first_col = s.cols[col_ix(s, row, chunk)]

        if cols < first_col || cols > (first_col + L-1)
            continue
        end

        n = cols-first_col+1
        return s.vals[val_ix(s, row, chunk, n)]

    end

    zero(T)
end

_nrows(s::SplineSparse) = Int(length(s.vals) / (n_chunks(s)*chunk_len(s)))
Base.size(s::SplineSparse) = (_nrows(s), s.ncol)
Base.size(s::SplineSparse, i::Integer) = i == 1 ? _nrows(s) :
                                         i == 2 ? s.ncol :
                                         1

@generated function row_kron(
        s1::SplineSparse{T1,I1,N1,len1}, s2::SplineSparse{T2,I2,N2,len2}
    ) where {T1,I1,N1,len1,T2,I2,N2,len2}

    # new number of chunks the length chunks times the number of chunks in
    # first matrix times number of chunks in second matrix
    N = len1*N1*N2

    # new chunk length is the chunk length from the second matrix
    len = len2

    # T and I are  easy...
    T = promote_type(T1, T2)
    I = promote_type(I1, I2)

    quote

    nrow = _nrows(s1)
    _nrows(s2) == nrow || error("s1 and s2 must have same number of rows")

    cols = Array{$I}(undef, $N*nrow)
    vals = Array{$T}(undef, $N*$len*nrow)

    ix = 0
    c_ix = 0
    @inbounds for row in 1:nrow
        for chunk1 in 1:N1
            first_col1 = s1.cols[col_ix(s1, row, chunk1)]

            for n1 in 1:len1
                v1 = s1.vals[val_ix(s1, row, chunk1, n1)]

                for chunk2 in 1:N2
                    first_col2 = s2.cols[col_ix(s2, row, chunk2)]

                    cols[c_ix+=1] = (first_col1+n1-2)*s2.ncol+first_col2

                    @simd for n2 in 1:len2
                        vals[ix+=1] = v1*s2.vals[val_ix(s2, row, chunk2, n2)]
                    end
                end
            end
        end
    end

    SplineSparse{$T,$I,$N,$len}(s1.ncol*s2.ncol, vals, cols)
    end
end

function LinearAlgebra.mul!(out::AbstractVector{Tout},
                       s::SplineSparse{T,I,N,L},
                       v::AbstractVector) where {T,I,N,L,Tout}
    @inbounds for row in eachindex(out)
        val = zero(Tout)
        for chunk in 1:N
            first_col = s.cols[col_ix(s, row, chunk)]

            @simd for n in 1:L
                val += s.vals[val_ix(s, row, chunk, n)] * v[first_col+(n-1)]
            end
        end

        out[row] = val
    end
    out
end

function *(s::SplineSparse{T}, v::AbstractVector{T2}) where {T,T2}
    size(s, 2) == size(v, 1) || throw(DimensionMismatch())

    out_T = promote_type(T, T2)
    out = Array{out_T}(undef, size(s, 1))
    mul!(out, s, v)
end

# TODO: define method A_mul_B!(ss::SplineSparse, csc::SparseMatrixCSC)
# TODO: define method A_mul_B!(ss::SplineSparse, ss2::SplineSparse)

function *(s::SplineSparse{T,I,N,L}, m::AbstractMatrix{T2}) where {T,I,N,L,T2}
    size(s, 2) == size(m, 1) || throw(DimensionMismatch())

    out_T = promote_type(T, T2)
    out = zeros(out_T, size(s, 1), size(m, 2))

    @inbounds for row in 1:size(s,1)
        for chunk in 1:N
            first_col = s.cols[col_ix(s, row, chunk)]

            for n in 1:L
                s_val = s.vals[val_ix(s, row, chunk, n)]
                s_col = first_col+(n-1)

                for m_col in 1:size(m, 2)
                    out[row, m_col] += s_val * m[s_col, m_col]
                end

            end
        end
    end

    out

end

# TODO: implement me for real to avoid conversion to CSC
function \(s::SplineSparse, x::Union{AbstractVector,AbstractMatrix})
    sparse(s) \ x
end

function tensor_prod(t::Type{T}, syms, inds, lens, add_index) where T<:AbstractArray
    if length(syms) == 0
        subs = [:($(Symbol("i_", j)) + $(inds[j])) for j in length(inds):-1:1]
        out = Expr(:ref, :c, subs...)
        T <: AbstractMatrix && push!(out.args, :col)
        return out
    else
        exprs = []
        for i in 1:lens[1]
            e = Expr(
                :call,
                :(*),
                Symbol(syms[1], "_", i),
                tensor_prod(
                    t,
                    syms[2:end],
                    cat(inds, [i-1], dims=1),
                    lens[2:end],
                    add_index
                )
            )
            push!(exprs, e)
        end
        return Expr(:call, :(+), exprs...)
    end
end

_ncol(x::AbstractArray) = size(x, 2)

shape_c_expr(::Type{T}) where {T<:AbstractVector} = :(reshape(_c, reverse(_ncol.(rk.B))))
shape_c_expr(::Type{T}) where {T<:AbstractMatrix} = :(reshape(_c, reverse(_ncol.(rk.B))..., size(_c, 2)))

const RKSS = RowKron{<:Tuple{Vararg{<:SplineSparse}}}

# if we have all `SplineSparse`s we can special case out = rk*c
@generated function LinearAlgebra.mul!(out::StridedVecOrMat,
                                  rk::RKSS,
                                  _c::StridedVecOrMat)
    N = length(rk.parameters[1].parameters)
    first_v = Expr(:(=), Symbol("v_", N), :(val_ix(rk.B[$N], row, 1, 1)))
    Ls = Int[chunk_len(i) for i in rk.parameters[1].parameters]

    unpack_Bs_args = []
    for i in 1:N
        L = Ls[i]
        iv_sym = Symbol("iv_", i)
        new_ex = [Expr(
            :(=),
            Symbol("v", i, "_", j),
            :(rk.B[$i].vals[$iv_sym + $(j-1)])
        ) for j in 1:L]
        push!(unpack_Bs_args, new_ex...)
    end
    unpack_Bs = Expr(:block, unpack_Bs_args...)

    syms = [Symbol("v", i) for i in 1:N]
    prod = tensor_prod(_c, syms, [], Ls, false)

    code = quote
        for B in rk.B
            n_chunks(B) == 1 && continue
            msg = "Only supported for combining univariate bases for now..."
            throw(ArgumentError(msg))
        end

        Base.@boundscheck begin
            if size(out, 1) != size(rk, 1)
                msg = "out should have $(size(rk, 1)) columns"
                throw(DimensionMismatch(msg))
            end
        end

        c = $(shape_c_expr(_c))

        @inbounds for col in 1:size(out, 2)
            for row in 1:size(out, 1)
                @nexprs $N j -> i_j = rk.B[j].cols[col_ix(rk.B[j], row, 1)]
                @nexprs $N j -> iv_j = val_ix(rk.B[j], row, 1, 1)
                $(unpack_Bs)
                out[row, col] = $(prod)
            end
        end
    end

    code

end


text:
######### Stan program example  ###########

using Mamba, Stan

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

bernoulli = "
data {
  int<lower=0> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
    y ~ bernoulli(theta);
}
"

bernoullidata = [
  Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
  Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 0, 1, 0, 0, 1]),
  Dict("N" => 10, "y" => [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]),
  Dict("N" => 10, "y" => [0, 0, 0, 1, 0, 0, 0, 1, 0, 1])
]

stanmodel = Stanmodel(Variational(), name="bernoulli", model=bernoulli);

sim1 = stan(stanmodel, bernoullidata, ProjDir, CmdStanDir=CMDSTAN_HOME);
monitor = ["theta"]
sim = sim1[1:size(sim1, 1), monitor, 1:size(sim1, 3)]
Mamba.describe(sim)
println()

## Plotting
p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
draw(p, ncol=1, nrow=4, filename="$(stanmodel.name)-variationalplot", fmt=:svg)
draw(p, ncol=1, nrow=4, filename="$(stanmodel.name)-variationalplot", fmt=:pdf)

# Below will only work on OSX, please adjust for your environment.
# JULIA_SVG_BROWSER is set from environment variable JULIA_SVG_BROWSER
@static is_apple() ? if isdefined(Main, :JULIA_SVG_BROWSER) && length(JULIA_SVG_BROWSER) > 0
        for i in 1:4
          isfile("$(stanmodel.name)-variationalplot-$(i).svg") &&
            run(`open -a $(JULIA_SVG_BROWSER) "$(stanmodel.name)-variationalplot-$(i).svg"`)
        end
      end : println()

end # cd