# # Systems with uncertain coefficients
#
# The worksheet *Polynomials with uncertain coefficients* enclosed the roots of a family of
# univariate polynomials. We do the same for systems: the coefficients lie in intervals,
# and we look for boxes that hold, for *every* system of the family, exactly one
# solution. The three steps are the same:
#
# 1. solve exactly an instance of the family, with rational coefficients (worksheet
#    *Real solutions of a zero-dimensional system*);
# 2. at each real solution, prove with the Kantorovich test that every system of the
#    family has one solution, and only one, in a ball around it;
# 3. shrink these balls by interval Newton.
#
# A last section shows what this approach cannot see.

#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 PACE, AbstractAlgebra
using PACE.DynamicPolynomials: @polyvar

# ## 1. The exact instance
#
# The system of the previous worksheet, the cubic $y^2 = x^3 - x$ and the circle
# $x^2 + y^2 = 3$:

R, (x, y) = polynomial_ring(QQ, ["x", "y"])
sys0 = [y^2 - x^3 + x, x^2 + y^2 - 3]
sols = real_solutions(sys0; prec = 64)

# ## 2. The family
#
# The coefficient $1$ of $x$ in the first equation and the radius $3$ in the second are now
# known up to $\pm e$. A system with interval coefficients is written with
# DynamicPolynomials:

@polyvar X Y
uncertain(q, e) = BigInterval(q; precision = 128) + BigInterval(e; precision = 128) * BigInterval(-1, 1)
family(e) = [Y^2 - X^3 + uncertain(1, e) * X, X^2 + Y^2 - uncertain(3, e)]
F = family(1//100)

# ## 3. Kantorovich, then Newton
#
# At the center `x0` of each exact solution, the Kantorovich test gives a radius $\rho$ and
# a number $\mu$; when $\mu < 1$, every system of the family has a unique solution in the
# ball of radius $\rho$ around `x0`. Interval Newton started from this ball shrinks it:

function enclose(F, sol; flags = 0)
    x0 = [BigFloat(mid(c); precision = 128) for c in sol]
    k = kantorovich_test(F, x0, 128, 64; flags = flags)
    mu, rho = right(k[2]), right(k[1])
    mu < 1 || return (mu = Float64(mu), box = nothing)
    _, box, _ = interval_newton(F, [BigInterval(c - rho, c + rho) for c in x0], 128, 128, 64, 30;
                                verbose = Cint(0), flags = flags)
    return (mu = Float64(mu), box = box)
end

enclosures = [enclose(F, s) for s in sols]

# Both solutions are certified: every system of the family has a solution in each box.

# ## 4. A check
#
# The four systems at the corners of the family, $a = 1 \pm e$ and $c = 3 \pm e$, are
# solved exactly; each of their real solutions lies in one of the boxes:

e = 1//100
inbox(s, B) = all(is_inside(s[j], B[j]) for j in eachindex(B))
[(a, c, [any(inbox(s, E.box) for E in enclosures)
         for s in real_solutions([y^2 - x^3 + a * x, x^2 + y^2 - c]; prec = 64)])
 for a in (1 - e, 1 + e) for c in (3 - e, 3 + e)]

# As $e$ grows, the boxes widen, and the certificate fails when $\mu$ reaches 1:

[(e, [E.mu for E in (enclose(family(e), s) for s in sols)]) for e in (1//10^6, 1//100, 1//10, 1//5)]

# ## 5. More variables: Katsura-6
#
# The Katsura-6 system has 32 real solutions (worksheet *Real solutions of a
# zero-dimensional system*). We give *every* coefficient a relative uncertainty $e$: the
# coefficient $c$ becomes $[c - e|c|, c + e|c|]$.

function katsura(n)
    K, u = polynomial_ring(QQ, ["u$i" for i in 0:n])
    U(i) = abs(i) <= n ? u[abs(i)+1] : zero(K)
    eqs = [sum(U(l) * U(m - l) for l in -n:n) - U(m) for m in 0:n-1]
    push!(eqs, sum(U(l) for l in -n:n) - 1)
    return eqs
end
k6 = katsura(6)
k6_sols = real_solutions(k6; prec = 64)
length(k6_sols)

# The conversion of a polynomial of AbstractAlgebra into a polynomial with interval
# coefficients:

@polyvar v[1:7]
relative(p, e) = sum(uncertain(c, e * abs(c)) * prod(v[i]^m[i] for i in eachindex(m))
                     for (c, m) in zip(coefficients(p), exponent_vectors(p)))
relative(k6[end], 1//10^6)

# Two boxes that overlap could hold the same solution of a system of the family: to count
# distinct solutions, the boxes must also be pairwise disjoint.

disjoint(A, B) = any(right(A[j]) < left(B[j]) || right(B[j]) < left(A[j]) for j in eachindex(A))

function certify(e; flags = 0)
    F = [relative(p, e) for p in k6]
    boxes = [E.box for E in (enclose(F, s; flags = flags) for s in k6_sols) if E.box !== nothing]
    overlaps = count(!disjoint(boxes[i], boxes[j]) for i in eachindex(boxes) for j in i+1:length(boxes))
    width = isempty(boxes) ? 0.0 : maximum(Float64(diam(c)) for B in boxes for c in B)
    return (e = e, certified = length(boxes), overlaps = overlaps, max_width = width)
end

[certify(e) for e in (1//10^12, 1//10^6, 1//10^5, 1//10^4, 1//10^3)]

# With $e = 10^{-12}$ the 32 solutions are certified, but the widest box is already
# about $10^{-8}$ wide: some solutions amplify the uncertainty of the data by $10^4$. The
# number of certified solutions then drops: at $e = 10^{-3}$ a single one remains. A
# solution that is not certified may still exist for every system of the family; the
# test does not decide it.

# ## 6. Preconditioning options
#
# `kantorovich_test` and `interval_newton` take a keyword `flags`, a sum of:
#
# - `PACE.KANTO_HESSIAN`: a bound independent of a diagonal scaling of the equations;
# - `PACE.KANTO_PRECOND`: the test is run on the system multiplied by $C$, an approximate
#   inverse of the Jacobian at `x0` (the inverse of its midpoint, in floating point at the
#   working precision); the bound no longer depends on how the equations are written;
# - `PACE.KANTO_SCHULZ`, with `PACE.KANTO_PRECOND`: when $C$ is poor, it is improved by
#   Schulz iterations $C \leftarrow C(2I - J C)$.
#
# `flags = 0` is the classical test, which returns $[2\eta, \mu]$. Otherwise the test
# returns $[2\eta, \mu, \mu_1, \mu_2, \delta]$: $\mu$ the best of the bounds computed,
# $\mu_1$ the one of `KANTO_HESSIAN`, $\mu_2$ the one of `KANTO_PRECOND` (`Inf` when not
# computed) and $\delta = \|I - CJ\|$, the defect of the preconditioner ($-1$ without it).
# In `interval_newton`, `flags` sets the Kantorovich test run at the starting point, whose
# constants are its third result; the Newton steps themselves do not depend on it.

# ### Scaling the equations
#
# Multiplying an equation by a constant changes neither the solutions nor the Newton
# iterates, but it changes the classical Kantorovich bound. With the second equation of
# section 2 multiplied by $s$:

scaled(s) = [F[1], s * F[2]]
[(s, enclose(scaled(s), sols[1]).mu) for s in (1, 1000, 10^6)]

# At $s = 1000$ the solution is no longer certified, although nothing has changed. With
# the preconditioned test:

[(s, enclose(scaled(s), sols[1]; flags = PACE.KANTO_PRECOND).mu) for s in (1, 1000, 10^6)]

# The bound no longer depends on $s$, and is smaller than the classical one even for
# $s = 1$.

# ### An ill-conditioned system, from a point
#
# Two curves that are almost the same near $(1, 1)$: their Jacobian there has a condition
# number of about $10^{18}$.

t = QQ(1, big(10)^18)
near = [x^2 + x + y - 3, x^2 + x + (1 + t) * y - 3 - t]

# `interval_newton` is started from a point close to the solution, at a working precision
# of 64 bits, with each option:

x64 = [BigFloat(1; precision = 64) + BigFloat(2)^-30, BigFloat(1; precision = 64) - BigFloat(2)^-31]
options = ["classical" => 0, "HESSIAN" => PACE.KANTO_HESSIAN, "PRECOND" => PACE.KANTO_PRECOND,
           "PRECOND + SCHULZ" => PACE.KANTO_PRECOND | PACE.KANTO_SCHULZ]
function from_point(sys, x0, prec, flags)
    status, _, k = interval_newton(sys, x0, prec, prec, 100, 64, 30; verbose = Cint(0), flags = flags)
    return (status = status, mu = Float64(right(k[2])),
            delta = length(k) == 5 ? Float64(right(k[5])) : nothing)
end
[(name, from_point(near, x64, 64, f)...) for (name, f) in options]

# The classical bound is about $10^{28}$; with the preconditioner it is below 1: the
# solution is certified, unique in the ball of radius $2\eta$ around the starting point.
# Here $\delta \approx 0.4 > 0.1$, so the Schulz iterations are run, and lower $\mu$ a little.
#
# The Newton iteration itself fails: status 0, a singular interval Jacobian. Its steps
# solve the interval linear system by Gaussian elimination, without preconditioning, and at
# 64 bits the interval Jacobian on the box contains singular matrices. It needs more
# precision to conclude on its own:

[(prec, from_point(near, [BigFloat(c; precision = prec) for c in x64], prec,
                   PACE.KANTO_PRECOND | PACE.KANTO_SCHULZ).status) for prec in (64, 80, 128, 200)]

# At 200 bits Newton encloses the solution (status 1). At 64 bits the certificate is the
# ball of the preconditioned Kantorovich test.

# ### Katsura-6
#
# With uncertain coefficients, the starting balls of section 5 come from the Kantorovich
# test, and the preconditioned test certifies many more solutions when the uncertainty
# grows:

[(e, certify(e).certified, certify(e; flags = PACE.KANTO_PRECOND).certified)
 for e in (1//10^5, 1//10^4, 1//10^3)]

# At $e = 10^{-3}$, 14 solutions are certified instead of one, in boxes that are still
# pairwise disjoint:

certify(1//10^3; flags = PACE.KANTO_PRECOND)

# Some of them are wide (up to about 2): for these solutions an uncertainty of $10^{-3}$
# on the data leaves a large uncertainty on the result, but each box still holds exactly
# one solution of every system of the family.

# ## 7. What this approach cannot see
#
# The method follows the real solutions of the exact instance. It proves nothing about
# the solutions the instance does not have. The line $y = c$ and the unit circle, with
# $c = 1 + 1/1000$, do not meet:

c0 = 1 + 1//1000
real_solutions([x^2 + y^2 - 1, y - c0]; prec = 64)

# Nothing to enclose. Yet with an uncertainty $e = 1/100$ on $c$, the systems with
# $c \le 1$ have real solutions: at the corner $c = c_0 - e$,

real_solutions([x^2 + y^2 - 1, y - (c0 - 1//100)]; prec = 64)

# Two complex solutions of the instance became real. To describe *all* the systems of a
# family, the uncertain coefficients must be treated as parameters, and the parameter
# space split into regions where the number of real solutions is constant: this is the
# subject of the worksheet on systems depending on parameters.