diff --git a/Project.toml b/Project.toml index 1d92704..fac42ed 100644 --- a/Project.toml +++ b/Project.toml @@ -33,4 +33,4 @@ REPL = "1.6" StaticArrays = "1" Test = "1.6" julia = "1.10" -msolve_jll = "0.900.502" +msolve_jll = "0.900.502" \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index 9c761f2..563d01f 100644 --- a/docs/make.jl +++ b/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" ] ) diff --git a/docs/src/connectivity.md b/docs/src/connectivity.md new file mode 100644 index 0000000..b6b20fc --- /dev/null +++ b/docs/src/connectivity.md @@ -0,0 +1,138 @@ +```@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) +``` + +# Algorithms for connectivity queries on real algebraic curve + +## Introduction + +AlgebraicSolving allows to compute a graph whose associated (3D) piecewise linear curve is semi-algebraically homeomorphic to a real algebraic curve given as input. +This input can be: + * an `Ideal` given by generators over the rationals, whose real trace is the curve; + * a one-dimensionnal parametrization given as a `RationalCurveParametrization`. We refer to the `Types` page for such a structure. + +It is also possible to compute the arrangement of several curves by simply providing the corresponding vector of Ideals. + +In both cases, user can, in addition, provide control points that will be correctly identified in the final graph, from a topological point of view. + +The underlying engines are + * internal subresultants computations for paramtrization of critical points; + * `msolve` for univariate real root isolation. + +## Computing curve graphs + +```@docs + curve_graph( + I::Ideal{P}, args...; + generic=true, + outf=true, + v=0, + kwargs... + ) where {P <: QQMPolyRingElem} +``` + +```@docs + curve_arrangement_graph( + curves::Vector{Ideal{P}}; + generic=true, + outf=true, + v=0, + kwargs... + ) where (P <: MPolyRingElem) +``` + +## Computing with curve graphs + +Once such a graph is computed, it is encoded in a CurveGraph object, whose description can be found in the Types section of this documentation. +These are nothing but sets of verteices and nodes. +One can perform different operations on such strutures. + +```@docs + connected_components( + G::CurveGraph{T} + ) where T +``` + +```@docs + group_by_components( + G::CurveGraph + ) +``` + +```@docs + number_of_connected_components( + G::CurveGraph + ) +``` + +```@docs + merge_graphs( + graphs::Vector{CurveGraph} + ) +``` + +## Displaying curve graph + +As e.g. Plots.jl is not part of AlgebraicSolving.jl, it cannot provide direct functions for plotting curve graphs. However, functions for exporting curve graphs in a format adapted to most ploting libraries, encoded in a `GraphPlotData` object. We refer to the Types section for more information on this structuree. + +```@docs + build_graph_data( + G::CurveGraph{T}; + width=3.0, + vert=true, + color="#FF0000" + ) where T <: Union{Float64,QQFieldElem} +``` + +```@docs + build_graph_data( + CG::Vector{CurveGraph{T}}; + width=3.0, + vert=true + ) where T <: Union{Float64,QQFieldElem} +``` \ No newline at end of file diff --git a/docs/src/types.md b/docs/src/types.md index b60edce..af668a5 100644 --- a/docs/src/types.md +++ b/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,9 +41,7 @@ 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 @@ -51,3 +49,244 @@ 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 intersects 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. + +# Graph Data Structures & Visualization + +## Curve Graphs + +The `CurveGraph{T}` structure represents a planar straight-line graph +that is homeomorphic to a real algebraic curve. It is the output of +algorithms [compute_graph](@ref), which compute the planar projection +of curves while resolving apparent singularities. + +These real algebraic curves are typically output of roadmap algorithms +presented above. Hence `CurveGraph{T}` object can provide a +discretized skeleton of real algebraic sets that preserves +connectivity properties. In particular they share the same number of +connected components (respectively as graphs and as semi-algebraic +sets) + + +The type `CurveGraph{T}` (where `T` is typically `Float64` or `QQFieldElem`) caches the following attributes: + + * `vertices::Vector{Tuple{T, T}}`: a list of 2D coordinates representing the nodes of the graph (critical points, regular routing points, and user-defined control points). + + * `edges::Vector{Tuple{Int, Int}}`: a list of index pairs defining the undirected edges connecting the vertices in `vertices`. + + * `control_nodes::Dict{Int, Vector{Int}}`: a dictionary mapping user-defined control point IDs to their corresponding local vertex indices in the `vertices` array. These IDs typically refer to other `CurveGraph{T}` and the associated vertices: their mutual intersection. This is useful when computing arrangement of curves using [merge_graphs](@ref). + +```julia +struct CurveGraph{T} + vertices::Vector{Tuple{T, T}} + edges::Vector{Tuple{Int, Int}} + control_nodes::Dict{Int, Vector{Int}} +end +``` + +## Graph Plotting Data + +To seamlessly visualize `CurveGraph` objects using standard plotting libraries (such as `Plots.jl`), AlgebraicSolving provides an intermediate geometric representation. + +The primary structure is `GraphPlotData`, which decomposes a mathematical graph into purely visual components: + +* `edges::EdgeGroup`: contains the geometric line segments representing the graph's connections, alongside styling data (color and width). + +* `points::Vector{PointGroup}`: separates standard vertices and special control nodes into distinct scatter groups with specific markers (e.g., `:x` for standard vertices, `:+` for control nodes). + +```julia +struct EdgeGroup + edges::Vector{Tuple{Vector{Float64}, Vector{Float64}}} + color::String + width::Float64 +end + +struct PointGroup + vertices::Tuple{Vector{Float64}, Vector{Float64}} + color::String + marker::Symbol +end + +struct GraphPlotData + edge_group::EdgeGroup + point_groups::Vector{PointGroup} +end +``` + +### Core Builders + +The functions [build_graph_data](@ref) allow to compute such `GraphPlotData` from `CurveGraph{T}` objects. These former objects can be plotted with `Plot.jl` as follows. Below, `G` is a `CurveGraph{T}`. + +```julia +using Plots + +function plot_graph(P) + plot(legend=false) + E = P.edge_group + plot!(E.edges, color=E.color, width = E.width) + [scatter!(V.vertices, color=V.color, marker = V.marker) for V in P.point_groups] + gui() +end + +plot_graph(build_graphs_data(G)) +``` + +If one has `CG` a list of `CurveGraph{T}` (typically the connected +components of `G` obtained with [connected_components](@ref)). + +```julia +# How to plot: +using Plots + +function plot_graphs(CP) + plot(legend=false) + for P in CP + E = P.edge_group + plot!(E.edges, color=E.color, width = E.width) + [ scatter!(V.vertices, color=V.color, marker = V.marker) for V in P.point_groups ] + end + gui() +end + +plot_graphs(build_graph_data(CG)) +``` + +This will plot graphs in CG with distinct and distinguishable colors. +See the documentation of the [build_graph_data](@ref) function for further details. \ No newline at end of file diff --git a/src/AlgebraicSolving.jl b/src/AlgebraicSolving.jl index f48e3d3..50ac27c 100644 --- a/src/AlgebraicSolving.jl +++ b/src/AlgebraicSolving.jl @@ -22,8 +22,11 @@ include("siggb/siggb.jl") include("progress/main.jl") #= interp =# include("interp/main.jl") +#= connectivity =# +include("connectivity/connectcurves.jl") #= examples =# include("examples/katsura.jl") include("examples/cyclic.jl") + end # module AlgebraicSolving diff --git a/src/algorithms/param-curve.jl b/src/algorithms/param-curve.jl index 1e553f2..55512e6 100644 --- a/src/algorithms/param-curve.jl +++ b/src/algorithms/param-curve.jl @@ -103,7 +103,7 @@ function rational_curve_parametrization( # Compute parametrization of each evaluation Lr = Vector{RationalParametrization}(undef, length(free_ind)) for j in 1:length(free_ind) - info_level>0 && print("Evaluated parametrizations: $(j+DEG+2-length(free_ind))/$(DEG+2)", "\r") + info_level>0 && print("Evaluated parametrizations: $(j)/$(length(free_ind))", "\r") Lr[j] = rational_parametrization(LFeval[j], nr_thrds=nr_thrds) end info_level>0 && println() @@ -118,8 +118,8 @@ function rational_curve_parametrization( fact = lc / leading_coefficient(Lr[j].elim) rr = [ p*fact for p in vcat(Lr[j].elim, Lr[j].denom, Lr[j].param) ] end - PARAM[j] = rr - _values[j] = curr_values[j] + PARAM[free_ind[j]] = rr + _values[free_ind[j]] = curr_values[j] used_ind[j] = true else info_level>0 && println("bad specialization: ", i+j-1) @@ -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, diff --git a/src/algorithms/solvers.jl b/src/algorithms/solvers.jl index 05805a2..7646ee5 100644 --- a/src/algorithms/solvers.jl +++ b/src/algorithms/solvers.jl @@ -272,7 +272,7 @@ end rational_solutions(I::Ideal{T} where T <: MPolyRingElem, ) 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. diff --git a/src/connectivity/arbtools.jl b/src/connectivity/arbtools.jl new file mode 100644 index 0000000..d31d312 --- /dev/null +++ b/src/connectivity/arbtools.jl @@ -0,0 +1,76 @@ +""" + arb_to_rat(x::ArbFieldElem) -> Tuple{Rational, Rational} + +Return a rational interval enclosing the Arb ball. +""" +function arb_to_rat(x::ArbFieldElem) + m, r = midpoint(x), 2*radius(x) + for _ in 2:10 + lo = m - r + hi = m + r + + I = map(simplest_rational_inside, [lo, hi]) + I[1] <= x <= I[2] && return I + r *= 10 + end + error("Problem in boxes computations with Arb (arb_to_rat). Try increasing precision.") +end + + +""" + rat_to_arb(interval::Vector{QQFieldElem}, field::ArbField) -> ArbFieldElem + +Convert rational interval (x1, x2) to Arb ball. +""" +function rat_to_arb(interval::Vector{QQFieldElem}, field::ArbField) + x1, x2 = interval + @assert x1 <= x2 "Invalid interval: x1 > x2" + new_field = field + for i in 2:5 + mid = field((x1 + x2) / 2) + rad = field((x2 - x1) / 2) + I = ball(mid, rad) + all(contains(I, x) for x in interval) && return I + new_field = ArbField(i * field.prec) + end + error("Problem in boxes computations with Arb (rat_to_arb). Try increasing precision.") +end + + +""" + evaluate_arb(f, g, x::Vector{QQFieldElem}) + +Evaluate rational function f/g at x. +Return a rational interval enclosing the result, computed using Arb. +""" +function evaluate_arb(f, g, x::Vector{QQFieldElem}, field::ArbField) + @assert !is_zero(g) "Denominator must be non-zero" + + is_zero(f) && return zero(parent(x)) + + cf = field.(coefficients_of_univariate(f)) + cg = field.(coefficients_of_univariate(g)) + + x_arb = rat_to_arb(x, field) + num = evalpoly(x_arb, cf) + den = evalpoly(x_arb, cg) + + return arb_to_rat(num / den) +end + + +""" + arb_eval(f, B::Vector{Tuple}, prec::Int) + +Evaluate polynomial f over interval box B. +Return an Arb ball containing the image of B by f. +""" +function arb_eval(f, B::Vector{Tuple}, field::ArbField) + BA = Vector{ArbFieldElem}(undef, length(B)) + for i in eachindex(B) + BA[i] = rat_to_arb(B[i], field) + end + + f_arb = change_base_ring(field, f) + return evaluate(f_arb, BA) +end \ No newline at end of file diff --git a/src/connectivity/buildpoly.jl b/src/connectivity/buildpoly.jl new file mode 100644 index 0000000..08f8a84 --- /dev/null +++ b/src/connectivity/buildpoly.jl @@ -0,0 +1,145 @@ +""" + change_ringvar(F, new_S::Vector{Symbol}) + +Inject polynomial(s) into a new multivariate polynomial ring defined by `new_S`. +Uses the direct ring constructor `R(coeffs, exps)` to automatically handle +monomial ordering, coalescing, and zero-term removal. +""" +function change_ringvar(F::AbstractVector{<:Union{PolyRingElem, MPolyRingElem}}, new_S::Vector{Symbol}) + isempty(F) && return typeof(first(F))[] + + R = parent(first(F)) + old_S = symbols(R) + newR, _ = polynomial_ring(base_ring(R), new_S) + + # Pre-calculate index mapping: old_index -> new_index (0 if removed) + idx_map = [isnothing(idx) ? 0 : idx for idx in (findfirst(==(v), new_S) for v in old_S)] + + T = typeof(zero(base_ring(R))) + n_vars = length(new_S) + + return map(F) do f + coeffs = T[] + exps = Vector{Int}[] + + if f isa MPolyRingElem + for (e, c) in zip(exponent_vectors(f), coefficients(f)) + e_new = zeros(Int, n_vars) + for (i, power) in enumerate(e) + if power > 0 + @assert idx_map[i] != 0 "Occurrence of removed variable '$(old_S[i])' found!" + e_new[idx_map[i]] = power + end + end + push!(coeffs, c) + push!(exps, e_new) + end + else # PolyRingElem + @assert idx_map[1] != 0 "Occurrence of removed variable '$(old_S[1])' found!" + for (deg_plus_1, c) in enumerate(coefficients(f)) + if !iszero(c) + e_new = zeros(Int, n_vars) + if deg_plus_1 > 1 + e_new[idx_map[1]] = deg_plus_1 - 1 + end + push!(coeffs, c) + push!(exps, e_new) + end + end + end + + # Create the new polynomial associated to f + return newR(coeffs, exps) + end +end + +# Dispatch for a single polynomial +change_ringvar(f::Union{PolyRingElem, MPolyRingElem}, new_S::Vector{Symbol}) = change_ringvar([f], new_S)[1] + + """ + change_ringvar(F) + +Automatically detect all unique variables present in the polynomials and inject them into a ring containing only those variables. +""" +function change_ringvar(F::AbstractVector{<:MPolyRingElem}) + used_vars = unique(reduce(vcat, [Symbol.(vars(f)) for f in F])) + return change_ringvar(F, used_vars) +end + +change_ringvar(f::MPolyRingElem) = change_ringvar([f])[1] + + +""" + MPolyBuild(F::Vector{Vector{RingElem}}, new_S::Vector{Symbol}, idx::Int) + +Construct multivariate polynomials in a single variable. +`F` is a list of coefficient lists in degree-increasing order. +The polynomial will use the variable at index `idx` in `new_S`. +""" +function MPolyBuild(F::AbstractVector{<:AbstractVector{<:RingElement}}, new_S::Vector{Symbol}, idx::Int) + isempty(F) && return [] + + A = parent(first(first(F))) + R, _ = polynomial_ring(A, new_S) + + T = typeof(zero(A)) + n_vars = length(new_S) + + return map(F) do f_coeffs + coeffs = T[] + exps = Vector{Int}[] + + # Construct coeff/exp lists + for (deg_plus_1, c) in enumerate(f_coeffs) + if !iszero(c) + e_new = zeros(Int, n_vars) + e_new[idx] = deg_plus_1 - 1 + push!(coeffs, c) + push!(exps, e_new) + end + end + + # Create the polynomial associated to f_coeffs + return R(coeffs, exps) + end +end + +# Dispatch for a single coefficient vector +MPolyBuild(f::AbstractVector{<:RingElement}, new_S::Vector{Symbol}, idx::Int) = MPolyBuild([f], new_S, idx)[1] + +""" + to_univariate(P::MPolyRingElem, var_idx::Int) + +Converts a bivariate polynomial in `K[x, y]` into a univariate polynomial in `(K[x])[y]`. +""" +function to_univariate(P::MPolyRingElem, var_idx::Int) + R = parent(P) + K = base_ring(R) + coeff_idx = 3 - var_idx # The variable to be pushed to the coefficients + + Rx, x = polynomial_ring(K, symbols(R)[coeff_idx]) + Rxy, y = polynomial_ring(Rx, symbols(R)[var_idx]) + + res = zero(Rxy) + for (c, exp) in zip(coefficients(P), exponent_vectors(P)) + # Build the coefficient in Rx, then attach to y in Rxy + res += (c * x^exp[coeff_idx]) * y^exp[var_idx] + end + return res +end + + +""" + to_bivariate(P::PolyRingElem, R_orig::MPolyRing) + +Maps a univariate polynomial in `(K[x])[y]` back to the original bivariate ring `K[x, y]`. +""" +function to_bivariate(P::PolyRingElem, R_orig::MPolyRing) + ctx = MPolyBuildCtx(R_orig) + for (deg_y, c_y) in enumerate(coefficients(P)) + for (deg_x, c_x) in enumerate(coefficients(c_y)) + push_term!(ctx, c_x, [deg_x - 1, deg_y - 1]) + end + end + return finish(ctx) +end \ No newline at end of file diff --git a/src/connectivity/connectcurves.jl b/src/connectivity/connectcurves.jl new file mode 100644 index 0000000..82ca459 --- /dev/null +++ b/src/connectivity/connectcurves.jl @@ -0,0 +1,349 @@ +include("tools.jl") +include("datastruct.jl") +include("buildpoly.jl") +include("paramcrit.jl") +include("arbtools.jl") +include("isolateboxes.jl") +include("graph.jl") +include("plots.jl") + +@doc Markdown.doc""" + curve_arrangement_graph(curves::Vector{Ideal}; generic=true, outf=true, kwargs...) + curve_arrangement_graph(curves::Vector{Ideal}, C::Vector; kwargs...) + +Computes the combined planar graph of an arrangement of multiple space curves. +Automatically computes the mutual intersections between all curves and guarantees +they are projected using a shared, unified linear form. +""" +function curve_arrangement_graph(curves::Vector{Ideal{P}}; generic=true, outf=true, v=0, kwargs...) where {P <: QQMPolyRingElem} + N = length(curves) + @assert N > 0 "Must provide at least one curve." + + R = parent(curves[1]) + n = nvars(R) + typeout = outf ? Float64 : QQFieldElem + + + # 1. Establish the Shared Projection Context + # We MUST generate the linear forms here and force all curves/intersections to use them. + lfs = nothing + if !generic + make_vec(i) = [j <= n ? ZZRingElem(rand(-100:100)) : (j == n+i ? one(ZZRingElem) : zero(ZZRingElem)) for j in 1:n+3] + lfs = make_vec.(1:3) + u_lfs = lfs[2][1:end-3] # for zero-dim param of intersect pts + end + + # 3. Compute Intersections and Build Individual Graphs + graphs = Vector{CurveGraph{typeout}}(undef, N) + p_inter = Dict{Set{Int}, RationalParametrization}() + for i in 1:N + v > 0 && println("Compute graph of curve number $i/$N...A,") + p_I = lfs === nothing ? rational_curve_parametrization(curves[i]) : + rational_curve_parametrization(curves[i], cfs_lfs=lfs) + new_RS = symbols(parent(p_I.elim)) + + # Control pts are intersections with all other curves; use Dict to avoid re-computation + C_i = Dict{Int, RationalParametrization}( j => p_inter[Set((i,j))] for j in 1:i-1) + for j in i+1:N + I_ij = vcat(curves[i].gens, curves[j].gens) |> Ideal + C_i[j] = isnothing(lfs) ? rational_parametrization(I_ij) : param_use_lfs(I_ij, u_lfs, new_RS[end-1]) + p_inter[Set((i,j))] = + RationalParametrization(C_i[j].vars, C_i[j].cfs_lf, C_i[j].elim, C_i[j].denom, C_i[j].param) + end + @show C_i + @show graphs + graphs[i] = curve_graph(p_I, C_i; outf=outf, v=v-1, kwargs...) + end + + # 4. Merge the graph according to their intersections + return merge_graphs(graphs) +end + +@doc Markdown.doc""" + curve_graph(I::Ideal; generic=true, outf=true, kwargs...) + curve_graph(P::RationalCurveParametrization; kwargs...) + + # Both functions accept optional control points C in various formats: + curve_graph(I, C::Vector{P}; kwargs...) + curve_graph(I, C::Vector{Vector{P}}; kwargs...) + curve_graph(I, C::Dict{Int, Vector{P}}; kwargs...) + +Computes a planar straight-line graph that is homeomorphic to a real algebraic (space) curve. + +### Core Workflow & Pre-processing +1. **Parametrization:** If given an `Ideal`, it first computes a `RationalCurveParametrization`. + If `generic=false`, it applies a random integer shear transformation to place the curve into generic position. +2. **Coefficient Extraction:** Pulls the planar projection `f(x,y) = 0` and the vertical lift `z = g(x,y) / df/dy(x,y)` from the parametrization. +3. **Graph Construction:** Computes bounding boxes for critical points, routes connections, + and identifies "apparent singularities" (2D crossings that do not intersect in 3D). + +### Output Data Structure +Returns a `CurveGraph{T}` object (where `T` is determined by the `outf` flag), containing: +* `Vert::Vector{Tuple{T, T}}`: 2D coordinates of the graph vertices (critical points, routing nodes, control points). +* `Edg::Vector{Tuple{Int, Int}}`: Index pairs defining undirected edges between `Vert` indices. +* `Vcon::Dict{Int, Vector{Int}}`: A mapping from the original index of a control point keys its vertex index in `Vert`. + +### Arguments +* **`I`** (`Ideal`): The algebraic ideal defining the curve. +* **`P`** (`RationalCurveParametrization`): A pre-computed rational parametrization. +* **`C`** (Optional): User-defined plane points on the curve (control points). Given either as Ideal or RationalParametrization. +* **`generic`** (`Bool`, default `true`): If `false`, applies a random shear transformation. +* **`precx`** (`Int`, default `150`): Base numerical precision for real root isolation. +* **`v`** (`Int`, default `0`): Verbosity level. +* **`force_app`** (`Bool`, default `false`): Skips 3D intersection checks, treating all 2D nodes as apparent singularities. +* **`outf`** (`Bool`, default `true`): Output coordinates as `Float64`. If `false`, outputs exact `QQFieldElem`. + +### Example +```jldoctest +julia> using AlgebraicSolving + +julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]); + +julia> f = -42*x^2 + 101*x*y - 8*x*z - 88*x + 53*y^2 + 71*y*z + 154*y + 53*z^2 + 2*z - 32; + g = -7*x^2 - 144*x*y - 126*x*z - 7*x + 30*y^2 - 34*y*z + 72*y + 71*z^2 + 46*z + 5; + +julia> I = Ideal([f, g]); + +julia> G = curve_graph(I); + +julia> number_of_connected_components(G) +3 +``` +""" +function curve_graph(I::Ideal{P}, args...; generic=true, outf=true, v=0, kwargs...) where {P <: QQMPolyRingElem} + R = parent(I) + n = nvars(R) + typeout = outf ? Float64 : QQFieldElem + + # Base Case: 1D Ideal (Points) + if n == 1 + sols = real_solutions(I) + Vert = [(typeout(xi[1]), zero(typeout)) for xi in sols] + # If C is a Dict, instantiate keys. If it's a Vector or absent, return empty Dict + C_keys = length(args) > 0 && args[1] isa Dict ? keys(args[1]) : Int[] + Vcon = Dict{Int, Vector{Int}}(k => Int[] for k in C_keys) + return CurveGraph{typeout}(Vert, Tuple{Int, Int}[], Vcon) + end + + # Pre-processing: Generic Position Shear + # We explicitly define linear forms if !generic to share them with C + v > 0 && println("Compute curve rational parametrization...") + lfs = nothing + if !generic + make_vec(i) = [j <= n ? ZZRingElem(rand(-100:100)) : (j == n+i ? one(ZZRingElem) : zero(ZZRingElem)) for j in 1:n+3] + lfs = make_vec.(1:3) + u_lfs = lfs[2][1:end-3] # for zero-dim param of control pts + p_I = rational_curve_parametrization(I, cfs_lfs = lfs) + else + p_I = rational_curve_parametrization(I) + end + + # 2. Process Control Ideals (C) if provided + if length(args) > 0 + C = args[1] + new_RS = symbols(parent(p_I.elim)) + + # Maps Ideal to the new ring and parameterize it with aligned linear forms + C_param = [ is_nothing(lfs) ? rational_parametrization(c) : param_use_lfs(c, u_lfs, new_RS[end-1]) for c in C ] + + # Map C based on its input structure + if C_input isa AbstractDict + C_param = Dict(k => param_C(v) for (k, v) in C_input) + else + C_input isa AbstractVector || error("Control points C must be a Vector or Dict of Ideals.") + end + + return curve_graph(p_I, C_param; outf=outf, v=v, kwargs...) + end + + # Delegate to the parametrization function, blindly passing any `C` format down + return curve_graph(p_I, args...; outf=outf, v=v, kwargs...) +end + +# ========================================================================= +# MULTIPLE DISPATCH WRAPPERS (INTERNAL) +# ========================================================================= + +# Small helper +function prepare_param(P,Q) + R = parent(Q) + param = is_constant(P.elim) ? [P.elim,P.elim,P.elim] : [P.elim, P.param[end], P.denom] + + return [ evaluate(p, gens(R)[1]) for p in param ] +end + +# Case 1: No control points +curve_graph(p::RationalCurveParametrization; generic=true, kwargs...) = + _compute_graph_core(p.elim, length(p.param) > 0 ? p.param[end] : zero(parent(p.elim)), + Dict{Int, Vector{typeof(p.elim)}}(); kwargs...) + +# Case 2: C is a Vector of Parametrizations +curve_graph(p::RationalCurveParametrization, C::Vector{RationalParametrization}; generic=true, kwargs...) = + _compute_graph_core(p.elim, length(p.param) > 0 ? p.param[end] : zero(parent(p.elim)), + Dict( i => prepare_param(c, p.elim) for (i,c) in enumerate(C)); kwargs...) + +# Case 3: C is a Dictionary of Parametrizations +curve_graph(p::RationalCurveParametrization, C::Dict{Int, RationalParametrization}; generic=true, kwargs...) = + _compute_graph_core(p.elim, length(p.param) > 0 ? p.param[end] : zero(parent(p.elim)), + Dict( i => prepare_param(c, p.elim) for (i,c) in C); kwargs...) + +# ========================================================================= +# CORE IMPLEMENTATION +# ========================================================================= + +function _compute_graph_core(f::P, g::P, C::Dict{Int, Vector{P}}; + precx=150, v=0, force_app=false, outf=true) where {P <: MPolyRingElem} + + @assert !iszero(f) "Input does not define a curve" + + R = parent(f) + x, y = gens(R) + typeout = outf ? Float64 : QQFieldElem + + # Empty set + total_degree(f) == 0 && + return CurveGraph( + Tuple{typeout, typeout}[], + Tuple{Int, Int}[], + Dict{Int, Vector{Int}}() + ) + + # Pre-processing the input + f, g = int_coeffs([f, g]) + precx = max(2, precx) + v > 2 && println("f = $f \n g = $g\n $C") + + # Zero-dim param conditions + d = total_degree(f) + @assert degree(f, 2) == d && degree(f, 1) == d && degree(g, 2) < d && + is_squarefree(f) "Curve not in generic position. Try with generic=false." + + v > 0 && println("Compute parametrization of critical pts...") + @iftime (v > 0) params = param_crit_split(f, g, v=v-1, force_app=force_app) + keys_C = collect(keys(C)) + for (i, k) in enumerate(keys_C) + params[-i] = [ [int_coeffs(C[k][1])], int_coeffs(C[k][2]), int_coeffs(C[k][3]) ] + end + + v > 0 && println("Computing insulating critical boxes") + @iftime (v > 0) LBcrit, Lprecx = insulate_crit_boxes(f, params, precx, v=v-1) + + v > 0 && println("Compute intersections with critical boxes..") + @iftime (v > 0) LPCside, LnPCside = intersect_vertical_boxes(f, params, LBcrit, Lprecx, v=v-1) + + # Critical values and their order + xcrit = Dict(i => [LBcrit[i][j][1] for j in eachindex(LBcrit[i])] for i in keys(LBcrit)) + xcritpermut = order_permut2d(xcrit) + + # Graph data + Vert = Tuple{QQFieldElem, QQFieldElem}[] # List of points (x,y) + Edg = Tuple{Int, Int}[] # List of tuples (idx, idy) + Vcon = Dict{Int, Vector{Int}}(k => Int[] for k in keys(C)) # Index of control vertices + + # Empty or unbounded real curves + if isempty(xcrit) + Vert = reduce(vcat, [ (typeout(e), typeout(yy[1])) for e=0:1 for yy in isolate_eval(f,1,1-2*e)]) + nV = length(Vert) / 2 + Edg = [ (i, i + nV) for i=1:nV ] + return CurveGraph{typeout}(Vert, Edg, Vcon) + end + + Lapp_isolated = Tuple{Int, Int}[] + Lapp_nodes = Tuple{Int, Int}[] + + # Keep track of connection between all computed points when increasing x + Corr = Dict(m => [CriticalNodeConnections() for _ in xcrit[m]] for m in keys(xcrit)) + + for (ind, (i, j)) in enumerate(xcritpermut) + i1, j1 = ind > 1 ? xcritpermut[ind - 1] : (0, 0) # previous + i2, j2 = ind < length(xcritpermut) ? xcritpermut[ind + 1] : (0, 0) # next + + # Current cylinder above xcrit[i][j] + PCside = LPCside[i][j] + I2L_points = ind < length(xcritpermut) ? LPCside[i2][j2].left.points : nothing + + # Midpoints of the critical box itself + xcmid = sum(LBcrit[i][j][1]) / 2 + ycmid = sum(LBcrit[i][j][2]) / 2 + + nI_left = PCside.left.indices_inside + nI_right = PCside.right.indices_inside + ymincrit = minimum(vcat(nI_left, nI_right, [length(PCside.left.points) + 1])) + + # 1. Left side vertices + # one-to-one connection with previous right side + if ind > 1 + for k in 1:length(PCside.left.points) + push!(Corr[i][j].left, Corr[i1][j1].right[k]) + end + else + for k in 1:length(PCside.left.points) + xval = xcrit[i][j][1] + yval = sum(PCside.left.points[k]) / 2 + push!(Vert, (xval, yval)) + push!(Corr[i][j].left, length(Vert)) + end + end + + # 2. Right side vertices + if ind < length(xcritpermut) + for k in 1:length(PCside.right.points) + xval = (xcrit[i][j][2] + xcrit[i2][j2][1]) / 2 + yval = sum(PCside.right.points[k] + I2L_points[k]) / 4 + push!(Vert, (xval, yval)) + push!(Corr[i][j].right, length(Vert)) + end + else + for k in 1:length(PCside.right.points) + xval = xcrit[i][j][2] + yval = sum(PCside.right.points[k]) / 2 + push!(Vert, (xval, yval)) + push!(Corr[i][j].right, length(Vert)) + end + end + + # 3. Below the critical point + for k in 1:ymincrit-1 + yval = sum(PCside.left.points[k] + PCside.right.points[k]) / 4 + push!(Vert, (xcmid, yval)) + push!(Corr[i][j].bottom, length(Vert)) + push!(Edg, (Corr[i][j].left[k], length(Vert))) + push!(Edg, (length(Vert), Corr[i][j].right[k])) + end + + # 4. The critical point + if i == 0 && isempty(nI_left) && isempty(nI_right) + push!(Lapp_isolated, (i, j)) + elseif i == 0 + push!(Edg, (Corr[i][j].left[nI_left[1]], Corr[i][j].right[nI_right[2]])) + push!(Edg, (Corr[i][j].left[nI_left[2]], Corr[i][j].right[nI_right[1]])) + push!(Lapp_nodes, (i, j)) + else + push!(Vert, (xcmid, ycmid)) + push!(Corr[i][j].center, length(Vert)) + + for k in nI_left; push!(Edg, (Corr[i][j].left[k], length(Vert))); end + for k in nI_right; push!(Edg, (length(Vert), Corr[i][j].right[k])); end + + if i < 0 # Control point + push!(Vcon[keys_C[-i]], length(Vert)) + end + end + + # 5. Above the critical point + for k = (length(PCside.left.points) - length(nI_left) - ymincrit + 1):-1:1 + yval = sum(PCside.left.points[end-k+1] + PCside.right.points[end-k+1]) / 4 + push!(Vert, (xcmid, yval)) + push!(Corr[i][j].top, length(Vert)) + push!(Edg, (Corr[i][j].left[end-k+1], length(Vert))) + push!(Edg, (length(Vert), Corr[i][j].right[end-k+1])) + end + end + + if outf + Vert = [ (Float64(x), Float64(y)) for (x,y) in Vert ] + end + + # Return the unified data structure + return CurveGraph{typeout}(Vert, Edg, Vcon) +end \ No newline at end of file diff --git a/src/connectivity/datastruct.jl b/src/connectivity/datastruct.jl new file mode 100644 index 0000000..e0349a0 --- /dev/null +++ b/src/connectivity/datastruct.jl @@ -0,0 +1,22 @@ +# Represents the isolated points and valid connections on a single box edge +struct BoxEdge{T} + points::Vector{Vector{T}} + indices_inside::Vector{Int} +end + +# Represents a complete critical box intersection profile +struct BoxIntersections{T} + bottom::BoxEdge{T} + top::BoxEdge{T} + left::BoxEdge{T} + right::BoxEdge{T} +end + +# Tracks the vertex IDs attached to different parts of a critical point box +Base.@kwdef struct CriticalNodeConnections + left::Vector{Int} = Int[] + bottom::Vector{Int} = Int[] + center::Vector{Int} = Int[] + top::Vector{Int} = Int[] + right::Vector{Int} = Int[] +end \ No newline at end of file diff --git a/src/connectivity/graph.jl b/src/connectivity/graph.jl new file mode 100644 index 0000000..1505beb --- /dev/null +++ b/src/connectivity/graph.jl @@ -0,0 +1,225 @@ +# ========================================================================= +# GRAPH DECOMPOSITION +# ========================================================================= + +@doc Markdown.doc""" + connected_components(G::CurveGraph{T}) where T + +Decomposes a graph into its topologically connected components using a Depth-First Search (DFS). + +This function is particularly useful when a single algebraic curve resolves into multiple +disconnected real branches in space. It guarantees that isolated vertices are preserved as +their own zero-edge components. + +Crucially, this function preserves the `control_nodes` mappings. If a global control node +belonged to a specific vertex, it is correctly remapped to that vertex's new local index +within its specific disconnected subgraph. + +# Inputs +- `G::CurveGraph{T}`: The unified curve graph to be decomposed. + +# Outputs +- `Vector{CurveGraph{T}}`: A list of disjoint subgraphs. Each subgraph retains + the mapped coordinates, edges, and control nodes belonging to its component. + +""" +function connected_components(G::CurveGraph{T}) where T + nv = length(G.vertices) + + # Work with vertices to avoid missing isolated vertices + adj_list = [Int[] for _ in 1:nv] + for (u, v) in G.edges + push!(adj_list[u], v) + push!(adj_list[v], u) + end + + visited = falses(nv) + components = Vector{Vector{Int}}() + + # Iterative DFS + for i in 1:nv + if !visited[i] + comp = Int[] + stack = [i] + visited[i] = true + while !isempty(stack) + curr = pop!(stack) + push!(comp, curr) + for neighbor in adj_list[curr] + if !visited[neighbor] + visited[neighbor] = true + push!(stack, neighbor) + end + end + end + push!(components, comp) + end + end + + # Mapping global vertex ID -> (component ID, local vertex ID) + vertex_map = Vector{Tuple{Int, Int}}(undef, nv) + for (cid, comp) in enumerate(components) + for (local_id, global_id) in enumerate(comp) + vertex_map[global_id] = (cid, local_id) + end + end + + # Group edges + comp_edges = [Tuple{Int, Int}[] for _ in 1:length(components)] + for (u, v) in G.edges + cid_u, local_u = vertex_map[u] + cid_v, local_v = vertex_map[v] + push!(comp_edges[cid_u], (local_u, local_v)) + end + + # Reconstruct subgraphs into CurveGraphs + subgraphs = CurveGraph{T}[] + for (cid, comp) in enumerate(components) + new_verts = G.vertices[comp] + new_edges = comp_edges[cid] + new_control = Dict{Int, Vector{Int}}() + + # Keep track of control points falling into this specific component + for (k, v_list) in G.control_nodes + local_v_list = Int[] + for global_v in v_list + c, local_v = vertex_map[global_v] + if c == cid + push!(local_v_list, local_v) + end + end + if !isempty(local_v_list) + new_control[k] = local_v_list + end + end + + push!(subgraphs, CurveGraph{T}(new_verts, new_edges, new_control)) + end + + return subgraphs +end + +@doc Markdown.doc""" + group_by_component(G::CurveGraph) + +Extracts the control_nodes dictionaries partitioned by the graph's connected components. +This is primarily used to determine which external connections (or user-defined points) +belong to which isolated branch of the curve. + +# Inputs +- `G::CurveGraph`: The parent graph structure. + +# Outputs +- `Vector{Dict{Int, Vector{Int}}}`: A list of dictionaries, one for each connected component, +detailing the local vertex indices mapped to specific control IDs. + +""" +function group_by_component(G::CurveGraph) + subgraphs = connected_components(G) + return [sg.control_nodes for sg in subgraphs if !isempty(sg.vertices)] +end + +@doc Markdown.doc""" + number_of_connected_components(G::CurveGraph) + +Computes the total number of topologically disconnected branches/components in the graph. +This is a convenience wrapper for topological validation. + +# Inputs +- `G::CurveGraph`: The parent graph structure. + +# Outputs +- `Int`: The count of isolated subgraphs. + +""" +number_of_connected_components(G::CurveGraph) = length(connected_components(G)) + + +# ========================================================================= +# GRAPH MERGING +# ========================================================================= + +""" + merge_graphs(graphs::Vector{CurveGraph{T}}) where T + +Merges a collection of disjoint graphs into a single graph, mathematically fusing +vertices based on their control_nodes dictionnaries that describe a one-to-one +correspondance between subset of vertices of each graph. + + +# Fusion Mechanism: +This work as follows. Given V1 = [u1,u2] and V2 = [w1,w2,w3] the vertices of two graphs G1 and G2. +Let C1= {2 => [1]} and C2 = {1 => [2]} be their respective control_nodes disctionnary. +Then it merges the first vertex w1 of G2 with the second vertex u2 of G1. +The merged graph has vertices [u1,u2,w2,w3]. + +Note that control nodes that refer to external bounds (e.g., indices k <= 0 or k > length(graphs)) +are preserved in the final merged graph to allow for subsequent higher-level merging or tracking. + +# Inputs +- `graphs::Vector{CurveGraph{T}}`: A list of graphs. `graphs[i].control_nodes[k]` + should contain the local vertex indices in graph `i` that connect to graph `k`. + This vertex correspondance must be one-to-one. + +# Outputs +- `CurveGraph{T}`: A unified graph struct. +""" +function merge_graphs(graphs::Vector{CurveGraph{T}}) where T + isempty(graphs) && return CurveGraph{T}(Tuple{T,T}[], Tuple{Int,Int}[], Dict{Int, Vector{Int}}()) + length(graphs) == 1 && return deepcopy(graphs[1]) + + Vtot = copy(graphs[1].vertices) + Etot = copy(graphs[1].edges) + + # Dictionary to track: (graph_id, local_vertex_id) -> merged_global_id + global_map = Dict{Tuple{Int, Int}, Int}() + for i in 1:length(graphs[1].vertices) + global_map[(1, i)] = i + end + + for i in 2:length(graphs) + G_i = graphs[i] + + # Step 1: Pre-map common control nodes connected to previously merged graphs + for (k, common_indices_in_i) in G_i.control_nodes + if k < i && !isempty(common_indices_in_i) + @assert haskey(graphs[k].control_nodes, i) "Correspondance is not one-to-one" + common_indices_in_k = graphs[k].control_nodes[i] + @assert length(common_indices_in_i) == length(common_indices_in_k) "Correspondance is not one-to-one" + # Zip ensures we pair corresponding control nodes precisely + for (idx_i, idx_k) in zip(common_indices_in_i, common_indices_in_k) + global_map[(i, idx_i)] = global_map[(k, idx_k)] + end + end + end + + # Step 2: Add completely new (unmapped) vertices from G_i to Vtot + for idx_i in 1:length(G_i.vertices) + if !haskey(global_map, (i, idx_i)) + push!(Vtot, G_i.vertices[idx_i]) + global_map[(i, idx_i)] = length(Vtot) + end + end + + # Step 3: Append the edges using the newly mapped global indices + for (u, v) in G_i.edges + push!(Etot, (global_map[(i, u)], global_map[(i, v)])) + end + end + + # Step 4: Preserve external control nodes pointing OUTSIDE the merged structure + merged_control = Dict{Int, Vector{Int}}() + num_graphs = length(graphs) + + for i in 1:num_graphs + for (k, v_list) in graphs[i].control_nodes + # If `k` points to a graph outside of the ones we just merged + if k > num_graphs || k <= 0 + mapped_v = [global_map[(i, v)] for v in v_list] + append!(get!(merged_control, k, Int[]), mapped_v) + end + end + end + + return CurveGraph{T}(Vtot, Etot, merged_control) +end \ No newline at end of file diff --git a/src/connectivity/isolateboxes.jl b/src/connectivity/isolateboxes.jl new file mode 100644 index 0000000..8e7edbd --- /dev/null +++ b/src/connectivity/isolateboxes.jl @@ -0,0 +1,322 @@ +in_inter(I, J) = J[1] <= I[1] && I[2] <= J[2] +overlap_inter(I, J) = max(I[1], J[1]) <= min(I[2], J[2]) + +""" + isolate(f; prec = 32) + +Isolate real roots of a multivariate polynomial with a single variable +""" +function isolate(f; prec = 32) + total_degree(f) == 0 && return [] + @assert is_univariate(f) "Not univariate polynomial" + + sols = real_solutions(Ideal([change_ringvar(f)]), precision=prec, interval=true, info_level=0) + return getindex.(sols, 1) +end + +""" + isolate_eval(f, ivar, val; prec=64) + +Isolate real roots of a bivariate polynomial whose ivar-th variable is +evaluated at val +""" +function isolate_eval(f, ivar, val; prec=64) + fev = evaluate(f, [ivar], [val]) + fev = fev / content(fev) + return isolate(fev, prec=prec) +end + +""" + _expand_degenerate_intervals!(xnew, prec) + +Expand intervals where lower == upper to avoid zero-width boxes. +""" +function _expand_degenerate_intervals!(xnew, prec) + d = 1 // (ZZ(1) << prec) + for i in eachindex(xnew) + if xnew[i][1] == xnew[i][2] + x = xnew[i][1] + xnew[i] = [x - d, x + d] + end + end + return xnew +end + +""" + _compute_boxes_for_index!(i, params, LBcrit, precx) + +Compute isolating boxes for a single index i. +""" +function _compute_boxes_for_index!(i, params, LBcrit, precx) + xvals = reduce(vcat, (isolate(pp, prec=precx) for pp in params[i][1])) + _expand_degenerate_intervals!(xvals, precx) + + arbField = ArbField(precx) + yvals = [ + evaluate_arb(params[i][2], params[i][3], xc, arbField) + for xc in xvals + ] + _expand_degenerate_intervals!(yvals, precx) + + @inbounds LBcrit[i] = [[xvals[j], yvals[j]] for j in eachindex(xvals)] +end + + +""" + _needs_refinement(i, f, LBcrit, precx) + +Check if boxes at index i require refinement. +""" +function _needs_refinement(i, Lfyk, LBcrit, precx) + arbField = ArbField(precx) + m = i <= 1 ? 2 : i + for box in LBcrit[i] + pcrit = [rat_to_arb(c, arbField) for c in box] + + if contains_zero(evaluate(Lfyk[m+1], pcrit)) + return true + end + end + return false +end + +""" + _diff_list(p, v, n) + +Compute the list of the n-th first derivative of p w.r.t v +""" +function _diff_list(p, v, n) + L = Vector{typeof(p)}(undef, n + 1) + L[1] = p + @inbounds for j in 2:n+1 + L[j] = derivative(L[j-1], v) + end + return L +end + +""" + compute_crit_and_singular_boxes(f, params, precx; max_attempts=5, v=0) + +Compute insulating boxes for critical points, refining precision if +necessary. Insulating means that each box contains a single critical +point and all branches of the curve intersecting the box are incident +to this critical point. +""" +function insulate_crit_boxes(f, params, precx; max_attempts=5, v=0) + LBcrit, Lprecx = Dict(), Dict() + Lfyk = _diff_list(f, 2, maximum(keys(params), init=2)) + + for i in keys(params) + attempts, precxi = 0, precx + while attempts <= max_attempts + try + attempts > 0 && v > 0 && + println("Refining index $i -> prec = $precxi") + + _compute_boxes_for_index!(i, params, LBcrit, precxi) + + # singularity check + if !_needs_refinement(i, Lfyk, LBcrit, precxi) + break + end + + # refine ONLY this index + precxi *= 2 + attempts += 1 + + catch e + precxi *= 2 + attempts += 1 + v > 1 && println("Error at index $i: ", e) + end + end + attempts > max_attempts && + error("Failed to isolate boxes for index $i") + + Lprecx[i] = precxi + end + + return LBcrit, Lprecx +end + +""" + _isolate_box_edge(f, fixed_dim, fixed_val, target_interval, initial_prec, v; max_retries=5) + +Function that computes the intersection of a curve with a vertical +(fixed_dim=1)/horizontal (fixed_dim=2) line, and identify those in +target interval. +Returns isolating intervals and indices of the interval in target_interval +""" +function _isolate_box_edge(f, fixed_dim, fixed_val, target_interval, initial_prec, v; max_retries=5) + prec = initial_prec + + for _ in 1:max_retries + roots = isolate_eval(f, fixed_dim, fixed_val, prec=prec) + indices_inside = Int[] + overlap_found = false + + for (j, l) in pairs(roots) + if in_inter(l, target_interval) + push!(indices_inside, j) + elseif overlap_inter(l, target_interval) + prec *= 2 + v > 0 && println("IntersectBox: increase precision to $prec bits") + overlap_found = true + break + end + end + # All roots without overlap on this edge + if !overlap_found + T = QQFieldElem + return BoxEdge{T}(roots, indices_inside) + end + end + error("Problem when isolating on one side of a box after $max_retries attempts\n + Consider increasing the initial precision or genericity position.") +end + +""" + intersect_box(f, B; prec=100, v=0) + +Computes the intersection of a curve with a box by +isolating roots on the edges. +""" +function intersect_box(f, B; prec=100, v=0) + # Evaluate horizontal edges (fix y to B[2][1] and B[2][2], target x-interval B[1]) + edge_y1 = _isolate_box_edge(f, 2, B[2][1], B[1], prec, v) + edge_y2 = _isolate_box_edge(f, 2, B[2][2], B[1], prec, v) + # Evaluate vertical edges (fix x to B[1][1] and B[1][2], target y-interval B[2]) + edge_x1 = _isolate_box_edge(f, 1, B[1][1], B[2], prec, v) + edge_x2 = _isolate_box_edge(f, 1, B[1][2], B[2], prec, v) + + # Return clean struct + T = QQFieldElem + return BoxIntersections{T}(edge_y1, edge_y2, edge_x1, edge_x2) +end + +""" + _update_LB_first_axis!(LB, xnew) + +Replace first coordinate of each box in LB with corresponding xnew interval. +""" +function _update_LB_first_axis!(LB, xnew) + @inbounds for i in eachindex(LB) + LB[i][1] = xnew[i] + end + return LB +end + +""" + refine_xboxes(f, LB, prec) + +Refine LB along first axis using roots of polynomials in F +The order in F is assumed to match the one in LB +""" +function refine_xboxes(F::Vector{T} where T<:Union{PolyRingElem, MPolyRingElem}, LB, prec) + # Concatenate all roots from isolating each polynomial in F + xnew = reduce(vcat, [isolate(f, prec=prec) for f in F]) + + _expand_degenerate_intervals!(xnew, prec) + return _update_LB_first_axis!(LB, xnew) +end + +""" + intersect_vertical_boxes(f, params, LBcrit, Lprecx; max_attempts=5, v=0) + +Computes the intersection of the vertical sides of the critical boxes +in LBcrit with the curve defined by f. + +If an intersection on top/bottom side is detected, refine the width +using params and Lprecx until this does not happen or max_attemps is +reached. +""" +function intersect_vertical_boxes(f, params, LBcrit, Lprecx; max_attempts=5, v=0) + T = QQFieldElem + LPCside = Dict{Int, Vector{BoxIntersections{T}}}() + LnPCside = Dict{Int, Vector{NTuple{4, Int}}}() + + for (i, boxes) in pairs(LBcrit) + prec = Lprecx[i] + n = length(boxes) + sides_vec = Vector{BoxIntersections{T}}(undef, n) + nsides_vec = Vector{NTuple{4, Int}}(undef,n) + + for attempt in 1:max_attempts + refined = false + + @inbounds for j in 1:n + v >= 0 && print("i=$i ($j/$n)\r") + pc = intersect_box(f, boxes[j], prec=prec, v=v) + + c1 = length(pc.bottom.indices_inside) + c2 = length(pc.top.indices_inside) + c3 = length(pc.left.indices_inside) + c4 = length(pc.right.indices_inside) + + if (i == 1 && (c1 + c2 + c3 + c4 > 2)) || + (i != 1 && (c1 != 0 || c2 != 0)) + prec *= 2 + v > 0 && println("\nRefine -> prec = $prec") + refine_xboxes(params[i][1], boxes, prec) + refined = true + break + end + + sides_vec[j] = pc + nsides_vec[j] = (c1,c2,c3,c4) + end + + refined || break + attempt == max_attempts && + error("Problem in computing intersections (i=$i)") + end + v > 0 && print("\n") + LPCside[i] = sides_vec + LnPCside[i] = nsides_vec + Lprecx[i] = prec + end + + # ---- Extremal case (i == 1) we vertically enlarge boxes manually ----- + if haskey(LBcrit, 1) + boxes, pcs, npcs = LBcrit[1], LPCside[1], LnPCside[1] + + @inbounds for j in eachindex(boxes) + c1, c2, c3, c4 = npcs[j] + + (c1 == 0 && c2 == 0) && continue # Nothing to do + + pc = pcs[j] + ycrit = boxes[j][2] + + # Orientation of the curve left(s=1)/right(s=2) + s = length(pc.left.points) >= length(pc.right.points) ? 1 : 2 + side_pts = s == 1 ? pc.left.points : pc.right.points + side_indices = s == 1 ? pc.left.indices_inside : pc.right.indices_inside + + # ---- intersection on bottom ---- + if c1 == 1 + # Attach to the max intersection below the box on the correct side + ind_yinf = maximum(i for (i, yy) in pairs(side_pts) if yy[1] < ycrit[1]) + push!(side_indices, ind_yinf) + empty!(pc.bottom.indices_inside) + npcs[j] = (0, c2, s == 1 ? c3 + 1 : c3, s == 2 ? c4 + 1 : c4) + # Update also LBcrit for consistency purpose only + boxes[j][2][1] = side_pts[ind_yinf][1] + end + + c1, c2, c3, c4 = npcs[j] # In case this has changed + # ---- intersection on top ---- + if c2 == 1 + # Attach to the min intersection above the bow on the correct side + ind_ymax = minimum(i for (i, yy) in pairs(side_pts) if yy[2] > ycrit[2]) + push!(side_indices, ind_ymax) + empty!(pc.top.indices_inside) + npcs[j] = (c1, 0, s == 1 ? c3 + 1 : c3, s == 2 ? c4 + 1 : c4) + # Update also LBcrit for consistency purpose only + boxes[j][2][2] = side_pts[ind_ymax][2] + end + end + end + + return LPCside, LnPCside +end \ No newline at end of file diff --git a/src/connectivity/paramcrit.jl b/src/connectivity/paramcrit.jl new file mode 100644 index 0000000..043fda3 --- /dev/null +++ b/src/connectivity/paramcrit.jl @@ -0,0 +1,292 @@ +# ========================================================================= +# SUBRESULTANTS +# ========================================================================= + +# Univariate subresultants +function subresultants(P::Union{PolyRingElem{T}, FqPolyRingElem}, Q::Union{PolyRingElem{T}, FqPolyRingElem}, is_Fq::Bool=false) where T <: RingElem + degree(P) < degree(Q) && ((P, Q) = (Q, P)) + + S = typeof(P)[Q] + s = leading_coefficient(Q)^(degree(P) - degree(Q)) # Initial scaling factor + # rem is faster when working with finite fields (relies on FLINT) + A = Q + B = is_Fq ? (-leading_coefficient(Q))^(degree(P) - degree(Q) + 1) * rem(P, Q) : pseudorem(P, -Q) + + while !iszero(B) + d, e = degree(A), degree(B) + push!(S, B) # current remainder + delta = d - e # degree drop + + if delta > 1 + if length(S) > 2 + last_add, prev_add = S[end], S[end-1] + n = degree(prev_add) - degree(last_add) - 1 # stabilize scaling + if n == 0 + C = last_add + else + x, y = leading_coefficient(last_add), leading_coefficient(prev_add) + a = 1 << (length(bits(ZZ(n))) - 1) + c, n = x, n - a + while a > 1 + a >>= 1 + c = divexact(c^2, y) + if n >= a + c = c * divexact(x, y) + n -= a + end + end + C = c * divexact(last_add, y) + end + else + # First step: fallback normalization + C = divexact(leading_coefficient(B)^(delta - 1) * B, s^(delta - 1)) + end + push!(S, C) + else + C = B # no degree drop + end + + # Termination: constant polynomial + e == 0 && return reverse!(S) + + # --- Next pseudo-remainder --- + B = is_Fq ? (-leading_coefficient(B))^(d - e + 1) * rem(A, B) : pseudorem(A, -B) + B = divexact(B, s^delta * leading_coefficient(A)) + + A = C + s = leading_coefficient(A) + end + return reverse!(S) +end + +# Bivariate subresultants seen in (K[x])[y] +function subresultants(P::MPolyRingElem, Q::MPolyRingElem, var_idx::Int; list=false) + UP, UQ = to_univariate(P, var_idx), to_univariate(Q, var_idx) + sr_uni = subresultants(UP, UQ) + + R_orig = parent(P) + v_orig = symbols(R_orig) + + if list + # Extract coefficients as bivariate polynomials + return [change_ringvar(collect(coefficients(sr)), v_orig) for sr in sr_uni] + else + return [to_bivariate(sr, R_orig) for sr in sr_uni] + end +end + +# ========================================================================= +# MODULAR ARITHMETIC & INTERSECTION +# ========================================================================= + +""" + num_biv_rat_mod(A::MPolyRingElem, P::Vector) + +Computes the numerator of `A(x, a(x)/b(x)) mod q(x)` efficiently using Residue Rings. +`P = [q, a, b]`. +""" +function num_biv_rat_mod(A::MPolyRingElem, P::Vector) + q, a, b = P[1], P[2], P[3] + Rx = parent(q) + + # Work modulo q and compute y = a/b in the quotient ring + U, mod_q = residue_ring(Rx, q) + y_u = mod_q(a) / mod_q(b) + + # Build coefficients of A as polynomials in y with coefficients in U + degA_y = maximum(e[2] for e in exponent_vectors(A)) + coeff_A_u = [ MPolyBuildCtx(Rx) for _ in 0:degA_y ] + for (c, exp) in zip(coefficients(A), exponent_vectors(A)) + push_term!(coeff_A_u[exp[2] + 1], c, [exp[1], 0]) + end + coeff_A_u = [ mod_q(finish(a)) for a in coeff_A_u ] + + # Evaluate the univariate polynomial at y_u in the quotient ring + U_y,_ = polynomial_ring(U) + res_u = evaluate(U_y(coeff_A_u), y_u) + + return lift(res_u) +end + +""" + intersect_biv(P::Vector, A::MPolyRingElem) + +P = [q, a, b] encodes a finite set (x,a(x)/b(x)) where q(x)=0 +Computes the divisor of `q` that encodes the intersection with `A(x,y)=0`. +Uses a modular CRT loop. +""" +function intersect_biv(P::Vector, A::MPolyRingElem) + iszero(A) && return P[1] + + RSA = symbols(parent(A)) + dA_Q_prev, dA_Q_final = QQFieldElem[], QQFieldElem[] + dA_Z_curr = ZZRingElem[] + pprod, p = ZZ(1), ZZ(1) << 60 + compt = 0 + + while compt < 12 + p = next_prime(p) + Fp = GF(p) + # Prime check + lcA, lcP = ZZ(leading_coefficient(A)), ZZ(leading_coefficient(P[1])) + iszero(rem(gcd(lcA, lcP), p)) && continue + + # Map polynomials natively to Fp + Pp = [map_coefficients(Fp, poly) for poly in P] + Ap = map_coefficients(Fp, A) + + Apev = num_biv_rat_mod(Ap, Pp) + dAp = gcd(Pp[1], Apev) + + # Lift coefficients back to integers for CRT + dA_Z_next = [lift(ZZ, c) for c in coefficients_of_univariate(dAp)] + + if compt > 0 + dA_Z_next = [crt([d1, d2], [pprod, p]) for (d1, d2) in zip(dA_Z_curr, dA_Z_next)] + end + + # Trivial GCD : reduced gcd has larger degree + # Then, A(x,y) does not vanish on the points defined by q(x)=0 + all(iszero(c) for c in dA_Z_next[2:end]) && return one(parent(A)) + + pprod *= p + try + dA_Q_final = [reconstruct(c, pprod) for c in dA_Z_next] + if (compt > 0 && dA_Q_final == dA_Q_prev) + fact = MPolyBuild(dA_Q_final, RSA, 1) + divides(P[1], fact)[1] && return fact + end + catch + end + + dA_Z_curr = dA_Z_next + dA_Q_prev = dA_Q_final + compt += 1 + end + error("Failed multi-modular computation of apparent sing: $compt primes used") +end + +# ========================================================================= + +@doc Markdown.doc""" + fact_gcd(delta::T, LP::Vector{T}) where T <: Union{PolyRingElem, MPolyRingElem} + +Decomposes the polynomial `delta` by separating it into factors based on its greatest common divisors with the sequence of polynomials in `LP`. +The $i$-th entry corresponds to the exact factor of `delta` that divides $LP[1], \dots, LP[i-1]$, but is coprime to $LP[i]$. Factors of degree $0$ are discarded. +""" +function fact_gcd(delta::T, LP::Vector{T}) where {T <: PolyRingElem} + Ldelta = Dict{Int, T}() + curr_phi = delta + + for (i, p) in pairs(LP) + degree(curr_phi) == 0 && break + + next_phi = gcd(curr_phi, p) + q = divexact(curr_phi, next_phi) + + if degree(q) > 0 + Ldelta[i] = q + end + + curr_phi = next_phi + end + + return Ldelta +end + +# wrapper for the above function +function fact_gcd(delta::T, LP::Vector{T}) where (T <:MPolyRingElem) + @assert is_univariate(delta) && all(is_univariate.(LP)) "Not univariate polynomial" + R = parent(delta) + K, RS = coefficient_ring(R), symbols(R) + A, = polynomial_ring(K, RS[1]) + + # Convert to standard univariate ring + d_uni = A(coefficients_of_univariate(delta)) + LP_uni = [A(coefficients_of_univariate(p)) for p in LP] + + out = fact_gcd(d_uni, LP_uni) + return Dict(i => change_ringvar(f, [first(RS)]) for (i,f) in out) +end + + +# ========================================================================= + +""" + param_crit_split(f::MPolyRingElem, g::MPolyRingElem; v=0, force_app=false) + +Computes the critical points (grouped by multiplicity) and apparent singularities +of the space curve defined by: + `f(x, y) = 0` and `(df/dy)(x,y) * z = g(x,y)` + +This relies on subresultant computations and apparent singularity criterion from: +A.Poteaux, N.Islam, R.Prébet - Algorithm for Connectivity Queries on Real Algebraic Curves - ISSAC'23 + +If force_app, it assumes that all nodes are apparent singularities. +""" +function param_crit_split(f::MPolyRingElem, g::MPolyRingElem; v=0, force_app=false) + v > 0 && println("Compute subresultant sequence") + f_y = derivative(f, 2) + @iftime v>0 sr = subresultants(f, f_y, 2, list=true) + (isempty(sr) || total_degree(sr[1][1]) == 0) && return Dict{Int, Vector}() + + v > 1 && println("Factorization") + sqr_factors = factor(sr[1][1]) + + # Filter non-zero multiplicity and sort + sqr = sort([(fac, mult) for (fac, mult) in sqr_factors if mult > 0], by = x -> x[2]) + isempty(sqr) && return Dict{Int, Vector}() + + sqrmult = unique([s[2] for s in sqr]) + group_sqr = Dict(m => [s[1] for s in sqr if s[2] == m] for m in sqrmult) + + param_crit = Dict{Int, Vector}() + + # Helper function to prevent out-of-bounds access on subresultant lists + get_sr = (p, idx)-> (length(sr) >= p && length(sr[p]) >= abs(idx) + 1) ? sr[p][end+idx] : zero(parent(f)) + + # Critical points : multiplicity 1 in res + (1 in sqrmult) && (param_crit[1] = [group_sqr[1], -get_sr(2, -1), get_sr(2, 0)]) + + # Singular points + singmult = filter(p -> p*(p-1) <= sqrmult[end], 2:sqrmult[end]) + isempty(singmult) && return filter(p -> length(p.second[1]) > 0, param_crit) + + @assert length(sr) >= singmult[end] "Curve not in generic position. Try with generic=false." + for p in singmult + @assert length(sr[p]) >= 2 "Curve not in generic position. Try with generic=false." + param_crit[p] = [MPolyRingElem[], -get_sr(p, -1), (p-1)*get_sr(p, 0)] + end + + # Nodes : multiplicity 2 in res + if 2 in sqrmult + if !(force_app) + v > 0 && println("Compute apparent singularities") + f_x = derivative(f, 1) + A = derivative(f_y, 2) * derivative(g, 1) - derivative(f_x, 2) * derivative(g, 2) + + dA = [intersect_biv([q, -get_sr(2, -1), get_sr(2, 0)], A) for q in group_sqr[2]] + + param_crit[0] = [group_sqr[2] ./ dA, -get_sr(2, -1), get_sr(2, 0)] + append!(param_crit[2][1], dA) + else + param_crit[0] = [group_sqr[2], -get_sr(2, -1), get_sr(2, 0)] + end + end + + # Other singularities + filter!(m -> !(m in [1, 2]), sqrmult) + lsr = [get_sr(p, 0) for p in 2:(singmult[end]+1)] + + v > 0 && println("Compute gcd decomposition") + for m in sqrmult + for q in group_sqr[m] + for (i, dji) in fact_gcd(q, lsr) + @assert haskey(param_crit, i+1) "Curve not in generic position. Try with generic=false." + push!(param_crit[i+1][1], dji) + end + end + end + + return filter(p -> length(p.second[1]) > 0, param_crit) +end \ No newline at end of file diff --git a/src/connectivity/plots.jl b/src/connectivity/plots.jl new file mode 100644 index 0000000..cec0ce7 --- /dev/null +++ b/src/connectivity/plots.jl @@ -0,0 +1,133 @@ +""" + distinguishable_colors(n; s=0.8, v=0.9) + +Return `n` visually distinct HEX colors using manual HSV->RGB conversion. +No external libraries used. +""" +function distinguishable_colors(n::Int; s::Float64=0.6, v::Float64=0.9) + n <= 0 && return String[] + + colors = Vector{String}(undef, n) + + for i in 0:n-1 + h = i / n # [0,1) + + c = v * s + h6 = h * 6 + x = c * (1 - abs(mod(h6, 2) - 1)) + m = v - c + + r1 = g1 = b1 = 0.0 + + if 0 ≤ h6 < 1 + r1, g1, b1 = c, x, 0 + elseif 1 ≤ h6 < 2 + r1, g1, b1 = x, c, 0 + elseif 2 ≤ h6 < 3 + r1, g1, b1 = 0, c, x + elseif 3 ≤ h6 < 4 + r1, g1, b1 = 0, x, c + elseif 4 ≤ h6 < 5 + r1, g1, b1 = x, 0, c + else + r1, g1, b1 = c, 0, x + end + + r = Int(round((r1 + m) * 255)) + g = Int(round((g1 + m) * 255)) + b = Int(round((b1 + m) * 255)) + + colors[i+1] = @sprintf("#%02X%02X%02X", r, g, b) + end + + return colors +end + + +# ---------- CORE BUILDERS ---------- + +@doc Markdown.doc""" + build_graph_data(G::CurveGraph{T}; width=3.0, vert=true, color="#FF0000") where T <: Union{Float64,QQFieldElem} + +Construct plotting data from a `CurveGraph`. + +This function converts the piecewise linear structure encoded in `G` into a +`GraphPlotData` object, separating edges and point groups (vertices and control points), +ready for visualization with plotting libraries such as `Plots.jl` (see documentation +of `GraphPlotData`) + +# Arguments +- `G::CurveGraph{T}`: Input graph containing vertices, control nodes, and edges. +- `width::Real=3.0`: Line width used for rendering edges. +- `vert::Bool=true`: If `true`, include vertices in the graph. +- `color::AbstractString="#FF0000"`: Color used for both edges and points. + +# Behavior +- Edges are converted into line segments between vertex coordinates. +- Vertices (if enabled) are displayed with marker `:x`. +- Control nodes are grouped and displayed with marker `:+`. + +# Returns +- `GraphPlotData`: A structure containing: + - one `EdgeGroup` for all edges + - multiple `PointGroup`s for vertices and control nodes + +# Notes +- All coordinates are converted to `Float64` for compatibility with plotting backends. +""" +function build_graph_data(G::CurveGraph{T}; width=3.0, vert=true, color="#FF0000") where T <: Union{Float64,QQFieldElem} + V, Vcon, E = G.vertices, G.control_nodes, G.edges + + # edges + edges = [ + begin + v1, v2 = [map(Float64, V[idx]) for idx in e] + ([v1[1], v2[1]], [v1[2], v2[2]]) + end + for e in E + ] + edge_group = EdgeGroup(edges, color, width) + + #points + point_groups = PointGroup[] + + if vert + vx = map(v -> Float64(v[1]), V) + vy = map(v -> Float64(v[2]), V) + push!(point_groups, PointGroup((vx, vy), color, :x)) + end + + for group in Vcon + hx = map(j -> Float64(V[j][1]), group) + hy = map(j -> Float64(V[j][2]), group) + push!(point_groups, PointGroup((hx, hy), color, :+)) + end + + return GraphPlotData(edge_group, point_groups) +end + +""" + build_graph_data(CG::Vector{CurveGraph{T}}; width=3.0, vert=true) where T <: Union{Float64,QQFieldElem} + +Construct plotting data for multiple `CurveGraph`s with automatically assigned colors. + +# Arguments +- `CG::Vector{CurveGraph{T}}`: Collection of curve graphs to visualize. +- `width::Real=3.0`: Line width used for rendering edges. +- `vert::Bool=true`: If `true`, include vertices in each graph. + +# Behavior +- Colors are generated using `distinguishable_colors` to ensure visual contrast. +- Each graph is assigned a unique color consistently across its edges and points. + +# Returns +- `Vector{GraphPlotData}`: One plotting data object per input graph. + +# Notes +- The number of colors scales with `length(CG)`. +""" +function build_graph_data(CG::Vector{CurveGraph{T}}; width=3.0, vert=true) where T <: Union{Float64,QQFieldElem} + c = distinguishable_colors(length(CG)) + return [build_graph_data(G; width=width, vert=vert, color=c[i]) + for (i, G) in enumerate(CG)] +end \ No newline at end of file diff --git a/src/connectivity/polar.jl b/src/connectivity/polar.jl new file mode 100644 index 0000000..ad0805c --- /dev/null +++ b/src/connectivity/polar.jl @@ -0,0 +1,112 @@ +@doc Markdown.doc""" + computepolar(J::Union{Vector{Int},UnitRange{Int}}, V::Ideal{P} where P <: MPolyRingElem; ) + +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 \ No newline at end of file diff --git a/src/connectivity/roadmap.jl b/src/connectivity/roadmap.jl new file mode 100644 index 0000000..ae3da33 --- /dev/null +++ b/src/connectivity/roadmap.jl @@ -0,0 +1,226 @@ +@doc Markdown.doc""" + roadmap(I::Ideal{T} where T <: MPolyRingElem, ) + +Given a **radical** ideal `I` with solution set X, that is smooth and +whose real trace XR is bounded, return a roadmap of XR + +The output is given as a Roadmap structure, encoding the recursive structure +of roadmaps. It is encoded as a chained list, whose root is containing the equations defining X +and each node representing a curve component, that is defined by additional polar equation and base point< +Moreover it is linked to fibers, that share the same base point. + +# Arguments +- `I::Ideal{T} where T <: QQMPolyRingElem`: input generators. +- `C::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[]`: query points with rational coefficients +- `info_level::Int=0`: verbosity level +- `checks::Bool=false`: whether perform checks (dimension, regularity, etc.) +- 'generic_change=false": whether it performs a prior random linear change of variables (TODO) +) + +# Examples +```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> RM = roadmap(I, checks=true) +Vector{QQFieldElem}[[], [-3], [-3, -2], [-2], [-2, -1], [-2, 0], [-2, 1], [3], [3, 2]] + +julia> nb_nodes(RM) +9 + +julia> all_eqs(RM) +9-element Vector{Ideal{QQMPolyRingElem}}: + 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^2*x3 - 4*x2^2*x3 - 4*x3^3 - 4*x3*x4^2 + 40*x3, -4*x1^2*x4 - 4*x2^2*x4 - 4*x3^2*x4 - 4*x4^3 - 32*x4] + 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^2*x4 - 4*x2^2*x4 - 4*x3^2*x4 - 4*x4^3 - 32*x4, x1 + 3] + 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, x1 + 3, x2 + 2] + 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^2*x4 - 4*x2^2*x4 - 4*x3^2*x4 - 4*x4^3 - 32*x4, x1 + 2] + 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, x1 + 2, x2 + 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, x1 + 2, x2] + 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, x1 + 2, x2 - 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, -4*x1^2*x4 - 4*x2^2*x4 - 4*x3^2*x4 - 4*x4^3 - 32*x4, x1 - 3] + 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, x1 - 3, x2 - 2] +``` +""" +function roadmap( + I::Ideal{P}; # input ideal + C::Vector{Vector{Vector{QQFieldElem}}}=Vector{Vector{QQFieldElem}}[], # query points: interval with rational coefficients + info_level::Int=0, # verbosity level + checks::Bool=false # perform checks (dimension, regularity, etc.) +) where (P <: QQMPolyRingElem) + return _roadmap_rec(I, QQFieldElem[], C, info_level, checks) +end + +@doc Markdown.doc""" + roadmap(I::Ideal{T} where T <: MPolyRingElem, I::Ideal{P}, C::Ideal{P}; info_level::Int=0, checks::Bool=false) +``` +""" +function roadmap( + I::Ideal{P}, # input ideal + C::Ideal{P}; # ideal defining query points + info_level::Int=0, # verbosity level + checks::Bool=false # perform checks (dimension, regularity, etc.) +) where (P <: QQMPolyRingElem) + @assert(parent(I)==parent(C), "Equations for variety and query points must live the same ring") + CI = real_solutions(C, info_level=max(info_level-1,0), nr_thrds=Threads.nthreads()) + return _roadmap_rec(I, QFieldElem[], CI, info_level, checks) +end + +function _roadmap_rec( + I::Ideal{T} where T <: QQMPolyRingElem, # input ideal + q::Vector{QQFieldElem}, # single base point with rational coefficients + C::Vector{Vector{Vector{QQFieldElem}}}, # query points with rational coefficients + info_level::Int, # verbosity level + checks::Bool # perform checks (dimension, regularity, etc.) +) + # Some base cases + if nvars(parent(I))<=2 + return [I.gens] + end + + # Some preprocessng + if isnothing(I.dim) + lucky_prime = _generate_lucky_primes(I.gens, one(ZZ)<<30, one(ZZ)<<31-1, 1) |> first + local INEW = Ideal(change_base_ring.(Ref(GF(lucky_prime)), I.gens)) + dimension(INEW) + I.dim = INEW.dim + end + e = length(q) + + ## Fq ## + # Genericity assumption (can be checked) + if checks + local Fnew = _fbr(I,q).gens + new_lucky_prime = _generate_lucky_primes(Fnew, one(ZZ)<<30, one(ZZ)<<31-1, 1) |> first + local INEW = Ideal(change_base_ring.(Ref(GF(new_lucky_prime)), Fnew)) + @assert(dimension(INEW) == I.dim - e, "Non-generic coordinates") + end + + # Terminal case (dim <=1) + if I.dim - e <= 1 + return RMnode(q, [], RMnode[]) + end + + ## sing(Fq) ## + if checks + info_level>0 && println("Check smoothness") + local FNEW = _fbr(computepolar(1:e, I)|> Ideal, q).gens + new_lucky_prime = _generate_lucky_primes(FNEW, one(ZZ)<<30, one(ZZ)<<31-1, 1) |> first + local INEW = Ideal(change_base_ring.(Ref(GF(new_lucky_prime)), FNEW)) + @assert(dimension(INEW) == -1, "Non-empty sing locus!") + end + + ## K(pi_1,Fq) ## + info_level>0 && println("Compute x1-critical points: K1") + K1Fq = computepolar(1:e+1, I) |> Ideal + K1Fq = real_solutions(_fbr(K1Fq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) + + ## K(pi_2, Fq) ## + info_level>0 && println("Compute (x1,x2)-polar variety: W") + K2Fqmins = computepolar(1:e+2, I, only_mins=true) + K2Fq = vcat(I.gens, K2Fqmins) |> Ideal + if checks + local FNEW = _fbr(K2Fq, q).gens + new_lucky_prime = _generate_lucky_primes(FNEW, one(ZZ)<<30, one(ZZ)<<31-1, 1) |> first + local INEW = Ideal(change_base_ring.(Ref(GF(new_lucky_prime)), FNEW)) + @assert(dimension(INEW) == 1, "Non-generic polar variety") + else + K2Fq.dim = e + 1 + end + RM = RMnode(q, K2Fqmins, RMnode[]) + + ## Points with vertical tg in K(pi_2, Fq) ## + info_level>0 && println("Compute W-critical points with vertical tangent: K1W") + K1WmFq = computepolar(1:e+2, K2Fq, dimproj=e) |> Ideal + K1WmFq = real_solutions(_fbr(K1WmFq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) + + ## New base and query points ## + Cq = isempty(q) ? C : [ c for c in C if q[e] in c[e]] + K1W = vcat(K1Fq, K1WmFq) + # Heuristic to be proven (Reeb's th) + #K1W = K1W[2:end-1] + ########## + K1WRat = _mid_rational_points_inter(getindex.(K1W,e+1), unique(getindex.(Cq, e+1))) + newQ = vcat.(Ref(q), K1WRat) + + # Recursively compute roadmap of possible fibers + if !isempty(newQ) + for newq in newQ + RMFq = _roadmap_rec(I, newq, Cq, info_level, checks) + push!(RM.children, RMFq) + end + end + + if e == 0 + return Roadmap(I, RM) + else + return RM + end +end + +function _mid_rational_points_inter(S::Vector{Vector{T}}, Q::Vector{Vector{T}} = Vector{T}[]) where {T <: QQFieldElem} + # * S is a list of [ [l_1,r_1], ..., [l_n, r_n] ] + # such that the [l_i, r_i] are rational and disjoint open intervals. + # * Same assumptions on Q + # * Intervals in S and Q do not intersect as well + # + # It orders the [l_i,r_i], and compute a list ratioP + # intersecting all intervals of Q and such that + # *strictly* between each of the [l_i,r_i] there is: + # - either at least one element inside an interval of Q + # - or the simplest rational number + isempty(S) && return Q + + S1, Q1 = sort(S, lt=(x, y) -> x[2] <= y[1]), sort(Q, lt=(x, y) -> x[2] <= y[1]) + ratioP = T[] + qidx = 1 + qlen = length(Q1) + + # Handle left gap before first interval + while qidx <= qlen && Q1[qidx][2] < S1[1][1] + ql, qr = Q1[qidx] + push!(ratioP, _open_simplest_between(ql, qr, abs(qr - ql)//1000)) + qidx += 1 + end + + # Loop through gaps between sorted disjoint intervals + for i in 1:(length(S1) - 1) + ri, li1 = S1[i][2], S1[i+1][1] + @assert ri < li1 "Intervals are not disjoint." + inserted = false + while qidx <= qlen && Q1[qidx][2] < li1 + ql, qr = Q1[qidx] + @assert(ql > ri, "A query point might be singular") + push!(ratioP, _open_simplest_between(ql, qr, abs(qr - ql)//1000)) + inserted = true + qidx += 1 + end + @assert qidx > qlen || Q1[qidx][1] > S1[i+1][2] "A query point might be singular" + # If there's already rational between no need to add new + !inserted && push!(ratioP, _open_simplest_between(ri, li1, abs(li1 - ri)//1000)) + end + + # Append remaining right-side Q points + while qidx <= qlen + ql, qr = Q1[qidx] + push!(ratioP, _open_simplest_between(ql, qr, abs(qr - ql)//1000)) + qidx += 1 + end + + return ratioP +end + +function _open_simplest_between(a::QQFieldElem, b::QQFieldElem, eps::QQFieldElem) + # We choose the simplest in absolute value + if -a > b # this means a is negative and the largest in absolute value + return -simplest_between(-a - eps, -b + eps) + else + return simplest_between( a + eps, b - eps) + end +end + + diff --git a/src/connectivity/tools.jl b/src/connectivity/tools.jl new file mode 100644 index 0000000..6ceb7ae --- /dev/null +++ b/src/connectivity/tools.jl @@ -0,0 +1,71 @@ +""" + order_permut2d(L) + +Return the indices that sort a 2D collection by its values. + +Given a nested array `L` (e.g. a matrix or vector of vectors), this function +flattens all entries, sorts them by value, and returns the corresponding +2D indices. +""" +function order_permut2d(L) + # Create a list of tuples with elements and their corresponding indices + LL = [(L[i][j], (i, j)) for i in eachindex(L) for j in eachindex(L[i])] + # Sort the enumerated list based on the values + sorted_LL = sort(LL, by = x -> x[1]) + # Extract the sorted values and their corresponding indices + sorted_ind = [pair[2] for pair in sorted_LL] + + return sorted_ind +end + +""" + trimat_rand(A, n; up=true, range=-100:100) + +Generate a random triangular matrix over a given coefficient domain. +""" +function trimat_rand(A, n; up=true, range=-100:100) + if up + return [ i==j ? one(A) : (ij ? A(rand(range)) : zero(A)) for i in 1:n, j in 1:n ] + end +end + +""" + int_coeffs(F::Vector{P}) where P <: Union{QQMPolyRingElem, QQPolyRingElem} + int_coeffs(f::Union{QQMPolyRingElem, QQPolyRingElem}) + +Clear denominators of polynomial coefficients. + +This function rescales each polynomial so that all coefficients become integers +(by multiplying with the least common multiple of denominators). +""" +function int_coeffs(F::Vector{P} where P <: Union{QQMPolyRingElem, QQPolyRingElem}) + CD = [ iszero(f) ? f : lcm(map(denominator, collect(coefficients(f)))) for f in F ] + return (F .* CD) +end + +int_coeffs(f::Union{QQMPolyRingElem, QQPolyRingElem}) = int_coeffs([f])[1] + +macro iftime(v, ex) + quote + if $(esc(v)) + @time $(esc(ex)) + else + $(esc(ex)) + end + end +end + +# We compute a zero-dim param of I w.r.t. a new generic variable v +# given by Sv (Symbol) and defined by the linear form in cfs_lf +function param_use_lfs(I::Ideal, cfs_lf::Vector{T} where T<: RingElem, Sv::Symbol) + R = parent(I) + new_RS = vcat(symbols(R), Sv) + _, new_V = polynomial_ring(base_ring(R), new_RS) + + new_gens = [change_ringvar(f, new_RS) for f in I.gens] + I_new = Ideal(vcat(new_gens, [new_V[end] + transpose(cfs_lf) * new_V[1:end-1]])) + + return rational_parametrization(I_new) +end \ No newline at end of file diff --git a/src/exports.jl b/src/exports.jl index 5e3d57e..3ac8301 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -5,4 +5,7 @@ export polynomial_ring, MPolyRing, GFElem, MPolyRingElem, prime_field, sig_groebner_basis, cyclic, leading_coefficient, equidimensional_decomposition, homogenize, dimension, FqMPolyRingElem, hilbert_series, hilbert_dimension, hilbert_degree, hilbert_polynomial, - rational_curve_parametrization + rational_curve_parametrization, roadmap, computepolar, all_eqs, + all_base_pts, nb_nodes, connected_components, number_of_connected_components, + group_by_component, merge_graphs, build_graph_data,curve_graph, + curve_arrangement_graph \ No newline at end of file diff --git a/src/types.jl b/src/types.jl index 20896d6..47cbb17 100644 --- a/src/types.jl +++ b/src/types.jl @@ -74,4 +74,78 @@ Base.parent(I::Ideal) = Nemo.parent(I.gens[1]) Base.show(io::IO, I::Ideal) = print(io, I.gens) -Base.getindex(I::Ideal, idx::Int) = I.gens[idx] +Base.getindex(I::Ideal, idx::Union{Int, UnitRange}) = I.gens[idx] + +Base.lastindex(I::Ideal) = lastindex(I.gens) + +# ---------- ROADMAPS ---------- + +mutable struct RMnode + base_pt::Vector{QQFieldElem} + polar_eqs::Vector{QQMPolyRingElem} + children::Vector{RMnode} +end + +mutable struct Roadmap + initial_ideal::Ideal{QQMPolyRingElem} + root::RMnode +end + +function _collect_roadmap(RMn::RMnode, F) + data = [F(RMn)] + for child in RMn.children + append!(data, _collect_roadmap(child, F)) + end + return data +end + +function _fbr(I::Ideal{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) + @assert(!isempty(I.gens), "Empty polynomial vector") + vars = gens(parent(first(I.gens))) + return Ideal(vcat(I.gens, [vars[i] - Q[i] for i in 1:min(length(vars),length(Q))])) +end + +function all_eqs(RM::Roadmap) + func(s) = _fbr(vcat(RM.initial_ideal.gens, s.polar_eqs) |> Ideal, s.base_pt) + return _collect_roadmap(RM.root, func) +end + +function all_base_pts(RM::Roadmap) + return _collect_roadmap(RM.root, s->s.base_pt) +end + +function nb_nodes(RM::Roadmap) + return length(_collect_roadmap(RM.root, s -> true)) +end + +Base.show(io::IO, RM::Roadmap) = print(io, all_base_pts(RM)) +Base.getindex(RM::Roadmap, idx::Union{Int, UnitRange}) = all_eqs(RM)[idx] +Base.lastindex(RM::Roadmap) = nb_nodes(RM) +Base.length(RM::Roadmap) = nb_nodes(RM) + +# ---------- CURVE CONNECTIVITY ---------- + +struct CurveGraph{T} + vertices::Vector{Tuple{T, T}} + edges::Vector{Tuple{Int, Int}} + control_nodes::Dict{Int, Vector{Int}} +end + +# ----- PLOT ---- + +struct EdgeGroup + edges::Vector{Tuple{Vector{Float64}, Vector{Float64}}} + color::String + width::Float64 +end + +struct PointGroup + vertices::Tuple{Vector{Float64}, Vector{Float64}} + color::String + marker::Symbol +end + +struct GraphPlotData + edge_group::EdgeGroup + point_groups::Vector{PointGroup} +end \ No newline at end of file diff --git a/test/connectivity/curve-connectivity.jl b/test/connectivity/curve-connectivity.jl new file mode 100644 index 0000000..edfb610 --- /dev/null +++ b/test/connectivity/curve-connectivity.jl @@ -0,0 +1,11 @@ +@testset "Algorithms -> Curve Graph" begin + R, (x,y) = polynomial_ring(QQ, [:x,:y]) + + # param2 (generic) + f = -49303382*x^4 + 9395599*x^3*y + 67366686*x^3 - 27407214*x^2*y^2 - 71298014*x^2*y - 24150320*x^2 - 12067817*x*y^3 + 6797370*x*y^2 + 20838420*x*y - 18118613*x + 2357706*y^4 - 13736522*y^3 - 37604516*y^2 - 13868221*y + 6980802 + df = 9395599*x^3 - 54814428*x^2*y - 71298014*x^2 - 36203451*x*y^2 + 13594740*x*y + 20838420*x + 9430824*y^3 - 41209566*y^2 - 75209032*y - 13868221 + g = 32091920*x^4 - 97772598*x^3*y - 245256584*x^3 + 99093410*x^2*y^2 + 162335273*x^2*y + 239310556*x^2 + 28995368*x*y^3 - 59544597*x*y^2 - 20499914*x*y + 98243402*x + 22738226*y^3 + 111105506*y^2 + 86287748*y - 62439535 + P = AlgebraicSolving.RationalCurveParametrization([:z,:x,:y], Vector{Vector{ZZRingElem}}(), f, df, [g]) + G = curve_graph(P) + @test number_of_connected_components(G) == 3 +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index f6cda36..c1c4b78 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,8 +11,9 @@ include("algorithms/hilbert.jl") include("algorithms/decomposition.jl") include("algorithms/param-curves.jl") include("examples/katsura.jl") -include("interp/thiele.jl") +include("interp/thiele.jl") include("interp/newton.jl") include("interp/cuyt_lee.jl") include("interp/resultant.jl") +include("connectivity/curve-connectivity.jl") end