# # Structural stability of a 2-D system: from a polydisc to real solutions
#
# A discrete 2-D system with transfer function $N/D$ is structurally stable when
# $D(z_1,z_2) \neq 0$ on the closed unit bidisc $\overline{\mathbb{U}}^2$. This is a
# question about infinitely many points. We turn it, exactly, into the question
# *does a system of two polynomial equations have a real solution?*, and answer it
# with certified tools.
#
# Running example:
# $$D(z_1,z_2) = (2z_1^2+10z_1+12) + (z_1^2+5z_1+6)\,z_2 = (z_1+2)(z_1+3)(z_2+2),$$
# stable, which lets us check each step. The algorithms do not factor $D$.

#nb # Notebook only: Nemo and AbstractAlgebra objects are shown as plain text
#nb # (their HTML and LaTeX output throws on purpose, which IJulia 1.34 no longer catches).
#nb filter!(m -> !(m == MIME("text/latex") || (m isa AbstractVector && MIME("text/html") in m)), IJulia.ijulia_mime_types)

using AbstractAlgebra, PACE, PaceControl

S, (z1, z2) = polynomial_ring(QQ, ["z1", "z2"])
D = (2z1^2 + 10z1 + 12) + (z1^2 + 5z1 + 6) * z2

# ## Step 1: two 1-D tests (DeCarlo, Murray, Saeks)
#
# $D$ is stable if and only if $D(z_1,1)$ and $D(1,z_2)$ have no root in the closed
# unit disc and $D$ has no zero on the torus $\mathbb{T}^2$.
#
# A univariate $p$ of degree $d$ has all its roots outside the closed disc if and only
# if $(1+s)^d\,p\bigl(\tfrac{1-s}{1+s}\bigr)$ is Hurwitz: the map $s\mapsto\frac{1-s}{1+s}$
# sends the open left half-plane onto the exterior of the disc. The exact Hurwitz test
# of [the first worksheet](01_one_question_four_paradigms.md) then decides.

U, s = polynomial_ring(QQ, "s")

function no_root_in_disc(p)
    d = degree(p)
    h = sum(coeff(p, j) * (1 - s)^j * (1 + s)^(d - j) for j in 0:d)
    hurwitz_stable(Rational{BigInt}[coeff(h, k) for k in degree(h):-1:0])
end

p1 = AbstractAlgebra.evaluate(D, [gen(U), one(U)])      ## D(z1, 1)
p2 = AbstractAlgebra.evaluate(D, [one(U), gen(U)])      ## D(1, z2)
(p1, no_root_in_disc(p1)), (p2, no_root_in_disc(p2))

# ## Step 2: from the torus to real points (Möbius)
#
# $z = \frac{x-i}{x+i}$ maps the real line onto $\mathbb{T}\setminus\{1\}$. Substituting
# it for $z_1$ and $z_2$ and clearing denominators, the numerator of $D$ is
# $\mathcal{R}(x_1,x_2) + i\,\mathcal{I}(x_1,x_2)$ with rational $\mathcal{R}$,
# $\mathcal{I}$, and $D$ has no zero on $(\mathbb{T}\setminus\{1\})^2$ if and only if
# $\mathcal{R}=\mathcal{I}=0$ has no real solution. (The points with a coordinate
# equal to $1$ are covered by step 1.)
#
# Complex polynomials are handled as pairs (real part, imaginary part):

R, (x1, x2) = polynomial_ring(QQ, ["x1", "x2"])
cmul((a, b), (c, d)) = (a * c - b * d, a * d + b * c)
cpow(p, k) = k == 0 ? (one(R), zero(R)) : cmul(p, cpow(p, k - 1))

function moebius(D)
    d1, d2 = degree(D, 1), degree(D, 2)
    re, im = zero(R), zero(R)
    for (c, (j, k)) in zip(coefficients(D), exponent_vectors(D))
        t = cmul(cmul(cpow((x1, -one(R)), j), cpow((x1, one(R)), d1 - j)),
                 cmul(cpow((x2, -one(R)), k), cpow((x2, one(R)), d2 - k)))
        re += c * t[1]
        im += c * t[2]
    end
    re, im
end

Re, Im = moebius(D)

# ## Step 3: real solutions of a bivariate system
#
# $\mathcal{R}$ and $\mathcal{I}$ have no common factor, so the system has finitely
# many complex solutions. For two equations in two unknowns, LACE has a dedicated
# solver, which uses subresultants over the integers, computed modulo many primes,
# and no Gröbner basis. It returns a *rational univariate representation* (RUR): a
# polynomial `rext` in a new variable $T$ (here $T = x_2 - x_1$), and two rational
# functions of $T$ that give $x_1$ and $x_2$ at each of its roots.

gcd(Re, Im)

#-

comps = biv_parameterization([Re, Im]; verbose = false)
comps[1].rext

# This is $324T^4 + 369T^2 + 100 = (36T^2+25)(9T^2+4)$: no real root. The certified
# isolation of the real solutions finds none:

# The bivariate RUR of LACE (no Gröbner basis), then the certified real solutions.
# A constant resultant means no solution, reported as a `PACEError` of code -1.
function biv_real_solutions(a, b)
    comps = try
        biv_parameterization([a, b]; verbose = false)
    catch e
        e isa PACEError && e.code == -1 || rethrow()
        []
    end
    vcat([real_solutions(c.rext, c.rden, [c.num_Y, c.num_X]; prec = 30) for c in comps]...)
end

real_torus_zeros(a, b) = length(biv_real_solutions(a, b))
real_torus_zeros(Re, Im)

# $D$ is structurally stable. This is a proof: every step is exact or certified.

# ## Near the boundary: a tolerance cannot decide
#
# Take $D_\varepsilon = 2 + z_1 + z_2 + \varepsilon$. For $\varepsilon > 0$ it is stable
# ($|z_1+z_2| \le 2$ on the bidisc); for $\varepsilon < 0$ it vanishes at
# $z_1 = z_2 = -1-\varepsilon/2$, inside the bidisc. With $\varepsilon = \pm 10^{-12}$,
# a floating-point test samples $|D_\varepsilon|$ on the torus:

function min_on_torus(ε, N)
    m = Inf
    for j in 0:N-1, k in 0:N-1
        w1, w2 = cis(2π * j / N), cis(2π * k / N)
        m = min(m, abs(2 + w1 + w2 + ε))
    end
    m
end

[(ε, min_on_torus(ε, 1000)) for ε in (1e-12, -1e-12)]

# Both minima are positive and equal (the grid contains $z_1=z_2=-1$): a sampling (or any test with a
# tolerance above $10^{-12}$) cannot tell the stable case from the unstable one. The
# exact test does:

for ε in (1 // 10^12, -1 // 10^12)
    Dε = 2 + z1 + z2 + ε
    a, b = moebius(Dε)
    q1 = AbstractAlgebra.evaluate(Dε, [gen(U), one(U)])
    q2 = AbstractAlgebra.evaluate(Dε, [one(U), gen(U)])
    println("ε = ", ε, ": disc tests ", (no_root_in_disc(q1), no_root_in_disc(q2)),
            ", real zeros on the torus: ", real_torus_zeros(a, b))
end

# For $\varepsilon = -10^{-12}$, the two real solutions are the two zeros of $D_\varepsilon$
# on the torus, at distance about $10^{-6}$ from $z_1 = z_2 = -1$: the system is not
# stable.

# ## Summary
#
# | Step | Question | Tool |
# |---|---|---|
# | 1 | two univariate disc tests | exact Liénard-Chipart test after a Möbius map (PaceControl) |
# | 2 | zeros on the torus ⟺ real zeros of $\{\mathcal{R},\mathcal{I}\}$ | exact substitution |
# | 3 | does $\{\mathcal{R}=\mathcal{I}=0\}$ have a real solution? | bivariate RUR by subresultants + certified isolation (LACE) |
#
# Sampling the torus gives no proof, and near the boundary no tolerance separates the
# two answers.