working kbase (but broken differentiate_output)

src/parametrizations.jl

@@ -222,7 +222,6 @@ $sat_string
         tagged_mqs,
         ordering = ord,
         homogenize = :no,
-        loglevel = _groebner_loglevel[],
     )
     # Relations between tags in K[T]
     relations_between_tags = filter(

src/power_series_utils.jl

@@ -52,7 +52,7 @@ function ps_matrix_inv(
     power_series_ring = base_ring(parent(M))
     result = map(a -> power_series_ring(a), Base.inv(const_term))
     cur_prec = 1
-    while cur_prec < prec
+    while cur_prec <= prec
         result = _matrix_inv_newton_iteration(M, result)
         cur_prec *= 2
     end
@@ -161,7 +161,7 @@ function ps_matrix_homlinear_de(
         (one(parent(A)) + gen(ps_ring) * truncate_matrix(A, cur_prec)) *
         map(e -> truncate(ps_ring(e), cur_prec), Y0)
     Z = map(e -> truncate(ps_ring(e), cur_prec), Base.inv(Y0))
-    while cur_prec < prec
+    while cur_prec <= prec
         matrix_set_precision!(Y, 2 * cur_prec)
         matrix_set_precision!(Z, cur_prec)
         Y, Z = _matrix_homlinear_de_newton_iteration(A, Y, Z, cur_prec)
@@ -260,9 +260,9 @@ function ps_ode_solution(
     end
 
     cur_prec = 1
-    while cur_prec < prec
+    while cur_prec <= prec
         @debug "\t Computing power series solution, currently at precision $cur_prec"
-        new_prec = min(prec, 2 * cur_prec)
+        new_prec = min(prec + 1, 2 * cur_prec)
         for i in 1:length(x_vars)
             set_precision!(solution[x_vars[i]], new_prec)
             set_precision!(solution[x_dot_vars[i]], new_prec)

src/primality_check.jl

@@ -1,8 +1,44 @@
+# ------------------------------------------------------------------------------
+# adapted from https://gitlab.inria.fr/newrur/code/-/blob/main/Julia/RationalUnivariateRepresentation.jl/src/RationalUnivariateRepresentation.jl?ref_type=heads#L180
+# thanks to Alexander Demin
+"""
+    quotient_basis(J::Array{QQMPolyRingElem, 1})
+
+Takes as input a Groebner basis J of a zero-dimensional ideal and
+returns a monomial basis of the quotient ring
+(more precisely, the list of standard monomials)
+"""
+function quotient_basis(J::Array{QQMPolyRingElem, 1})
+    if !Groebner.isgroebner(J)
+        throw(DomainError("Input is not a Groebner basis"))
+    end
+    n = length(gens(parent(first(J))))
+    leading_exponents = [first(Nemo.exponent_vectors(Nemo.leading_monomial(p))) for p in J]
+    exponents_to_check = [[0 for _ in 1:n]]
+    exponents_checked = []
+    basis_exponents = []
+    while length(exponents_to_check) > 0
+        e = popfirst!(exponents_to_check)
+        push!(exponents_checked, e)
+        if !any(map(le -> all(e .>= le), leading_exponents))
+            push!(basis_exponents, e)
+            for i in 1:n
+                next_e = copy(e)
+                next_e[i] += 1
+                if !(next_e in exponents_checked) && !(next_e in exponents_to_check)
+                    push!(exponents_to_check, next_e)
+                end
+            end
+        end
+    end
+    return [prod(gens(parent(first(J))) .^ e) for e in basis_exponents]
+end
+
 # ------------------------------------------------------------------------------
 
 function check_primality_zerodim(J::Array{QQMPolyRingElem, 1})
     J = Groebner.groebner(J)
-    basis = Groebner.kbase(J)
+    basis = quotient_basis(J)
     dim = length(basis)
     S = Nemo.matrix_space(Nemo.QQ, dim, dim)
     matrices = []

test/differentiate_output.jl

@@ -123,6 +123,57 @@ end
         ),
     )
 
+    ode = @ODEmodel(
+        x'(t) = x(t)^2 + 2 * x(t) * y(t) - 3 * a * y(t),
+        y'(t) = x(t)^2 + a * b - b^2 + 4 * b * x(t),
+        y1(t) = a * x(t),
+        y2(t) = b * y(t)^2 - y(t)
+    )
+    push!(
+        test_cases,
+        Dict(
+            :ODE => ode,
+            :ic => Dict(x => Nemo.QQ(rand(1:10)), y => Nemo.QQ(rand(1:10))),
+            :param_vals => Dict(a => Nemo.QQ(rand(1:10)), b => Nemo.QQ(rand(1:10))),
+            :inputs => Dict{P, Array{QQFieldElem, 1}}(),
+            :prec => 7,
+        ),
+    )
+
+    ode = @ODEmodel(
+        x'(t) = x(t)^2 + 2 * x(t) * y(t) - 3 * a * y(t),
+        y'(t) = x(t)^2 + a * b - b^2 + 4 * b * x(t),
+        y1(t) = a * x(t),
+        y2(t) = b * y(t)^2 - y(t)
+    )
+    push!(
+        test_cases,
+        Dict(
+            :ODE => ode,
+            :ic => Dict(x => Nemo.QQ(rand(1:10)), y => Nemo.QQ(rand(1:10))),
+            :param_vals => Dict(a => Nemo.QQ(rand(1:10)), b => Nemo.QQ(rand(1:10))),
+            :inputs => Dict{P, Array{QQFieldElem, 1}}(),
+            :prec => 9,
+        ),
+    )
+
+    ode = @ODEmodel(
+        x'(t) = x(t)^2 + 2 * x(t) * y(t) - 3 * a * y(t),
+        y'(t) = x(t)^2 + a * b - b^2 + 4 * b * x(t),
+        y1(t) = a * x(t),
+        y2(t) = b * y(t)^2 - y(t)
+    )
+    push!(
+        test_cases,
+        Dict(
+            :ODE => ode,
+            :ic => Dict(x => Nemo.QQ(rand(1:10)), y => Nemo.QQ(rand(1:10))),
+            :param_vals => Dict(a => Nemo.QQ(rand(1:10)), b => Nemo.QQ(rand(1:10))),
+            :inputs => Dict{P, Array{QQFieldElem, 1}}(),
+            :prec => 3,
+        ),
+    )
+
     ode = @ODEmodel(x'(t) = u(t) + a, y(t) = x(t))
     push!(
         test_cases,
@@ -201,6 +252,40 @@ end
             :prec => 4,
         ),
     )
+    push!(
+        test_cases,
+        Dict(
+            :ODE => ODE{P}(
+                Dict{P, DType}(
+                    vars[i] => rand_poly(1, vars[1:5]) // (vars[1] + vars[3]) for i in 1:2
+                ),
+                Dict{P, DType}(vars[i] => rand_poly(1, vars[1:5]) for i in 6:7),
+                [vars[5]],
+            ),
+            :ic => Dict(vars[i] => F(rand(1:50)) for i in 1:2),
+            :param_vals => Dict(vars[i + 2] => F(rand(1:50)) for i in 1:2),
+            :inputs => Dict(vars[5] => [F(rand(-30:30)) for i in 1:4]),
+            :prec => 3,
+        ),
+    )
+    push!(
+        test_cases,
+        Dict(
+            :ODE => ODE{P}(
+                Dict{P, DType}(
+                    vars[i] => rand_poly(1, vars[1:5]) // (vars[1] + vars[3]) for i in 1:2
+                ),
+                Dict{P, DType}(vars[i] => rand_poly(1, vars[1:5]) for i in 6:7),
+                [vars[5]],
+            ),
+            :ic => Dict(vars[i] => F(rand(1:50)) for i in 1:2),
+            :param_vals => Dict(vars[i + 2] => F(rand(1:50)) for i in 1:2),
+            :inputs => Dict(vars[5] => [F(rand(-30:30)) for i in 1:4]),
+            :prec => 5,
+        ),
+    )
+
+
 
     for case in test_cases
         ode, prec = case[:ODE], case[:prec]

test/runtests.jl

@@ -65,7 +65,8 @@ using StructuralIdentifiability:
     y_vars,
     x_equations,
     y_equations,
-    inputs
+    inputs,
+    quotient_basis
 
 const GROUP = get(ENV, "GROUP", "All")