diff --git a/src/parametrizations.jl b/src/parametrizations.jl index 0118db5c..6b873ac5 100644 --- a/src/parametrizations.jl +++ b/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( diff --git a/src/power_series_utils.jl b/src/power_series_utils.jl index ac442924..3054e503 100644 --- a/src/power_series_utils.jl +++ b/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) diff --git a/src/primality_check.jl b/src/primality_check.jl index 53d16ed0..9fa4764e 100644 --- a/src/primality_check.jl +++ b/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 = [] diff --git a/test/differentiate_output.jl b/test/differentiate_output.jl index 615ce15a..2630097c 100644 --- a/test/differentiate_output.jl +++ b/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] diff --git a/test/runtests.jl b/test/runtests.jl index 50936512..40324819 100644 --- a/test/runtests.jl +++ b/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")