Roadmap algorithm

docs/make.jl

@@ -1,6 +1,9 @@
+#using Pkg
+#Pkg.activate("../../AlgebraicSolving.jl/")
 using AlgebraicSolving
 using Documenter
 
+
 push!(LOAD_PATH, "../src")
 
 DocMeta.setdocmeta!(AlgebraicSolving, :DocTestSetup, :(using AlgebraicSolving); recursive=true)
@@ -24,7 +27,8 @@ makedocs(
                          "hilbert.md",
                          "dimension.md",
                          "solvers.md",
-                         "decomposition.md"],
+                         "decomposition.md",
+                         "connectivity.md"],
         "Examples" => "katsura.md"
         ]
 )

docs/src/connectivity.md

@@ -0,0 +1,48 @@
+```@meta
+CurrentModule = AlgebraicSolving
+DocTestSetup = quote
+  using AlgebraicSolving
+end
+```
+
+```@setup algebraicsolving
+using AlgebraicSolving
+```
+
+```@contents
+Pages = ["connectivity.md"]
+```
+
+# Algorithms for connectivity queries on smooth bounded algebraic sets
+
+## Introduction
+
+AlgebraicSolving allows to compute a roadmap for the real trace of the zero-set
+of the ideal spanned by given input generators over the rationals.
+
+It assumes that the underlying algebraic set is **smooth**, and its real trace is **bounded**.
+
+The underlying engine is provided by msolve.
+
+## Functionality
+
+```@docs
+    roadmap(
+      I::Ideal{P} where P <: QQMPolyRingElem;
+      C::Vector{Vector{Vector{QQFieldElem}}}=Vector{Vector{QQFieldElem}}[],
+      info_level::Int=0,
+      checks::Bool=false
+    )
+```
+
+In addition, AlgebraicSolving can compute equations definition critical loci of polynomial maps over the given algebraic set.
+
+```@docs
+   computepolar(
+        J::Union{Vector{Int},UnitRange{Int}},
+        V::Ideal{P};
+        phi::Vector{P} = P[],
+        dimproj = length(J)-1,
+        only_mins = false
+    ) where (P <: MPolyRingElem)
+```

docs/src/types.md

@@ -17,21 +17,21 @@ Pages = ["types.md"]
 
 ## Introduction
 
-AlgebraicSolving handles ideals in multivariate polynomial rings over a prime 
+AlgebraicSolving handles ideals in multivariate polynomial rings over a prime
 field of characteristic $0$ or $p$ where $p$ is a prime number $<2^{31}$.
 
 ## Polynomial Rings
 
-We use [Nemo](https://www.nemocas.org/)'s multivariate polynomial 
+We use [Nemo](https://www.nemocas.org/)'s multivariate polynomial
 ring structures:
 
 ```@repl
 using AlgebraicSolving
 R, (x,y,z) = polynomial_ring(QQ, ["x", "y", "z"], internal_ordering=:degrevlex)
 ```
-The above example defines a multivariate polynomial ring in three variables `x`, 
-`y`, and `z` over the rationals using the dgree reverse lexicographical ordering 
-for printing polynomials in the following. One can also define polynomial rings 
+The above example defines a multivariate polynomial ring in three variables `x`,
+`y`, and `z` over the rationals using the dgree reverse lexicographical ordering
+for printing polynomials in the following. One can also define polynomial rings
 over finite fields:
 
 ```@repl
@@ -41,13 +41,149 @@ R, (x,y,z) = polynomial_ring(GF(101), ["x", "y", "z"], internal_ordering=:degrev
 
 ## Ideals
 
-Ideals can be constructed by giving an array of generators. Ideals cache varies 
-data structures connected to ideals in order to make computational algebra more 
-effective:
+Ideals can be constructed by giving an array of generators:
 
 ```@repl
 using AlgebraicSolving
 R, (x,y,z) = polynomial_ring(QQ, ["x", "y", "z"], internal_ordering=:degrevlex)
 I = Ideal([x+y+1, y*z^2-13*y^2])
 ```
 
+Ideals cache various data structures connected to ideals in order to make
+computational algebra more effective. Assume that `T <: MPolyRingElem`, then
+  * `gens::Vector{T}`: generators of the Ideal provded by the user;
+  * `deg::Int`: degree of the ideal (`Nothing` if unknown);
+  * `dim::Int`: Krull dimension of the Ideal (`Nothing` if unknown);
+  * `gb::Dict{Int, Vector{T}}`: Gröbner bases of the Ideal where the key
+    correspond to the number of variables eliminated (starting with the first);
+  * `inter_sols::Vector{Vector{Vector{QQFieldElem}}}`: intervals with rational
+    endpoints, containing the real solutions of the system (if `dim=0`);
+  * `real_sols::Vector{Vector{QQFieldElem}}`: midpoints of the above intervals;
+  * `rat_sols::Vector{Vector{QQFieldElem}}`: rational solutions of the system
+    (if `dim=0`).
+  * `rat_param::Union{RationalParametrization, RationalCurveParametrization}`:
+    rational parametrization encoding the complex solution set of the Ideal (see
+    below).
+
+
+## Rational Parametrizations
+Rational Parametrizations are special data structure for representing equidimensionnal
+algebraic sets of dimension 0 and 1. They are typical outputs of algorithm for solving
+polynomial systems encoding such algebraic sets thanks to their nice behaviour with
+subsequent computations.
+
+### Zero-dimensionnal parametrizations: finite set of points (see also [there](https://msolve.lip6.fr/downloads/msolve-tutorial.pdf#section.6))
+
+A *zero-dimensional parametrization* $\mathscr{P}$ with
+coefficients in a field $\mathbb{Q}$ consists of:
+* polynomials $(w, w', v_1, \ldots, v_n)$ in $\mathbb{K}[t]$ where $t$ is
+  a new variable and such that
+  *  $w$ is a monic square-free polynomial;
+  *  $w'=1$ when $\mathbb{K}$ is a prime field, $w'=\frac{dw}{dt}$ else;
+  *  $\deg(v_i) < \deg(w)$.
+* a linear form $l$  in the variables $x_1,\dotsc, x_n$,
+such that
+$l(\rho_1,\dotsc,\rho_n) = t \mod w$.
+
+Such a data-structure encodes the following finite set of points
+*  $\left\{\left(\frac{\rho_1(\vartheta)}{w'(\vartheta)}, \ldots, \frac{\rho_n(\vartheta)}{w'(\vartheta)} \right) \middle|\, w(\vartheta) = 0\right\}$.
+According to this definition, the roots of $w$ are exactly the
+values taken by $l$ on this set.
+
+The type `RationalParametrization` therefore caches the following attributes:
+  * `vars::Vector{Symbol}`: variables used for the parametrization the last one playing the role of the above $t$
+  (hence, maybe with one more variable than the ones given as input);
+  * `cfs_lf::Vector{ZZRingElem}`: coefficients on the linear form $l$ (when variables are added);
+  * `elim::QQPolyRingElem`: elimination polynomial $w$ ;
+  * `denom::QQPolyRingElem`: denominator polynomial (usually $w'$);
+  * `param::Vector{QQPolyRingElem}`: numerators $\rho_i$'s.
+
+See the documentation of the [`rational_parametrization`](@ref) function for further details.
+
+### One-dimensionnal parametrizations: algebraic curves
+
+A *one-dimensional parametrization* $\mathscr{C}$ with
+coefficients in a field $\mathbb{C}$ consists of:
+* polynomials $(w, w', v_1, \ldots, v_n)$ in $\mathbb{K}[t,s]$ where $t,s$ are
+  new variables and such that
+    *  $w$ is a monic square-free polynomial and $w'=\frac{\partial w}{\partial t}$;
+    *  $\deg(v_i) < \deg(w)$.
+*  two linear forms $(l,l')$  in the variables $x_1,\dotsc, x_n$,
+such that
+$l(\rho_1,\dotsc,\rho_n) = t\, w'\mod w \quad \text{and} \quad l'(\rho_1,\dotsc,\rho_n) = s\,w' \mod w$
+
+Such a data-structure encodes the curve defined as the Zariski closure of the following set
+*  $\left\{\left(\frac{\rho_1(\vartheta,\eta)}{w'(\vartheta,\eta)}, \ldots, \frac{\rho_n(\vartheta,\eta)}{w'(\vartheta,\eta)} \right) \middle|\, w(\vartheta,\eta) = 0,\, w'(\vartheta, \eta) \neq 0\right\}$.
+According to this definition, the roots of $w$ are exactly the
+values taken by $l$ on this set.
+
+The type `RationalCurveParametrization` therefore caches the following
+attributes:
+  * `vars::Vector{Symbol}`: variables used for the parametrization the last two
+  ones playing the role of the above $(t,s)$ (hence, maybe with up to two more
+  variable than the ones given as input);;
+  * `cfs_lfs::Vector{Vector{ZZRingElem}}`: coefficients on the linear form $l$
+    (when variables are added);
+  * `elim::QQMPolyRingElem`: elimination polynomial $w$ ;
+  * `denom::QQMPolyRingElem`: denominator polynomial (usually $w'$);
+  * `param::Vector{QQMPolyRingElem}`: numerators $\rho_i$'s.
+
+See the documentation of the [`rational_curve_parametrization`](@ref) function for
+further details.
+
+## Roadmaps
+
+Consider an algebraic set $V\subset \mathbb{C}^n$ and a finite set of query
+points $\mathcal{P} \subset V$, both defined by polynomials with coefficients in
+$\mathbb{Q}$. A *roadmap* $\mathcal{R}$ associated to $(V,\mathcal{P})$ is an
+algebraic curve such that
+  *  $\mathcal{P} \subset \mathcal{R} \subset V$;
+  *  $C \cap \mathcal{R}$ is non-empty and connected, for each connected
+    component $C$ of $V\cap \mathbb{R}^n$.
+
+Roadmaps are algebraic curves capturing the connectivity properties of an
+algebraic set containing it so that points in $\mathcal{P}$ are equivalently
+connected both on $\mathcal{R}$ and $V\cap \mathbb{R}^n$.
+
+The effective construction of roadmaps relies on connectivity statements which
+allow one to construct real algebraic subsets of $V\cap \mathbb{R}^n$, of
+smaller dimension, having a connected intersection with the connected components
+of $V\cap \mathbb{R}^n$ Therefore, the output is union of several curves,
+organized in a tree due to the recursive nature of the roadmap algorithms.
+
+```julia
+mutable struct Roadmap
+    initial_ideal::Ideal{QQMPolyRingElem}
+    root::RMnode
+end
+```
+
+Hence, it is composed of a main node, containing the equations of the initial
+algebraic set $V$ and a reference to the root node of the tree.
+
+```julia
+mutable struct RMnode
+    base_pt::Vector{QQFieldElem}
+    polar_eqs::Vector{QQMPolyRingElem}
+    children::Vector{RMnode}
+end
+```
+
+Each roadmap node (including the root) correspond to a curve component of the
+roadmap. More precisely, they are defined by as an algebraic subset $W$ (called
+polar variety) of a fiber $F$ of the initial variety $V$ such that:
+  * if `base_pt` contains $\mathbf{q}=(q_1,\dotsc,q_e) \in \mathbb{Q}^e$ then
+    *  $F = V \cap \left\{x_1=q_1,\,\dotsc,\, x_e=q_e\right\}$;
+  * if `polar_eqs` contains the polynomials $g_1,\dotsc,g_s \in   \mathbb{Q}[x_1,\dotsc,x_n]$ then 
+    *  $W = F \cap \{g_1=0,\, \ldots,\, g_s = 0\}$.
+
+Moreover, `children` rassembles all the tree nodes that contains curves
+component for which their attribute `base_pt` contains a point $\mathbf{q}'$
+that shares the same first $e$ coordinates with $\mathbf{q}$.
+
+The reason of this tree arrangement is mainly because of the following property:
+*two roadmap components intersect if, and only if, one's node is a descendant
+of the other's*. In many cases, one can even replace "descendant" by "parent".
+
+
+See the documentation of the [roadmap](@ref) function for further details.

src/AlgebraicSolving.jl

@@ -18,6 +18,9 @@ include("algorithms/param-curve.jl")
 include("algorithms/hilbert.jl")
 #= siggb =#
 include("siggb/siggb.jl")
+#= connectivity =#
+include("connectivity/polar.jl")
+include("connectivity/roadmap.jl")
 #= examples =#
 include("examples/katsura.jl")
 include("examples/cyclic.jl")

src/algorithms/param-curve.jl

@@ -232,7 +232,8 @@ function _evalvar(
     return LFeval
 end
 
-# Generate N primes > start that do not divide any numerator/denominator
+# Generate N random primes between low and up
+# that do not divide any numerator/denominator
 # of any coefficient in polynomials from LP
 function _generate_lucky_primes(
     LF::Vector{P} where P<:MPolyRingElem,

src/algorithms/solvers.jl

@@ -272,7 +272,7 @@ end
     rational_solutions(I::Ideal{T} where T <: MPolyRingElem, <keyword arguments>)
 
 Given an ideal `I` with a finite solution set over the complex numbers, return
-the rational roots of the ideal. 
+the rational roots of the ideal.
 
 # Arguments
 - `I::Ideal{T} where T <: MPolyRingElem`: input generators.

src/connectivity/polar.jl

@@ -0,0 +1,112 @@
+@doc Markdown.doc"""
+    computepolar(J::Union{Vector{Int},UnitRange{Int}}, V::Ideal{P} where P <: MPolyRingElem; <keyword arguments>)
+
+Compute the **polar variety** associated with the map whose components are
+`[phi_1, …, phi_p, x_{p+1}, …, x_n]` indexed by `J`, for an algebraic variety defined
+by the ideal `V`.
+
+More precisely, this function computes the set of points `x` in `V` such that,
+if `psi` denotes the above map, then the image of the tangent space `Tₓ(V)` under
+the differential `psi` has dimension strictly less than `dimproj`.
+
+This is a key geometric construction in the computation of polar varieties,
+used in the roadmap algorithm and critical point methods.
+
+# Arguments
+- `J::Union{Vector{Int},UnitRange{Int}}`:
+  Indices of the selected coordinate functions
+- `V::Ideal{P} where P <: MPolyRingElem`:
+  Input ideal defining the variety V on which critical loci is computed
+- `phi::Vector{P}=P[]`:
+  Polynomial map possibly completed by projection coordinates to have a map of total length `n`.
+- `dimproj::Int=length(J)-1`:
+  Expected maximal dimension of the image of the tangent space.
+  Typically equals the number of projection coordinates minus one.
+- `only_mins::Bool=false`:
+  If `true`, only the computed minors of the Jacobian are returned;
+  otherwise, the output includes both the generators of `V` and those minors.
+
+# Returns
+- If `only_mins=false` (default):
+  A `Vector{P}` containing the union of the generators of `V` and the computed minors.
+  This defines the ideal of the polar variety.
+- If `only_mins=true`:
+  A `Vector{P}` containing only the minors (without the original equations of `V`).
+
+# Example
+```jldoctest
+julia> using AlgebraicSolving
+
+julia> R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"])
+(Multivariate polynomial ring in 4 variables over QQ, QQMPolyRingElem[x1, x2, x3, x4])
+
+julia> I = Ideal([(x1^2+x2^2+x3^2+x4^2+9-1)^2-4*9*(x1^2+x2^2+x3^2) + 1])
+QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65]
+
+julia> computepolar(1:3, I, dimproj=1, phi=[x1^2+x2^2+x3^2+x4^2])
+5-element Vector{QQMPolyRingElem}:
+ x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65
+ -4*x1^3 - 4*x1*x2^2 - 4*x1*x3^2 - 4*x1*x4^2 + 40*x1
+ -4*x1^2*x4 - 4*x2^2*x4 - 4*x3^2*x4 - 4*x4^3 - 32*x4
+ -2*x1
+ -2*x4
+```
+"""
+function computepolar(
+        J::Union{Vector{Int},UnitRange{Int}},   # coordinate images of [phi,proj] (as above)
+        V::Ideal{P};                            # input ideal
+        phi::Vector{P} = P[],                   # polynomial map in consideration (completed by sufficiently many projections)
+        dimproj = length(J)-1,                  # maximum dimension of tangent space of phi
+        only_mins = false                       # return only minors without eqns of V
+    ) where (P <: MPolyRingElem)
+
+    R = parent(V)
+    n = nvars(R)
+    nphi = length(phi)
+    @assert all([1<=j<=n for j in J])
+
+    isnothing(V.dim) && dimension(V)
+    c = n - V.dim
+
+    ##
+    Jphi = [ j for j in J if j <= nphi ]
+    Jproj = setdiff(J, Jphi)
+
+    # Construct the truncated Jacobian matrix
+    psi = vcat(V.gens, phi[Jphi])
+    JW = matrix(R, QQMPolyRingElem[ derivative(f, k) for f in psi, k in setdiff(1:n, Jproj)])
+    # Compute the minors
+    sizeminors = c + length(Jphi) + min(dimproj, length(J)-1) - (length(J)-1)
+    minors = _compute_minors(sizeminors, JW)
+
+    if only_mins
+        return minors
+    else
+        return vcat(V.gens, minors)
+    end
+end
+
+function _compute_minors(p, A)
+    # Computes the p-minors of a matrix A
+    rowsmins = _combinations(1:nrows(A), p)
+    colsmins = _combinations(1:ncols(A), p)
+    # We use charpoly for a division-free determinant method
+    return [ coeff(charpoly(A[rows, cols]), 0) for rows in rowsmins for cols in colsmins ]
+end
+
+function _combinations(v::UnitRange{Int}, k::Int)
+    # Compute the k-subsets of v
+    n = length(v)
+    ans = Vector{Int}[]
+    k > n && return ans
+    _combinations_dfs!(ans, Vector{Int}(undef, k), v, n, k)
+    return ans
+end
+
+function _combinations_dfs!(ans::Vector{Vector{Int}}, comb::Vector{Int}, v::UnitRange{Int}, n::Int, k::Int)
+    k < 1 && (pushfirst!(ans, comb[:]); return)
+    for m in n:-1:k
+        comb[k] = v[m]
+        _combinations_dfs!(ans, comb, v, m - 1, k - 1)
+    end
+end