closes #16

benchmark/biomodels_groebner.jl

@@ -0,0 +1,42 @@
+
+if !isdefined(Main, :Groebner)
+    import Groebner
+end
+
+include("biomodels/parser.jl")
+
+import AbstractAlgebra
+using BenchmarkTools
+using Logging
+global_logger(ConsoleLogger(stderr, Logging.Error))
+
+BenchmarkTools.DEFAULT_PARAMETERS.samples = 3
+
+function benchmark_system_my(system)
+    system = Groebner.change_ordering(system, :degrevlex)
+    Groebner.groebner([system[1]])
+
+    # gb = Groebner.groebner(system)
+    # println("length = $(length(gb))")
+    # println("degree = $(maximum(AbstractAlgebra.total_degree, gb))")
+
+    @btime gb = Groebner.groebner($system)
+    # gb = Groebner.groebner(system, reduced=false)
+end
+
+function run_f4_biomodels_benchmarks()
+    systems = read_BIOMDs(5:20)
+
+    println()
+    for (name, system) in systems
+        println("$name")
+        # benchmark_system_my(system)
+        benchmark_system_my(system)
+    end
+end
+
+run_f4_biomodels_benchmarks()
+
+#=
+
+=#

benchmark/run_biomodels_groebner.jl → benchmark/biomodels_groebner_vs_singular.jl

@@ -5,6 +5,7 @@ end
 
 include("biomodels/parser.jl")
 
+import Singular
 import AbstractAlgebra
 using BenchmarkTools
 using Logging
@@ -44,18 +45,20 @@ function benchmark_system_singular(system)
     @btime Singular.std($ideal_s)
 end
 
-function run_f4_MQ_benchmarks()
+function run_biomodels_benchmarks()
     systems = read_BIOMDs(5:20)
 
     println()
     for (name, system) in systems
         println("$name")
-        # benchmark_system_my(system)
+        println("my")
         benchmark_system_my(system)
+        println("singular")
+        benchmark_system_singular(system)
     end
 end
 
-run_f4_MQ_benchmarks()
+run_biomodels_benchmarks()
 
 #=
 

src/Groebner.jl

@@ -48,7 +48,7 @@ end
 # most of these are not needed
 # TODO
 import AbstractAlgebra
-import AbstractAlgebra.Generic: MPoly, GFElem, MPolyRing
+import AbstractAlgebra.Generic: MPoly, GFElem, MPolyRing, Poly
 import AbstractAlgebra: leading_term, QQ, PolynomialRing, terms,
                         coeff, divides, base_ring, elem_type,
                         rref, isconstant, leading_coefficient,

src/io.jl

@@ -70,13 +70,13 @@ function convert_to_internal(
 end
 
 function extract_ring(R::T) where {T}
-    @assert hasmethod(nvars, Tuple{typeof(R)})
-    @assert hasmethod(ordering, Tuple{typeof(R)})
-    @assert hasmethod(characteristic, Tuple{typeof(R)})
+    @assert hasmethod(nvars, Tuple{T})
+    # @assert hasmethod(ordering, Tuple{typeof(R)})
+    @assert hasmethod(characteristic, Tuple{T})
 
     nv     = nvars(R)
     explen = nv + 1
-    ord    = ordering(R)
+    ord = hasmethod(ordering, Tuple{T}) ? ordering(R) : :lex
     ch     = characteristic(R)
 
     @assert nv + 1 == explen
@@ -94,11 +94,27 @@ function extract_coeffs(ring::PolyRing, orig_polys::Vector{T}) where {T}
     end
 end
 
+function extract_coeffs_ff(ring::PolyRing, poly::Poly)
+    reverse(map(CoeffFF ∘ data, filter(!iszero, collect(coefficients(poly)))))
+end
+
+function extract_coeffs_qq(ring::PolyRing, poly::Poly)
+    reverse(map(Rational, filter(!iszero, collect(coefficients(poly)))))
+end
+
+function extract_coeffs_ff(ring::PolyRing, poly::MPoly)
+    map(CoeffFF ∘ data, coefficients(poly))
+end
+
+function extract_coeffs_qq(ring::PolyRing, poly::MPoly)
+    map(Rational, coefficients(poly))
+end
+
 function extract_coeffs_ff(ring::PolyRing, orig_polys::Vector{T}) where {T}
     npolys = length(orig_polys)
     coeffs = Vector{Vector{CoeffFF}}(undef, npolys)
     for i in 1:npolys
-        coeffs[i] = map(CoeffFF ∘ data, coefficients(orig_polys[i]))
+        coeffs[i] = extract_coeffs_ff(ring, orig_polys[i])
     end
     coeffs
 end
@@ -107,22 +123,39 @@ function extract_coeffs_qq(ring::PolyRing, orig_polys::Vector{T}) where {T}
     npolys = length(orig_polys)
     coeffs = Vector{Vector{CoeffQQ}}(undef, npolys)
     for i in 1:npolys
-        coeffs[i] = map(Rational, coefficients(orig_polys[i]))
+        coeffs[i] = extract_coeffs_qq(ring, orig_polys[i])
     end
     coeffs
 end
 
+function extract_exponents(ring, poly::MPoly)
+    exps = Vector{ExponentVector}(undef, length(poly))
+    @inbounds for j in 1:length(poly)
+        exps[j] = ExponentVector(undef, ring.explen)
+        exps[j][1:ring.nvars] .= exponent_vector(poly, j)
+        exps[j][end] = sum(exps[i][j][k] for k in 1:ring.nvars)
+    end
+    exps
+end
+
+function extract_exponents(ring, poly::Poly)
+    # Why define length of univeriate polynomial to be the dense length??
+    exps = Vector{ExponentVector}(undef, 0)
+    @inbounds while !iszero(poly)
+        push!(exps, ExponentVector(undef, ring.explen))
+        exps[end][1] = degree(poly)
+        exps[end][end] = degree(poly)
+        poly = AbstractAlgebra.tail(poly)
+    end
+    exps
+end
+
 function extract_exponents(ring::PolyRing, orig_polys::Vector{T}) where {T}
     npolys = length(orig_polys)
     exps = Vector{Vector{ExponentVector}}(undef, npolys)
     for i in 1:npolys
         poly = orig_polys[i]
-        exps[i] = Vector{ExponentVector}(undef, length(poly))
-        @inbounds for j in 1:length(poly)
-            exps[i][j] = ExponentVector(undef, ring.explen)
-            exps[i][j][1:ring.nvars] .= exponent_vector(poly, j)
-            exps[i][j][end] = sum(exps[i][j][k] for k in 1:ring.nvars)
-        end
+        exps[i] = extract_exponents(ring, poly)
     end
     exps
 end
@@ -400,12 +433,14 @@ function convert_to_output(
     @assert hasmethod(base_ring, Tuple{typeof(origring)})
 
     # TODO: hardcoded
-    (metainfo.targetord != ordering(origring)) && @warn "Unknown polynomial type. Computed basis is in $(metainfo.targetord), but terms are ordered in $(ordering(origring)) in output"
+    if hasmethod(ordering, Tuple{M})
+        (metainfo.targetord != ordering(origring)) && @warn "Unknown polynomial type. Computed basis is in $(metainfo.targetord), but terms are ordered in $(ordering(origring)) in output"
+    end
 
     etype = elem_type(base_ring(origring))
     # rather weak but okay for now
     if etype <: Integer
-        coeffs_zz = scale_denominators!(gbcoeffs)
+        coeffs_zz = scale_denominators(gbcoeffs)
         convert_to_output(origring, gbexps, coeffs_zz, metainfo)
     else
         convert_to_output(origring, gbexps, gbcoeffs, metainfo)
@@ -416,7 +451,26 @@ function convert_to_output(
             origring::M,
             gbexps::Vector{Vector{ExponentVector}},
             gbcoeffs::Vector{Vector{I}},
-            metainfo::GroebnerMetainfo) where {M, I}
+            metainfo::GroebnerMetainfo) where {M<:AbstractAlgebra.Generic.PolyRing, I}
+
+    ground   = base_ring(origring)
+    exported = Vector{elem_type(origring)}(undef, length(gbexps))
+    @inbounds for i in 1:length(gbexps)
+        cfs = Vector{}
+        cfs    = zeros(ground, gbexps[i][1][1] + 1)
+        for (idx, j) in enumerate(gbexps[i])
+            cfs[j[1] + 1] = ground(gbcoeffs[i][idx])
+        end
+        exported[i] = origring(cfs)
+    end
+    exported
+end
+
+function convert_to_output(
+            origring::M,
+            gbexps::Vector{Vector{ExponentVector}},
+            gbcoeffs::Vector{Vector{I}},
+            metainfo::GroebnerMetainfo) where {M<:AbstractAlgebra.Generic.MPolyRing, I}
 
     ground   = base_ring(origring)
     exported = Vector{elem_type(origring)}(undef, length(gbexps))

test/runtests.jl

@@ -24,9 +24,10 @@ end
 
 @testset "All tests" begin
 
-    @includetests ["commons", "io", "io_consistency"]
+    @includetests ["commons", "io", "io_consistency",
+                    "univariate_aa"]
 
-    @includetests ["crt_reconstruction"] # "mod_reconstruction"]#, 
+    @includetests ["crt_reconstruction"] # "mod_reconstruction"]#,
 
     @includetests ["core_f4_reduce", "core_f4_stress",
                     "core_f4", "rational_f4", "groebner_certify"]

test/univariate_aa.jl

@@ -0,0 +1,13 @@
+
+using AbstractAlgebra
+
+@testset "AbstractAlgebra io consistency" begin
+    R, x = PolynomialRing(GF(2^31 - 1), "x")
+    @test Groebner.groebner([x^2 - 4, x + 2]) == [x + 2]
+
+    R, x = PolynomialRing(QQ, "x")
+    @test Groebner.groebner([x^2 - 4, x + 2]) == [x + 2]
+
+    R, x = PolynomialRing(ZZ, "x")
+    @test Groebner.groebner([x^2 - 4, x + 2]) == [x + 2]
+end