From ed11f7572239b127ddac69f27a12bd11b33aca07 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Fri, 16 Feb 2024 01:13:05 +0100 Subject: [PATCH 1/4] remember the lumping matrix --- src/Reduction.jl | 12 ++++++++---- src/parametrization.jl | 2 +- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/Reduction.jl b/src/Reduction.jl index b69e51a7..257389b9 100644 --- a/src/Reduction.jl +++ b/src/Reduction.jl @@ -38,14 +38,18 @@ struct Reduction{A, B, C} new_system::ODE{A} new_vars::Dict{A, B} old_system::ODE{C} + lumping_matrix::Matrix{Any} function Reduction{A, B, C}(new_system::ODE{A}, new_vars::Dict{A, B}, old_system::ODE{C}) where {A, B, C} - return new{A, B, C}(new_system, new_vars, old_system) + lumping_vars = [new_vars[x] for x in gens(parent(old_system))] + lumping_combs = [sum(collect(coefficients(comb)) .* exponent_vectors(comb)) for comb in lumping_vars] + lumping_matrix = mapreduce(permutedims, vcat, lumping_combs) + return new{A, B, C}(new_system, new_vars, old_system, lumping_matrix) end function Reduction{P}(old_ode::ODE{P}, subspace, parameter_strategy=:inhertiance) where {P} - (transformation, new_equations) = perform_change_of_variables(equations_extended(old_ode), subspace) - + (lumping_matrix, transformation, new_equations) = perform_change_of_variables(equations_extended(old_ode), subspace) + # (!) assumes the correct order in new_equations, i.e., # ∂y[1] ~ new_equations[1], # ∂y[2] ~ new_equations[2], @@ -82,7 +86,7 @@ struct Reduction{A, B, C} @warn "Unknown parameter handling strategy $parameter_strategy, using none" end - return new{elem_type(parent(new_ode)), elem_type(parent(first(transformation))), elem_type(parent(old_ode))}(new_ode, new_vars, old_ode) + return new{elem_type(parent(new_ode)), elem_type(parent(first(transformation))), elem_type(parent(old_ode))}(new_ode, new_vars, old_ode, lumping_matrix) end end diff --git a/src/parametrization.jl b/src/parametrization.jl index d025c633..cc58e47c 100644 --- a/src/parametrization.jl +++ b/src/parametrization.jl @@ -55,7 +55,7 @@ function perform_change_of_variables(system, invariants; new_vars_name="y") newsystem = [Nemo.evaluate(p, substitutions) for p in newsystem] - return (transform, newsystem) + return (Matrix(transform_matrix), transform, newsystem) end """ From 23a4390881d05e5a35cca0ecfacbd134565c3662 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Fri, 16 Feb 2024 01:27:33 +0100 Subject: [PATCH 2/4] up --- src/Reduction.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Reduction.jl b/src/Reduction.jl index 257389b9..7228cc7e 100644 --- a/src/Reduction.jl +++ b/src/Reduction.jl @@ -41,7 +41,7 @@ struct Reduction{A, B, C} lumping_matrix::Matrix{Any} function Reduction{A, B, C}(new_system::ODE{A}, new_vars::Dict{A, B}, old_system::ODE{C}) where {A, B, C} - lumping_vars = [new_vars[x] for x in gens(parent(old_system))] + lumping_vars = [new_vars[x] for x in gens(parent(new_system))] lumping_combs = [sum(collect(coefficients(comb)) .* exponent_vectors(comb)) for comb in lumping_vars] lumping_matrix = mapreduce(permutedims, vcat, lumping_combs) return new{A, B, C}(new_system, new_vars, old_system, lumping_matrix) From 632a93e267887182881ae04cd2a5c9edf92bf328 Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Fri, 16 Feb 2024 23:49:36 +0300 Subject: [PATCH 3/4] add find_conservation_laws --- src/ExactODEReduction.jl | 22 ++++++++++++++++++++++ src/linalg_utils.jl | 21 +++++++++++++++++++++ 2 files changed, 43 insertions(+) diff --git a/src/ExactODEReduction.jl b/src/ExactODEReduction.jl index 2f707b39..7c01cd0b 100644 --- a/src/ExactODEReduction.jl +++ b/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 diff --git a/src/linalg_utils.jl b/src/linalg_utils.jl index bbb9ba72..ec4c7443 100644 --- a/src/linalg_utils.jl +++ b/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 From 311f5fbe10a53cf37e4888a27f4e9489ad916acd Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Sat, 17 Feb 2024 00:07:03 +0300 Subject: [PATCH 4/4] fix tests --- src/Reduction.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/Reduction.jl b/src/Reduction.jl index 7228cc7e..e897dcda 100644 --- a/src/Reduction.jl +++ b/src/Reduction.jl @@ -38,12 +38,15 @@ struct Reduction{A, B, C} new_system::ODE{A} new_vars::Dict{A, B} old_system::ODE{C} - lumping_matrix::Matrix{Any} + lumping_matrix::Union{Nothing, Matrix{Any}} function Reduction{A, B, C}(new_system::ODE{A}, new_vars::Dict{A, B}, old_system::ODE{C}) where {A, B, C} - lumping_vars = [new_vars[x] for x in gens(parent(new_system))] - lumping_combs = [sum(collect(coefficients(comb)) .* exponent_vectors(comb)) for comb in lumping_vars] - lumping_matrix = mapreduce(permutedims, vcat, lumping_combs) + lumping_matrix = nothing + if !isempty(new_vars) + lumping_vars = [new_vars[x] for x in gens(parent(new_system))] + lumping_combs = [sum(collect(coefficients(comb)) .* exponent_vectors(comb)) for comb in lumping_vars] + lumping_matrix = mapreduce(permutedims, vcat, lumping_combs) + end return new{A, B, C}(new_system, new_vars, old_system, lumping_matrix) end