More flexible choice of type, algo, solver (#38)

* More flexible choice of type, algo, solver

* Add NoAlgorithm

* Fix

* Fix format

Author: blegat

src/groebner.jl

@@ -106,8 +106,15 @@ end
 # TODO sugar and double sugar selections
 
 default_gröbner_basis_algorithm(p) = Buchberger()
+
 abstract type AbstractGröbnerBasisAlgorithm end
 
+function promote_for end
+
+struct NoAlgorithm <: AbstractGröbnerBasisAlgorithm end
+
+promote_for(::Type{T}, ::Type{NoAlgorithm}) where {T} = T
+
 struct Buchberger <: AbstractGröbnerBasisAlgorithm
     ztol::Float64
     pre!::Function
@@ -118,6 +125,15 @@ function Buchberger(ztol = Base.rtoldefault(Float64))
     return Buchberger(ztol, presort!, normal_selection)
 end
 
+promote_for_division(::Type{T}) where {T} = T
+promote_for_division(::Type{T}) where {T<:Integer} = Rational{big(T)}
+promote_for_division(::Type{Rational{T}}) where {T} = Rational{big(T)}
+function promote_for_division(::Type{Complex{T}}) where {T}
+    return Complex{promote_for_division(T)}
+end
+
+promote_for(::Type{T}, ::Type{Buchberger}) where {T} = promote_for_division(T)
+
 # Taken from Theorem 9 of
 # Ideals, Varieties, and Algorithms
 # Cox, Little and O'Shea, Fourth edition

src/ideal.jl

@@ -8,13 +8,6 @@ function ideal(p::FullSpace, _ = default_gröbner_basis_algorithm(p))
 end
 Base.rem(p::AbstractPolynomialLike, ::EmptyPolynomialIdeal) = p
 
-promote_for_division(::Type{T}) where {T} = T
-promote_for_division(::Type{T}) where {T<:Integer} = Rational{big(T)}
-promote_for_division(::Type{Rational{T}}) where {T} = Rational{big(T)}
-function promote_for_division(::Type{Complex{T}}) where {T}
-    return Complex{promote_for_division(T)}
-end
-
 mutable struct PolynomialIdeal{T,PT<:APL{T},A<:AbstractGröbnerBasisAlgorithm} <:
                AbstractPolynomialIdeal
     p::Vector{PT}
@@ -28,7 +21,7 @@ function PolynomialIdeal{T,PT}(
     return PolynomialIdeal{T,PT,A}(p, false, algo)
 end
 function PolynomialIdeal(p::Vector{PT}, algo) where {T,PT<:APL{T}}
-    S = promote_for_division(T)
+    S = promote_for(T, typeof(algo))
     return PolynomialIdeal{S,polynomial_type(PT, S)}(
         AbstractVector{polynomial_type(PT, S)}(p),
         algo,

src/solve.jl

@@ -7,6 +7,10 @@ Solver of algebraic equations.
 """
 abstract type AbstractAlgebraicSolver end
 
+function default_gröbner_basis_algorithm(p, ::AbstractAlgebraicSolver)
+    return default_gröbner_basis_algorithm(p)
+end
+
 """
     solvealgebraicequations(V::AbstractAlgebraicSet, algo::AbstractAlgebraicSolver)::Union{Nothing, Vector{<:Vector}}}
 
@@ -146,6 +150,13 @@ function algebraic_solver(
     return SolverUsingMultiplicationMatrices(algo, solver)
 end
 
+function promote_for(
+    ::Type{T},
+    ::Type{<:SolverUsingMultiplicationMatrices},
+) where {T}
+    return float(T)
+end
+
 function default_multiplication_matrices_algorithm(p)
     return GröbnerBasisMultiplicationMatricesAlgorithm()
 end
@@ -179,3 +190,7 @@ function default_algebraic_solver(
         solver,
     )
 end
+
+function default_algebraic_solver(p, ::AbstractGröbnerBasisAlgorithm)
+    return default_algebraic_solver(p)
+end

src/variety.jl

@@ -5,6 +5,10 @@ struct DefaultAlgebraicSetLibrary{S<:AbstractAlgebraicSolver}
     solver::S
 end
 
+function default_gröbner_basis_algorithm(p, lib::DefaultAlgebraicSetLibrary)
+    return default_gröbner_basis_algorithm(p, lib.solver)
+end
+
 function default_algebraic_set_library(
     ::Vector{<:APL},
     solver::AbstractAlgebraicSolver,
@@ -34,7 +38,7 @@ function AlgebraicSet{T,PT,A,S,U}(
     return AlgebraicSet{T,PT,A,S,U}(I, false, Vector{U}[], false, false, solver)
 end
 function AlgebraicSet(I::PolynomialIdeal{T,PT,A}, solver::S) where {T,PT,A,S}
-    return AlgebraicSet{T,PT,A,S,float(T)}(I, solver)
+    return AlgebraicSet{T,PT,A,S,promote_for(T, S)}(I, solver)
 end
 function AlgebraicSet{T,PT}() where {T,PT}
     return AlgebraicSet(PolynomialIdeal{T,PT}(), default_algebraic_solver(T))
@@ -69,52 +73,46 @@ function Base.convert(
     )
 end
 
-function algebraic_set(p::Vector, lib::DefaultAlgebraicSetLibrary)
-    return AlgebraicSet(p, default_gröbner_basis_algorithm(p), lib.solver)
-end
 function algebraic_set(
     p::Vector,
     algo::AbstractGröbnerBasisAlgorithm = default_gröbner_basis_algorithm(p),
-    lib::DefaultAlgebraicSetLibrary = default_algebraic_set_library(p),
+    lib::DefaultAlgebraicSetLibrary = default_algebraic_set_library(p, algo);
+    projective::Bool = false,
 )
-    return AlgebraicSet(p, algo, lib.solver)
-end
-function algebraic_set(p::Vector, solver)
-    return algebraic_set(p, default_algebraic_set_library(p, solver))
-end
-function algebraic_set(p::Vector, algo::AbstractGröbnerBasisAlgorithm, solver)
-    return algebraic_set(p, algo, default_algebraic_set_library(p, solver))
+    set = AlgebraicSet(p, algo, lib.solver)
+    set.projective = projective
+    return set
 end
 
-function projective_algebraic_set(p::Vector, lib::DefaultAlgebraicSetLibrary)
-    return projective_algebraic_set(
+function algebraic_set(p::Vector, lib::DefaultAlgebraicSetLibrary; kws...)
+    return algebraic_set(
         p,
-        default_gröbner_basis_algorithm(p),
-        lib.solver,
+        default_gröbner_basis_algorithm(p, lib),
+        lib.solver;
+        kws...,
     )
 end
-function projective_algebraic_set(
+
+function algebraic_set(p::Vector, solver; kws...)
+    return algebraic_set(p, default_algebraic_set_library(p, solver); kws...)
+end
+
+function algebraic_set(
     p::Vector,
-    algo::AbstractGröbnerBasisAlgorithm = default_gröbner_basis_algorithm(p),
-    lib::DefaultAlgebraicSetLibrary = default_algebraic_set_library(p),
+    algo::AbstractGröbnerBasisAlgorithm,
+    solver;
+    kws...,
 )
-    V = AlgebraicSet(p, algo, lib.solver)
-    V.projective = true
-    return V
-end
-function projective_algebraic_set(p::Vector, algo, solver)
-    return projective_algebraic_set(
+    return algebraic_set(
         p,
         algo,
-        default_algebraic_set_library(p, solver),
+        default_algebraic_set_library(p, solver);
+        kws...,
     )
 end
-function projective_algebraic_set(p::Vector, solver)
-    return projective_algebraic_set(
-        p,
-        default_gröbner_basis_algorithm(p),
-        default_algebraic_set_library(p, solver),
-    )
+
+function projective_algebraic_set(p::Vector, args...)
+    return algebraic_set(p, args...; projective = true)
 end
 
 ideal(V::AlgebraicSet) = V.I

test/macro.jl

@@ -1,7 +1,17 @@
+module TestMacro
+
+using Test
+using SemialgebraicSets
+using MultivariatePolynomials
+
 struct DummySolver <: SemialgebraicSets.AbstractAlgebraicSolver end
+function SemialgebraicSets.default_gröbner_basis_algorithm(p, ::DummySolver)
+    return SemialgebraicSets.default_gröbner_basis_algorithm(p)
+end
+SemialgebraicSets.promote_for(::Type{T}, ::Type{DummySolver}) where {T} = T
 
-@testset "Basic semialgebraic set" begin
-    Mod.@polyvar x y
+function runtests()
+    Main.Mod.@polyvar x y
     @test isa(FullSpace(), FullSpace)
     V = @set x * y == 1
     @test V isa AlgebraicSet{Rational{BigInt}}
@@ -161,3 +171,7 @@ Algebraic Set defined by no equality
         end
     end
 end
+
+end
+
+TestMacro.runtests()

test/variety.jl

@@ -1,11 +1,12 @@
 using Test
 
-using MultivariatePolynomials
-const MP = MultivariatePolynomials
+import MultivariatePolynomials as MP
 
 using SemialgebraicSets
 
 struct DummyAlgo <: SemialgebraicSets.AbstractGröbnerBasisAlgorithm end
+SemialgebraicSets.default_algebraic_solver(::Any, ::DummyAlgo) = DummySolver()
+SemialgebraicSets.promote_for(::Type, ::Type{DummyAlgo}) = Int
 function SemialgebraicSets.gröbner_basis!(
     ::AbstractVector{<:MP.APL},
     ::DummyAlgo,
@@ -14,6 +15,10 @@ function SemialgebraicSets.gröbner_basis!(
 end
 
 struct DummySolver <: SemialgebraicSets.AbstractAlgebraicSolver end
+function SemialgebraicSets.default_gröbner_basis_algorithm(::Any, ::DummySolver)
+    return SemialgebraicSets.NoAlgorithm()
+end
+SemialgebraicSets.promote_for(::Type{T}, ::Type{DummySolver}) where {T} = T
 function SemialgebraicSets.solve(
     ::SemialgebraicSets.AlgebraicSet,
     ::DummySolver,
@@ -24,17 +29,30 @@ end
 function algo_and_solver_test(func)
     Mod.@polyvar x y
     p = [y^2 - x * y, x^2 + x * y]
-    V = func(p, DummyAlgo())
+    V = func(p, DummyAlgo(), SemialgebraicSets.default_algebraic_set_library(p))
+    @test V.I.algo isa DummyAlgo
     @test_throws ErrorException("Dummy algo") p[1] in V
     @test_throws ErrorException("Dummy algo") rem(p[1], V.I)
+    V = func(p, DummyAlgo())
+    @test V.I.algo isa DummyAlgo
+    @test V.solver isa DummySolver
+    @test eltype(V) == Vector{Int}
+    @test_throws ErrorException("Dummy solver") p[1] in V
+    @test_throws ErrorException("Dummy algo") rem(p[1], V.I)
     V = func(p, DummySolver())
+    @test V.I.algo isa SemialgebraicSets.NoAlgorithm
+    @test V.solver isa DummySolver
+    @test eltype(V) == Vector{Int}
     @test_throws ErrorException("Dummy solver") collect(V)
     V = func(p, DummyAlgo(), DummySolver())
+    @test V.I.algo isa DummyAlgo
+    @test V.solver isa DummySolver
+    @test eltype(V) == Vector{Int}
     @test_throws ErrorException("Dummy solver") collect(V)
     @test_throws ErrorException("Dummy algo") rem(p[1], V.I)
 end
 
 @testset "Algo and solver" begin
     algo_and_solver_test(algebraic_set)
-    algo_and_solver_test(projective_algebraic_set)
+    #algo_and_solver_test(projective_algebraic_set)
 end