# # Polynomials with uncertain coefficients
#
# The coefficients of a polynomial often come from measurements: they are known up to an
# error, that is, they lie in intervals. Such a polynomial is a *family* of polynomials,
# and its "roots" are the set of the roots of all its members. We show how to enclose
# this set, with a proof, by combining an exact computation and interval arithmetic:
#
# 1. solve exactly an *instance* of the family, with small rational coefficients;
# 2. prove with the Kantorovich test that each polynomial of the family has one root, and
#    only one, in a ball around the exact root;
# 3. shrink this ball by interval Newton to an enclosure of the roots of the family.
#
# The example is $x^3 - 2x - 5$, the polynomial on which Newton presented his method,
# with an uncertain constant term $-5 \pm e$.

#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. An exact instance
#
# The instance with $e = 0$ has rational coefficients: its real roots are isolated and
# refined exactly (see the worksheet *Real roots of a univariate polynomial*).

R, x = polynomial_ring(QQ, "x")
coeffs(p) = Rational{BigInt}.(collect(coefficients(p)))

function real_roots(p, prec)
    iso = isolation(coeffs(p))
    PACE.refine!(iso, prec)
    isolate(iso)
end

f0 = x^3 - 2x - 5
real_roots(f0, 64)

# One real root. Its approximation `x0` is a floating-point number, at 64 bits:

x0 = BigFloat(PACE.lo(real_roots(f0, 64)[1]); precision = 128)

# ## 2. The family
#
# A polynomial with interval coefficients is written with DynamicPolynomials, the
# coefficients being intervals of MPFI (`BigInterval`). With $e = 1/1000$:

@polyvar y
family(e) = y^3 - 2y + (BigInterval(-5; precision = 128) +
                        BigInterval(e; precision = 128) * BigInterval(-1, 1))
F = family(1//1000)

# ## 3. The Kantorovich test
#
# `kantorovich_test` computes, in interval arithmetic, two numbers at `x0`: a radius $\rho$
# and a number $\mu$. When $\mu < 1$, Newton's method converges from `x0` to a root, unique
# in the ball of radius $\rho$ around `x0`. The computation being done on the
# intervals of the coefficients, the bounds hold for *every* polynomial of the family.

k = kantorovich_test([F], [x0], 128, 64)

# Here $\mu \approx 2 \cdot 10^{-4} < 1$: each polynomial of the family has exactly one root
# in $[x_0 - \rho, x_0 + \rho]$.

rho = right(k[1])
X0 = BigInterval(x0 - rho, x0 + rho)

# ## 4. Interval Newton
#
# `interval_newton` starts from this box. Each step intersects the Newton image with the
# current box, so that the roots in the box are never lost:

status, sol, _ = interval_newton([F], [X0], 128, 128, 64, 30; verbose = Cint(0))
status, sol[1]

# The status 2 says that the box no longer shrinks: with uncertain coefficients the
# roots fill an interval of positive width, and the box cannot be smaller than it. As the
# starting box was proved to hold the roots, the result is an enclosure of them.
#
# Here the exact set of roots is known: the root increases when the constant term
# decreases, so the roots of the family go from the root of $x^3 - 2x - 5 + e$ to that
# of $x^3 - 2x - 5 - e$. Both are isolated exactly:

e = 1//1000
a = real_roots(f0 + e, 80)[1]
b = real_roots(f0 - e, 80)[1]
Float64(PACE.lo(a)), Float64(PACE.hi(b))

# and the enclosure holds them, with little overestimation:

left(sol[1]) <= PACE.lo(a) && PACE.hi(b) <= right(sol[1])

# ## 5. Larger uncertainties
#
# The same three steps for several widths:

function enclose(F, x0)
    k = kantorovich_test([F], [x0], 128, 64)
    mu, rho = right(k[2]), right(k[1])
    mu < 1 || return (Float64(mu), "not certified")
    _, sol, _ = interval_newton([F], [BigInterval(x0 - rho, x0 + rho)], 128, 128, 64, 30;
                                verbose = Cint(0))
    return (Float64(mu), sol[1])
end

[(e, enclose(family(e), x0)...) for e in (1//10^6, 1//10^3, 1//10, 1//2, 1)]

# The certificate holds up to $e = 1$: all the polynomials of the family still have a
# single real root, and the enclosure $[1.994, 2.196]$ holds the exact set of these roots,
# $[2, 2.18]$. The overestimation grows with the width of the coefficients, as always with
# interval arithmetic.

# ## 6. When the uncertainty is too large
#
# $y^2 - 2y + c$ with $c = 1 - 10^{-4}$ has two real roots, $1 \pm 1/100$. With an uncertainty
# $e$ on $c$, the family contains $(y-1)^2$ as soon as $e \ge 10^{-4}$: two roots merge, and
# beyond, some members have no real root.

c = 1 - 1//10^4
close_family(e) = y^2 - 2y + (BigInterval(c; precision = 128) +
                              BigInterval(e; precision = 128) * BigInterval(-1, 1))
starts = [BigFloat(99//100; precision = 128), BigFloat(101//100; precision = 128)]
[(e, [enclose(close_family(e), s) for s in starts]...) for e in (1//10^6, 5//10^5, 1//10^4, 2//10^4)]

# Up to $e = 5 \cdot 10^{-5}$ both roots are enclosed, in disjoint intervals: every polynomial
# of the family has two simple real roots, one in each. At $e = 10^{-4}$ the test fails
# ($\mu = 1$), as it must: no ball can hold one root of each member, since one of them has
# a double root. The failure of the certificate is the information: the uncertainty on
# the data is too large to tell the two roots apart.