# Systems depending on parameters
#
# Maple version of the worksheet 06_parametric_systems.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 06_parametric_systems.mpl
#
# with UDX_MPL_LIB set, and UDX_SRC, PACE_SRC, JULIA, JULIA_PROJECT, UDX_LIB set or replaced
# below by their paths. The drawings are SVG files, written in the current directory.
#
# The coefficients of a system often depend on parameters. For each value of the parameters
# the system has finitely many complex solutions, but the number of *real* ones changes with
# the parameters. We describe this for *all* the values at once:
#
# 1. the *discriminant variety* is a hypersurface of the parameter space outside of which the
#    number of real solutions is locally constant: it is constant on each connected component
#    of its complement;
# 2. one sample point in each component, and one exact solving there, give this number
#    everywhere.
#
# This answers the question left open by the worksheet "Systems with uncertain coefficients":
# when does a family of systems keep the real solutions of one instance, and no other?

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

svg := name -> cat(currentdir(), "/", name, ".svg"):
field := (h, f) -> rhs(select(e -> lhs(e) = f, fields(h))[1]):
# the number of real solutions of a zero-dimensional system
nreal := sys -> nops(real_solutions(sys, prec = 32)):
# the Jacobian determinant of F with respect to vars
jacobian_det := (F, vars) -> expand(LinearAlgebra:-Determinant(
                                 Matrix([seq([seq(diff(f, v), v in vars)], f in F)]))):

# The polynomials of a call are sent in a common ring whose variables are sorted by name,
# unless ring = [...] is given: the unknowns and the parameters are passed explicitly to
# discriminant_variety, and for a curve the first variable is the one of the projection.

# ## 1. One parameter
#
# The unit circle and the horizontal line y = c. The unknowns are x, y, the parameter is c:

line_circle := [x^2 + y^2 - 1, y - c]:
printf("## 1. %a\ndiscriminant variety: %a\n", line_circle,
       discriminant_variety(line_circle, [x, y], [c], ring = [x, y, c]));

# The discriminant variety is {c = -1} union {c = 1}: the line is tangent to the circle. Its
# complement in R has three components, c < -1, -1 < c < 1 and c > 1; one sample point in
# each, and the system solved there:

printf("[c, number of real solutions]: %a\n",
       [seq([cv, nreal(eval(line_circle, c = cv))], cv in [-2, 0, 2])]);

# Two real solutions for -1 < c < 1, none otherwise. In the worksheet on uncertain
# coefficients, c = 1 + 1/1000 +- 1/100: this interval contains the point c = 1 of the
# discriminant variety, so that the number of real solutions is not the same for all the
# systems of the family. This is what the exact instance c = 1 + 1/1000 could not see.

# ## 2. Two parameters
#
# The family of the worksheet on uncertain coefficients, with its two coefficients now
# parameters a and c:

family := [y^2 - x^3 + a * x, x^2 + y^2 - c]:
dv := discriminant_variety(family, [x, y], [a, c], ring = [x, y, a, c]):
printf("\n## 2. %a\ndiscriminant variety: %a\n", family, dv);

# Three components: the line a = c, the line c = 0 and a cubic curve. Together they form a
# plane curve, whose topology and complement are computed as in the worksheet "Topology of a
# real plane curve", a being the first variable:

Dp := expand(mul(mul(f, f in comp), comp in dv)):
T := curve_topology(Dp, ring = [a, c]):
K := curve_complement(T):
curve_summary(T, K);

# Eight components. In each, the system is solved at the sample point:

count_at := (av, cv) -> nreal(eval(family, [a = av, c = cv])):
printf("[component, sample (a, c), number of real solutions]:\n");
for k to nops(K) do
    s := field(K[k], sample);
    printf("  %a\n", [k, s, count_at(op(s))]);
end do:

# The plane of the parameters, each component in its colour (numbered as above):

printf("written: %s\n", write_curve_svg(svg("parameters"), T,
       curve_drawing(T, window = [-3, 3, -3, 4], pixels = [600, 600]), K,
       title = "discriminant variety in the (a, c) plane"));

# The two bounded components, #4 and #5, are small; around the origin:

printf("written: %s\n", write_curve_svg(svg("parameters-zoom"), T,
       curve_drawing(T, window = [-1/2, 1/4, -1/16, 1/16], pixels = [600, 400]), K,
       title = "around the origin"));

# The systems have 0, 2 or 6 real solutions according to the component. The number is the
# same at every point of a component: the sample point stands for all of them.

# ## 3. Back to the uncertain family
#
# The worksheet on uncertain coefficients took a = 1 +- e and c = 3 +- e. If this box does
# not meet the discriminant variety, it lies in one component of its complement, and all the
# systems of the family have the same number of real solutions as the instance a = 1, c = 3,
# that is 2.
#
# A sufficient condition: no polynomial of the discriminant variety vanishes on the box,
# which interval arithmetic proves when the enclosure of its values does not contain 0:

DVQ := [seq(op(comp), comp in dv)]:
uncertain := (q, e) -> INTERVAL(q - e .. q + e):
meets := e -> map(has_zero, evaluate(DVQ, [uncertain(1, e), uncertain(3, e)], precision = 64,
                                    ring = [a, c])):
printf("\n## 3. The box a = 1 +- e, c = 3 +- e\npolynomials: %a\n", DVQ);
for e in [1/100, 1/10, 1/2, 1] do printf("  e = %-5a meets: %a\n", e, meets(e)) end do:

# Up to e = 1/2, no component vanishes on the box: every system of the family has exactly two
# real solutions. With the worksheet on uncertain coefficients, which enclosed the two
# solutions of every system for e <= 1/10, the description is complete: two real solutions,
# one in each box, and no other.
#
# At e = 1 the test fails for a - c and for the cubic. An interval enclosure may contain 0
# without the polynomial vanishing; here the box does meet them: a = c = 2 is in it, and the
# cubic changes sign between two corners:

cub := op(select(f -> degree(f) = 3, DVQ)):
printf("signs of the cubic at (0, 2) and (2, 2): %a\n",
       [seq(signum(eval(cub, [a = pt[1], c = pt[2]])), pt in [[0, 2], [2, 2]])]);

# ## 4. More parameters: a partial CAD
#
# With more than two parameters, the components of the complement are no longer given by a
# curve. partial_CAD gives sample points of the cells of a *cylindrical algebraic
# decomposition* adapted to the polynomials of the discriminant variety, in any number of
# variables; each component of the complement contains at least one cell. Here, in the plane
# of (a, c):

points := partial_CAD(DVQ, ring = [a, c]):
printf("\n## 4. partial_CAD: %d cells\n", nops(points));
printf("[sample, number of real solutions]:\n");
for pt in points do printf("  %a\n", [pt, count_at(op(pt))]) end do:

# 23 cells for 8 components: the decomposition cuts the plane by vertical lines, and several
# cells lie in the same component. The numbers of real solutions found are again 0, 2 and 6.

# ## 5. Square systems: the discriminant variety from LACE
#
# discriminant_variety handles any system that is generically zero-dimensional in the
# unknowns, and splits its result into irreducible components (with DiscriminantVariety.jl).
# When the system has as many equations as unknowns, the discriminant variety is the union of
# two sets that the C library of LACE computes directly:
#
# - W_c, the *critical values*: the projection on the parameters of the solutions where the
#   Jacobian determinant with respect to the unknowns vanishes, that is the elimination ideal
#   of the system augmented by this determinant (eliminate);
# - W_infinity, the parameters where solutions go to infinity (ideals_at_infinity).
#
# No factorization is needed. For the family of section 2:

Wc := eliminate([op(family), jacobian_det(family, [x, y])], [x, y], ring = [x, y, a, c]):
printf("\n## 5. W_c = %a\n", factor(Wc));
printf("ideals_at_infinity: %a\n", ideals_at_infinity(family, [x, y], ring = [x, y, a, c]));

# W_c is the polynomial of section 2, the product of the three components, and W_infinity is
# empty ([1] for each unknown): no solution goes to infinity.

# ### Solutions at infinity
#
# For a x^2 + y^2 = 1 and x = y, that is (a + 1) x^2 = 1:

F5 := [a * x^2 + y^2 - 1, x - y]:
printf("\n%a\nW_c = %a, at infinity: %a\n", F5,
       eliminate([op(F5), jacobian_det(F5, [x, y])], [x, y], ring = [x, y, a]),
       ideals_at_infinity(F5, [x, y], ring = [x, y, a]));

# W_c is empty: the solutions are never double. But at a = -1 the two solutions
# +-1/sqrt(a + 1) go to infinity, and W_infinity gives a + 1. On each side:

printf("[a, number of real solutions]: %a\n", [seq([av, nreal(eval(F5, a = av))], av in [-2, 0])]);

# ### A larger system
#
# Three equations in three unknowns u, v, w, two parameters p, q:

F6 := [u^2 + p * v^2 + w - 1, v^2 - q * u * w + u - 2, w^2 + u * v + p * q * w - 3]:
R6 := [u, v, w, p, q]:
t := time[real]():
Wc6 := eliminate([op(F6), jacobian_det(F6, [u, v, w])], [u, v, w], ring = R6):
zi := ideals_at_infinity(F6, [u, v, w], ring = R6):
t_lace := time[real]() - t:
Winf6 := zi[2]:
printf("\n%a\ndegrees of W_c: %a, at infinity: %a, %.2f s\n", F6, map(degree, Wc6), Winf6, t_lace);

# W_c is a curve of degree 45, and solutions go to infinity on p^3 q^4 + 1 = 0. The same
# discriminant variety with DiscriminantVariety.jl alone (set_dv_backend(julia)), after a
# first call on the small system above, so that its compilation is not counted:

set_dv_backend(julia):
discriminant_variety(F5, [x, y], [a], ring = [x, y, a]):
t := time[real]():
discriminant_variety(F6, [u, v, w], [p, q], ring = R6):
t_dv := time[real]() - t:
set_dv_backend(lace):
printf("LACE: %.2f s, DiscriminantVariety.jl: %.2f s\n", t_lace, t_dv);

# The regions of the (p, q) plane, and the number of real solutions in each:

D6 := expand(mul(f, f in remove(f -> f = 1, [op(Wc6), seq(op(w), w in Winf6)]))):
T6 := curve_topology(D6, ring = [p, q]):
K6 := curve_complement(T6):
printf("%d regions\n", nops(K6));
count6 := (pv, qv) -> nreal(eval(F6, [p = pv, q = qv])):
counts6 := [seq(count6(op(field(k, sample))), k in K6)]:
printf("[number of real solutions, number of regions]: %a\n",
       [seq([n, numboccur(counts6, n)], n in sort([op({op(counts6)})]))]);

# 42 regions, with 0 to 8 real solutions. They spread from p ~ -150 to |q| ~ 10^9; the
# window [-6, 2] x [-2, 2] of the (p, q) plane:

printf("written: %s\n", write_curve_svg(svg("pq"), T6,
       curve_drawing(T6, window = [-6, 2, -2, 2], pixels = [600, 600]), K6,
       title = "discriminant variety in the (p, q) plane"));

# All of them, up to a homeomorphism of the plane:

printf("written: %s\n", write_curve_svg_topo(svg("pq-topo"), T6, K6));

# ## 6. A constraint g <> 0: a discriminant variety by hand
#
# The system of section 2, with its second equation multiplied by x, and the constraint
# x <> 0:

F7 := [y^2 - x^3 + a * x, x * (x^2 + y^2 - c)]:
R7 := [x, y, a, c]:

# The point (0, 0) is a solution for every value of the parameters, and the Jacobian
# determinant vanishes there: the critical values fill the whole plane, and the elimination
# ideal is the zero ideal.

J7 := jacobian_det(F7, [x, y]):
printf("\n## 6. %a, x <> 0\nelimination ideal of F7 and J7: %a\n", F7,
       eliminate([op(F7), J7], [x, y], ring = R7));

# The constraint removes this solution. saturate(F, g, vars) computes
# (<F> : g^infinity) inter Q[params], the elimination ideal without the solutions where g
# vanishes; with no variable to eliminate, it gives the saturated system itself, whose
# solutions are those of F with x <> 0:

S7 := saturate(F7, x, [], ring = R7):
printf("saturated system: %a\n", S7);

# The discriminant variety of the system under the constraint is assembled from three parts,
# each computed by LACE:
#
# - W_c, the critical values, away from x = 0;
# - W_g, the parameters where a solution reaches x = 0: the elimination ideal of the
#   saturated system augmented by x;
# - W_infinity, the parameters where solutions go to infinity.

Wc7 := saturate([op(F7), J7], x, [x, y], ring = R7):
Wg7 := eliminate([op(S7), x], [x, y], ring = R7):
printf("W_c = %a\nW_g = %a\nat infinity: %a\n", map(factor, Wc7), Wg7,
       ideals_at_infinity(F7, [x, y], ring = R7));

# W_c is (a - c) times the cubic of section 2, W_g is the line c = 0, and no solution goes to
# infinity. Their union is the curve of section 2:

D7 := expand(mul(f, f in remove(f -> f = 1, [op(Wc7), op(Wg7)]))):
printf("the same curve as in section 2: %a\n", evalb(expand(D7 - Dp) = 0 or expand(D7 + Dp) = 0));

# As expected: with x <> 0, the system is the one of section 2. In each component of the
# complement, the solutions with x <> 0 are those of the saturated system; the system itself
# has one more, the point (0, 0):

printf("[component, solutions with x <> 0, all solutions]:\n");
for k to nops(K) do
    s := field(K[k], sample);
    printf("  %a\n", [k, nreal(eval(S7, [a = s[1], c = s[2]])), nreal(eval(F7, [a = s[1], c = s[2]]))]);
end do: