From 2daef114dffeccd9d279fdbd5c4801de785115d8 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Thu, 18 Apr 2024 07:22:55 +0200 Subject: [PATCH 1/7] Update versions, fix issues, remove Oscar and msolve --- src/ParameterEstimation.jl | 7 +- src/estimation/fixed_degree.jl | 14 ++-- src/estimation/solve.jl | 67 ++++++++--------- src/filtering.jl | 6 +- src/identifiability/check_identifiability.jl | 62 ++++++++-------- src/identifiability/identifiability_data.jl | 34 ++++----- src/identifiability/preprocessor.jl | 61 ++++++++-------- src/identifiability/transcendence_basis.jl | 22 +++--- src/identifiability/utils.jl | 77 ++++++++++---------- src/rational_interpolation/utils.jl | 10 +-- src/result.jl | 2 +- src/utils.jl | 6 +- test/msolve_estimation_tests.jl | 53 +++++++------- test/nemo_hc_conversion.jl | 14 ++-- test/runtests.jl | 2 +- test/test_alg_indep.jl | 8 +- 16 files changed, 226 insertions(+), 219 deletions(-) diff --git a/src/ParameterEstimation.jl b/src/ParameterEstimation.jl index 0034d649..292fd484 100644 --- a/src/ParameterEstimation.jl +++ b/src/ParameterEstimation.jl @@ -15,9 +15,10 @@ using Suppressor using ProgressMeter, Logging, Printf using ModelingToolkit, LinearSolve, LinearAlgebra -using SIAN, HomotopyContinuation, Groebner, Oscar +using SIAN, HomotopyContinuation, Groebner, Nemo using .ReturnCode -import StructuralIdentifiability: eval_at_nemo, ODE +import StructuralIdentifiability +import StructuralIdentifiability: ODE using BaryRational using ForwardDiff using ArbNumerics @@ -53,10 +54,10 @@ export check_identifiability, estimate, filter_solutions p_true = [2.0, 3.0, 4.0, 5.0] time_interval = [-4.0, 4.0] datasize = 9 + model = complete(model) data_sample = ParameterEstimation.sample_data(model, measured_quantities, time_interval, p_true, ic, datasize;) res = ParameterEstimation.estimate(model, measured_quantities, data_sample) end end end - diff --git a/src/estimation/fixed_degree.jl b/src/estimation/fixed_degree.jl index d1e1e9ee..e0b3f695 100644 --- a/src/estimation/fixed_degree.jl +++ b/src/estimation/fixed_degree.jl @@ -2,7 +2,7 @@ backsolve_initial_conditions(model, E, report_time, inputs::Vector{Equation}, data_sample; solver = Vern9(), abstol = 1e-14, reltol = 1e-14) - initial_conditions = [E[s] for s in ModelingToolkit.states(model)] + initial_conditions = [E[s] for s in ModelingToolkit.unknowns(model)] parameter_values = [E[p] for p in ModelingToolkit.parameters(model)] tspan = (E.at_time, report_time) @@ -10,7 +10,7 @@ """ function backsolve_initial_conditions(model, E, report_time, inputs::Vector{Equation}, data_sample; solver = Vern9(), abstol = 1e-14, reltol = 1e-14) - initial_conditions = [E[s] for s in ModelingToolkit.states(model)] + initial_conditions = [E[s] for s in ModelingToolkit.unknowns(model)] parameter_values = [E[p] for p in ModelingToolkit.parameters(model)] tspan = (E.at_time, report_time) #this is almost always backwards! @@ -18,20 +18,20 @@ function backsolve_initial_conditions(model, E, report_time, inputs::Vector{Equa ode_equations = substitute(ode_equations, Dict(each.lhs => Num(each.rhs) for each in inputs)) t = ModelingToolkit.get_iv(model) - @named new_model = ODESystem(ode_equations, t, ModelingToolkit.states(model), + @named new_model = ODESystem(ode_equations, t, ModelingToolkit.unknowns(model), ModelingToolkit.parameters(model)) - prob = ODEProblem(new_model, initial_conditions, tspan, parameter_values) + prob = ModelingToolkit.complete(ODEProblem(new_model, initial_conditions, tspan, parameter_values)) saveat = range(tspan[1], tspan[2], length = length(data_sample["t"])) ode_solution = ModelingToolkit.solve(prob, solver, p = parameter_values, saveat = saveat, abstol = abstol, reltol = reltol) state_param_map = (Dict(x => replace(string(x), "(t)" => "") - for x in ModelingToolkit.states(model))) + for x in ModelingToolkit.unknowns(model))) newstates = copy(E.states) - for s in ModelingToolkit.states(model) + for s in ModelingToolkit.unknowns(model) temp = ode_solution[Symbol(state_param_map[s])][end] newstates[s] = temp end @@ -94,7 +94,7 @@ function estimate_single_interpolator(model::ModelingToolkit.ODESystem, check_inputs(measured_quantities, data_sample) #TODO(orebas): I took out checking the degree. Do we want to check the interpolator otherwise? datasize = length(first(values(data_sample))) parameters = ModelingToolkit.parameters(model) - states = ModelingToolkit.states(model) + states = ModelingToolkit.unknowns(model) num_parameters = length(parameters) + length(states) @debug "Interpolating sample data via interpolation method $(interpolator.first)" if !haskey(data_sample, "t") diff --git a/src/estimation/solve.jl b/src/estimation/solve.jl index ab40a79a..d6aabd26 100644 --- a/src/estimation/solve.jl +++ b/src/estimation/solve.jl @@ -3,7 +3,7 @@ function solve_via_homotopy(identifiability_result, model; real_tol = 1e-12) polynomial_system = identifiability_result["polynomial_system_to_solve"] state_param_map = merge(Dict(replace(string(x), "(t)" => "") => x - for x in ModelingToolkit.states(model)), + for x in ModelingToolkit.unknowns(model)), Dict(string(x) => x for x in ModelingToolkit.parameters(model))) results = HomotopyContinuation.solve(polynomial_system; show_progress = false) all_solutions = HomotopyContinuation.real_solutions(results) @@ -34,38 +34,39 @@ function solve_via_homotopy(identifiability_result, model; real_tol = 1e-12) end function solve_via_msolve(identifiability_result, model; real_tol = 1e-12) - polynomial_system = identifiability_result["polynomial_system_to_solve"] - state_param_map = merge(Dict(replace(string(x), "(t)" => "") => x - for x in ModelingToolkit.states(model)), - Dict(string(x) => x for x in ModelingToolkit.parameters(model))) + throw("Not implemented, sorry.") + # polynomial_system = identifiability_result["polynomial_system_to_solve"] + # state_param_map = merge(Dict(replace(string(x), "(t)" => "") => x + # for x in ModelingToolkit.states(model)), + # Dict(string(x) => x for x in ModelingToolkit.parameters(model))) - all_vars = reduce((x, y) -> union(x, y), vars(poly) for poly in polynomial_system) - R, _ = Oscar.PolynomialRing(QQ, string.(all_vars); ordering = :degrevlex) - ps = [SIAN.parent_ring_change(poly, R) for poly in polynomial_system] - i = Oscar.ideal(R, ps) - all_solutions_ = Vector{Dict}([]) + # all_vars = reduce((x, y) -> union(x, y), vars(poly) for poly in polynomial_system) + # R, _ = Nemo.polynomial_ring(QQ, string.(all_vars); ordering = :degrevlex) + # ps = [SIAN.parent_ring_change(poly, R) for poly in polynomial_system] + # i = Oscar.ideal(R, ps) + # all_solutions_ = Vector{Dict}([]) - if Oscar.VERSION_NUMBER == v"0.11.3" - solutions, rat_param = Oscar.real_solutions(i) - elseif Oscar.VERSION_NUMBER == v"0.10.0" - rat_param, solutions = Oscar.msolve(i) - else - solutions, rat_param = Oscar.real_solutions(i) - end - for sol in solutions - tmp = Dict() - sol = map(each -> Float64(each), sol) - for (idx, v) in enumerate(all_vars) - if endswith(string(v), "_0") - tmp[state_param_map[string(v)[1:(end-2)]]] = sol[idx] - end - end - for (key, val) in identifiability_result.transcendence_basis_subs - if endswith(string(key), "_0") - tmp[state_param_map[string(key)[1:(end-2)]]] = Int(Meta.parse("$val")) - end - end - push!(all_solutions_, tmp) - end - return all_solutions_ + # if Oscar.VERSION_NUMBER == v"0.11.3" + # solutions, rat_param = Oscar.real_solutions(i) + # elseif Oscar.VERSION_NUMBER == v"0.10.0" + # rat_param, solutions = Oscar.msolve(i) + # else + # solutions, rat_param = Oscar.real_solutions(i) + # end + # for sol in solutions + # tmp = Dict() + # sol = map(each -> Float64(each), sol) + # for (idx, v) in enumerate(all_vars) + # if endswith(string(v), "_0") + # tmp[state_param_map[string(v)[1:(end-2)]]] = sol[idx] + # end + # end + # for (key, val) in identifiability_result.transcendence_basis_subs + # if endswith(string(key), "_0") + # tmp[state_param_map[string(key)[1:(end-2)]]] = Int(Meta.parse("$val")) + # end + # end + # push!(all_solutions_, tmp) + # end + # return all_solutions_ end diff --git a/src/filtering.jl b/src/filtering.jl index 2fdbd58d..dfa4e45d 100644 --- a/src/filtering.jl +++ b/src/filtering.jl @@ -20,16 +20,16 @@ Compute the error between the solution and the data sample. The error is recorde """ function solve_ode(model, estimate::EstimationResult, inputs::Vector{Equation}, data_sample; solver = Tsit5(), return_ode = false, abstol = 1e-12, reltol = 1e-12) - initial_conditions = [estimate[s] for s in ModelingToolkit.states(model)] + initial_conditions = [estimate[s] for s in ModelingToolkit.unknowns(model)] parameter_values = [estimate[p] for p in ModelingToolkit.parameters(model)] tspan = (estimate.at_time, data_sample["t"][end] + estimate.at_time) ode_equations = ModelingToolkit.equations(model) ode_equations = substitute(ode_equations, Dict(each.lhs => Num(each.rhs) for each in inputs)) t = ModelingToolkit.get_iv(model) - @named new_model = ODESystem(ode_equations, t, ModelingToolkit.states(model), + @named new_model = ODESystem(ode_equations, t, ModelingToolkit.unknowns(model), ModelingToolkit.parameters(model)) - prob = ODEProblem(new_model, initial_conditions, tspan, parameter_values) + prob = ODEProblem(ModelingToolkit.complete(new_model), initial_conditions, tspan, parameter_values) ode_solution = ModelingToolkit.solve(prob, solver, saveat = range(tspan[1], tspan[2], length = length(data_sample["t"])), abstol = 1e-12, reltol = 1e-12) diff --git a/src/identifiability/check_identifiability.jl b/src/identifiability/check_identifiability.jl index f67b86bc..c0b74fa1 100644 --- a/src/identifiability/check_identifiability.jl +++ b/src/identifiability/check_identifiability.jl @@ -38,7 +38,7 @@ object that contains the results of the identifiability analysis. end end ode_prep, input_syms, gens_ = preprocess_ode(ode, measured_quantities, inputs) - t = ModelingToolkit.arguments(ModelingToolkit.states(ode)[1])[1] + t = ModelingToolkit.arguments(ModelingToolkit.unknowns(ode)[1])[1] params_to_assess_ = SIAN.get_parameters(ode_prep) nemo2mtk = Dict(gens_ .=> input_syms) @@ -113,8 +113,8 @@ function identifiability_ode(ode, params_to_assess; p = 0.99, p_mod = 0, infolev @info "Truncating" # (a) ----------------------- - d0 = BigInt(maximum(vcat([SIAN.Nemo.total_degree(SIAN.unpack_fraction(Q * eq[2])[1]) - for eq in eqs], SIAN.Nemo.total_degree(Q)))) + d0 = BigInt(maximum(vcat([Nemo.total_degree(SIAN.unpack_fraction(Q * eq[2])[1]) + for eq in eqs], Nemo.total_degree(Q)))) # (b) ----------------------- D1 = floor(BigInt, @@ -130,19 +130,19 @@ function identifiability_ode(ode, params_to_assess; p = 0.99, p_mod = 0, infolev # (e) ------------------ alpha = [1 for i in 1:n] beta = [0 for i in 1:m] - Et = Array{SIAN.Nemo.fmpq_mpoly}(undef, 0) + Et = Array{Nemo.QQMPolyRingElem}(undef, 0) x_theta_vars = all_params prolongation_possible = [1 for i in 1:m] # (f) ------------------ all_x_theta_vars_subs = SIAN.insert_zeros_to_vals(all_subs[1], all_subs[2]) - eqs_i_old = Array{SIAN.Nemo.fmpq_mpoly}(undef, 0) - evl_old = Array{SIAN.Nemo.fmpq_mpoly}(undef, 0) + eqs_i_old = Array{Nemo.QQMPolyRingElem}(undef, 0) + evl_old = Array{Nemo.QQMPolyRingElem}(undef, 0) while sum(prolongation_possible) > 0 for i in 1:m if prolongation_possible[i] == 1 eqs_i = vcat(Et, Y[i][beta[i] + 1]) - evl = [SIAN.Nemo.evaluate(eq, vcat(u_hat[1], y_hat[1]), + evl = [Nemo.evaluate(eq, vcat(u_hat[1], y_hat[1]), vcat(u_hat[2], y_hat[2])) for eq in eqs_i if !(eq in eqs_i_old)] evl_old = vcat(evl_old, evl) @@ -154,20 +154,20 @@ function identifiability_ode(ode, params_to_assess; p = 0.99, p_mod = 0, infolev # adding necessary X-equations polys_to_process = vcat(Et, [Y[k][beta[k] + 1] for k in 1:m]) while length(polys_to_process) != 0 - new_to_process = Array{SIAN.Nemo.fmpq_mpoly}(undef, 0) - vrs = Set{SIAN.Nemo.fmpq_mpoly}() + new_to_process = Array{Nemo.QQMPolyRingElem}(undef, 0) + vrs = Set{Nemo.QQMPolyRingElem}() for poly in polys_to_process vrs = union(vrs, [v - for v in SIAN.Nemo.vars(poly) if v in x_variables]) + for v in Nemo.vars(poly) if v in x_variables]) end - vars_to_add = Set{SIAN.Nemo.fmpq_mpoly}(v + vars_to_add = Set{Nemo.QQMPolyRingElem}(v for v in vrs if !(v in x_theta_vars)) for v in vars_to_add x_theta_vars = vcat(x_theta_vars, v) ord_var = SIAN.get_order_var2(v, all_indets, n + m + u, s) - var_idx = SIAN.Nemo.var_index(ord_var[1]) + var_idx = Nemo.var_index(ord_var[1]) poly = X[var_idx][ord_var[2]] Et = vcat(Et, poly) new_to_process = vcat(new_to_process, poly) @@ -200,15 +200,15 @@ function identifiability_ode(ode, params_to_assess; p = 0.99, p_mod = 0, infolev end end - theta_l = Array{SIAN.Nemo.fmpq_mpoly}(undef, 0) + theta_l = Array{Nemo.QQMPolyRingElem}(undef, 0) params_to_assess_ = [SIAN.add_to_var(param, Rjet, 0) for param in params_to_assess] - Et_eval_base = [SIAN.Nemo.evaluate(e, vcat(u_hat[1], y_hat[1]), + Et_eval_base = [Nemo.evaluate(e, vcat(u_hat[1], y_hat[1]), vcat(u_hat[2], y_hat[2])) for e in Et] for param_0 in params_to_assess_ other_params = [v for v in x_theta_vars if v != param_0] - Et_subs = [SIAN.Nemo.evaluate(e, [param_0], - [SIAN.Nemo.evaluate(param_0, all_x_theta_vars_subs)]) + Et_subs = [Nemo.evaluate(e, [param_0], + [Nemo.evaluate(param_0, all_x_theta_vars_subs)]) for e in Et_eval_base] JacX = SIAN.jacobi_matrix(Et_subs, other_params, all_x_theta_vars_subs) if LinearAlgebra.rank(JacX) != max_rank @@ -251,7 +251,7 @@ function identifiability_ode(ode, params_to_assess; p = 0.99, p_mod = 0, infolev # 3. Randomize. @info "Randomizing" # (a) ------------ - deg_variety = foldl(*, [BigInt(SIAN.Nemo.total_degree(e)) for e in Et]) + deg_variety = foldl(*, [BigInt(Nemo.total_degree(e)) for e in Et]) D2 = floor(BigInt, 3 / 4 * 6 * length(theta_l) * deg_variety * (1 + 2 * d0 * maximum(beta)) / @@ -264,22 +264,22 @@ function identifiability_ode(ode, params_to_assess; p = 0.99, p_mod = 0, infolev theta_hat = sample[3] # (d) ------------ - Et_hat = [SIAN.Nemo.evaluate(e, vcat(y_hat[1], u_hat[1]), vcat(y_hat[2], u_hat[2])) + Et_hat = [Nemo.evaluate(e, vcat(y_hat[1], u_hat[1]), vcat(y_hat[2], u_hat[2])) for e in Et] - transcendence_substitutions = Array{SIAN.Nemo.fmpq}(undef, 0) + transcendence_substitutions = Array{Nemo.QQFieldElem}(undef, 0) for (idx, var) in enumerate(theta_hat[1]) if var in alg_indep push!(transcendence_substitutions, theta_hat[2][idx]) end end - Et_hat = [SIAN.Nemo.evaluate(e, alg_indep, transcendence_substitutions) + Et_hat = [Nemo.evaluate(e, alg_indep, transcendence_substitutions) for e in Et_hat] - Et_x_vars = Set{SIAN.Nemo.fmpq_mpoly}() + Et_x_vars = Set{Nemo.QQMPolyRingElem}() for poly in Et_hat - Et_x_vars = union(Et_x_vars, Set(SIAN.Nemo.vars(poly))) + Et_x_vars = union(Et_x_vars, Set(Nemo.vars(poly))) end Et_x_vars = setdiff(Et_x_vars, not_int_cond_params) - Q_hat = SIAN.Nemo.evaluate(Q, u_hat[1], u_hat[2]) + Q_hat = Nemo.evaluate(Q, u_hat[1], u_hat[2]) vrs_sorted = vcat(sort([e for e in Et_x_vars], lt = (x, y) -> SIAN.compare_diff_var(x, y, all_indets, n + m + u, s)), z_aux, @@ -288,14 +288,14 @@ function identifiability_ode(ode, params_to_assess; p = 0.99, p_mod = 0, infolev # assign weights to variables if weighted_ordering for i in eachindex(Et_hat) - for _var in Set(SIAN.Nemo.vars(Et_hat[i])) + for _var in Set(Nemo.vars(Et_hat[i])) _var_non_jet, _var_order = SIAN.get_order_var(_var, non_jet_ring) Et_hat[i] = SIAN.make_substitution(Et_hat[i], _var, _var^get(weights, _var_non_jet, 1), parent(_var)(1)) end end - for _var in Set(SIAN.Nemo.vars(Q_hat)) + for _var in Set(Nemo.vars(Q_hat)) _var_non_jet, _var_order = SIAN.get_order_var(_var, non_jet_ring) Q_hat = SIAN.make_substitution(Q_hat, _var, _var^get(weights, _var_non_jet, 1), @@ -309,13 +309,13 @@ function identifiability_ode(ode, params_to_assess; p = 0.99, p_mod = 0, infolev Et_hat = [SIAN._reduce_poly_mod_p(e, p_mod) for e in Et_hat] z_aux = SIAN._reduce_poly_mod_p(z_aux, p_mod) Q_hat = SIAN._reduce_poly_mod_p(Q_hat, p_mod) - Rjet_new, vrs_sorted = SIAN.Nemo.PolynomialRing(SIAN.Nemo.GF(p_mod), + Rjet_new, vrs_sorted = Nemo.polynomial_ring(Nemo.Native.GF(p_mod), [string(v) for v in vrs_sorted], - ordering = :degrevlex) + internal_ordering = :degrevlex) else - Rjet_new, vrs_sorted = SIAN.Nemo.PolynomialRing(SIAN.Nemo.QQ, + Rjet_new, vrs_sorted = Nemo.polynomial_ring(Nemo.QQ, [string(v) for v in vrs_sorted], - ordering = :degrevlex) + internal_ordering = :degrevlex) end Et_hat = [SIAN.parent_ring_change(e, Rjet_new) for e in Et_hat] @@ -382,10 +382,10 @@ function identifiability_ode(ode, params_to_assess; p = 0.99, p_mod = 0, infolev u_derivative_dict[each] = order end Et_ = [Et[idx] for idx in Et_ids] - full_result = Dict("full_polynomial_system" => [SIAN.Nemo.evaluate(e, alg_indep, + full_result = Dict("full_polynomial_system" => [Nemo.evaluate(e, alg_indep, transcendence_substitutions) for e in Et], - "polynomial_system" => [SIAN.Nemo.evaluate(e, alg_indep, + "polynomial_system" => [Nemo.evaluate(e, alg_indep, transcendence_substitutions) for e in Et_], "polynomial_system_to_solve" => HomotopyContinuation.System([]), diff --git a/src/identifiability/identifiability_data.jl b/src/identifiability/identifiability_data.jl index d4eca3a2..6d745f6b 100644 --- a/src/identifiability/identifiability_data.jl +++ b/src/identifiability/identifiability_data.jl @@ -1,4 +1,4 @@ -PolySystem = Union{HomotopyContinuation.ModelKit.System, Vector{SIAN.Nemo.fmpq_mpoly}} +PolySystem = Union{HomotopyContinuation.ModelKit.System, Vector{Nemo.QQMPolyRingElem}} """ IdentifiabilityData @@ -7,36 +7,36 @@ A struct that contains the data from identifiability analysis. This is used for parameter estimation. # Fields -- `polynomial_system::Vector{SIAN.Nemo.fmpq_mpoly}`: The polynomial system. +- `polynomial_system::Vector{Nemo.QQMPolyRingElem}`: The polynomial system. - `polynomial_system_to_solve::PolySystem`: The polynomial system with derivatives substitutited and ready to be solved. -- `denominator::SIAN.Nemo.fmpq_mpoly`: The denominator of the polynomial system. -- `variables::Vector{SIAN.Nemo.fmpq_mpoly}`: The variables of the polynomial system. +- `denominator::Nemo.QQMPolyRingElem`: The denominator of the polynomial system. +- `variables::Vector{Nemo.QQMPolyRingElem}`: The variables of the polynomial system. - `substitutions::Vector{Vector}`: The substitutions used to assess identifiability. - `identifiability_nemo::Any`: The identifiability data from SIAN in Nemo data type. - `identifiability::Dict`: The identifiability data from SIAN in HomotopyContinuation compatible data type. -- `basis::Vector{SIAN.Nemo.fmpq_mpoly}`: The transcendence basis of the polynomial system. -- `transcendence_basis_subs::Vector{SIAN.Nemo.RingElem}`: The transcendence basis substitutions of the polynomial system. -- `weights::Dict{SIAN.Nemo.fmpq_mpoly, Int64}`: The weights of the variables used by SIAN to assess GroebnerBasis. +- `basis::Vector{Nemo.QQMPolyRingElem}`: The transcendence basis of the polynomial system. +- `transcendence_basis_subs::Vector{Nemo.RingElem}`: The transcendence basis substitutions of the polynomial system. +- `weights::Dict{Nemo.QQMPolyRingElem, Int64}`: The weights of the variables used by SIAN to assess GroebnerBasis. """ mutable struct IdentifiabilityData - full_polynomial_system::Vector{SIAN.Nemo.fmpq_mpoly} - polynomial_system::Vector{SIAN.Nemo.fmpq_mpoly} + full_polynomial_system::Vector{Nemo.QQMPolyRingElem} + polynomial_system::Vector{Nemo.QQMPolyRingElem} polynomial_system_to_solve::PolySystem - denomiantor::SIAN.Nemo.fmpq_mpoly - variables::Vector{SIAN.Nemo.fmpq_mpoly} + denomiantor::Nemo.QQMPolyRingElem + variables::Vector{Nemo.QQMPolyRingElem} substitutions::Vector{Vector} identifiability_nemo::Any identifiability::AbstractDict transcendence_basis_subs::AbstractDict Y_eq::AbstractDict u_variables::AbstractDict - basis::Vector{SIAN.Nemo.fmpq_mpoly} - weights::AbstractDict{SIAN.Nemo.fmpq_mpoly, Int64} - non_jet_ring::SIAN.Nemo.FmpqMPolyRing + basis::Vector{Nemo.QQMPolyRingElem} + weights::AbstractDict{Nemo.QQMPolyRingElem, Int64} + non_jet_ring::Nemo.QQMPolyRing nemo_mtk::AbstractDict - solution_counts::AbstractDict + solution_counts::AbstractDict # TODO: temporarily not meaningful function IdentifiabilityData(input::AbstractDict) - solution_counts = count_solutions(input) + # solution_counts = count_solutions(input) return new(input["full_polynomial_system"], input["polynomial_system"], input["polynomial_system_to_solve"], input["denominator"], input["vars"], @@ -44,6 +44,6 @@ mutable struct IdentifiabilityData input["transcendence_basis_subs"], input["Y_eq"], input["u_variables"], input["basis"], input["weights"], input["non_jet_ring"], input["nemo_mtk"], - solution_counts) + Dict())#solution_counts) end end diff --git a/src/identifiability/preprocessor.jl b/src/identifiability/preprocessor.jl index 389ab40f..484347ed 100644 --- a/src/identifiability/preprocessor.jl +++ b/src/identifiability/preprocessor.jl @@ -9,57 +9,60 @@ Output: - `ODE` object containing required data for identifiability assessment """ function preprocess_ode(de::ModelingToolkit.ODESystem, - measured_quantities::Array{ModelingToolkit.Equation}, - inputs = Vector{Num}()) + measured_quantities::Array{ModelingToolkit.Equation}, + inputs = Vector{Num}()) @info "Preproccessing `ModelingToolkit.ODESystem` object" diff_eqs = filter(eq -> !(ModelingToolkit.isoutput(eq.lhs)), - ModelingToolkit.equations(de)) + ModelingToolkit.equations(de)) y_functions = [each.lhs for each in measured_quantities] state_vars = filter(s -> !(ModelingToolkit.isinput(s) || ModelingToolkit.isoutput(s)), - ModelingToolkit.states(de)) + ModelingToolkit.unknowns(de)) params = ModelingToolkit.parameters(de) t = ModelingToolkit.arguments(measured_quantities[1].lhs)[1] - params_from_measured_quantities = ModelingToolkit.parameters(ModelingToolkit.ODESystem(measured_quantities, - t, - name = :DataSeries)) + params_from_measured_quantities = ModelingToolkit.parameters(ModelingToolkit.ODESystem( + measured_quantities, + t, + name = :DataSeries)) params = union(params, params_from_measured_quantities) input_symbols = vcat(state_vars, y_functions, inputs, params) generators = string.(input_symbols) generators = map(g -> replace(g, "(t)" => ""), generators) - R, gens_ = Nemo.PolynomialRing(Nemo.QQ, generators) - state_eqn_dict = Dict{Oscar.fmpq_mpoly, - Union{Oscar.fmpq_mpoly, - Oscar.Generic.Frac{ - Oscar.fmpq_mpoly - }}}() - out_eqn_dict = Dict{Oscar.fmpq_mpoly, - Union{Oscar.fmpq_mpoly, - Oscar.Generic.Frac{ - Oscar.fmpq_mpoly - }}}() + R, gens_ = Nemo.polynomial_ring(Nemo.QQ, generators) + state_eqn_dict = Dict{Nemo.QQMPolyRingElem, + Union{Nemo.QQMPolyRingElem, + Nemo.Generic.FracFieldElem{ + Nemo.QQMPolyRingElem + }}}() + out_eqn_dict = Dict{Nemo.QQMPolyRingElem, + Union{Nemo.QQMPolyRingElem, + Nemo.Generic.FracFieldElem{ + Nemo.QQMPolyRingElem + }}}() for i in eachindex(diff_eqs) if !(typeof(diff_eqs[i].rhs) <: Number) - state_eqn_dict[substitute(state_vars[i], input_symbols .=> gens_)] = eval_at_nemo(diff_eqs[i].rhs, - Dict(input_symbols .=> - gens_)) + state_eqn_dict[substitute(state_vars[i], input_symbols .=> gens_)] = StructuralIdentifiability.eval_at_nemo( + diff_eqs[i].rhs, + Dict(input_symbols .=> + gens_)) else state_eqn_dict[substitute(state_vars[i], input_symbols .=> gens_)] = R(diff_eqs[i].rhs) end end for i in 1:length(measured_quantities) - out_eqn_dict[substitute(y_functions[i], input_symbols .=> gens_)] = eval_at_nemo(measured_quantities[i].rhs, - Dict(input_symbols .=> - gens_)) + out_eqn_dict[substitute(y_functions[i], input_symbols .=> gens_)] = StructuralIdentifiability.eval_at_nemo( + measured_quantities[i].rhs, + Dict(input_symbols .=> + gens_)) end inputs_ = [substitute(each, input_symbols .=> gens_) for each in inputs] if isequal(length(inputs_), 0) - inputs_ = Vector{Oscar.fmpq_mpoly}() + inputs_ = Vector{Nemo.QQMPolyRingElem}() end - return (ODE{Oscar.fmpq_mpoly}(state_eqn_dict, - out_eqn_dict, - inputs_), - input_symbols, gens_) + return (ODE{Nemo.QQMPolyRingElem}(state_eqn_dict, + out_eqn_dict, + inputs_), + input_symbols, gens_) end diff --git a/src/identifiability/transcendence_basis.jl b/src/identifiability/transcendence_basis.jl index d45cadd0..e7ad6992 100644 --- a/src/identifiability/transcendence_basis.jl +++ b/src/identifiability/transcendence_basis.jl @@ -1,22 +1,22 @@ """ - algebraic_independence(Et::Vector{SIAN.Nemo.fmpq_mpoly}, - indets::Vector{SIAN.Nemo.fmpq_mpoly}, + algebraic_independence(Et::Vector{Nemo.QQMPolyRingElem}, + indets::Vector{Nemo.QQMPolyRingElem}, vals) Returns the indices of the equations in Et to be used for polynomial solving and the variables that form a transcendence basis. # Arguments -- `Et::Vector{SIAN.Nemo.fmpq_mpoly}`: The equations to be solved (must come from identifiability check). -- `indets::Vector{SIAN.Nemo.fmpq_mpoly}`: The indeterminates. -- `vals::Vector{SIAN.Nemo.fmpq_mpoly}`: The values of the indeterminates sampled by identifiability algorithm. +- `Et::Vector{Nemo.QQMPolyRingElem}`: The equations to be solved (must come from identifiability check). +- `indets::Vector{Nemo.QQMPolyRingElem}`: The indeterminates. +- `vals::Vector{Nemo.QQMPolyRingElem}`: The values of the indeterminates sampled by identifiability algorithm. """ -function algebraic_independence(Et::Vector{SIAN.Nemo.fmpq_mpoly}, - indets::Vector{SIAN.Nemo.fmpq_mpoly}, - vals) - pivots = Vector{SIAN.Nemo.fmpq_mpoly}() +function algebraic_independence(Et::Vector{Nemo.QQMPolyRingElem}, + indets::Vector{Nemo.QQMPolyRingElem}, + vals) + pivots = Vector{Nemo.QQMPolyRingElem}() Jacobian = SIAN.jacobi_matrix(Et, indets, vals) - U = SIAN.Nemo.lu(Jacobian)[end] + U = Nemo.lu(Jacobian)[end] #find pivot columns in u for row_idx in 1:size(U, 1) row = U[row_idx, :] @@ -31,7 +31,7 @@ function algebraic_independence(Et::Vector{SIAN.Nemo.fmpq_mpoly}, output_ids = [1] for current_idx in 2:length(Et) current = [output_rows; Jacobian[current_idx, :]] - if SIAN.Nemo.rank(current) > current_rank + if Nemo.rank(current) > current_rank output_rows = current push!(output_ids, current_idx) current_rank += 1 diff --git a/src/identifiability/utils.jl b/src/identifiability/utils.jl index 085ba14a..72e6bea7 100644 --- a/src/identifiability/utils.jl +++ b/src/identifiability/utils.jl @@ -15,12 +15,12 @@ function Base.getindex(identifiability_result::IdentifiabilityData, key::String) end function Base.setindex!(identifiability_result::IdentifiabilityData, - value::PolySystem, key::String) + value::PolySystem, key::String) return setproperty!(identifiability_result, Symbol(key), value) end function Base.setindex!(identifiability_result::IdentifiabilityData, - value::PolySystem, key::Symbol) + value::PolySystem, key::Symbol) return setproperty!(identifiability_result, key, value) end @@ -36,41 +36,42 @@ Non-identifable parameters are assigned 0. Globally identifable parameters are a - `identifiability_result::IdentifiabilityData`: Identifiability result from `identifiability_analysis`. """ function count_solutions(identifiability_result) - globally_id = identifiability_result["identifiability_nemo"]["globally"] - locally_not_globally_id = identifiability_result["identifiability_nemo"]["locally_not_globally"] - non_id = identifiability_result["identifiability_nemo"]["nonidentifiable"] - weights_table = identifiability_result["weights"] - non_jet_ring = identifiability_result["non_jet_ring"] - n2m = identifiability_result["nemo_mtk"] - solutions_table = Dict{Any, Int}() - for param in globally_id - v = SIAN.get_order_var(param, non_jet_ring)[1] - solutions_table[n2m[v]] = 1 - end - for param in non_id - v = SIAN.get_order_var(param, non_jet_ring)[1] - solutions_table[n2m[v]] = 0 - end - # basis_lex = Groebner.fglm(identifiability_result.basis) - R = parent(identifiability_result["basis"][1]) - OR, svars = Oscar.PolynomialRing(Oscar.QQ, string.(gens(R)); ordering = :lex) - basis_sing = [SIAN.parent_ring_change(identifiability_result["basis"][i], OR) - for i in 1:length(identifiability_result["basis"])] - locally_not_globally_id_sing = [SIAN.parent_ring_change(locally_not_globally_id[i], - OR) - for i in eachindex(locally_not_globally_id)] - OIdeal = Oscar.ideal(OR, basis_sing) + throw("Not implemented, sorry.") + # globally_id = identifiability_result["identifiability_nemo"]["globally"] + # locally_not_globally_id = identifiability_result["identifiability_nemo"]["locally_not_globally"] + # non_id = identifiability_result["identifiability_nemo"]["nonidentifiable"] + # weights_table = identifiability_result["weights"] + # non_jet_ring = identifiability_result["non_jet_ring"] + # n2m = identifiability_result["nemo_mtk"] + # solutions_table = Dict{Any, Int}() + # for param in globally_id + # v = SIAN.get_order_var(param, non_jet_ring)[1] + # solutions_table[n2m[v]] = 1 + # end + # for param in non_id + # v = SIAN.get_order_var(param, non_jet_ring)[1] + # solutions_table[n2m[v]] = 0 + # end + # # basis_lex = Groebner.fglm(identifiability_result.basis) + # R = parent(identifiability_result["basis"][1]) + # OR, svars = Oscar.PolynomialRing(Oscar.QQ, string.(gens(R)); ordering = :lex) + # basis_sing = [SIAN.parent_ring_change(identifiability_result["basis"][i], OR) + # for i in 1:length(identifiability_result["basis"])] + # locally_not_globally_id_sing = [SIAN.parent_ring_change(locally_not_globally_id[i], + # OR) + # for i in eachindex(locally_not_globally_id)] + # OIdeal = Oscar.ideal(OR, basis_sing) - for (idx, param) in enumerate(locally_not_globally_id_sing) - others = setdiff(svars, [param]) - elim = Oscar.gens(Oscar.eliminate(OIdeal, others)) - polynomials = filter(x -> param in Oscar.vars(x), - elim)[1] - param_idx = findfirst(x -> x == param, svars) - v = SIAN.get_order_var(locally_not_globally_id[idx], non_jet_ring)[1] - dgrs = Oscar.degrees(polynomials)[param_idx] / - get(weights_table, v, 1) - solutions_table[n2m[v]] = Int(dgrs) - end - return solutions_table + # for (idx, param) in enumerate(locally_not_globally_id_sing) + # others = setdiff(svars, [param]) + # elim = Oscar.gens(Oscar.eliminate(OIdeal, others)) + # polynomials = filter(x -> param in Oscar.vars(x), + # elim)[1] + # param_idx = findfirst(x -> x == param, svars) + # v = SIAN.get_order_var(locally_not_globally_id[idx], non_jet_ring)[1] + # dgrs = Oscar.degrees(polynomials)[param_idx] / + # get(weights_table, v, 1) + # solutions_table[n2m[v]] = Int(dgrs) + # end + # return solutions_table end diff --git a/src/rational_interpolation/utils.jl b/src/rational_interpolation/utils.jl index fbf42378..a7cb5549 100644 --- a/src/rational_interpolation/utils.jl +++ b/src/rational_interpolation/utils.jl @@ -38,10 +38,10 @@ function eval_derivs(polynomial_system, interpolant::Interpolant, end end elseif isequal(method, :msolve) - y_derivs = Vector{SIAN.Nemo.fmpq_mpoly}() - y_vals = Vector{SIAN.Nemo.fmpq}() - u_derivs = Vector{SIAN.Nemo.fmpq_mpoly}() - u_vals = Vector{SIAN.Nemo.fmpq}() + y_derivs = Vector{Nemo.QQMPolyRingElem}() + y_vals = Vector{Nemo.fmpq}() + u_derivs = Vector{Nemo.QQMPolyRingElem}() + u_vals = Vector{Nemo.fmpq}() for (y_func, y_deriv_order) in pairs(identifiability_result["Y_eq"]) if occursin(y_function_name, string(y_func)) push!(y_derivs, y_func) @@ -62,7 +62,7 @@ function eval_derivs(polynomial_system, interpolant::Interpolant, end end end - polynomial_system = [SIAN.Nemo.evaluate(poly, y_derivs, y_vals) + polynomial_system = [Nemo.evaluate(poly, y_derivs, y_vals) for poly in polynomial_system] end return polynomial_system diff --git a/src/result.jl b/src/result.jl index 9f3f3c8f..b3b5277c 100644 --- a/src/result.jl +++ b/src/result.jl @@ -37,7 +37,7 @@ struct EstimationResult for p in ModelingToolkit.parameters(model) parameters[ModelingToolkit.Num(p)] = get(poly_sol, p, nothing) end - for s in ModelingToolkit.states(model) + for s in ModelingToolkit.unknowns(model) states[ModelingToolkit.Num(s)] = get(poly_sol, s, nothing) end new(parameters, states, degree, at_time, nothing, interpolants, return_code, diff --git a/src/utils.jl b/src/utils.jl index 8e548ed3..f5204919 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -23,7 +23,7 @@ function nemo2hc(expr_tree::Union{Expr, Symbol}) end end -function nemo2hc(expr_tree::fmpq_mpoly) +function nemo2hc(expr_tree::QQMPolyRingElem) # println(expr_tree) return nemo2hc(Meta.parse(string(expr_tree))) end @@ -32,8 +32,8 @@ function nemo2hc(expr_tree::Number) return expr_tree end -function nemo2hc(expr_tree::Oscar.Generic.Frac) - numer, denom = Oscar.numerator(expr_tree), Oscar.denominator(expr_tree) +function nemo2hc(expr_tree::Nemo.Generic.FracFieldElem) + numer, denom = Nemo.numerator(expr_tree), Nemo.denominator(expr_tree) return nemo2hc(numer) / nemo2hc(denom) end diff --git a/test/msolve_estimation_tests.jl b/test/msolve_estimation_tests.jl index 98bc85df..d3a2c301 100644 --- a/test/msolve_estimation_tests.jl +++ b/test/msolve_estimation_tests.jl @@ -1,30 +1,31 @@ +# TODO: msolve functionality is broken. @testset "Run small estimation code to test correctness" begin - using ParameterEstimation - using ModelingToolkit # ODE definitions + # using ParameterEstimation + # using ModelingToolkit # ODE definitions - # define toy model - @parameters mu - @variables t x1(t) y1(t) - D = Differential(t) - @named model = ODESystem([D(x1) ~ -mu * x1], - t, [x1], [mu]) - outs = [y1 ~ x1 + x1^2] - time = [0.0, 1.0] # sampling interval - data = Dict{Any, Vector{Float64}}("t" => [0.0, 1 / 3, 2 / 3, 1.0], - x1 + x1^2 => [2.0, 1.56301, - 1.22995, 0.97441]) - res = ParameterEstimation.estimate(model, outs, data; method = :msolve) - # @test length(res) == 1 - @test Symbol(res[1].return_code) == :Success - @test isapprox(res[1].parameters[mu], 0.5, atol = 1e-3) - @test isapprox(res[1].states[x1], 1.0, atol = 1e-3) + # # define toy model + # @parameters mu + # @variables t x1(t) y1(t) + # D = Differential(t) + # @named model = ODESystem([D(x1) ~ -mu * x1], + # t, [x1], [mu]) + # outs = [y1 ~ x1 + x1^2] + # time = [0.0, 1.0] # sampling interval + # data = Dict{Any, Vector{Float64}}("t" => [0.0, 1 / 3, 2 / 3, 1.0], + # x1 + x1^2 => [2.0, 1.56301, + # 1.22995, 0.97441]) + # res = ParameterEstimation.estimate(model, outs, data; method = :msolve) + # # @test length(res) == 1 + # @test Symbol(res[1].return_code) == :Success + # @test isapprox(res[1].parameters[mu], 0.5, atol = 1e-3) + # @test isapprox(res[1].states[x1], 1.0, atol = 1e-3) - # uneven time sample - data = Dict{Any, Vector{Float64}}(x1^2 + x1 => [2.0, 1.98506, 1.17611, 0.97441], - "t" => [0.0, 0.01, 0.73, 1.0]) - res = ParameterEstimation.estimate(model, outs, data; method = :msolve) - #@test length(res) == 1 - @test Symbol(res[1].return_code) == :Success - @test isapprox(res[1].parameters[mu], 0.5, atol = 1e-3) - @test isapprox(res[1].states[x1], 1.0, atol = 1e-3) + # # uneven time sample + # data = Dict{Any, Vector{Float64}}(x1^2 + x1 => [2.0, 1.98506, 1.17611, 0.97441], + # "t" => [0.0, 0.01, 0.73, 1.0]) + # res = ParameterEstimation.estimate(model, outs, data; method = :msolve) + # #@test length(res) == 1 + # @test Symbol(res[1].return_code) == :Success + # @test isapprox(res[1].parameters[mu], 0.5, atol = 1e-3) + # @test isapprox(res[1].states[x1], 1.0, atol = 1e-3) end diff --git a/test/nemo_hc_conversion.jl b/test/nemo_hc_conversion.jl index 7391e402..c4c9ea53 100644 --- a/test/nemo_hc_conversion.jl +++ b/test/nemo_hc_conversion.jl @@ -1,5 +1,5 @@ @testset "Convert Nemo polynomial systems into HomotopyContinuation type" begin - R, (X, Y, Z) = Oscar.PolynomialRing(Oscar.QQ, ["x", "y", "z"]) + R, (X, Y, Z) = Nemo.polynomial_ring(Nemo.QQ, ["x", "y", "z"]) @var x y z Et = [X + Y + Z - 1, X + Y - 2, Z * X - 1] E = HomotopyContinuation.System([ParameterEstimation.nemo2hc(e) for e in Et]) @@ -11,8 +11,8 @@ Et = [X^2 + 1 // 3 * Y + Z * X - 1, -X + Y^3 * X^2 - 2, Z * X - 1] E = HomotopyContinuation.System([ParameterEstimation.nemo2hc(e) for e in Et]) E_native = HomotopyContinuation.System([x^2 + 1 // 3 * y + z * x - 1, - -x + y^3 * x^2 - 2, - x * z - 1]) + -x + y^3 * x^2 - 2, + x * z - 1]) @info "Source System: $Et" @info "Converted System: $E" @test E_native == E @@ -20,10 +20,10 @@ Et = [X * Y * Z * X // 10 - 1, -X - Y - 2, Z * X - 1] E = HomotopyContinuation.System([ParameterEstimation.nemo2hc(e) for e in Et]) E_native = HomotopyContinuation.System([ - x * y * z * x // 10 - 1, - -x - y - 2, - z * x - 1, - ]) + x * y * z * x // 10 - 1, + -x - y - 2, + z * x - 1 + ]) @info "Source System: $Et" @info "Converted System: $E" @test E_native == E diff --git a/test/runtests.jl b/test/runtests.jl index 3bdaaea4..cf060bb5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,6 @@ using Test, TestSetExtensions using ModelingToolkit, SIAN, HomotopyContinuation -using ParameterEstimation, Oscar +using ParameterEstimation @info "Testing started" diff --git a/test/test_alg_indep.jl b/test/test_alg_indep.jl index fd791c7a..32020b9e 100644 --- a/test/test_alg_indep.jl +++ b/test/test_alg_indep.jl @@ -1,6 +1,6 @@ @testset "Check size of transcendence_basis of polynomial systems" begin - R, (a, b, c) = SIAN.Nemo.PolynomialRing(SIAN.Nemo.QQ, ["a", "b", "c"]; - ordering = :degrevlex) + R, (a, b, c) = Nemo.polynomial_ring(Nemo.QQ, ["a", "b", "c"]; + internal_ordering = :degrevlex) Et = [a + b, a + b + c] vals = [1, 1, 1] _, tr_basis = ParameterEstimation.algebraic_independence(Et, [a, b, c], vals) @@ -21,8 +21,8 @@ _, tr_basis = ParameterEstimation.algebraic_independence(Et, [a, b, c], vals) @test length(tr_basis) == 0 - R, (a, b, c, d) = SIAN.Nemo.PolynomialRing(SIAN.Nemo.QQ, ["a", "b", "c", "d"]; - ordering = :degrevlex) + R, (a, b, c, d) = Nemo.polynomial_ring(Nemo.QQ, ["a", "b", "c", "d"]; + internal_ordering = :degrevlex) vals = [1, 1, 1, 1] Et = [a + d, b + c, a + a * b * c + a^2 * d] _, tr_basis = ParameterEstimation.algebraic_independence(Et, [a, b, d, c], vals) From 2f612eee68d0363f507f28fdd42a458a1cd409b5 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Thu, 18 Apr 2024 07:24:31 +0200 Subject: [PATCH 2/7] bump versions, now for real --- Project.toml | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index a1526ed0..51475e65 100644 --- a/Project.toml +++ b/Project.toml @@ -14,8 +14,8 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -Oscar = "f1435218-dba5-11e9-1e4d-f1a5fab5fc13" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" @@ -33,19 +33,26 @@ ArbNumerics = "1" BaryRational = "1" DifferentialEquations = "7" ForwardDiff = "0.10" -Groebner = "0.2, 0.3, 0.4 - 0.5" +Groebner = "0.6, 0.7" HomotopyContinuation = "2" LinearSolve = "2" -ModelingToolkit = "8" +ModelingToolkit = "9" OrderedCollections = "1" -Oscar = "0.10 - 0.20" PrecompileTools = "1" ProgressMeter = "1" -SIAN = "1 - 1.5.10" -StructuralIdentifiability = "0.4.9 - 0.5.30" +SIAN = "1.5" +StructuralIdentifiability = "0.5" Suppressor = "0.2" Symbolics = "4, 5" TaylorDiff = "0.2" TaylorSeries = "0.12 - 0.30" +Test = "1.9" TestSetExtensions = "2" julia = "1.9" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04" + +[targets] +test = ["Test", "TestSetExtensions"] From 48c0dd7a63c33851075d85cf2e8aa078cd051b34 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Fri, 19 Apr 2024 17:04:47 +0200 Subject: [PATCH 3/7] fix a bug in algebraic_independence --- Project.toml | 2 +- src/estimation/fixed_degree.jl | 9 ++++----- src/identifiability/check_identifiability.jl | 2 +- src/identifiability/preprocessor.jl | 3 +-- src/identifiability/transcendence_basis.jl | 6 +++--- 5 files changed, 10 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index 51475e65..9452f87d 100644 --- a/Project.toml +++ b/Project.toml @@ -41,7 +41,7 @@ OrderedCollections = "1" PrecompileTools = "1" ProgressMeter = "1" SIAN = "1.5" -StructuralIdentifiability = "0.5" +StructuralIdentifiability = "0.5.7" Suppressor = "0.2" Symbolics = "4, 5" TaylorDiff = "0.2" diff --git a/src/estimation/fixed_degree.jl b/src/estimation/fixed_degree.jl index e0b3f695..f0d8c395 100644 --- a/src/estimation/fixed_degree.jl +++ b/src/estimation/fixed_degree.jl @@ -13,14 +13,14 @@ function backsolve_initial_conditions(model, E, report_time, inputs::Vector{Equa initial_conditions = [E[s] for s in ModelingToolkit.unknowns(model)] parameter_values = [E[p] for p in ModelingToolkit.parameters(model)] tspan = (E.at_time, report_time) #this is almost always backwards! - + ode_equations = ModelingToolkit.equations(model) ode_equations = substitute(ode_equations, Dict(each.lhs => Num(each.rhs) for each in inputs)) t = ModelingToolkit.get_iv(model) @named new_model = ODESystem(ode_equations, t, ModelingToolkit.unknowns(model), ModelingToolkit.parameters(model)) - prob = ModelingToolkit.complete(ODEProblem(new_model, initial_conditions, tspan, parameter_values)) + prob = ODEProblem(ModelingToolkit.complete(new_model), initial_conditions, tspan, parameter_values) saveat = range(tspan[1], tspan[2], length = length(data_sample["t"])) ode_solution = ModelingToolkit.solve(prob, solver, p = parameter_values, saveat = saveat, abstol = abstol, reltol = reltol) @@ -107,7 +107,7 @@ function estimate_single_interpolator(model::ModelingToolkit.ODESystem, diff_order = num_parameters + 1, #todo(OREBAS): is this always forcing num_parameters + 1 derivatives? at_t = at_time, method = method) - + if method == :homotopy all_solutions = solve_via_homotopy(identifiability_result, model; real_tol = real_tol) @@ -117,8 +117,7 @@ function estimate_single_interpolator(model::ModelingToolkit.ODESystem, else throw(ArgumentError("Method $method not supported")) end - #println(all_solutions) - #println("HERE") + all_solutions = [EstimationResult(model, each, interpolator.first, at_time, interpolants, ReturnCode.Success, datasize, report_time) for each in all_solutions] diff --git a/src/identifiability/check_identifiability.jl b/src/identifiability/check_identifiability.jl index c0b74fa1..102c68ac 100644 --- a/src/identifiability/check_identifiability.jl +++ b/src/identifiability/check_identifiability.jl @@ -216,7 +216,7 @@ function identifiability_ode(ode, params_to_assess; p = 0.99, p_mod = 0, infolev end end x_theta_vars_reorder = vcat(theta_l, - reverse([x for x in x_theta_vars if !(x in theta_l)])) + reverse([x for x in x_theta_vars if !(x in theta_l)])) Et_ids, alg_indep = algebraic_independence(Et_eval_base, x_theta_vars_reorder, all_x_theta_vars_subs) @info "Found Pivots: [$(join(alg_indep, ", "))]" diff --git a/src/identifiability/preprocessor.jl b/src/identifiability/preprocessor.jl index 484347ed..8556ced2 100644 --- a/src/identifiability/preprocessor.jl +++ b/src/identifiability/preprocessor.jl @@ -15,8 +15,7 @@ function preprocess_ode(de::ModelingToolkit.ODESystem, diff_eqs = filter(eq -> !(ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(de)) y_functions = [each.lhs for each in measured_quantities] - state_vars = filter(s -> !(ModelingToolkit.isinput(s) || ModelingToolkit.isoutput(s)), - ModelingToolkit.unknowns(de)) + state_vars = ModelingToolkit.unknowns(de) params = ModelingToolkit.parameters(de) t = ModelingToolkit.arguments(measured_quantities[1].lhs)[1] params_from_measured_quantities = ModelingToolkit.parameters(ModelingToolkit.ODESystem( diff --git a/src/identifiability/transcendence_basis.jl b/src/identifiability/transcendence_basis.jl index e7ad6992..c00dc58a 100644 --- a/src/identifiability/transcendence_basis.jl +++ b/src/identifiability/transcendence_basis.jl @@ -22,15 +22,15 @@ function algebraic_independence(Et::Vector{Nemo.QQMPolyRingElem}, row = U[row_idx, :] if !all(row .== 0) pivot_col = findfirst(row .!= 0) - push!(pivots, indets[pivot_col[2]]) + push!(pivots, indets[pivot_col]) end end current_idx = 1 - output_rows = Jacobian[current_idx, :] + output_rows = Jacobian[[current_idx], :] current_rank = 1 output_ids = [1] for current_idx in 2:length(Et) - current = [output_rows; Jacobian[current_idx, :]] + current = vcat(output_rows, Jacobian[[current_idx], :]) if Nemo.rank(current) > current_rank output_rows = current push!(output_ids, current_idx) From 7eafde6453f8d5a0bc8553dafc197d6f2defa273 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Fri, 19 Apr 2024 17:06:08 +0200 Subject: [PATCH 4/7] update Julia version --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 9452f87d..21f36acd 100644 --- a/Project.toml +++ b/Project.toml @@ -46,9 +46,9 @@ Suppressor = "0.2" Symbolics = "4, 5" TaylorDiff = "0.2" TaylorSeries = "0.12 - 0.30" -Test = "1.9" +Test = "1.10" TestSetExtensions = "2" -julia = "1.9" +julia = "1.10" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 6480723af2eccea5146eff5302bb4e7be452253e Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Fri, 19 Apr 2024 17:52:50 +0200 Subject: [PATCH 5/7] replace Nemo.fmpq -> Nemo.QQFieldElem --- src/rational_interpolation/utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/rational_interpolation/utils.jl b/src/rational_interpolation/utils.jl index a7cb5549..7f02cff0 100644 --- a/src/rational_interpolation/utils.jl +++ b/src/rational_interpolation/utils.jl @@ -39,9 +39,9 @@ function eval_derivs(polynomial_system, interpolant::Interpolant, end elseif isequal(method, :msolve) y_derivs = Vector{Nemo.QQMPolyRingElem}() - y_vals = Vector{Nemo.fmpq}() + y_vals = Vector{Nemo.QQFieldElem}() u_derivs = Vector{Nemo.QQMPolyRingElem}() - u_vals = Vector{Nemo.fmpq}() + u_vals = Vector{Nemo.QQFieldElem}() for (y_func, y_deriv_order) in pairs(identifiability_result["Y_eq"]) if occursin(y_function_name, string(y_func)) push!(y_derivs, y_func) From 37fb3cb3c3ed612b2e70b737225ebc976941a374 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Sat, 27 Apr 2024 16:10:36 +0200 Subject: [PATCH 6/7] fix tests (note that threaded code seems broken) --- src/estimation/serial.jl | 2 +- src/estimation/threaded.jl | 17 +++++++++-------- src/estimation/utils.jl | 4 ++-- src/utils.jl | 2 +- test/runtests.jl | 1 + 5 files changed, 14 insertions(+), 12 deletions(-) diff --git a/src/estimation/serial.jl b/src/estimation/serial.jl index e7070378..0e3cd2e1 100644 --- a/src/estimation/serial.jl +++ b/src/estimation/serial.jl @@ -38,5 +38,5 @@ function estimate_serial(model::ModelingToolkit.ODESystem, ]) end end - return post_process(estimates, filtermode, parameter_constraints, ic_constraints) + return post_process(estimates, filtermode, parameter_constraints, ic_constraints; threaded=false) end diff --git a/src/estimation/threaded.jl b/src/estimation/threaded.jl index 0cb4ec52..e0ca1af8 100644 --- a/src/estimation/threaded.jl +++ b/src/estimation/threaded.jl @@ -8,7 +8,7 @@ function estimate_threaded(model::ModelingToolkit.ODESystem, check_inputs(measured_quantities, data_sample) datasize = length(first(values(data_sample))) - if isnothing(interpolators) === nothing + if isnothing(interpolators) interpolators = default_interpolator(datasize) end @@ -20,12 +20,13 @@ estimates = Vector{Vector{ParameterEstimation.EstimationResult}}() @info "Estimating via the interpolators: $(keys(interpolators))" n_threads = Threads.nthreads() N = length(interpolators) + interpolators = collect(interpolators) estimates = Vector{Any}(nothing, n_threads) Threads.@threads :static for t in 1:N - interp = interp[t] - id = Threads.threadid() - if isnothing(estimates[id]) - estimates[id] = [] + interp = interpolators[t] + tid = Threads.threadid() + if isnothing(estimates[tid]) + estimates[tid] = [] end unfiltered = estimate_single_interpolator(model, measured_quantities, inputs, data_sample; identifiability_result = id, @@ -35,10 +36,10 @@ estimates = Vector{Vector{ParameterEstimation.EstimationResult}}() if length(unfiltered) > 0 filtered = filter_solutions(unfiltered, id, model, inputs, data_sample; solver = solver, filtermode) - push!(estimates[id], filtered) + push!(estimates[tid], filtered) else - push!(estimates[id], + push!(estimates[tid], [ EstimationResult(model, Dict(), interp, at_time, Dict{Any, Interpolant}(), ReturnCode.Failure, @@ -46,5 +47,5 @@ estimates = Vector{Vector{ParameterEstimation.EstimationResult}}() ]) end end - return post_process(estimates,filtermode,parameter_constraints,ic_constraints) + return post_process(estimates,filtermode,parameter_constraints,ic_constraints,threaded=true) end diff --git a/src/estimation/utils.jl b/src/estimation/utils.jl index ce9b0c67..5634b726 100644 --- a/src/estimation/utils.jl +++ b/src/estimation/utils.jl @@ -68,8 +68,8 @@ function new_clustering(estimates, tol = 1e-3) return new_estimates end -function post_process(estimates, filtermode = :new, parameter_constraints = nothing, ic_constraints = nothing) - if Threads.nthreads() > 1 +function post_process(estimates, filtermode = :new, parameter_constraints = nothing, ic_constraints = nothing; threaded=false) + if threaded estimates = filter(x -> !isnothing(x), estimates) estimates = vcat(estimates...) end diff --git a/src/utils.jl b/src/utils.jl index f5204919..1c55d1e9 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -107,7 +107,7 @@ function sample_data(model::ModelingToolkit.ODESystem, else sampling_times = range(time_interval[1], time_interval[2], length = num_points) end - problem = ODEProblem(model, u0, time_interval, p_true) + problem = ODEProblem(ModelingToolkit.complete(model), u0, time_interval, p_true) solution_true = ModelingToolkit.solve(problem, solver, p = p_true, saveat = sampling_times; abstol, reltol) diff --git a/test/runtests.jl b/test/runtests.jl index cf060bb5..a177ae8f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,6 @@ using Test, TestSetExtensions using ModelingToolkit, SIAN, HomotopyContinuation +using Nemo using ParameterEstimation @info "Testing started" From 99f4bd59d9c6cedcfeb644672168aaa3f0088984 Mon Sep 17 00:00:00 2001 From: Sasha Demin Date: Mon, 29 Apr 2024 03:03:15 +0200 Subject: [PATCH 7/7] Use Dict instead of Vector to communicate parameters to ODEProblem --- src/estimation/fixed_degree.jl | 8 ++++++-- src/filtering.jl | 6 +++++- src/utils.jl | 4 ++-- 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/estimation/fixed_degree.jl b/src/estimation/fixed_degree.jl index f0d8c395..784ba01e 100644 --- a/src/estimation/fixed_degree.jl +++ b/src/estimation/fixed_degree.jl @@ -20,10 +20,14 @@ function backsolve_initial_conditions(model, E, report_time, inputs::Vector{Equa t = ModelingToolkit.get_iv(model) @named new_model = ODESystem(ode_equations, t, ModelingToolkit.unknowns(model), ModelingToolkit.parameters(model)) - prob = ODEProblem(ModelingToolkit.complete(new_model), initial_conditions, tspan, parameter_values) + prob = ODEProblem( + ModelingToolkit.complete(new_model), + initial_conditions, + tspan, + Dict(ModelingToolkit.parameters(model) .=> parameter_values)) saveat = range(tspan[1], tspan[2], length = length(data_sample["t"])) - ode_solution = ModelingToolkit.solve(prob, solver, p = parameter_values, saveat = saveat, abstol = abstol, reltol = reltol) + ode_solution = ModelingToolkit.solve(prob, solver, saveat = saveat, abstol = abstol, reltol = reltol) state_param_map = (Dict(x => replace(string(x), "(t)" => "") for x in ModelingToolkit.unknowns(model))) diff --git a/src/filtering.jl b/src/filtering.jl index dfa4e45d..deff8803 100644 --- a/src/filtering.jl +++ b/src/filtering.jl @@ -29,7 +29,11 @@ function solve_ode(model, estimate::EstimationResult, inputs::Vector{Equation}, t = ModelingToolkit.get_iv(model) @named new_model = ODESystem(ode_equations, t, ModelingToolkit.unknowns(model), ModelingToolkit.parameters(model)) - prob = ODEProblem(ModelingToolkit.complete(new_model), initial_conditions, tspan, parameter_values) + prob = ODEProblem( + ModelingToolkit.complete(new_model), + initial_conditions, tspan, + Dict(ModelingToolkit.parameters(new_model) .=> parameter_values) + ) ode_solution = ModelingToolkit.solve(prob, solver, saveat = range(tspan[1], tspan[2], length = length(data_sample["t"])), abstol = 1e-12, reltol = 1e-12) diff --git a/src/utils.jl b/src/utils.jl index 1c55d1e9..26bfa53d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -107,8 +107,8 @@ function sample_data(model::ModelingToolkit.ODESystem, else sampling_times = range(time_interval[1], time_interval[2], length = num_points) end - problem = ODEProblem(ModelingToolkit.complete(model), u0, time_interval, p_true) - solution_true = ModelingToolkit.solve(problem, solver, p = p_true, + problem = ODEProblem(ModelingToolkit.complete(model), u0, time_interval, Dict(ModelingToolkit.parameters(model) .=> p_true)) + solution_true = ModelingToolkit.solve(problem, solver, saveat = sampling_times; abstol, reltol) data_sample = OrderedDict{Any, Vector{T}}(Num(v.rhs) => solution_true[Num(v.rhs)]