Merge pull request #374 from SciML/bump_groebner

Bump groebner

Author: pogudingleb

.github/workflows/CI.yml

@@ -15,8 +15,8 @@ jobs:
           - Core
           - ModelingToolkitSIExt
         version:
-          - '<1.10.3 || >=1.10.4'
-          - '1.6'
+          - '1'
+          - '1.10'
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2

Project.toml

@@ -29,31 +29,31 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 ModelingToolkitSIExt = ["ModelingToolkit", "SymbolicUtils", "Symbolics"]
 
 [compat]
-AbstractAlgebra = "0.40, 0.41, 0.42, 0.43"
+AbstractAlgebra = "0.42, 0.43, 0.44"
 Aqua = "0.8"
 CPUSummary = "0.2"
 Combinatorics = "1"
 DataStructures = "0.18"
-Dates = "1.6, 1.7"
-Groebner = "0.7.3"
+Dates = "1.10, 1.11"
+Groebner = "0.8.1"
 IterTools = "1"
-LinearAlgebra = "1.6, 1.7"
-Logging = "1.6, 1.7"
+LinearAlgebra = "1.10, 1.11"
+Logging = "1.10, 1.11"
 MacroTools = "0.5"
 ModelingToolkit = "9.33"
-Nemo = "0.43, 0.44, 0.45, 0.46, 0.47"
-ParamPunPam = "0.4"
-Pkg = "1.6, 1.7"
+Nemo = "0.46, 0.47, 0.48"
+ParamPunPam = "0.5"
+Pkg = "1.10, 1.11"
 PrecompileTools = "1.2"
 Primes = "0.5"
-Random = "1.6, 1.7"
+Random = "1.10, 1.11"
 SpecialFunctions = "2"
-SymbolicUtils = "3.2"
-Symbolics = "6.2"
-Test = "1.6, 1.7"
+SymbolicUtils = "3.7"
+Symbolics = "6.16"
+Test = "1.10, 1.11"
 TestSetExtensions = "2"
 TimerOutputs = "0.5"
-julia = "1.6 - 1.10.2, 1.10.4, 1.11"
+julia = "1.10.4, 1.11.2"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"

benchmarking/IdentifiableFunctions/experiments.jl

@@ -1478,10 +1478,10 @@ mqs_spec = StructuralIdentifiability.ParamPunPam.specialize_mod_p(mqs, point);
 
 Groebner.logging_enabled() = false
 
-@time graph, gb = Groebner.groebner_learn(mqs_spec, loglevel = 0, sweep = true);
-@time Groebner.groebner_apply!(graph, mqs_spec, loglevel = 0, sweep = true);
+@time graph, gb = Groebner.groebner_learn(mqs_spec, sweep = true);
+@time Groebner.groebner_apply!(graph, mqs_spec, sweep = true);
 
-@benchmark Groebner.groebner_apply!($graph, $mqs_spec, loglevel = 0, sweep = true)
+@benchmark Groebner.groebner_apply!($graph, $mqs_spec, sweep = true)
 
 # Results for covid
 #=
@@ -1523,4 +1523,4 @@ BenchmarkTools.Trial: 180 samples with 1 evaluation.
  Memory estimate: 775.77 KiB, allocs estimate: 3662.
 =#
 
-@my_profview_allocs Groebner.groebner_apply!(graph, mqs_spec, loglevel = 0, sweep = true);
+@my_profview_allocs Groebner.groebner_apply!(graph, mqs_spec, sweep = true);

benchmarking/IdentifiableFunctions/homogenization.jl

@@ -31,7 +31,7 @@ ideal_spec = StructuralIdentifiability.specialize_mod_p(mqs, point)
 @time gb = groebner(ideal_spec, ordering = Groebner.DegRevLex());
 
 # There is an existent possibility that this would not finish in two and a half lifetimes
-# @time gb = groebner(ideal_spec, ordering = Groebner.Lex(), loglevel = -3);
+# @time gb = groebner(ideal_spec, ordering = Groebner.Lex());
 
 homogeneous_ideal_spec = StructuralIdentifiability.homogenize(ideal_spec);
 # 100 ms

ext/ModelingToolkitSIExt.jl

@@ -210,9 +210,19 @@ function __mtk_to_si(
     all_funcs = collect(Set(clean_calls(ModelingToolkit.unknowns(de))))
     inputs = filter(s -> !ModelingToolkit.isoutput(s), setdiff(all_funcs, state_vars))
     params = ModelingToolkit.parameters(de)
-    t = ModelingToolkit.arguments(diff_eqs[1].lhs)[1]
+    t = ModelingToolkit.arguments(clean_calls([diff_eqs[1].lhs])[1])[1]
+    # very long if in order to avoid duplication
     params_from_measured_quantities = union(
-        [filter(s -> !iscall(s), get_variables(y[2])) for y in measured_quantities]...,
+        [
+            filter(
+                s ->
+                    !iscall(s) &&
+                        !(string(s) in string.(state_vars)) &&
+                        !(string(s) * "(t)" in string.(state_vars)) &&
+                        (string(s) != string(t)),
+                get_variables(y[2]),
+            ) for y in measured_quantities
+        ]...,
     )
     params = union(params, params_from_measured_quantities)
 

src/RationalFunctionFields/RationalFunctionField.jl

@@ -149,8 +149,8 @@ Output:
     mqs_ratfuncs = specialize(IdealMQS(ratfuncs), point; saturated = false)
     @assert parent(first(mqs_specialized)) == parent(first(mqs_ratfuncs))
     @debug "Starting the groebner basis computation"
-    gb = groebner(mqs_specialized, loglevel = _groebner_loglevel[])
-    result = map(iszero, normalform(gb, mqs_ratfuncs, loglevel = _groebner_loglevel[]))
+    gb = groebner(mqs_specialized)
+    result = map(iszero, normalform(gb, mqs_ratfuncs))
     return result
 end
 
@@ -256,8 +256,8 @@ end
     polys_specialized =
         ParamPunPam.specialize_mod_p(mqs_tobereduced, point, saturated = false)
     @assert parent(first(gens_specialized)) == parent(first(polys_specialized))
-    gb = groebner(gens_specialized, loglevel = _groebner_loglevel[])
-    nf = normalform(gb, polys_specialized, loglevel = _groebner_loglevel[])
+    gb = groebner(gens_specialized)
+    nf = normalform(gb, polys_specialized)
     result = map(iszero, nf)
     return result
 end

src/RationalFunctionFields/normalforms.jl

@@ -56,7 +56,7 @@ function local_normal_forms(
     @assert parent(first(point)) == finite_field
     point_ff_ext = append_at_index(point, mqs.sat_var_index, one(finite_field))
     gens_ff_spec = specialize_mod_p(mqs, point)
-    gb_ff_spec = Groebner.groebner(gens_ff_spec, loglevel = _groebner_loglevel[])
+    gb_ff_spec = Groebner.groebner(gens_ff_spec)
     ring_ff = parent(gb_ff_spec[1])
     xs_ff = gens(ring_ff)
     normal_forms_ff = Vector{elem_type(ring_ff)}(undef, 0)

src/logging.jl

@@ -29,8 +29,6 @@ else
     Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(stderr, Logging.Info))
 end
 
-const _groebner_loglevel = Ref{Int}(0)
-
 function restart_logging(; loglevel = Logging.Info)
     @assert loglevel isa Base.CoreLogging.LogLevel
     _si_logger[] = @static if VERSION >= v"1.7.0"
@@ -41,13 +39,6 @@ function restart_logging(; loglevel = Logging.Info)
     for r in _runtime_rubrics
         _runtime_logger[r] = 0
     end
-    if loglevel < Logging.Info
-        _groebner_loglevel[] = 0
-    elseif loglevel < Logging.Warn
-        _groebner_loglevel[] = 0
-    else
-        _groebner_loglevel[] = 10
-    end
     return nothing
 end
 

src/parametrizations.jl

@@ -218,12 +218,7 @@ $sat_string
     $tagged_mqs
     Monom ordering:
     $(ord)"""
-    tagged_mqs_gb = groebner(
-        tagged_mqs,
-        ordering = ord,
-        homogenize = :no,
-        loglevel = _groebner_loglevel[],
-    )
+    tagged_mqs_gb = groebner(tagged_mqs, ordering = ord, homogenize = :no)
     # Relations between tags in K[T]
     relations_between_tags = filter(
         poly -> isempty(intersect(vars(poly), vcat(sat_var, orig_vars))),

src/power_series_utils.jl

@@ -239,8 +239,7 @@ function ps_ode_solution(
     Svconst = AbstractAlgebra.matrix_space(base_ring(ring), n, 1)
     eqs = Sv(equations)
 
-    x_vars = filter(v -> ("$(v)_dot" in map(string, gens(ring))), gens(ring))
-    x_vars = [x for x in x_vars]
+    x_vars = collect(filter(v -> ("$(v)_dot" in map(var_to_str, gens(ring))), gens(ring)))
     x_dot_vars = [str_to_var(var_to_str(x) * "_dot", ring) for x in x_vars]
 
     Jac_dots = S([derivative(p, xd) for p in equations, xd in x_dot_vars])
@@ -267,7 +266,7 @@ function ps_ode_solution(
             set_precision!(solution[x_vars[i]], new_prec)
             set_precision!(solution[x_dot_vars[i]], new_prec)
         end
-        eval_point = [solution[v] for v in gens(ring)]
+        eval_point = [copy(solution[v]) for v in gens(ring)]
         map(ps -> set_precision!(ps, 2 * cur_prec), eval_point)
         eqs_eval = map(p -> evaluate(p, eval_point), eqs)
         J_eval = map(p -> evaluate(p, eval_point), Jac_xs)

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{<:MPolyRingElem, 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, loglevel = _groebner_loglevel[])
-    basis = Groebner.kbase(J, loglevel = _groebner_loglevel[])
+    J = Groebner.groebner(J)
+    basis = quotient_basis(J)
     dim = length(basis)
     S = Nemo.matrix_space(Nemo.QQ, dim, dim)
     matrices = []
@@ -11,7 +47,7 @@ function check_primality_zerodim(J::Array{QQMPolyRingElem, 1})
     for v in gens(parent(first(J)))
         M = zero(S)
         for (i, vec) in enumerate(basis)
-            image = Groebner.normalform(J, v * vec, loglevel = _groebner_loglevel[])
+            image = Groebner.normalform(J, v * vec)
             for (j, base_vec) in enumerate(basis)
                 M[i, j] = Nemo.QQ(coeff(image, base_vec))
             end

src/wronskian.jl

@@ -237,7 +237,7 @@ Computes the Wronskians of io_equations
 
     @debug "Constructing Wronskians"
     result = []
-    for (i, tlist) in enumerate(termlists)
+    for tlist in termlists
         n = length(tlist)
         evaled = massive_eval(tlist, ps_ext)
         S = Nemo.matrix_space(F, n, n)

test/differentiate_output.jl

@@ -64,7 +64,7 @@ function diff_sol_Lie_derivatives(ode::ODE, params, ic, inputs, prec::Int)
                 push!(
                     result[y][v],
                     eval_at_dict(
-                        derivative(Lie_derivatives[y][j], str_to_var("$v", new_ring)),
+                        derivative(Lie_derivatives[y][j], switch_ring(v, new_ring)),
                         eval_point,
                     ),
                 )
@@ -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,
@@ -161,6 +212,7 @@ end
         ),
     )
 
+    t = copy(test_cases)
     varnames = vcat(
         ["x_$i" for i in 1:3],
         ["p_$i" for i in 1:3],
@@ -201,8 +253,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:5]),
+            :prec => 5,
+        ),
+    )
 
-    for case in test_cases
+    for case in t
         ode, prec = case[:ODE], case[:prec]
         @time sol1 =
             differentiate_output(ode, case[:param_vals], case[:ic], case[:inputs], 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")