Skip to content

Commit ad10edd

Browse files
wegankwegank
committed
Add CAD-based sample points
1 parent a3713b0 commit ad10edd

6 files changed

Lines changed: 269 additions & 1 deletion

File tree

‎src/AlgebraicSolving.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@ include("algorithms/dimension.jl")
1616
include("algorithms/decomposition.jl")
1717
include("algorithms/param-curve.jl")
1818
include("algorithms/hilbert.jl")
19+
include("algorithms/sampling.jl")
1920
#= siggb =#
2021
include("siggb/siggb.jl")
2122
#= progress =#

‎src/algorithms/sampling.jl‎

Lines changed: 228 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,228 @@
1+
@doc Markdown.doc"""
2+
_to_multivariate(R::QQMPolyRing, p::QQPolyRingElem)::QQMPolyRingElem
3+
4+
Convert a univariate polynomial `p` to a multivariate polynomial in the ring `R` with the same variable.
5+
6+
**Note**: This is an internal function.
7+
"""
8+
function _to_multivariate(R::QQMPolyRing, p::QQPolyRingElem)::QQMPolyRingElem
9+
@assert length(gens(R)) == 1
10+
M = MPolyBuildCtx(R)
11+
for i in 0:length(p)
12+
push_term!(M, coeff(p, i), [i])
13+
end
14+
finish(M)
15+
end
16+
17+
@doc Markdown.doc"""
18+
_sample_points(f::Vector{QQPolyRingElem}; worker_pool::AbstractWorkerPool=default_worker_pool())::Vector{QQFieldElem}
19+
20+
Given a vector of univariate polynomials, returns a vector of rational numbers such that each connected component of the non-vanishing set of the polynomials contains at least one of the returned points.
21+
22+
**Note**: This is an internal function.
23+
"""
24+
function _sample_points(f::Vector{QQPolyRingElem}; worker_pool::AbstractWorkerPool=default_worker_pool())::Vector{QQFieldElem}
25+
R, _ = polynomial_ring(QQ, [:x])
26+
fₓ = map(p -> _to_multivariate(R, p), f)
27+
factors = unique([p for fₓ′ in fₓ for (p, _) in factor(fₓ′)])
28+
# We map each factor to its roots, and each root to its factor
29+
roots_by_factor = Dict{QQMPolyRingElem,Vector{Vector{Vector{QQFieldElem}}}}()
30+
factors_by_root = Dict{Vector{Vector{QQFieldElem}},QQMPolyRingElem}()
31+
for factor in factors
32+
roots = real_solutions(Ideal([factor]); interval=true, worker_pool=worker_pool)
33+
roots_by_factor[factor] = roots
34+
for r in roots
35+
factors_by_root[r] = factor
36+
end
37+
end
38+
# We order intervals by their left endpoint
39+
roots = sort(collect(keys(factors_by_root)), by=x -> x[1][1])
40+
# If intervals are not ordered by their right endpoint,
41+
# we merge the offending factors
42+
while (!issorted(roots, by=x -> x[1][2]))
43+
i = findfirst(i -> roots[i][1][2] > roots[i+1][1][2], 1:length(roots)-1)
44+
bad_factors = [factors_by_root[roots[i]], factors_by_root[roots[i+1]]]
45+
# Remove old factors and roots
46+
for bad_factor in bad_factors
47+
for r in roots_by_factor[bad_factor]
48+
delete!(factors_by_root, r)
49+
end
50+
delete!(roots_by_factor, bad_factor)
51+
end
52+
# Replace with merged factor and its roots
53+
merged_factor = prod(bad_factors)
54+
merged_roots = real_solutions(Ideal([merged_factor]); interval=true, worker_pool=worker_pool)
55+
roots_by_factor[merged_factor] = merged_roots
56+
for r in merged_roots
57+
factors_by_root[r] = merged_factor
58+
end
59+
roots = sort(collect(keys(factors_by_root)), by=r -> r[1][1])
60+
end
61+
# Now the intervals are properly ordered, we can sample points
62+
points = Vector{QQFieldElem}()
63+
if length(roots) == 0
64+
push!(points, QQ(0))
65+
return points
66+
end
67+
push!(points, floor(roots[1][1][1]) - 1)
68+
for i in 1:length(roots)-1
69+
push!(points, (roots[i][1][2] + roots[i+1][1][1]) // 2)
70+
end
71+
push!(points, ceil(roots[end][1][2]) + 1)
72+
return points
73+
end
74+
75+
function _sample_points_0(f::Vector{QQMPolyRingElem})::Vector{Vector{QQFieldElem}}
76+
@assert all(map(is_constant, f))
77+
return [Vector{QQFieldElem}()]
78+
end
79+
80+
function _sample_points_1(f::Vector{QQMPolyRingElem}; worker_pool::AbstractWorkerPool=default_worker_pool())::Vector{Vector{QQFieldElem}}
81+
@assert all(map(is_univariate, f))
82+
R, _ = polynomial_ring(QQ, :x)
83+
return [[p] for p in _sample_points(map(p -> to_univariate(R, p), f); worker_pool=worker_pool)]
84+
end
85+
86+
function _sample_points_desc(n::Int)::String
87+
if n == 0
88+
return "Constant"
89+
elseif n == 1
90+
return "Univariate"
91+
elseif n == 2
92+
return "Bivariate"
93+
elseif n == 3
94+
return "Trivariate"
95+
else
96+
return "Multivariate"
97+
end
98+
end
99+
100+
function _sample_points_2(
101+
f::Vector{QQMPolyRingElem},
102+
xs::Vector{QQMPolyRingElem};
103+
nr_thrds::Int=1,
104+
worker_pool::AbstractWorkerPool=default_worker_pool(),
105+
show_progress::Bool=false,
106+
desc::String="$(_sample_points_desc(length(xs))) sample points"
107+
)::Vector{Vector{QQFieldElem}}
108+
@assert length(xs) >= 2
109+
x₁ = xs[1:end-1]
110+
x₂ = xs[end]
111+
factors = unique([p for fₓ in f for (p, _) in factor(fₓ)])
112+
v = Vector{QQMPolyRingElem}()
113+
for i in eachindex(factors)
114+
if !isempty(intersect(x₁, vars(factors[i])))
115+
push!(v, leading_coefficient(factors[i], length(xs)))
116+
end
117+
if x₂ in vars(factors[i])
118+
push!(v, Interpolation.discriminant(factors[i], length(xs); nr_thrds=nr_thrds))
119+
end
120+
end
121+
for i in 1:length(factors)-1
122+
for j in i+1:length(factors)
123+
if x₂ in vars(factors[i]) || x₂ in vars(factors[j])
124+
push!(v, Interpolation.resultant(factors[i], factors[j], length(xs); nr_thrds=nr_thrds))
125+
end
126+
end
127+
end
128+
p₁ = if length(xs) == 2
129+
_sample_points_1(v; worker_pool=worker_pool)
130+
else
131+
_sample_points_2(v, x₁; show_progress=show_progress, worker_pool=worker_pool)
132+
end
133+
prog = Progress.ProgressBar(total=length(p₁); desc=desc, enabled=show_progress)
134+
Progress.update!(prog, 0)
135+
function _points_chunk(p::Vector{Vector{QQFieldElem}})::Vector{Vector{QQFieldElem}}
136+
res_chunk = Vector{Vector{QQFieldElem}}()
137+
for i in eachindex(p)
138+
p₂ = _sample_points_1(map(p′ -> evaluate(p′, x₁, p[i]), f); worker_pool=worker_pool)
139+
for j in eachindex(p₂)
140+
push!(res_chunk, [p[i]; p₂[j][1]])
141+
end
142+
Progress.next!(prog)
143+
end
144+
res_chunk
145+
end
146+
chunk_size = ceil(Int, length(p₁) / nr_thrds)
147+
chunks = [p₁[i:min(i + chunk_size - 1, end)] for i in 1:chunk_size:length(p₁)]
148+
tasks = [Threads.@spawn _points_chunk(chunk) for chunk in chunks]
149+
res = sort(vcat(fetch.(tasks)...))
150+
Progress.finish!(prog)
151+
res
152+
end
153+
154+
@doc Markdown.doc"""
155+
_sample_points(f::Vector{QQMPolyRingElem}; nr_thrds::Int=1, worker_pool::AbstractWorkerPool=default_worker_pool(), show_progress::Bool=false)::Vector{Vector{QQFieldElem}}
156+
157+
Given a vector of multivariate polynomials, returns a vector of points such that each connected component of the non-vanishing set of the polynomials contains at least one of the returned points.
158+
159+
**Note**: This is an internal function.
160+
"""
161+
function _sample_points(
162+
f::Vector{QQMPolyRingElem};
163+
nr_thrds::Int=1,
164+
worker_pool::AbstractWorkerPool=default_worker_pool(),
165+
show_progress::Bool=false
166+
)::Vector{Vector{QQFieldElem}}
167+
if any(is_zero, f)
168+
error("Cannot sample points for polynomials containing zero polynomial")
169+
end
170+
p = Vector{Vector{QQFieldElem}}([])
171+
if length(gens(parent(f[1]))) == 0
172+
p = _sample_points_0(f)
173+
elseif length(gens(parent(f[1]))) == 1
174+
p = _sample_points_1(f; worker_pool=worker_pool)
175+
else
176+
p = _sample_points_2(f, gens(parent(f[1])); nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress)
177+
end
178+
return p
179+
end
180+
181+
@doc Markdown.doc"""
182+
points_per_components(eqs::Vector{QQMPolyRingElem}, pos::Vector{QQMPolyRingElem}, ineqs::Vector{QQMPolyRingElem}; nr_thrds::Int=1, worker_pool::AbstractWorkerPool=default_worker_pool(), info_level::Int=0)::Vector{Vector{Vector{QQFieldElem}}}
183+
184+
Given a semi-algebraic set defined by equations `eqs`, positivity constraints `pos`, and non-vanishing constraints `ineqs`, returns a vector of points meeting all connected components of the set.
185+
Each point is represented by an isolating box, i.e., a vector of intervals represented by pairs of rational numbers.
186+
187+
**NOTE**: Sampling points per components is currently only implemented for systems without equations, i.e., `eqs` must be empty.
188+
189+
# Arguments
190+
- `eqs::Vector{QQMPolyRingElem}`: equations defining the semi-algebraic set.
191+
- `pos::Vector{QQMPolyRingElem}`: positivity constraints defining the semi-algebraic set.
192+
- `ineqs::Vector{QQMPolyRingElem}`: non-vanishing constraints defining the semi-algebraic set.
193+
- `nr_thrds::Int=1`: the number of threads to use for parallel computations.
194+
- `worker_pool::AbstractWorkerPool=default_worker_pool()`: the worker pool to use for parallel computations.
195+
- `info_level::Int=0`: info level printout: off (`0`, default), summary (`1`), detailed (`2`).
196+
197+
# Examples
198+
199+
```jldoctest
200+
julia> using AlgebraicSolving
201+
202+
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
203+
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])
204+
205+
julia> length(points_per_components(QQMPolyRingElem[], [x, y], [x + y - 1])) >= 2
206+
true
207+
```
208+
"""
209+
function points_per_components(
210+
eqs::Vector{QQMPolyRingElem},
211+
pos::Vector{QQMPolyRingElem},
212+
ineqs::Vector{QQMPolyRingElem};
213+
nr_thrds::Int=1,
214+
worker_pool::AbstractWorkerPool=default_worker_pool(),
215+
info_level::Int=0
216+
)::Vector{Vector{Vector{QQFieldElem}}}
217+
if !isempty(eqs)
218+
throw(ArgumentError("Sampling points per components is not implemented for systems involving equations"))
219+
end
220+
points = _sample_points([pos; ineqs]; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=info_level > 0)
221+
res = Vector{Vector{Vector{QQFieldElem}}}()
222+
for point in points
223+
if all(evaluate(constraint, point) > 0 for constraint in pos)
224+
push!(res, map(x -> [x, x], point))
225+
end
226+
end
227+
return res
228+
end

‎src/exports.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ export polynomial_ring, MPolyRing, GFElem, MPolyRingElem,
33
QQ, vars, nvars, ngens, ZZRingElem, QQFieldElem,
44
QQMPolyRingElem, base_ring, coefficient_ring, evaluate,
55
prime_field, sig_groebner_basis, cyclic, leading_coefficient,
6+
points_per_components,
67
equidimensional_decomposition, homogenize, dimension, FqMPolyRingElem,
78
hilbert_series, hilbert_dimension, hilbert_degree, hilbert_polynomial,
89
rational_curve_parametrization

‎src/imports.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ using StaticArrays
77
import Random: MersenneTwister
88
import Logging: ConsoleLogger, with_logger, Warn, Info
99
import Printf: @sprintf, @printf
10-
import Distributed: AbstractWorkerPool, remotecall_fetch
10+
import Distributed: AbstractWorkerPool, default_worker_pool, remotecall_fetch
1111

1212
import Nemo:
1313
bell,

‎test/algorithms/sampling.jl‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
function points_per_components(ineqs::Vector{QQMPolyRingElem})::Vector{Vector{Vector{QQFieldElem}}}
2+
AlgebraicSolving.points_per_components(QQMPolyRingElem[], QQMPolyRingElem[], ineqs)
3+
end
4+
5+
@testset "Algorithms -> Univariate sample points" begin
6+
R, (x,) = polynomial_ring(QQ, ["x"])
7+
f = [one(R)]
8+
@test length(points_per_components(f)) == 1
9+
f = [x + 1, x + 4294967295 // 4294967296, x + 4294967297 // 4294967296]
10+
@test length(points_per_components(f)) == 4
11+
end
12+
13+
@testset "Algorithms -> Bivariate sample points" begin
14+
R, (x, y) = polynomial_ring(QQ, ["x", "y"])
15+
f = [one(R)]
16+
@test length(points_per_components(f)) == 1
17+
f = [x, x + 1]
18+
@test length(points_per_components(f)) == 3
19+
f = [y, y + 1]
20+
@test length(points_per_components(f)) == 3
21+
f = [x, x + 1, y - x, y + x - 1]
22+
@test length(points_per_components(f)) >= 9
23+
end
24+
25+
@testset "Algorithms -> Trivariate sample points" begin
26+
R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"])
27+
f = [one(R)]
28+
@test length(points_per_components(f)) == 1
29+
f = [x, x + 1]
30+
@test length(points_per_components(f)) == 3
31+
f = [y, y + 1]
32+
@test length(points_per_components(f)) == 3
33+
f = [x, x + 1, y - x, y + x - 1]
34+
@test length(points_per_components(f)) >= 9
35+
f = [x, y, z]
36+
@test length(points_per_components(f)) == 8
37+
end

‎test/runtests.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ include("algorithms/dimension.jl")
1010
include("algorithms/hilbert.jl")
1111
include("algorithms/decomposition.jl")
1212
include("algorithms/param-curves.jl")
13+
include("algorithms/sampling.jl")
1314
include("examples/katsura.jl")
1415
include("interp/thiele.jl")
1516
include("interp/newton.jl")

0 commit comments

Comments
 (0)