# Polynomials with uncertain coefficients
#
# Maple version of the worksheet 02_univariate_uncertain_coefficients.jl: the functions of
# PACE.jl are called from Maple through UDX (see "Using PACE.jl in Maple with UDX.jl"). Run
# in batch with
#
#     maple -q 02_univariate_uncertain_coefficients.mpl
#
# with UDX_MPL_LIB set, and UDX_SRC, PACE_SRC, JULIA, JULIA_PROJECT, UDX_LIB set or replaced
# below by their paths.
#
# 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 +- e.

interface(prettyprint = 0, warnlevel = 0):
read cat(getenv("UDX_SRC"), "/Maple/udx.mpl"):
read cat(getenv("PACE_SRC"), "/maple/PACE.mpl"):
UDX:-udx_set_julia(getenv("JULIA"), env = ["JULIA_PROJECT" = getenv("JULIA_PROJECT"),
                                           "UDX_LIB" = getenv("UDX_LIB")]):
with(PACE):

# the real roots of a polynomial with rational coefficients, refined to prec bits
real_roots := proc(p, prec)
    local iso;
    iso := isolation(p);
    refine(iso, prec);
    isolate(iso)
end proc:
# the upper bound of an interval INTERVAL(a .. b), exact
up := I -> rhs(op(1, I)):

# ## 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").

f0 := x^3 - 2*x - 5:
r0 := real_roots(f0, 64):
printf("## 1. The exact instance %a\n%d real root in %a\n", f0, nops(r0),
       evalf[20]([lo(r0[1]), hi(r0[1])]));

# One real root. Its approximation x0 is a floating-point number (40 digits, about 128 bits),
# and x0r the same number as an exact rational:

x0 := evalf[40](lo(r0[1])):
x0r := convert(x0, rational, exact):

# ## 2. The family
#
# A polynomial with interval coefficients is written with INTERVAL(a .. b), a and b exact
# rationals or floats; Julia receives it as a polynomial with MPFI coefficients (BigInterval),
# the bounds rounded outward. With e = 1/1000:

family := e -> y^3 - 2*y + INTERVAL(-5 - e .. -5 + e):
F := family(1/1000):
printf("\n## 2. The family\nF = %a\n", F);

# ## 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):
rho := up(k[1]):
printf("\n## 3. The Kantorovich test\nrho <= %.3e, mu <= %.3e\n", evalf(rho), evalf(up(k[2])));

# Here mu ~ 2e-4 < 1: each polynomial of the family has exactly one root in
# [x0 - rho, x0 + rho].

X0 := INTERVAL(x0r - rho .. x0r + 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. It returns the status, the box
# and the Kantorovich numbers.

res := interval_newton([F], [X0], 128, 128, 64, 30, verbose = 0):
sol := res[2][1]:
printf("\n## 4. Interval Newton\nstatus %a, enclosure %a\n", res[1], evalf[12](bounds(sol)));

# 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]:
printf("exact set of the roots: %a\n", evalf[12]([lo(a), hi(b)]));

# and the enclosure holds them, with little overestimation:

printf("held by the enclosure: %a\n", evalb(op([1, 1], sol) <= lo(a) and hi(b) <= up(sol)));

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

enclose := proc(F, x0)
    local k, mu, rho, x0r, res;
    k := kantorovich_test([F], [x0], 128, 64);
    mu, rho := up(k[2]), up(k[1]);
    if mu >= 1 then return [evalf[6](mu), "not certified"] end if;
    x0r := convert(x0, rational, exact);
    res := interval_newton([F], [INTERVAL(x0r - rho .. x0r + rho)], 128, 128, 64, 30,
                           verbose = 0);
    [evalf[6](mu), evalf[6](bounds(res[2][1]))]
end proc:

printf("\n## 5. Larger uncertainties\n   e        mu         enclosure\n");
for e in [1/10^6, 1/10^3, 1/10, 1/2, 1] do
    printf("  %-8a %a\n", e, enclose(family(e), x0));
end do:

# 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 +- 1/100. With an uncertainty e on
# c, the family contains (y - 1)^2 as soon as e >= 10^(-4): two roots merge, and beyond, some
# members have no real root.

c := 1 - 1/10^4:
close_family := e -> y^2 - 2*y + INTERVAL(c - e .. c + e):
starts := [evalf[40](99/100), evalf[40](101/100)]:
printf("\n## 6. Two roots 1 +- 1/100\n");
for e in [1/10^6, 5/10^5, 1/10^4, 2/10^4] do
    printf("  e = %a\n", e);
    for s in starts do printf("    from %.2f: %a\n", s, enclose(close_family(e), s)) end do:
end do:

# Up to e = 5e-5 both roots are enclosed, in disjoint intervals: every polynomial of the
# family has two simple real roots, one in each. At e = 1e-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.