add find_conservation_laws

src/ExactODEReduction.jl

@@ -405,6 +405,28 @@ function find_reductions(
     ChainOfReductions(results)
 end
 
+function find_conservation_laws(system::ODE{P}; seed=nothing, loglevel=Logging.Info) where {P}
+    prev_logger = Logging.global_logger(ConsoleLogger(stderr, loglevel))
+
+    isnothing(seed) && (seed = getnewrandomseed())
+    Random.seed!(seed)
+    @info "Global random seed: $seed"
+
+    eqs = equations_extended(system)
+
+    basis, kerdim, ker, relations = constant_linear_relations(eqs)
+    _relations_sym = ker * basis
+    relations_sym = [_relations_sym[i, :] for i in 1:kerdim]
+
+    Logging.global_logger(prev_logger)
+
+    laws = []
+    for i in 1:kerdim
+        push!(laws, (basis=basis, matrix=relations[i], law=relations_sym[i]))
+    end
+    laws
+end
+
 export ODE, @ODEsystem, equations, vars
 export states, parameters, initial_conditions, parameter_values
 export set_parameter_values, to_state, to_parameter

src/linalg_utils.jl

@@ -136,3 +136,24 @@ function construct_jacobians(system)
     ring = parent(first(keys(system)))
     construct_jacobians([system[x] for x in gens(ring)])
 end
+
+#------------------------------------------------------------------------------
+
+# Find the left kernel of a Macaulay matrix of `system`.
+function constant_linear_relations(system::AbstractArray{T}) where {T<:Nemo.MPolyElem}
+    domain = base_ring(first(system))
+    poly_ring = parent(first(system))
+    m = length(system)
+    labels = unique(sort(reduce(vcat, map(collect ∘ monomials, system)), rev=true))
+    n = length(labels)
+    S = Nemo.MatrixSpace(domain, m, n)
+    A = zero(S)
+    for (p_idx, poly) in enumerate(system)
+        for (l_idx, label) in enumerate(labels)
+            A[p_idx, l_idx] = coeff(poly, label)
+        end
+    end
+    kerdim, ker = Nemo.left_kernel(A)
+    relations = [ker[i, :] for i in 1:kerdim]
+    gens(poly_ring), kerdim, ker, relations
+end