# Systems with uncertain coefficients
#
# Maple version of the worksheet 04_systems_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 04_systems_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 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.

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 intervals INTERVAL(a .. b) are handled in Maple: bounds, midpoint, width, inclusion
lft := I -> lhs(op(1, I)):
rgt := I -> rhs(op(1, I)):
midp := I -> (lft(I) + rgt(I)) / 2:
# the width of MPFI (diam in Julia): relative when the interval does not hold 0
diameter := I -> `if`(lft(I) <= 0 and 0 <= rgt(I), rgt(I) - lft(I),
                      (rgt(I) - lft(I)) / abs(midp(I))):
fbox := b -> evalf[10](bounds(b)):

# ## 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:

sys0 := [y^2 - x^3 + x, x^2 + y^2 - 3]:
sols := real_solutions(sys0, prec = 64):
printf("## 1. %a\n", sys0);
for b in sols do printf("  %a\n", fbox(b)) end do:

# ## 2. The family
#
# The coefficient 1 of x in the first equation and the radius 3 in the second are now known
# up to +-e. A coefficient is an interval INTERVAL(a .. b):

uncertain := (q, e) -> INTERVAL(q - e .. q + e):
family := e -> [y^2 - x^3 + uncertain(1, e) * x, x^2 + y^2 - uncertain(3, e)]:
F := family(1/100):
printf("\n## 2. The family\n%a\n", F);

# ## 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. enclose returns
# [mu, box], the box being nothing when the test fails.

enclose := proc(F, sol, {flags := 0})
    local x0, k, mu, rho, res;
    x0 := map(midp, sol);
    k := kantorovich_test(F, evalf[40](x0), 128, 64, ':-flags' = flags);
    mu, rho := rgt(k[2]), rgt(k[1]);
    # mu is undefined (NaN) when the bound cannot be computed: the test fails as well
    if not (type(mu, numeric) and mu < 1) then return [evalf[6](mu), nothing] end if;
    res := interval_newton(F, [seq(INTERVAL(c - rho .. c + rho), c in x0)], 128, 128, 64, 30,
                           ':-verbose' = 0, ':-flags' = flags);
    [evalf[6](mu), res[2]]
end proc:

enclosures := [seq(enclose(F, s), s in sols)]:
printf("\n## 3. Kantorovich, then Newton\n");
for E in enclosures do printf("  mu = %a, box %a\n", E[1], fbox(E[2])) end do:

# 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 +- e and c = 3 +- e, are solved
# exactly; each of their real solutions lies in one of the boxes:

e := 1/100:
inside := (I, J) -> evalb(lft(J) <= lft(I) and rgt(I) <= rgt(J)):
inbox := (s, B) -> andmap(j -> inside(s[j], B[j]), [seq(1 .. nops(B))]):
printf("\n## 4. The corners of the family\n");
for a in [1 - e, 1 + e] do
    for c in [3 - e, 3 + e] do
        printf("  a = %a, c = %a: %a\n", a, c,
               [seq(ormap(E -> inbox(s, E[2]), enclosures),
                    s in real_solutions([y^2 - x^3 + a * x, x^2 + y^2 - c], prec = 64))]);
    end do:
end do:

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

printf("mu as e grows:\n");
for e in [1/10^6, 1/100, 1/10, 1/5] do
    printf("  e = %-9a %a\n", e, [seq(enclose(family(e), s)[1], s in sols)]);
end do:

# ## 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|].

katsura := proc(n)
    local U, eqs, m, l;
    U := i -> `if`(abs(i) <= n, u || (abs(i)), 0);
    eqs := [seq(add(U(l) * U(m - l), l = -n .. n) - U(m), m = 0 .. n - 1)];
    [op(eqs), add(U(l), l = -n .. n) - 1]
end proc:
k6 := katsura(6):
vars6 := [seq(u || i, i = 0 .. 6)]:
k6_sols := real_solutions(k6, prec = 64):
printf("\n## 5. Katsura-6: %d real solutions\n", nops(k6_sols));

# A polynomial with rational coefficients turned into one with interval coefficients:

relative := proc(p, e)
    local cs, ts;
    cs := [coeffs(expand(p), vars6, 'ts')];
    add(uncertain(cs[i], e * abs(cs[i])) * [ts][i], i = 1 .. nops(cs))
end proc:
printf("%a\n", relative(k6[-1], 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) -> ormap(j -> rgt(A[j]) < lft(B[j]) or rgt(B[j]) < lft(A[j]),
                            [seq(1 .. nops(A))]):

certify := proc(e, {flags := 0})
    local F, boxes, overlaps, i, j, w;
    F := map(relative, k6, e);
    boxes := [seq(`if`(E[2] = nothing, NULL, E[2]),
                  E in [seq(enclose(F, s, ':-flags' = flags), s in k6_sols)])];
    overlaps := add(add(`if`(disjoint(boxes[i], boxes[j]), 0, 1), j = i + 1 .. nops(boxes)),
                    i = 1 .. nops(boxes));
    w := `if`(boxes = [], 0, max(seq(seq(diameter(c), c in B), B in boxes)));
    [e, nops(boxes), overlaps, evalf[3](w)]
end proc:

printf("   e        certified  overlaps  max width\n");
for e in [1/10^12, 1/10^6, 1/10^5, 1/10^4, 1/10^3] do printf("  %a\n", certify(e)) end do:

# 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:
#
# - KANTO_HESSIAN: a bound independent of a diagonal scaling of the equations;
# - 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;
# - KANTO_SCHULZ, with KANTO_PRECOND: when C is poor, it is improved by Schulz iterations
#   C <- 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 (infinity 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:

F := family(1/100):
scaled := s -> [F[1], expand(s * F[2])]:
printf("\n## 6. Preconditioning\nscaling, classical test:\n");
for s in [1, 1000, 10^6] do printf("  s = %-8a mu = %a\n", s, enclose(scaled(s), sols[1])[1]) end do:

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

printf("scaling, KANTO_PRECOND:\n");
for s in [1, 1000, 10^6] do
    printf("  s = %-8a mu = %a\n", s, enclose(scaled(s), sols[1], flags = KANTO_PRECOND)[1]);
end do:

# 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 := 1/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 := [1 + 2^(-30), 1 - 2^(-31)]:
opts := ["classical" = 0, "HESSIAN" = KANTO_HESSIAN, "PRECOND" = KANTO_PRECOND,
            "PRECOND + SCHULZ" = KANTO_PRECOND + KANTO_SCHULZ]:
from_point := proc(sys, x0, prec, flags)
    local res, k;
    res := interval_newton(sys, evalf[30](x0), prec, prec, 100, 64, 30, ':-verbose' = 0,
                           ':-flags' = flags);
    k := res[3];
    [res[1], evalf[3](rgt(k[2])), `if`(nops(k) = 5, evalf[3](rgt(k[5])), nothing)]
end proc:
printf("from a point, 64 bits: [status, mu, delta]\n");
for o in opts do printf("  %-17s %a\n", lhs(o), from_point(near, x64, 64, rhs(o))) end do:

# 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 ~ 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:

printf("status of Newton by precision:\n");
for p in [64, 80, 128, 200] do
    printf("  %3d bits: %a\n", p, from_point(near, x64, p, KANTO_PRECOND + KANTO_SCHULZ)[1]);
end do:

# 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:

printf("Katsura-6, certified: [e, classical, KANTO_PRECOND]\n");
for e in [1/10^5, 1/10^4, 1/10^3] do
    printf("  %a\n", [e, certify(e)[2], certify(e, flags = KANTO_PRECOND)[2]]);
end do:

# At e = 10^(-3), 14 solutions are certified instead of one, in boxes that are still pairwise
# disjoint ([e, certified, overlaps, max width]):

printf("  %a\n", certify(1/10^3, flags = 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:
printf("\n## 7. y = %a and the unit circle: %a\n", c0,
       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 <= 1 have
# real solutions: at the corner c = c0 - e,

printf("y = %a: %a\n", c0 - 1/100,
       map(fbox, 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.