From 76f54e54e75dbc22a9128a9c7faf33287f9e2505 Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Mon, 22 Jun 2026 17:46:05 +0200 Subject: [PATCH 01/17] the cost of being simple --- compare.jl | 45 +++++++++++++++++++++++++++++++++++ src/global_identifiability.jl | 4 +++- src/identifiable_functions.jl | 8 ++++++- src/known_ic.jl | 3 +++ 4 files changed, 58 insertions(+), 2 deletions(-) create mode 100644 compare.jl diff --git a/compare.jl b/compare.jl new file mode 100644 index 000000000..8aa9dd1b6 --- /dev/null +++ b/compare.jl @@ -0,0 +1,45 @@ +using Revise, StructuralIdentifiability, ParamPunPam, RationalFunctionFields +import StructuralIdentifiability: parent_ring_change, lie_derivative +import RationalFunctionFields: fields_equal + +cmp = (f, g) -> begin + f = parent_ring_change(f, parent(ode)) + g = parent_ring_change(g, parent(ode)) + df = lie_derivative(f, ode) + dg = lie_derivative(g, ode) + default_cmp = RationalFunctionFields.rational_function_cmp + if iszero(df) && iszero(dg) + default_cmp(f, g) + else + default_cmp(df, dg) + end +end + +ode = @ODEmodel( + x1'(t) = + (p1 * x4(t)) - (p3 * x1(t)) - + p4 * ((x1(t)^2 / (p5 + x1(t))) * (1 + (p6 * u1(t) / (p7 + u1(t))))), + x2'(t) = + p8 - (p9 * x2(t)) - + p10 * ( + (x1(t) * x2(t) / (p11 + x2(t))) * (1 + (p12 * u1(t) / (p13 + u1(t)))) + ), + x3'(t) = + p14 - (p15 * x3(t)) - + p16 * x1(t) * x3(t) * (1 - p18 * u1(t)) / (p17 + x3(t)), + x4'(t) = + p20 - p21 * (1 - p24) * (1 - p25) / ((p22^4) + 1) - (p20 * x4(t)) + + (p21 * (x3(t)^4)) + + (1 + p23 * u1(t)) * (1 - p24 * x1(t)) * (1 - p25 * x2(t)) / + (p22^4 + x3(t)^4), + y1(t) = x1(t), + y2(t) = x2(t), + y3(t) = x3(t), + y4(t) = x4(t) + ) + +id1 = find_identifiable_functions(ode, with_states=true, cmp = cmp) + +id2 = find_identifiable_functions(ode, with_states=true, cmp = RationalFunctionFields.rational_function_cmp) + +@assert fields_equal(RationalFunctionField(id1), RationalFunctionField(id2), 0.99) diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl index 4f97226ff..8d0a174e3 100644 --- a/src/global_identifiability.jl +++ b/src/global_identifiability.jl @@ -91,6 +91,7 @@ The function returns a tuple containing the following: with_states::Bool = false, var_change_policy = :default, rational_interpolator = :VanDerHoevenLecerf, + cmp = RationalFunctionFields.rational_function_cmp ) where {T} @info "Computing IO-equations" ioeq_time = @elapsed io_equations = @@ -148,6 +149,7 @@ The function returns a tuple containing the following: seed = 42, simplify = :standard, rational_interpolator = rational_interpolator, + cmp = cmp ) dennums_simplified = fractions_to_dennums(no_states_simplified) # switch back the ring @@ -202,7 +204,7 @@ Output: a list L of booleans with L[i] being the identifiability status of the i known = known, prob_threshold = half_p, var_change_policy = var_change_policy, - with_states = states_needed, + with_states = states_needed ) funcs_to_check = Vector{Generic.FracFieldElem{P}}( diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index 02079564b..1447123c3 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -54,6 +54,7 @@ function find_identifiable_functions( simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, loglevel = Logging.Info, + cmp = RationalFunctionFields.rational_function_cmp, ) where {T <: MPolyRingElem{QQFieldElem}} restart_logging(loglevel = loglevel) reset_timings() @@ -66,6 +67,7 @@ function find_identifiable_functions( with_states = with_states, simplify = simplify, rational_interpolator = rational_interpolator, + cmp = cmp ) else id_funcs = _find_identifiable_functions_kic( @@ -75,6 +77,7 @@ function find_identifiable_functions( seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, + cmp = cmp ) # renaming variables from `x(t)` to `x(0)` return replace_with_ic(ode, id_funcs) @@ -89,6 +92,7 @@ function _find_identifiable_functions( with_states = false, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, + cmp = RationalFunctionFields.rational_function_cmp ) where {T <: MPolyRingElem{QQFieldElem}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) @@ -109,6 +113,7 @@ function _find_identifiable_functions( prob_threshold = half_p, with_states = with_states, rational_interpolator = rational_interpolator, + cmp = cmp ) # If simplification is needed if simplify !== :absent @@ -122,7 +127,8 @@ function _find_identifiable_functions( seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, - priority_variables = [parent_ring_change(p, bring) for p in ode.parameters], + # priority_variables = [parent_ring_change(p, bring) for p in ode.parameters], + cmp = cmp ) else id_funcs_fracs = RationalFunctionFields.dennums_to_fractions(id_funcs) diff --git a/src/known_ic.jl b/src/known_ic.jl index 83f359bce..bd0cccf08 100644 --- a/src/known_ic.jl +++ b/src/known_ic.jl @@ -33,6 +33,7 @@ function _find_identifiable_functions_kic( seed = 42, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, + cmp = RationalFunctionFields.rational_function_cmp ) where {T <: MPolyRingElem{Nemo.QQFieldElem}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) @@ -45,6 +46,7 @@ function _find_identifiable_functions_kic( simplify = :absent, rational_interpolator = rational_interpolator, seed = seed, + cmp = cmp ) id_funcs = vcat(id_funcs_general, [f // one(parent(ode)) for f in known_ic]) @@ -56,6 +58,7 @@ function _find_identifiable_functions_kic( seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, + cmp = cmp ) end From ca88bc06f161ed559f0c16a15ffed4259c93e49f Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Mon, 27 Jul 2026 15:35:37 +0200 Subject: [PATCH 02/17] experiment with comparator --- compare-units.jl | 70 ++++++++++++++++++++++++++++++++++++++++++++++++ compare.jl | 31 +++++++-------------- 2 files changed, 80 insertions(+), 21 deletions(-) create mode 100644 compare-units.jl diff --git a/compare-units.jl b/compare-units.jl new file mode 100644 index 000000000..67181f794 --- /dev/null +++ b/compare-units.jl @@ -0,0 +1,70 @@ +using Revise, StructuralIdentifiability, ParamPunPam, RationalFunctionFields, Unitful, Nemo, AbstractAlgebra +import StructuralIdentifiability: parent_ring_change, lie_derivative +import RationalFunctionFields: fields_equal + +function eval_at_dict(poly::P, d::Dict{P,S}) where {P<:MPolyRingElem,S} + R = parent(poly) + @assert R == parent(first(keys(d))) + xs = gens(parent(first(keys(d)))) + xs_sym = [get(d, x, 0) for x in xs if string(x) in map(string, gens(R))] + accum = nothing + try + for t in terms(poly) + cf = coeff(t, 1) + cf_ = BigInt(numerator(cf)) // BigInt(denominator(cf)) + ex = exponent_vector(t, 1) + if accum == nothing + accum = cf_ * prod(xs_sym .^ ex) + else + accum += cf_ * prod(xs_sym .^ ex) + end + end + catch + return "beda!" + end + return accum +end + +eval_at_dict(frac::AbstractAlgebra.Generic.FracFieldElem{T}, d) where {T} = + begin + num = eval_at_dict(numerator(frac), d) + den = eval_at_dict(denominator(frac), d) + "beda!" in (num, den) && return "beda!" + num // den + end + +ode = @ODEmodel( + S'(t) = -beta_W*S(t)*W(t) - beta_I*S(t)*I(t), + I'(t) = beta_W*S(t)*W(t) + beta_I*S(t)*I(t) - gamma*I(t), + W'(t) = alpha*I(t) - zeta*W(t), + y(t) = W(t) +) + +# Domain knowledge is that the units of S and I are the same. +units = Dict( + S => u"P", # meaning: "population" + I => u"P", + W => u"C", # meaning: "concentration" + gamma => u"1/s", # meaning: "1/second" + zeta => u"1/s", + beta_W => u"1/C * 1/s", + beta_I => u"1/P * 1/s", + alpha => u"C * 1/P * 1/s", +) + +cmp = (f, g) -> begin + f = parent_ring_change(f, parent(ode)) + g = parent_ring_change(g, parent(ode)) + uf = Unitful.unit(eval_at_dict(f, units)) + ug = Unitful.unit(eval_at_dict(g, units)) + length(string(uf)) < length(string(ug)) +end + +id1 = find_identifiable_functions(ode, with_states=true, cmp = cmp) + +id2 = find_identifiable_functions(ode, with_states=true, cmp = RationalFunctionFields.rational_function_cmp) + +@assert fields_equal(RationalFunctionField(id1), RationalFunctionField(id2), 0.99) + +u1 = Unitful.unit.(map(f -> eval_at_dict(f, units), id1)) +u2 = Unitful.unit.(map(f -> eval_at_dict(f, units), id2)) diff --git a/compare.jl b/compare.jl index 8aa9dd1b6..9e64b6ae8 100644 --- a/compare.jl +++ b/compare.jl @@ -16,27 +16,16 @@ cmp = (f, g) -> begin end ode = @ODEmodel( - x1'(t) = - (p1 * x4(t)) - (p3 * x1(t)) - - p4 * ((x1(t)^2 / (p5 + x1(t))) * (1 + (p6 * u1(t) / (p7 + u1(t))))), - x2'(t) = - p8 - (p9 * x2(t)) - - p10 * ( - (x1(t) * x2(t) / (p11 + x2(t))) * (1 + (p12 * u1(t) / (p13 + u1(t)))) - ), - x3'(t) = - p14 - (p15 * x3(t)) - - p16 * x1(t) * x3(t) * (1 - p18 * u1(t)) / (p17 + x3(t)), - x4'(t) = - p20 - p21 * (1 - p24) * (1 - p25) / ((p22^4) + 1) - (p20 * x4(t)) + - (p21 * (x3(t)^4)) + - (1 + p23 * u1(t)) * (1 - p24 * x1(t)) * (1 - p25 * x2(t)) / - (p22^4 + x3(t)^4), - y1(t) = x1(t), - y2(t) = x2(t), - y3(t) = x3(t), - y4(t) = x4(t) - ) + S'(t) = -beta_W*S(t)*W(t)-beta_I*S(t)*I(t), + I'(t) = beta_W*S(t)*W(t)+beta_I*S(t)*I(t) - gamma*I(t), + W'(t) = alpha*I(t)-zeta*W(t), + y(t) = W(t) +) + +# ode = @ODEmodel( +# x1'(t) = x1 + u(t), +# y(t) = x1 +# ) id1 = find_identifiable_functions(ode, with_states=true, cmp = cmp) From 5abbf8f92824005d7ee4bd5ce5477d3cb1a30190 Mon Sep 17 00:00:00 2001 From: Sasha Date: Mon, 27 Jul 2026 20:37:50 +0200 Subject: [PATCH 03/17] custom comparator for identifiable funcs --- Project.toml | 4 +- compare-units.jl | 70 --------------------------- compare.jl | 34 ------------- ext/ModelingToolkitSIExt.jl | 7 +++ src/identifiable_functions.jl | 6 ++- test/bodies/identifiable_functions.jl | 46 ++++++++++++++++++ 6 files changed, 60 insertions(+), 107 deletions(-) delete mode 100644 compare-units.jl delete mode 100644 compare.jl diff --git a/Project.toml b/Project.toml index a7ef97aa5..b530c5066 100644 --- a/Project.toml +++ b/Project.toml @@ -43,11 +43,11 @@ Logging = "1.10, 1.11" MacroTools = "0.5" ModelingToolkitBase = "1.12" Nemo = "0.46, 0.47, 0.48, 0.49, 0.50, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56" -ParamPunPam = "0.5.5" +ParamPunPam = "0.5" PrecompileTools = "1.2" Primes = "0.5" Random = "1.10, 1.11" -RationalFunctionFields = "0.3.1" +RationalFunctionFields = "0.3" SafeTestsets = "0.0.1, 0.1" SciMLTesting = "1" SpecialFunctions = "2" diff --git a/compare-units.jl b/compare-units.jl deleted file mode 100644 index 67181f794..000000000 --- a/compare-units.jl +++ /dev/null @@ -1,70 +0,0 @@ -using Revise, StructuralIdentifiability, ParamPunPam, RationalFunctionFields, Unitful, Nemo, AbstractAlgebra -import StructuralIdentifiability: parent_ring_change, lie_derivative -import RationalFunctionFields: fields_equal - -function eval_at_dict(poly::P, d::Dict{P,S}) where {P<:MPolyRingElem,S} - R = parent(poly) - @assert R == parent(first(keys(d))) - xs = gens(parent(first(keys(d)))) - xs_sym = [get(d, x, 0) for x in xs if string(x) in map(string, gens(R))] - accum = nothing - try - for t in terms(poly) - cf = coeff(t, 1) - cf_ = BigInt(numerator(cf)) // BigInt(denominator(cf)) - ex = exponent_vector(t, 1) - if accum == nothing - accum = cf_ * prod(xs_sym .^ ex) - else - accum += cf_ * prod(xs_sym .^ ex) - end - end - catch - return "beda!" - end - return accum -end - -eval_at_dict(frac::AbstractAlgebra.Generic.FracFieldElem{T}, d) where {T} = - begin - num = eval_at_dict(numerator(frac), d) - den = eval_at_dict(denominator(frac), d) - "beda!" in (num, den) && return "beda!" - num // den - end - -ode = @ODEmodel( - S'(t) = -beta_W*S(t)*W(t) - beta_I*S(t)*I(t), - I'(t) = beta_W*S(t)*W(t) + beta_I*S(t)*I(t) - gamma*I(t), - W'(t) = alpha*I(t) - zeta*W(t), - y(t) = W(t) -) - -# Domain knowledge is that the units of S and I are the same. -units = Dict( - S => u"P", # meaning: "population" - I => u"P", - W => u"C", # meaning: "concentration" - gamma => u"1/s", # meaning: "1/second" - zeta => u"1/s", - beta_W => u"1/C * 1/s", - beta_I => u"1/P * 1/s", - alpha => u"C * 1/P * 1/s", -) - -cmp = (f, g) -> begin - f = parent_ring_change(f, parent(ode)) - g = parent_ring_change(g, parent(ode)) - uf = Unitful.unit(eval_at_dict(f, units)) - ug = Unitful.unit(eval_at_dict(g, units)) - length(string(uf)) < length(string(ug)) -end - -id1 = find_identifiable_functions(ode, with_states=true, cmp = cmp) - -id2 = find_identifiable_functions(ode, with_states=true, cmp = RationalFunctionFields.rational_function_cmp) - -@assert fields_equal(RationalFunctionField(id1), RationalFunctionField(id2), 0.99) - -u1 = Unitful.unit.(map(f -> eval_at_dict(f, units), id1)) -u2 = Unitful.unit.(map(f -> eval_at_dict(f, units), id2)) diff --git a/compare.jl b/compare.jl deleted file mode 100644 index 9e64b6ae8..000000000 --- a/compare.jl +++ /dev/null @@ -1,34 +0,0 @@ -using Revise, StructuralIdentifiability, ParamPunPam, RationalFunctionFields -import StructuralIdentifiability: parent_ring_change, lie_derivative -import RationalFunctionFields: fields_equal - -cmp = (f, g) -> begin - f = parent_ring_change(f, parent(ode)) - g = parent_ring_change(g, parent(ode)) - df = lie_derivative(f, ode) - dg = lie_derivative(g, ode) - default_cmp = RationalFunctionFields.rational_function_cmp - if iszero(df) && iszero(dg) - default_cmp(f, g) - else - default_cmp(df, dg) - end -end - -ode = @ODEmodel( - S'(t) = -beta_W*S(t)*W(t)-beta_I*S(t)*I(t), - I'(t) = beta_W*S(t)*W(t)+beta_I*S(t)*I(t) - gamma*I(t), - W'(t) = alpha*I(t)-zeta*W(t), - y(t) = W(t) -) - -# ode = @ODEmodel( -# x1'(t) = x1 + u(t), -# y(t) = x1 -# ) - -id1 = find_identifiable_functions(ode, with_states=true, cmp = cmp) - -id2 = find_identifiable_functions(ode, with_states=true, cmp = RationalFunctionFields.rational_function_cmp) - -@assert fields_equal(RationalFunctionField(id1), RationalFunctionField(id2), 0.99) diff --git a/ext/ModelingToolkitSIExt.jl b/ext/ModelingToolkitSIExt.jl index 66f30474a..44022a99f 100644 --- a/ext/ModelingToolkitSIExt.jl +++ b/ext/ModelingToolkitSIExt.jl @@ -12,6 +12,8 @@ using TimerOutputs using ModelingToolkitBase +RationalFunctionFields = StructuralIdentifiability.RationalFunctionFields + # ------------------------------------------------------------------------------ # checking if it is a function of the form x(t), a bit dirty @@ -671,6 +673,7 @@ function StructuralIdentifiability.find_identifiable_functions( simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, loglevel = Logging.Info, + cmp = RationalFunctionFields.rational_function_cmp ) restart_logging(loglevel = loglevel) reset_timings() @@ -684,6 +687,7 @@ function StructuralIdentifiability.find_identifiable_functions( with_states = with_states, simplify = simplify, rational_interpolator = rational_interpolator, + cmp = cmp ) end end @@ -697,6 +701,7 @@ function _find_identifiable_functions( with_states = false, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, + cmp = RationalFunctionFields.rational_function_cmp ) Random.seed!(seed) ode, conversion = mtk_to_si(ode, measured_quantities) @@ -710,6 +715,7 @@ function _find_identifiable_functions( seed = seed, with_states = with_states, rational_interpolator = rational_interpolator, + cmp = cmp ) else result = StructuralIdentifiability._find_identifiable_functions_kic( @@ -719,6 +725,7 @@ function _find_identifiable_functions( prob_threshold = prob_threshold, seed = seed, rational_interpolator = rational_interpolator, + cmp = cmp ) end result = [parent_ring_change(f, ode.poly_ring) for f in result] diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index 1447123c3..fb990812f 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -22,6 +22,10 @@ This functions takes the following optional arguments: initial conditions, not states (this is an experimental functionality). - `prob_threshold`: A float in the range from 0 to 1, the probability of correctness. Default is `0.99`. +- `cmp`: Comparator for rational functions. + A custom comparator can be provided: for two rational functions `f1` and `f2`, + `cmp(f1,f2)` should be `true` if `f1` is simpler than `f2`. + Default is `RationalFunctionFields.rational_function_cmp`. - `seed`: The rng seed. Default value is `42`. - `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) @@ -39,7 +43,7 @@ ode = @ODEmodel( find_identifiable_functions(ode) # prints -3-element Vector{AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}: +2-element Vector{AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}: a12 + a01 + a21 a12*a01 ``` diff --git a/test/bodies/identifiable_functions.jl b/test/bodies/identifiable_functions.jl index 24cce99a2..dc04147bf 100644 --- a/test/bodies/identifiable_functions.jl +++ b/test/bodies/identifiable_functions.jl @@ -1197,5 +1197,51 @@ push!(test_cases, (ode = ode, ident_funcs = ident_funcs_states, with_states = tr end end +@testset "Identifiable functions with custom comparator" begin + cmp1 = (f, g) -> begin + length(string(f)) < length(string(g)) + end + + cmp2 = (f, g) -> begin + f = parent_ring_change(f, parent(ode)) + g = parent_ring_change(g, parent(ode)) + df = lie_derivative(f, ode) + dg = lie_derivative(g, ode) + default_cmp = RationalFunctionFields.rational_function_cmp + if iszero(df) && iszero(dg) + default_cmp(f, g) + else + default_cmp(df, dg) + end + end + + ode = @ODEmodel( + S'(t) = -beta_W * S(t) * W(t) - beta_I * S(t) * I(t), + I'(t) = beta_W * S(t) * W(t) + beta_I * S(t) * I(t) - gamma * I(t), + W'(t) = alpha * I(t) - zeta * W(t), + y(t) = W(t) + ) + + with_cmp1 = find_identifiable_functions(ode, with_states = true, cmp = cmp1) + + with_cmp2 = find_identifiable_functions(ode, with_states = true, cmp = cmp2) + + without_cmp = find_identifiable_functions(ode, with_states = true) + + @test fields_equal( + RationalFunctionField(with_cmp1), + RationalFunctionField(without_cmp), + 0.99, + ) + @test fields_equal( + RationalFunctionField(with_cmp2), + RationalFunctionField(without_cmp), + 0.99, + ) + + @test issorted(with_cmp1, lt = cmp1) + @test issorted(with_cmp2, lt = cmp2) +end + # restoring the default StructuralIdentifiability.RationalFunctionFields.SAT_FACTORIZATION_DEFAULT = :none From f77a975921a2d075e7b888e916658b03001235a9 Mon Sep 17 00:00:00 2001 From: Sasha Date: Mon, 27 Jul 2026 20:47:19 +0200 Subject: [PATCH 04/17] precise RFF version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index b530c5066..21f6ef1cd 100644 --- a/Project.toml +++ b/Project.toml @@ -47,7 +47,7 @@ ParamPunPam = "0.5" PrecompileTools = "1.2" Primes = "0.5" Random = "1.10, 1.11" -RationalFunctionFields = "0.3" +RationalFunctionFields = "0.3.2" SafeTestsets = "0.0.1, 0.1" SciMLTesting = "1" SpecialFunctions = "2" From c1990269c6bb1a7f34640940ce1f0092fa891361 Mon Sep 17 00:00:00 2001 From: Sasha Date: Mon, 27 Jul 2026 23:13:38 +0200 Subject: [PATCH 05/17] upd tests --- src/identifiable_functions.jl | 5 ++- src/util.jl | 35 ++++++++++++++++++++ test/bodies/identifiable_functions.jl | 46 +++++++++++++-------------- 3 files changed, 60 insertions(+), 26 deletions(-) diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index fb990812f..e6954629b 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -58,7 +58,7 @@ function find_identifiable_functions( simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, loglevel = Logging.Info, - cmp = RationalFunctionFields.rational_function_cmp, + cmp = cmp_prefer_params(ode, RationalFunctionFields.rational_function_cmp) ) where {T <: MPolyRingElem{QQFieldElem}} restart_logging(loglevel = loglevel) reset_timings() @@ -96,7 +96,7 @@ function _find_identifiable_functions( with_states = false, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, - cmp = RationalFunctionFields.rational_function_cmp + cmp = cmp_prefer_params(ode, RationalFunctionFields.rational_function_cmp) ) where {T <: MPolyRingElem{QQFieldElem}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) @@ -131,7 +131,6 @@ function _find_identifiable_functions( seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, - # priority_variables = [parent_ring_change(p, bring) for p in ode.parameters], cmp = cmp ) else diff --git a/src/util.jl b/src/util.jl index 81c85713b..6f1b8389c 100644 --- a/src/util.jl +++ b/src/util.jl @@ -302,3 +302,38 @@ function replace_with_ic(ode, funcs) end # ----------------------------------------------------------------------------- + +# Given an ODE and a comparator for rational functions returns a new comparator, +# which prefers functions containing only the parameters of the given ODE. +function cmp_prefer_params(ode, cmp) + some_vars = ode.parameters + new_cmp = (f, g) -> begin + isempty(some_vars) && return cmp(f, g) + some_vars = map(x -> parent_ring_change(x, parent(numerator(f))), some_vars) + prefer_f = issubset(vars(f), some_vars) + prefer_g = issubset(vars(g), some_vars) + (prefer_f && !prefer_g) && return true + (prefer_f == prefer_g) && return cmp(f, g) + false + end + new_cmp +end + +# Given an ODE and a comparator for rational functions returns a new comparator, +# which first compares the Lie derivatives when the latter are not zero. +function cmp_lie(ode, cmp) + new_cmp = (f, g) -> begin + f = parent_ring_change(f, parent(ode)) + g = parent_ring_change(g, parent(ode)) + df = lie_derivative(f, ode) + dg = lie_derivative(g, ode) + if iszero(df) && iszero(dg) + cmp(f, g) + else + cmp(df, dg) + end + end + new_cmp +end + +# ----------------------------------------------------------------------------- diff --git a/test/bodies/identifiable_functions.jl b/test/bodies/identifiable_functions.jl index dc04147bf..dc9febfa0 100644 --- a/test/bodies/identifiable_functions.jl +++ b/test/bodies/identifiable_functions.jl @@ -1113,10 +1113,10 @@ ident_funcs_states = [ k02 + k12, b1 * c, Km * c, - x1 * c, - x2 * c * k21, - (Vm * a2) // x2, Vm * a1 * c + Vm * a2 * c, + (a2 * k21) // (a1 + a2), + x1 * c, + x2 * c * k21 ] push!(test_cases, (ode = ode, ident_funcs = ident_funcs, with_states = false)) push!(test_cases, (ode = ode, ident_funcs = ident_funcs_states, with_states = true)) @@ -1198,23 +1198,6 @@ push!(test_cases, (ode = ode, ident_funcs = ident_funcs_states, with_states = tr end @testset "Identifiable functions with custom comparator" begin - cmp1 = (f, g) -> begin - length(string(f)) < length(string(g)) - end - - cmp2 = (f, g) -> begin - f = parent_ring_change(f, parent(ode)) - g = parent_ring_change(g, parent(ode)) - df = lie_derivative(f, ode) - dg = lie_derivative(g, ode) - default_cmp = RationalFunctionFields.rational_function_cmp - if iszero(df) && iszero(dg) - default_cmp(f, g) - else - default_cmp(df, dg) - end - end - ode = @ODEmodel( S'(t) = -beta_W * S(t) * W(t) - beta_I * S(t) * I(t), I'(t) = beta_W * S(t) * W(t) + beta_I * S(t) * I(t) - gamma * I(t), @@ -1222,25 +1205,42 @@ end y(t) = W(t) ) + cmp1 = (f, g) -> begin + length(string(f)) < length(string(g)) + end + + cmp2 = StructuralIdentifiability.cmp_lie( + ode, RationalFunctionFields.rational_function_cmp) + + cmp3 = StructuralIdentifiability.cmp_prefer_params(ode, cmp1) + with_cmp1 = find_identifiable_functions(ode, with_states = true, cmp = cmp1) with_cmp2 = find_identifiable_functions(ode, with_states = true, cmp = cmp2) - without_cmp = find_identifiable_functions(ode, with_states = true) + with_cmp3 = find_identifiable_functions(ode, with_states = true, cmp = cmp3) + + default_cmp = find_identifiable_functions(ode, with_states = true) @test fields_equal( RationalFunctionField(with_cmp1), - RationalFunctionField(without_cmp), + RationalFunctionField(default_cmp), 0.99, ) @test fields_equal( RationalFunctionField(with_cmp2), - RationalFunctionField(without_cmp), + RationalFunctionField(default_cmp), + 0.99, + ) + @test fields_equal( + RationalFunctionField(with_cmp3), + RationalFunctionField(default_cmp), 0.99, ) @test issorted(with_cmp1, lt = cmp1) @test issorted(with_cmp2, lt = cmp2) + @test issorted(with_cmp3, lt = cmp3) end # restoring the default From 3deb19fcce18daeec20cdce74d368a9166a51e47 Mon Sep 17 00:00:00 2001 From: Sasha Date: Mon, 27 Jul 2026 23:53:19 +0200 Subject: [PATCH 06/17] precise use of generators --- src/output_saturation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/output_saturation.jl b/src/output_saturation.jl index 1ccb739ba..7a987a58c 100644 --- a/src/output_saturation.jl +++ b/src/output_saturation.jl @@ -25,10 +25,10 @@ function saturate_outputs( sort!(lie_ders, lt = RationalFunctionFields.rational_function_cmp) for f in lie_ders if !first(RationalFunctionFields.check_algebraicity_modp(current_y, [f])) - current_y = RationalFunctionField(vcat(generators(current_y), [f])) + current_y = RationalFunctionField(vcat(RationalFunctionFields.generators(current_y), [f])) end end - new_outputs = generators(current_y)[(length(ode.y_vars) + length(ode.u_vars) + length(ode.parameters) + 1):end] + new_outputs = RationalFunctionFields.generators(current_y)[(length(ode.y_vars) + length(ode.u_vars) + length(ode.parameters) + 1):end] idx = 1 old_y_names = map(var_to_str, ode.y_vars) From e8677aeedb538bfdd495b3d1949c1191e923a29b Mon Sep 17 00:00:00 2001 From: Sasha Date: Tue, 28 Jul 2026 00:52:06 +0200 Subject: [PATCH 07/17] format --- src/ODE.jl | 5 +---- src/discrete.jl | 5 +---- src/elimination.jl | 12 ++++-------- src/global_identifiability.jl | 6 +++--- src/identifiable_functions.jl | 12 ++++++------ src/known_ic.jl | 13 +++++++++---- src/output_saturation.jl | 25 +++++++++++++++++++------ src/pb_representation.jl | 4 +++- src/primality_check.jl | 3 ++- src/states.jl | 10 ++++++++-- src/util.jl | 13 +++++++------ src/wronskian.jl | 4 +--- test/bodies/identifiable_functions.jl | 5 +++-- 13 files changed, 67 insertions(+), 50 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index dd36602e9..dcb1cb8ef 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -62,10 +62,7 @@ end #------------------------------------------------------------------------------ -function add_outputs( - ode::ODE{P}, - extra_y::Dict{String, <:RingElem}, - ) where {P <: MPolyRingElem} +function add_outputs(ode::ODE{P}, extra_y::Dict{String, <:RingElem}) where {P <: MPolyRingElem} new_var_names = vcat(collect(map(var_to_str, gens(ode.poly_ring))), collect(keys(extra_y))) new_ring, new_vars = Nemo.polynomial_ring(base_ring(ode.poly_ring), new_var_names) diff --git a/src/discrete.jl b/src/discrete.jl index 708d348e3..fd8e18551 100644 --- a/src/discrete.jl +++ b/src/discrete.jl @@ -57,10 +57,7 @@ end #------------------------------------------------------------------------------ # Some functions to transform DDS's -function add_outputs( - dds::DDS{P}, - extra_y::Dict{String, <:RingElem}, - ) where {P <: MPolyRingElem} +function add_outputs(dds::DDS{P}, extra_y::Dict{String, <:RingElem}) where {P <: MPolyRingElem} return DDS{P}(add_outputs(dds.ode, extra_y)) end diff --git a/src/elimination.jl b/src/elimination.jl index 0831b011c..bc631b92a 100644 --- a/src/elimination.jl +++ b/src/elimination.jl @@ -70,8 +70,8 @@ function Bezout_matrix(f::P, g::P, var_elim::P) where {P <: MPolyRingElem} for i in 1:n for j in 1:n M[i, j] = sum( - coeffs_f[j + k + 1] * coeffs_g[i - k] - - coeffs_g[j + k + 1] * coeffs_f[i - k] for k in 0:min(i - 1, n - j) + coeffs_f[j + k + 1] * coeffs_g[i - k] - coeffs_g[j + k + 1] * coeffs_f[i - k] for + k in 0:min(i - 1, n - j) ) end end @@ -188,10 +188,7 @@ mutable struct ODEPointGenerator{P} <: PointGenerator{P} cached_points::Array{Dict{P, <:FieldElem}, 1} number_type::Type - function ODEPointGenerator{P}( - ode::ODE{P}, - big_ring::MPolyRing, - ) where {P <: MPolyRingElem} + function ODEPointGenerator{P}(ode::ODE{P}, big_ring::MPolyRing) where {P <: MPolyRingElem} prec = length(ode.x_vars) + 1 number_type = typeof(one(base_ring(big_ring))) return new(ode, big_ring, prec, Array{Dict{P, number_type}}[], number_type) @@ -215,8 +212,7 @@ function Base.iterate( initial_conditions = Dict{P, Int}(x => rand(1:sample_max) for x in gpg.ode.x_vars) input_values = Dict{P, Array{Int, 1}}( - u => [rand(1:sample_max) for _ in 1:(gpg.precision)] for - u in gpg.ode.u_vars + u => [rand(1:sample_max) for _ in 1:(gpg.precision)] for u in gpg.ode.u_vars ) @debug "Computing a power series solution" ps_solution = undef diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl index 8d0a174e3..35ec1cf28 100644 --- a/src/global_identifiability.jl +++ b/src/global_identifiability.jl @@ -91,7 +91,7 @@ The function returns a tuple containing the following: with_states::Bool = false, var_change_policy = :default, rational_interpolator = :VanDerHoevenLecerf, - cmp = RationalFunctionFields.rational_function_cmp + cmp = RationalFunctionFields.rational_function_cmp, ) where {T} @info "Computing IO-equations" ioeq_time = @elapsed io_equations = @@ -149,7 +149,7 @@ The function returns a tuple containing the following: seed = 42, simplify = :standard, rational_interpolator = rational_interpolator, - cmp = cmp + cmp = cmp, ) dennums_simplified = fractions_to_dennums(no_states_simplified) # switch back the ring @@ -204,7 +204,7 @@ Output: a list L of booleans with L[i] being the identifiability status of the i known = known, prob_threshold = half_p, var_change_policy = var_change_policy, - with_states = states_needed + with_states = states_needed, ) funcs_to_check = Vector{Generic.FracFieldElem{P}}( diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index e6954629b..285ac0e65 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -58,7 +58,7 @@ function find_identifiable_functions( simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, loglevel = Logging.Info, - cmp = cmp_prefer_params(ode, RationalFunctionFields.rational_function_cmp) + cmp = cmp_prefer_params(ode, RationalFunctionFields.rational_function_cmp), ) where {T <: MPolyRingElem{QQFieldElem}} restart_logging(loglevel = loglevel) reset_timings() @@ -71,7 +71,7 @@ function find_identifiable_functions( with_states = with_states, simplify = simplify, rational_interpolator = rational_interpolator, - cmp = cmp + cmp = cmp, ) else id_funcs = _find_identifiable_functions_kic( @@ -81,7 +81,7 @@ function find_identifiable_functions( seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, - cmp = cmp + cmp = cmp, ) # renaming variables from `x(t)` to `x(0)` return replace_with_ic(ode, id_funcs) @@ -96,7 +96,7 @@ function _find_identifiable_functions( with_states = false, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, - cmp = cmp_prefer_params(ode, RationalFunctionFields.rational_function_cmp) + cmp = cmp_prefer_params(ode, RationalFunctionFields.rational_function_cmp), ) where {T <: MPolyRingElem{QQFieldElem}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) @@ -117,7 +117,7 @@ function _find_identifiable_functions( prob_threshold = half_p, with_states = with_states, rational_interpolator = rational_interpolator, - cmp = cmp + cmp = cmp, ) # If simplification is needed if simplify !== :absent @@ -131,7 +131,7 @@ function _find_identifiable_functions( seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, - cmp = cmp + cmp = cmp, ) else id_funcs_fracs = RationalFunctionFields.dennums_to_fractions(id_funcs) diff --git a/src/known_ic.jl b/src/known_ic.jl index bd0cccf08..6effd3f6d 100644 --- a/src/known_ic.jl +++ b/src/known_ic.jl @@ -33,7 +33,7 @@ function _find_identifiable_functions_kic( seed = 42, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, - cmp = RationalFunctionFields.rational_function_cmp + cmp = RationalFunctionFields.rational_function_cmp, ) where {T <: MPolyRingElem{Nemo.QQFieldElem}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) @@ -46,7 +46,7 @@ function _find_identifiable_functions_kic( simplify = :absent, rational_interpolator = rational_interpolator, seed = seed, - cmp = cmp + cmp = cmp, ) id_funcs = vcat(id_funcs_general, [f // one(parent(ode)) for f in known_ic]) @@ -58,7 +58,7 @@ function _find_identifiable_functions_kic( seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, - cmp = cmp + cmp = cmp, ) end @@ -94,7 +94,12 @@ function _assess_identifiability_kic( funcs_to_check = vcat(ode.x_vars, ode.parameters) end half_p = 0.5 + prob_threshold / 2 - id_funcs = _find_identifiable_functions_kic(ode, known_ic, prob_threshold = half_p, simplify = :absent) + id_funcs = _find_identifiable_functions_kic( + ode, + known_ic, + prob_threshold = half_p, + simplify = :absent, + ) #funcs_to_check = replace_with_ic(ode, funcs_to_check) result = OrderedDict(f => :globally for f in funcs_to_check) diff --git a/src/output_saturation.jl b/src/output_saturation.jl index 7a987a58c..0d0c277fb 100644 --- a/src/output_saturation.jl +++ b/src/output_saturation.jl @@ -17,18 +17,29 @@ function saturate_outputs( orders::Dict{P, Int}, max_deg = 5, ) where {P <: MPolyRingElem} - lie_ders = lie_derivatives_up_to(ode, Dict(ode.y_equations[y] => ord for (y, ord) in orders)) + lie_ders = + lie_derivatives_up_to(ode, Dict(ode.y_equations[y] => ord for (y, ord) in orders)) @info "Degrees: $([(total_degree(numerator(f)), total_degree(denominator(f))) for f in lie_ders])" - current_y = RationalFunctionField(vcat(collect(values(ode.y_equations)), ode.u_vars, ode.parameters)) + current_y = RationalFunctionField( + vcat(collect(values(ode.y_equations)), ode.u_vars, ode.parameters), + ) lie_ders = filter(f -> total_degree_frac(f) <= max_deg, lie_ders) sort!(lie_ders, lt = RationalFunctionFields.rational_function_cmp) for f in lie_ders if !first(RationalFunctionFields.check_algebraicity_modp(current_y, [f])) - current_y = RationalFunctionField(vcat(RationalFunctionFields.generators(current_y), [f])) + current_y = RationalFunctionField( + vcat(RationalFunctionFields.generators(current_y), [f]), + ) end end - new_outputs = RationalFunctionFields.generators(current_y)[(length(ode.y_vars) + length(ode.u_vars) + length(ode.parameters) + 1):end] + new_outputs = RationalFunctionFields.generators(current_y)[ + ( + length(ode.y_vars) + length( + ode.u_vars, + ) + length(ode.parameters) + 1 + ):end, + ] idx = 1 old_y_names = map(var_to_str, ode.y_vars) @@ -59,9 +70,11 @@ Inputs: Output: dictionary from outputs to orders """ function propose_orders(ode) - x_with_u = filter(x -> length(intersect(vars(ode.x_equations[x]), ode.u_vars)) > 0, ode.x_vars) + x_with_u = + filter(x -> length(intersect(vars(ode.x_equations[x]), ode.u_vars)) > 0, ode.x_vars) - graph = Dict(x => Set(intersect(vars(ode.x_equations[x]), ode.x_vars)) for x in ode.x_vars) + graph = + Dict(x => Set(intersect(vars(ode.x_equations[x]), ode.x_vars)) for x in ode.x_vars) result = Dict(y => 0 for y in ode.y_vars) for y in ode.y_vars diff --git a/src/pb_representation.jl b/src/pb_representation.jl index 96c2efa8d..756d11a10 100644 --- a/src/pb_representation.jl +++ b/src/pb_representation.jl @@ -109,7 +109,9 @@ function common_ring(poly::MPolyRingElem, pbr::PBRepresentation) "$(u)_$h" for h in 0:max( max_ords[u], - max_offset + max([difforder(p, u) for p in values(pbr.projections)]...), + max_offset + max( + [difforder(p, u) for p in values(pbr.projections)]..., + ), ) ], ) diff --git a/src/primality_check.jl b/src/primality_check.jl index 4bdaa7fb4..13d029366 100644 --- a/src/primality_check.jl +++ b/src/primality_check.jl @@ -26,7 +26,8 @@ function check_primality_zerodim(J::Array{QQMPolyRingElem, 1}) @debug "Prime is $p" F = Nemo.GF(p) S_F = Nemo.matrix_space(F, dim, dim) - generic_multiplication_modp = S_F([generic_multiplication[i, j] for i in 1:dim for j in 1:dim]) + generic_multiplication_modp = + S_F([generic_multiplication[i, j] for i in 1:dim for j in 1:dim]) R, t = Nemo.polynomial_ring(F, "t") chpoly = Nemo.charpoly(R, generic_multiplication_modp) @debug "Charpoly computed" diff --git a/src/states.jl b/src/states.jl index ae48e773d..2dd7d7775 100644 --- a/src/states.jl +++ b/src/states.jl @@ -69,7 +69,10 @@ end #------------------------------------------------------------------------------ -function lie_derivatives_up_to(ode::ODE, orders::Dict{P, Int}) where {P <: AbstractAlgebra.RingElem} +function lie_derivatives_up_to( + ode::ODE, + orders::Dict{P, Int}, + ) where {P <: AbstractAlgebra.RingElem} result = Array{Generic.FracFieldElem, 1}() for (f, ord) in orders curr = extract_coefficients_ratfunc(f, ode.u_vars) @@ -117,5 +120,8 @@ identifiable functions of parameters only end end - return lie_derivatives_up_to(ode, Dict(ode.y_equations[y] => ord for (y, ord) in y_to_ord)) + return lie_derivatives_up_to( + ode, + Dict(ode.y_equations[y] => ord for (y, ord) in y_to_ord), + ) end diff --git a/src/util.jl b/src/util.jl index 6f1b8389c..9490f83fa 100644 --- a/src/util.jl +++ b/src/util.jl @@ -303,11 +303,12 @@ end # ----------------------------------------------------------------------------- -# Given an ODE and a comparator for rational functions returns a new comparator, -# which prefers functions containing only the parameters of the given ODE. +# Given an ODE and a comparator for rational functions returns a new comparator, +# which prefers functions containing only the parameters of the given ODE. function cmp_prefer_params(ode, cmp) some_vars = ode.parameters - new_cmp = (f, g) -> begin + new_cmp = + (f, g) -> begin isempty(some_vars) && return cmp(f, g) some_vars = map(x -> parent_ring_change(x, parent(numerator(f))), some_vars) prefer_f = issubset(vars(f), some_vars) @@ -316,10 +317,10 @@ function cmp_prefer_params(ode, cmp) (prefer_f == prefer_g) && return cmp(f, g) false end - new_cmp + return new_cmp end -# Given an ODE and a comparator for rational functions returns a new comparator, +# Given an ODE and a comparator for rational functions returns a new comparator, # which first compares the Lie derivatives when the latter are not zero. function cmp_lie(ode, cmp) new_cmp = (f, g) -> begin @@ -333,7 +334,7 @@ function cmp_lie(ode, cmp) cmp(df, dg) end end - new_cmp + return new_cmp end # ----------------------------------------------------------------------------- diff --git a/src/wronskian.jl b/src/wronskian.jl index 244077783..b2e2a86e0 100644 --- a/src/wronskian.jl +++ b/src/wronskian.jl @@ -24,9 +24,7 @@ function monomial_compress(io_equation, params::Array{<:MPolyRingElem, 1}) params_xs = isempty(params) ? empty(params) : gens(parent(first(params))) # Pre-compute param string names for faster lookup param_names = Set(var_to_str(p, xs = params_xs) for p in params) - other_vars = [ - v for v in gens(parent(io_equation)) if !(var_to_str(v) in param_names) - ] + other_vars = [v for v in gens(parent(io_equation)) if !(var_to_str(v) in param_names)] coeffdict = extract_coefficients(io_equation, other_vars) expvect = collect(keys(coeffdict)) coeffs = collect(values(coeffdict)) diff --git a/test/bodies/identifiable_functions.jl b/test/bodies/identifiable_functions.jl index dc9febfa0..723abe4e6 100644 --- a/test/bodies/identifiable_functions.jl +++ b/test/bodies/identifiable_functions.jl @@ -1116,7 +1116,7 @@ ident_funcs_states = [ Vm * a1 * c + Vm * a2 * c, (a2 * k21) // (a1 + a2), x1 * c, - x2 * c * k21 + x2 * c * k21, ] push!(test_cases, (ode = ode, ident_funcs = ident_funcs, with_states = false)) push!(test_cases, (ode = ode, ident_funcs = ident_funcs_states, with_states = true)) @@ -1210,7 +1210,8 @@ end end cmp2 = StructuralIdentifiability.cmp_lie( - ode, RationalFunctionFields.rational_function_cmp) + ode, RationalFunctionFields.rational_function_cmp + ) cmp3 = StructuralIdentifiability.cmp_prefer_params(ode, cmp1) From f0b6a7b0d8d26a9cded6b7ee2df7f95cfb0c6778 Mon Sep 17 00:00:00 2001 From: Sasha Date: Wed, 29 Jul 2026 01:23:00 +0200 Subject: [PATCH 08/17] polish --- docs/src/tutorials/identifiable_functions.md | 11 +++++ src/StructuralIdentifiability.jl | 3 +- src/global_identifiability.jl | 2 +- src/identifiable_functions.jl | 13 +++--- src/known_ic.jl | 2 +- src/util.jl | 42 +++++++++++--------- test/bodies/identifiable_functions.jl | 6 ++- 7 files changed, 50 insertions(+), 29 deletions(-) diff --git a/docs/src/tutorials/identifiable_functions.md b/docs/src/tutorials/identifiable_functions.md index c161af2a5..75ddfa403 100644 --- a/docs/src/tutorials/identifiable_functions.md +++ b/docs/src/tutorials/identifiable_functions.md @@ -40,6 +40,17 @@ By default, `find_identifiable_functions` tries to simplify the output functions the degree of simplification. The default value is `:standard` but one could use `:strong` to try to simplify further (at the expense of heavier computation) or use `:weak` to simplify less (but compute faster). +When several generating sets describe the same field, the `cmp` keyword can be used to control which rational functions are considered simpler and end up in the final result. +It is possible to provide a custom function `cmp(f, g)` that returns `true` when `f` should be preferred to `g`. For example, the following comparator prefers functions with shorter printed representations: + +```@example funcs +find_identifiable_functions( + LLW1987, + with_states = true, + cmp = (f1, f2) -> length(string(f1)) < length(string(f2)), +) +``` + As `assess_identifiability` and `assess_local_identifiability`, `find_identifiable_functions` accepts an optional parameter `loglevel` (default: `Logging.Info`) to adjust the verbosity of logging. diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index 15d72355e..03c0daa20 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -33,7 +33,8 @@ using RationalFunctionFields: _reduce_mod_p, is_rational_func_const, fractions_to_dennums, - check_constructive_field_membership + check_constructive_field_membership, + rational_function_cmp # defining a model export ODE, @ODEmodel, @DDSmodel diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl index 35ec1cf28..f201386c0 100644 --- a/src/global_identifiability.jl +++ b/src/global_identifiability.jl @@ -91,7 +91,7 @@ The function returns a tuple containing the following: with_states::Bool = false, var_change_policy = :default, rational_interpolator = :VanDerHoevenLecerf, - cmp = RationalFunctionFields.rational_function_cmp, + cmp = default_cmp(ode), ) where {T} @info "Computing IO-equations" ioeq_time = @elapsed io_equations = diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index 285ac0e65..ed21486a9 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -6,7 +6,7 @@ system. ## Options -This functions takes the following optional arguments: +This function takes the following optional arguments: - `with_states`: When `true`, also reports the identifiabile functions in the ODE states. Default is `false`. - `simplify`: The extent to which the output functions are simplified. Stronger @@ -22,10 +22,9 @@ This functions takes the following optional arguments: initial conditions, not states (this is an experimental functionality). - `prob_threshold`: A float in the range from 0 to 1, the probability of correctness. Default is `0.99`. -- `cmp`: Comparator for rational functions. - A custom comparator can be provided: for two rational functions `f1` and `f2`, - `cmp(f1,f2)` should be `true` if `f1` is simpler than `f2`. - Default is `RationalFunctionFields.rational_function_cmp`. +- `cmp`: A comparator for rational functions. For two rational functions `f1` + and `f2`, `cmp(f1, f2)` should be `true` if `f1` is simpler than `f2`. + By default, `StructuralIdentifiability.default_cmp(ode)` is used. - `seed`: The rng seed. Default value is `42`. - `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) @@ -58,7 +57,7 @@ function find_identifiable_functions( simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, loglevel = Logging.Info, - cmp = cmp_prefer_params(ode, RationalFunctionFields.rational_function_cmp), + cmp = default_cmp(ode), ) where {T <: MPolyRingElem{QQFieldElem}} restart_logging(loglevel = loglevel) reset_timings() @@ -96,7 +95,7 @@ function _find_identifiable_functions( with_states = false, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, - cmp = cmp_prefer_params(ode, RationalFunctionFields.rational_function_cmp), + cmp = default_cmp(ode), ) where {T <: MPolyRingElem{QQFieldElem}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) diff --git a/src/known_ic.jl b/src/known_ic.jl index 6effd3f6d..30934ae1c 100644 --- a/src/known_ic.jl +++ b/src/known_ic.jl @@ -33,7 +33,7 @@ function _find_identifiable_functions_kic( seed = 42, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, - cmp = RationalFunctionFields.rational_function_cmp, + cmp = default_cmp(ode), ) where {T <: MPolyRingElem{Nemo.QQFieldElem}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) diff --git a/src/util.jl b/src/util.jl index 9490f83fa..1dbd80a59 100644 --- a/src/util.jl +++ b/src/util.jl @@ -303,38 +303,44 @@ end # ----------------------------------------------------------------------------- -# Given an ODE and a comparator for rational functions returns a new comparator, -# which prefers functions containing only the parameters of the given ODE. +default_cmp(ode) = cmp_prefer_params(ode, rational_function_cmp) + +""" + cmp_prefer_params(ode, cmp) + +Wrap a comparator of rational functions `cmp` so that functions containing +only parameters of the `ode` are preferred over functions containing states. +""" function cmp_prefer_params(ode, cmp) - some_vars = ode.parameters - new_cmp = - (f, g) -> begin - isempty(some_vars) && return cmp(f, g) - some_vars = map(x -> parent_ring_change(x, parent(numerator(f))), some_vars) - prefer_f = issubset(vars(f), some_vars) - prefer_g = issubset(vars(g), some_vars) - (prefer_f && !prefer_g) && return true - (prefer_f == prefer_g) && return cmp(f, g) - false + params = ode.parameters + return function (f, g) + isempty(params) && return cmp(f, g) + params = map(p -> parent_ring_change(p, parent(numerator(f))), params) + prefer_f = issubset(vars(f), params) + prefer_g = issubset(vars(g), params) + prefer_f != prefer_g && return prefer_f + return cmp(f, g) end - return new_cmp end -# Given an ODE and a comparator for rational functions returns a new comparator, -# which first compares the Lie derivatives when the latter are not zero. +""" + cmp_lie(ode, cmp) + +Wrap a comparator of rational functions `cmp` so Lie derivatives are compared +when they are nonzero. +""" function cmp_lie(ode, cmp) - new_cmp = (f, g) -> begin + return function (f, g) f = parent_ring_change(f, parent(ode)) g = parent_ring_change(g, parent(ode)) df = lie_derivative(f, ode) dg = lie_derivative(g, ode) - if iszero(df) && iszero(dg) + return if iszero(df) && iszero(dg) cmp(f, g) else cmp(df, dg) end end - return new_cmp end # ----------------------------------------------------------------------------- diff --git a/test/bodies/identifiable_functions.jl b/test/bodies/identifiable_functions.jl index 723abe4e6..237e4af7c 100644 --- a/test/bodies/identifiable_functions.jl +++ b/test/bodies/identifiable_functions.jl @@ -1206,7 +1206,7 @@ end ) cmp1 = (f, g) -> begin - length(string(f)) < length(string(g)) + return length(string(f)) < length(string(g)) end cmp2 = StructuralIdentifiability.cmp_lie( @@ -1215,6 +1215,10 @@ end cmp3 = StructuralIdentifiability.cmp_prefer_params(ode, cmp1) + o = one(parent(ode)) + @test cmp3(beta_W // o, S // o) + @test !cmp3(S // o, beta_W // o) + with_cmp1 = find_identifiable_functions(ode, with_states = true, cmp = cmp1) with_cmp2 = find_identifiable_functions(ode, with_states = true, cmp = cmp2) From 04a7cdf2bc29f41ba8982bac07b495d4947604c6 Mon Sep 17 00:00:00 2001 From: Sasha Date: Wed, 29 Jul 2026 01:26:53 +0200 Subject: [PATCH 09/17] polish ext and docs --- docs/src/tutorials/identifiable_functions.md | 3 +-- ext/ModelingToolkitSIExt.jl | 16 ++++++++-------- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/docs/src/tutorials/identifiable_functions.md b/docs/src/tutorials/identifiable_functions.md index 75ddfa403..b297e6336 100644 --- a/docs/src/tutorials/identifiable_functions.md +++ b/docs/src/tutorials/identifiable_functions.md @@ -40,8 +40,7 @@ By default, `find_identifiable_functions` tries to simplify the output functions the degree of simplification. The default value is `:standard` but one could use `:strong` to try to simplify further (at the expense of heavier computation) or use `:weak` to simplify less (but compute faster). -When several generating sets describe the same field, the `cmp` keyword can be used to control which rational functions are considered simpler and end up in the final result. -It is possible to provide a custom function `cmp(f, g)` that returns `true` when `f` should be preferred to `g`. For example, the following comparator prefers functions with shorter printed representations: +When several generating sets describe the same field, the `cmp` keyword can be used to control which rational functions are considered simpler and end up in the final result. It is possible to provide a custom function `cmp(f, g)` that returns `true` when `f` should be preferred to `g`. For example, the following comparator prefers functions with shorter printed representations: ```@example funcs find_identifiable_functions( diff --git a/ext/ModelingToolkitSIExt.jl b/ext/ModelingToolkitSIExt.jl index 44022a99f..3d0e4d003 100644 --- a/ext/ModelingToolkitSIExt.jl +++ b/ext/ModelingToolkitSIExt.jl @@ -12,8 +12,6 @@ using TimerOutputs using ModelingToolkitBase -RationalFunctionFields = StructuralIdentifiability.RationalFunctionFields - # ------------------------------------------------------------------------------ # checking if it is a function of the form x(t), a bit dirty @@ -631,11 +629,12 @@ system. ## Options -This functions takes the following optional arguments: +This function takes the following optional arguments: - `measured_quantities` - the output functions of the model. - `known_ic` - a list of functions whose initial conditions are assumed to be known, then the returned identifiable functions will be functions of parameters and initial conditions, not states (this is an experimental functionality). +- `cmp` - a comparator for rational functions. - `loglevel` - the verbosity of the logging (can be Logging.Error, Logging.Warn, Logging.Info, Logging.Debug) @@ -673,7 +672,7 @@ function StructuralIdentifiability.find_identifiable_functions( simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, loglevel = Logging.Info, - cmp = RationalFunctionFields.rational_function_cmp + cmp = nothing, ) restart_logging(loglevel = loglevel) reset_timings() @@ -687,7 +686,7 @@ function StructuralIdentifiability.find_identifiable_functions( with_states = with_states, simplify = simplify, rational_interpolator = rational_interpolator, - cmp = cmp + cmp = cmp, ) end end @@ -701,10 +700,11 @@ function _find_identifiable_functions( with_states = false, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, - cmp = RationalFunctionFields.rational_function_cmp + cmp = nothing, ) Random.seed!(seed) ode, conversion = mtk_to_si(ode, measured_quantities) + isnothing(cmp) && (cmp = StructuralIdentifiability.default_cmp(ode)) known_ic_ = [eval_at_nemo(each, conversion) for each in known_ic] result = nothing if isempty(known_ic) @@ -715,7 +715,7 @@ function _find_identifiable_functions( seed = seed, with_states = with_states, rational_interpolator = rational_interpolator, - cmp = cmp + cmp = cmp, ) else result = StructuralIdentifiability._find_identifiable_functions_kic( @@ -725,7 +725,7 @@ function _find_identifiable_functions( prob_threshold = prob_threshold, seed = seed, rational_interpolator = rational_interpolator, - cmp = cmp + cmp = cmp, ) end result = [parent_ring_change(f, ode.poly_ring) for f in result] From 8dd5208abf25fa93a2c32296925c24ecc538766f Mon Sep 17 00:00:00 2001 From: Sasha Date: Sun, 2 Aug 2026 02:44:18 +0200 Subject: [PATCH 10/17] clean benchmarking scripts --- benchmarking/IdentifiableFunctions/.gitignore | 1 + .../IdentifiableFunctions/benchmark_result.md | 156 +++---- .../IdentifiableFunctions/run_benchmarks.jl | 383 ++++++++---------- .../run_single_benchmark.jl | 133 +++--- benchmarking/IdentifiableFunctions/utils.jl | 131 ++---- benchmarking/benchmark_result.md | 54 +-- 6 files changed, 366 insertions(+), 492 deletions(-) diff --git a/benchmarking/IdentifiableFunctions/.gitignore b/benchmarking/IdentifiableFunctions/.gitignore index fbca22537..cd07ff2f3 100644 --- a/benchmarking/IdentifiableFunctions/.gitignore +++ b/benchmarking/IdentifiableFunctions/.gitignore @@ -1 +1,2 @@ results/ +logs/ diff --git a/benchmarking/IdentifiableFunctions/benchmark_result.md b/benchmarking/IdentifiableFunctions/benchmark_result.md index 9e46740c0..9a8cdbc9f 100644 --- a/benchmarking/IdentifiableFunctions/benchmark_result.md +++ b/benchmarking/IdentifiableFunctions/benchmark_result.md @@ -1,84 +1,92 @@ ## Benchmark results -2023-09-15T02:10:37.348 +2026-07-28T01:16:05.923 - - Benchmarked function: `find_identifiable_functions` - - Workers: 7 - - Timeout: 600 s +- Benchmarked function: `find_identifiable_functions` +- Workers: 30 +- Timeout: 7200 s +- Memory limit: 28192 MiB -**All timings in seconds.** - -| Model | (:normalforms, 2)_with_states / Runtime | (:normalforms, 2)_with_states / beautifulization | -|:----------------------------- |:--------------------------------------- |:------------------------------------------------ | -| Akt pathway | 2.44 | 0.28 | -| Bilirubin2_io | 4.08 | 0.76 | -| Biohydrogenation_io | 0.44 | 0.11 | -| Bruno2016 | 0.19 | 0.07 | -| CD8 T cell differentiation | 0.72 | 0.27 | -| CGV1990 | 2.25 | 0.31 | -| Chemical reaction network | 0.36 | 0.13 | -| Crauste_SI | 0.47 | 0.18 | -| Fujita | 2.50 | 0.30 | -| Goodwin oscillator | 0.29 | 0.09 | -| HIV | 0.78 | 0.19 | -| HIV2_io | 2.00 | 0.27 | -| HighDimNonLin | 20.55 | 0.68 | -| KD1999 | 0.97 | 0.14 | -| LLW1987_io | 0.25 | 0.08 | -| MAPK model (5 outputs) | 48.45 | 0.37 | -| MAPK model (6 outputs) | 7.99 | 0.42 | -| Modified LV for testing | 0.07 | 0.01 | -| PK1 | 0.74 | 0.18 | -| PK2 | 130.96 | 0.07 | -| Pharm | 152.26 | 0.09 | -| Phosphorylation | 0.45 | 0.11 | -| Pivastatin | 8.44 | 0.04 | -| QY | 49.21 | 0.23 | -| Ruminal lipolysis | 0.21 | 0.07 | -| SEAIJRC Covid model | 85.13 | 0.11 | -| SEIR 34 | 0.39 | 0.09 | -| SEIR 36 ref | 0.78 | 0.14 | -| SEIR2T | 0.23 | 0.09 | -| SEIRT | 0.36 | 0.08 | -| SEIR_1_io | 0.56 | 0.10 | -| SEUIR | - | - | -| SIR 19 | - | - | -| SIR 21 | 0.30 | 0.19 | -| SIR 24 | 0.36 | 0.07 | -| SIR 6 | 0.20 | 0.08 | -| SIRS forced | 10.62 | 0.40 | -| SIWR original | 13.72 | 0.07 | -| SIWR with extra output | 0.67 | 0.10 | -| SLIQR | 1.42 | 0.15 | -| St | 37.54 | 1.38 | -| Transfection_4State | 0.23 | 0.07 | -| Treatment_io | 0.44 | 0.14 | -| cLV1 (2o) | 0.80 | 0.20 | -| cholera | 0.66 | 0.09 | -| generalizedLoktaVolterra (1o) | 0.18 | 0.05 | -| generalizedLoktaVolterra (2o) | 0.09 | 0.07 | -| p53 | 1.21 | 0.34 | +|Model|funcs (cmp_default)|funcs (cmp_lie)| +|---|---|---| +|Akt pathway|[reaction_8_k1, reaction_7_k1, reaction_6_k1, reaction_5_k2, reaction_4_k1, reaction_3_k1, reaction_2_k2, reaction_1_k1 - reaction_1_k2 - reaction_9_k1, reaction_2_k1//reaction_5_k1, a3//reaction_5_k1, a2//reaction_5_k1, a1//reaction_5_k1, pS6(t)*reaction_5_k1, pAkt_S6(t)*reaction_5_k1, S6(t)*reaction_5_k1, pAkt(t)*reaction_5_k1, Akt(t)*reaction_5_k1, pEGFR_Akt(t)*reaction_5_k1, pEGFR(t)*reaction_5_k1, EGF_EGFR(t)*reaction_5_k1*reaction_9_k1]|[reaction_8_k1, reaction_7_k1, reaction_6_k1, reaction_5_k2, reaction_4_k1, reaction_3_k1, reaction_2_k2, reaction_1_k1 - reaction_1_k2 - reaction_9_k1, reaction_2_k1//reaction_5_k1, a3//reaction_5_k1, a2//reaction_5_k1, a1//reaction_5_k1, pS6(t)*reaction_5_k1, pAkt(t)*a2 + pAkt_S6(t)*a2, pEGFR(t)*a1 + pEGFR_Akt(t)*a1, EGF_EGFR(t)*reaction_5_k1*reaction_9_k1, S6(t)*reaction_5_k1, pAkt_S6(t)*reaction_5_k1, Akt(t)*reaction_5_k1, pEGFR_Akt(t)*reaction_5_k1]| +|Bilirubin2_io|[k01, k21 + k31 + k41, k12 + k13 + k14, k21*k31 + k21*k41 + k31*k41, k12*k13 + k12*k14 + k13*k14, k12*k31 + k12*k41 + k13*k21 + k13*k41 + k14*k21 + k14*k31, k21*k31*k41, k12*k13*k14, x1(t), x2(t) + x3(t) + x4(t), x2(t)*x3(t) + x2(t)*x4(t) + x3(t)*x4(t), x2(t)*k31 + x2(t)*k41 + x3(t)*k21 + x3(t)*k41 + x4(t)*k21 + x4(t)*k31, x2(t)*k13 + x2(t)*k14 + x3(t)*k12 + x3(t)*k14 + x4(t)*k12 + x4(t)*k13]|[k01, k21 + k31 + k41, k12 + k13 + k14, k21*k31 + k21*k41 + k31*k41, k12*k13 + k12*k14 + k13*k14, k12*k31 + k12*k41 + k13*k21 + k13*k41 + k14*k21 + k14*k31, k21*k31*k41, k12*k13*k14, x2(t) + x3(t) + x4(t), x1(t), x2(t)*k31 + x2(t)*k41 + x3(t)*k21 + x3(t)*k41 + x4(t)*k21 + x4(t)*k31, x2(t)*k13 + x2(t)*k14 + x3(t)*k12 + x3(t)*k14 + x4(t)*k12 + x4(t)*k13]| +|Biohydrogenation_io|[k7, k6, k5, k10 + 2*k8, k9^2, k10*k9, x5(t), x4(t), x6(t) + k8]|[k7, k6, k5, k10 + 2*k8, k9^2, k10*k9, x4(t), x5(t), x6(t) - 1//2*k10]| +|Bruno2016|[kbeta10, kbeta, kcryOH + kcrybeta, beta10(t), beta(t), cry(t)*kcryOH]|[kbeta10, kbeta, kcryOH + kcrybeta, beta(t), beta10(t), cry(t)*kcryOH]| +|CD8 T cell differentiation|[mu_PL, mu_PE, mu_P, mu_N, mu_M, mu_LL, mu_LE, mu_EE, delta_LM, delta_EL, rho_E//rho_P, delta_NE//rho_P, M(t), S(t), E(t), N(t), P(t)*rho_P]|[mu_PL, mu_PE, mu_P, mu_N, mu_M, mu_LL, mu_LE, mu_EE, delta_LM, delta_EL, rho_E//rho_P, delta_NE//rho_P, M(t), N(t), S(t), E(t), P(t)*rho_P]| +|CGV1990|[k7, k6, k4, k3, R*V36 + 1//25*S*V3, R*V3 + S*V36, R*V36*k5 + 1//5*S*V36*k5, (V3*k5 + 5*V36*k5)//V3, q7(t), q3(t), q1(t), q35(t) + q36(t), q35(t)*q36(t), q35(t)*S*V36 + q36(t)*R*V3]|[k7, k6, k4, k3, R*V36 + 1//25*S*V3, R*V3 + S*V36, R*V36*k5 + 1//5*S*V36*k5, (V3*k5 + 5*V36*k5)//V3, q7(t), q1(t), q35(t)*S*V36 + q36(t)*R*V3, q35(t)*R*V36 + 1//25*q36(t)*S*V3, q3(t)]| +|Chemical reaction network|[k6, k5, k4, k3, k2, k1, x6(t), x5(t), x4(t), x3(t), x2(t), x1(t)]|[k6, k5, k4, k3, k2, k1, x6(t), x5(t), x3(t), x1(t), x2(t), x4(t)]| +|Crauste_SI|[mu_PL, mu_PE, mu_P, mu_N, mu_M, mu_LL, mu_LE, mu_EE, delta_LM, delta_EL, rho_E//rho_P, delta_NE//rho_P, M(t), S(t), E(t), N(t), P(t)*rho_P]|[mu_PL, mu_PE, mu_P, mu_N, mu_M, mu_LL, mu_LE, mu_EE, delta_LM, delta_EL, rho_E//rho_P, delta_NE//rho_P, M(t), N(t), S(t), E(t), P(t)*rho_P]| +|Fujita|[reaction_8_k1, reaction_7_k1, reaction_6_k1, reaction_5_k2, reaction_4_k1, reaction_3_k1, reaction_2_k2, reaction_1_k1 - reaction_1_k2 - reaction_9_k1, reaction_2_k1//reaction_5_k1, a3//reaction_5_k1, a2//reaction_5_k1, a1//reaction_5_k1, pS6(t)*reaction_5_k1, pAkt_S6(t)*reaction_5_k1, S6(t)*reaction_5_k1, pAkt(t)*reaction_5_k1, Akt(t)*reaction_5_k1, pEGFR_Akt(t)*reaction_5_k1, pEGFR(t)*reaction_5_k1, EGF_EGFR(t)*reaction_5_k1*reaction_9_k1]|[reaction_8_k1, reaction_7_k1, reaction_6_k1, reaction_5_k2, reaction_4_k1, reaction_3_k1, reaction_2_k2, reaction_1_k1 - reaction_1_k2 - reaction_9_k1, reaction_2_k1//reaction_5_k1, a3//reaction_5_k1, a2//reaction_5_k1, a1//reaction_5_k1, pS6(t)*reaction_5_k1, pAkt(t)*a2 + pAkt_S6(t)*a2, pEGFR(t)*a1 + pEGFR_Akt(t)*a1, EGF_EGFR(t)*reaction_5_k1*reaction_9_k1, S6(t)*reaction_5_k1, pAkt_S6(t)*reaction_5_k1, Akt(t)*reaction_5_k1, pEGFR_Akt(t)*reaction_5_k1]| +|Goodwin oscillator|[sigma, c, b, beta + delta, beta*delta, x4(t), x1(t), (alpha*gama)//x3(t), (x2(t)*gama - x3(t)*delta)//x3(t)]|[sigma, c, b, beta + delta, beta*delta, x1(t), x4(t), (x2(t)*gama + 1//2*x3(t)*beta - 1//2*x3(t)*delta)//(alpha*gama), (alpha*gama)//x3(t)]| +|HIV|[u, h, d, b, a, lm//q, c*q^2, beta*k*q, z(t), w(t), v(t)*beta, y(t)//q, x(t)//q]|[u, h, d, b, a, lm//q, c*q^2, beta*k*q, v(t)*beta, z(t), y(t)//q, w(t), x(t)//q]| +|HIV2_io|[s, d, b, c + k1 + w1 + w2, k2*q2, c*k1 + c*w1 + c*w2 + k1*w2 + w1*w2, c*k1*w2 + c*w1*w2, k1*k2*q1 + k1*k2*q2 + k2*q2*w1, x4(t), x1(t), x3(t)*k2 - x4(t)*c, x2(t)*k1*k2 + x3(t)*k1*k2 + x3(t)*k2*w1 + x4(t)*k1*w2 + x4(t)*w1*w2]|[s, d, b, c + k1 + w1 + w2, k2*q2, c*k1 + c*w1 + c*w2 + k1*w2 + w1*w2, c*k1*w2 + c*w1*w2, k1*k2*q1 + k1*k2*q2 + k2*q2*w1, x4(t), x1(t), x3(t)*k2 - x4(t)*c, x2(t)*k1*k2 + x3(t)*k1*k2 + x3(t)*k2*w1 + x4(t)*k1*w2 + x4(t)*w1*w2]| +|HighDimNonLin|[vm, p9, p8, p7, p6, p5, p4, p3, p20, p2, p19, p18, p17, p16, p15, p14, p13, p12, p11, p10, p1, km, x20(t), x19(t), x18(t), x17(t), x16(t), x15(t), x14(t), x13(t), x12(t), x11(t), x10(t), x9(t), x8(t), x7(t), x6(t), x5(t), x4(t), x3(t), x2(t), x1(t)]|[vm, p9, p8, p7, p6, p5, p4, p3, p20, p2, p19, p18, p17, p16, p15, p14, p13, p12, p11, p10, p1, km, x20(t), x19(t), x18(t), x17(t), x16(t), x15(t), x14(t), x13(t), x12(t), x11(t), x10(t), x9(t), x8(t), x7(t), x6(t), x5(t), x4(t), x3(t), x2(t), x1(t)]| +|JAK-STAT 1|[t9, t8, t7, t6, t5, t4, t3, t20, t2, t19, t18, t16, t14, t13, t12, t10, t1, t17*t22, t15*t21, t11*t21, x10(t), x9(t), x7(t), x6(t), x5(t), x4(t), x3(t), x2(t), x1(t), x8(t)*t21]|[t9, t8, t7, t6, t5, t4, t3, t20, t2, t19, t18, t16, t14, t13, t12, t10, t1, t17*t22, t15*t21, t11*t21, x5(t), x2(t), x1(t) + x3(t) + x4(t), x9(t), (x8(t)*t17*t22)//t11, x4(t), x3(t), x7(t), x10(t), x6(t)]| +|KD1999|[Vh, V, Th, Ta, Ca0, DH//UA, E//R, (cph*roh)//UA, (cp*ro)//UA, Tj(t), T(t), Cb(t), Ca(t), Arr(t)*k0]|[Vh, V, Th, Ta, Ca0, DH//UA, E//R, (cph*roh)//UA, (cp*ro)//UA, Cb(t), Ca(t), Tj(t), T(t), Arr(t)*k0]| +|LLW1987_io|[p1 + p3, p2*p4, p1*p3, x3(t), x1(t)*x2(t), x1(t)*p4 + x2(t)*p2, x1(t)*p1*p4 + x2(t)*p2*p3]|[p1 + p3, p2*p4, p1*p3, x1(t)*p4 + x2(t)*p2, x3(t), x1(t)*x2(t), (x1(t)*p4 - x2(t)*p2)//(p1 - p3)]| +|LeukaemiaLeon2021|-|-| +|MAPK model (5 outputs bis)|-|-| +|MAPK model (5 outputs)|[gamma1110, gamma1101, gamma1100, gamma1000, gamma0100, c1011, c0111, c0011, c0010, c0001, beta11, beta10, beta01, b10, b01, b00, alpha11, alpha10, alpha01, a10, a01, a00, S11(t), S10(t), S01(t), S00(t), F(t), K(t), FS11(t), FS10(t), FS01(t), KS10(t), KS01(t), KS00(t)]|[gamma1110, gamma1101, gamma1100, gamma1000, gamma0100, c1011, c0111, c0011, c0010, c0001, beta11, beta10, beta01, b10, b01, b00, alpha11, alpha10, alpha01, a10, a01, a00, S00(t), F(t), K(t), FS11(t), S01(t), FS01(t), KS00(t), FS10(t), KS10(t), KS01(t), S11(t), S10(t)]| +|MAPK model (6 outputs)|[gamma1110, gamma1101, gamma1100, gamma1000, gamma0100, c1011, c0111, c0011, c0010, c0001, beta11, beta10, beta01, b10, b01, b00, alpha11, alpha10, alpha01, a10, a01, a00, S11(t), S10(t), S01(t), S00(t), F(t), K(t), FS11(t), FS10(t), FS01(t), KS10(t), KS01(t), KS00(t)]|[gamma1110, gamma1101, gamma1100, gamma1000, gamma0100, c1011, c0111, c0011, c0010, c0001, beta11, beta10, beta01, b10, b01, b00, alpha11, alpha10, alpha01, a10, a01, a00, S00(t), F(t), K(t), FS11(t), S01(t), FS01(t), KS00(t), FS10(t), KS10(t), KS01(t), S11(t), S10(t)]| +|Modified LV for testing|[d, a + b, a*b, x1(t), x2(t)*c]|[d, a + b, a*b, x1(t), x2(t)*c]| +|NFkB|-|-| +|PK1|[k6, k5, k4, k3 + k7, k1 + k2, k2*s3, k1*s2, k1*k3*s3, x1(t), x4(t)*s2, x3(t)*s3, x2(t)*s2]|[k6, k5, k4, k3 + k7, k1 + k2, k2*s3, k1*s2, k1*k3*s3, x1(t), x4(t)*s2, x3(t)*s3, x2(t)*s2]| +|PK2|[n, kc, ka, b2, b1, a2, a1, x3(t), x2(t), x1(t), x0(t)]|[n, kc, ka, b2, b1, a2, a1, x3(t), x1(t), x2(t), x0(t)]| +|Pharm|[n, kc, ka, b2, b1, a2, a1, x3(t), x2(t), x1(t), x0(t)]|[n, kc, ka, b2, b1, a2, a1, x3(t), x1(t), x2(t), x0(t)]| +|Phosphorylation|[k6, k5, k4, k3, k2, k1, x3(t), x1(t), x2(t), x4(t), x6(t), x5(t)]|[k6, k5, k4, k3, k2, k1, x4(t), x2(t), x1(t), x3(t), x6(t), x5(t)]| +|Pivastatin|[r3, r1, k4, k3, k2, T0*k1, T0*k, x3(t)*k1, x2(t)*k1, x1(t)*k1]|[r3, r1, k4, k3, k2, T0*k1, T0*k, x3(t)*k1, x2(t)*k1, x1(t)*k1]| +|QWWC|-|-| +|QY|[alpa, Mar, Ks, M + siga1 + siga2, phi//(Mar - siga1), M*siga1 + M*siga2 + siga1*siga2, M*siga1*siga2, (M*phi*siga2 - M*siga2)//beta_SA, (Ks*beta_SA + M*beta_SA - 1//3*Mar*beta_SA*phi + 1//3*Mar*beta_SA + 1//3*beta_SA*phi*siga2 + 2//3*beta_SA*siga1 + 2//3*beta_SA*siga2 - 1//3*beta_SI*phi*siga2 + 1//3*beta_SI*siga2)//beta_SA, P4(t), P3(t), P2(t), P1(t), P0(t), P5(t)*Mar + beta]|[alpa, Mar, Ks, M + siga1 + siga2, phi//(Mar - siga1), M*siga1 + M*siga2 + siga1*siga2, M*siga1*siga2, (M*phi*siga2 - M*siga2)//beta_SA, (Ks*beta_SA + M*beta_SA - 1//3*Mar*beta_SA*phi + 1//3*Mar*beta_SA + 1//3*beta_SA*phi*siga2 + 2//3*beta_SA*siga1 + 2//3*beta_SA*siga2 - 1//3*beta_SI*phi*siga2 + 1//3*beta_SI*siga2)//beta_SA, P3(t), P2(t), P1(t), P0(t), (P5(t)*Mar + beta)//Mar, P4(t)]| +|Ruminal lipolysis|[k4, k3, k2, x5(t), x4(t), x3(t), x2(t), x1(t)]|[k4, k3, k2, x5(t), x3(t), x1(t), x2(t), x4(t)]| +|SEAIJRC Covid model|[k, g2, g1, b, alpha, (q*r - q)//r, Ninv(t), C(t), J(t), I(t), A(t)*q, E(t)*r, S(t)*r]|[k, g2, g1, b, alpha, Ninv(t), (q*r - q)//r, C(t), J(t), I(t), A(t)*q, S(t)*r, E(t)*r]| +|SEIR 34|[mu, epsilon + gamma, epsilon*gamma, K*epsilon, beta*epsilon*r, S(t), N(t), A(t), E(t) + I(t), I(t)*gamma]|[mu, N(t), A(t), epsilon + gamma, epsilon*gamma, K*epsilon, beta*epsilon*r, I(t)*gamma, S(t), E(t) + I(t)]| +|SEIR 36 ref|[s_d, s, phi_e, phi, mu_i, mu_d, mu_0, gamma_d, gamma, beta_d, beta, F(t), Di(t), De(t), I(t), E(t), S(t), q(t), nu(t), N(t)]|[s_d, s, phi_e, phi, mu_i, mu_d, mu_0, gamma_d, gamma, beta_d, beta, q(t), nu(t), N(t), F(t), De(t), Di(t), I(t), S(t), E(t)]| +|SEIR2T|[nu, b, a, Cu(t), N(t), In(t), E(t), S(t)]|[nu, b, a, N(t), Cu(t), In(t), S(t), E(t)]| +|SEIRT|[beta, alpha + lambda, alpha*lambda, N(t), I(t), S(t)*alpha, E(t)*alpha - I(t)*lambda]|[beta, N(t), alpha + lambda, alpha*lambda, I(t), S(t)*alpha, (E(t) + I(t))//lambda]| +|SEIR_1_io|[gamma, beta//psi, gamma*psi - psi*v, psi*v - psi - v, Q(t), I(t)*psi, S(t)*psi*v, E(t)*psi*v + I(t)*psi*v]|[gamma, beta//psi, gamma*psi - psi*v, psi*v - psi - v, Q(t), I(t)*psi, S(t)*psi*v, E(t)*psi*v + I(t)*psi*v]| +|SEUIR|[z, d, (N*w)//beta, I(t), E(t)*w, S(t)*w, U(t)*w + I(t)*w]|[z, d, (N*w)//beta, I(t), U(t)*w + I(t)*w, S(t)*w, E(t)*w]| +|SIR 19|[r, q, pp, mu, beta, D(t), C(t), I(t), S(t), N(t)]|[r, q, pp, mu, beta, N(t), D(t), C(t), I(t), S(t)]| +|SIR 21|[r, q, pp, mu, beta, D(t), C(t), I(t), S(t), N(t)]|[r, q, pp, mu, beta, N(t), D(t), C(t), I(t), S(t)]| +|SIR 24|[K, gamma + mu, c*phi, I(t), S(t), A(t) - mu]|[K, gamma + mu, A(t) - mu, c*phi, I(t), S(t)]| +|SIR 6|[gamma, K//beta, N(t), I(t)*beta, S(t)*beta]|[gamma, N(t), K//beta, S(t)*beta, I(t)*beta]| +|SIRS forced|[nu, mu, g, b0, M^2, r(t), i(t), s(t), x1(t)*b1, x2(t)*M*b1]|[nu, mu, g, b0, M^2, (x2(t)*b1)//M, r(t), x1(t)*b1, i(t), s(t)]| +|SIWR original|[xi, mu, k, gam, bw, bi, a, R(t), W(t), I(t), S(t)]|[xi, mu, k, gam, bw, bi, a, W(t), R(t), I(t), S(t)]| +|SIWR with extra output|[xi, mu, k, gam, bw, bi, a, R(t), W(t), I(t), S(t)]|[xi, mu, k, gam, bw, bi, a, W(t), R(t), S(t) + I(t) + R(t), I(t)]| +|SLIQR|[s, b, Ninv, a + g, a*e*g, a*g + e*g*s - g*s, In(t), S(t)*a, Q(t)*a - Q(t)*s, L(t)*a - In(t)*g + Q(t)*s]|[s, b, Ninv, a + g, a*e*g, a*g + e*g*s - g*s, In(t), Q(t)*a - Q(t)*s, S(t)*a, L(t)*a - In(t)*g + Q(t)*s]| +|St|[T, Dd, T*a + T*d + T*dr + e + g - r - rR, (d*rR - dr*r)//(d - dr), (a*r - a*rR + d*e + d*g - dr*e - dr*g)//(d - dr), (a^2 + 2*a*d + d^2 + dr^2)//(a*dr + d*dr), (a*dr*r - a*dr*rR + d^2*g + d*dr*e - d*dr*g - dr^2*e)//(a*d - a*dr + d^2 - dr^2), S(t) + R(t), W(t)*a + W(t)*d + W(t)*dr + e + g - r - rR, S(t)*T*d - S(t)*r + R(t)*T*dr - R(t)*rR]|[T, Dd, T*a + T*d + T*dr + e + g - r - rR, (d*rR - dr*r)//(d - dr), (a*r - a*rR + d*e + d*g - dr*e - dr*g)//(d - dr), (a^2 + 2*a*d + d^2 + dr^2)//(a*dr + d*dr), (a*dr*r - a*dr*rR + d^2*g + d*dr*e - d*dr*g - dr^2*e)//(a*d - a*dr + d^2 - dr^2), S(t) + R(t), W(t)*a + W(t)*d + W(t)*dr + e + g - r - rR, (S(t)*d + R(t)*dr)//(a + d + dr)]| +|Transfection_4State|[d3, d1, b, d2//kTL, GFP(t), mRNAenz(t)*kTL, enz(t)*kTL, mRNA(t)*kTL]|[d3, d1, b, d2//kTL, GFP(t), mRNAenz(t)*kTL, enz(t)*kTL, mRNA(t)*kTL]| +|Treatment_io|[a + g + nu, b//g, d*g + nu, a*nu + g*nu, N(t), Tr(t), S(t)*g, In(t)*g - Tr(t)*nu]|[N(t), a + g + nu, b//g, d*g + nu, a*nu + g*nu, Tr(t), S(t)*g, In(t)*g - Tr(t)*nu]| +|TumorHu2019|-|-| +|TumorPillis2007|-|-| +|cLV1 (1o)|-|-| +|cLV1 (2o)|[g3, g2, g1, B31, B21, B11, A32, A31, A22, A21, A12, A11, A23//A33, A13//A33, pi2(t), pi1(t), pi3(t)*A33]|[g3, g2, g1, B31, B21, B11, A32, A31, A22, A21, A12, A11, A23//A33, A13//A33, pi2(t), pi1(t), pi3(t)*A33]| +|cholera|[xi, mu, k, gam, bw, bi, a, R(t), W(t), I(t), S(t)]|[xi, mu, k, gam, bw, bi, a, W(t), R(t), S(t) + I(t) + R(t), I(t)]| +|generalizedLoktaVolterra (1o)|[r2, r1, beta21, beta11, beta12//beta22, x1(t), x2(t)*beta22]|[r2, r1, beta21, beta11, beta12//beta22, x1(t), x2(t)*beta22]| +|generalizedLoktaVolterra (2o)|[r2, r1, beta22, beta21, beta12, beta11, x2(t), x1(t)]|[r2, r1, beta22, beta21, beta12, beta11, x2(t), x1(t)]| +|p53|[p9, p8, p7, p6, p5, p4, p3, p25, p24, p23, p21, p20, p18, p17, p16, p15, p14, p13, p12, p11, p10, p1, p22^4, x4(t), x3(t), x2(t), x1(t)]|[p9, p8, p7, p6, p5, p4, p3, p25, p24, p23, p21, p20, p18, p17, p16, p15, p14, p13, p12, p11, p10, p1, p22^4, x3(t), x2(t), x1(t), x4(t)]| *Benchmarking environment:* - - Total RAM (GiB): 15 - - Processor: Intel Xeon Processor (Icelake) - - Julia version: 1.9.3 +* Total RAM (GiB): 251 +* Processor: INTEL(R) XEON(R) GOLD 6538Y+ +* Julia version: 1.12.6 Versions of the dependencies: - - Primes : 0.5.4 - - BenchmarkTools : 1.3.2 - - IterTools : 1.8.0 - - PrecompileTools : 1.2.0 - - Symbolics : 5.5.3 - - Combinatorics : 1.0.2 - - SymbolicUtils : 1.3.0 - - DataStructures : 0.18.15 - - Groebner : 0.4.3 - - ParamPunPam : 0.0.3 - - ModelingToolkit : 8.68.0 - - AbstractAlgebra : 0.31.1 - - MacroTools : 0.5.11 - - Nemo : 0.35.3 - - SpecialFunctions : 2.3.1 +* AbstractAlgebra : 0.50.1 +* Combinatorics : 1.1.0 +* DataStructures : 0.19.6 +* Dates : 1.11.0 +* Groebner : 0.10.4 +* IterTools : 1.10.0 +* LinearAlgebra : 1.12.0 +* Logging : 1.11.0 +* MacroTools : 0.5.16 +* Nemo : 0.56.1 +* ParamPunPam : 0.5.8 +* PrecompileTools : 1.3.4 +* Primes : 0.5.7 +* Random : 1.11.0 +* RationalFunctionFields : 0.3.2 +* TimerOutputs : 0.5.29 diff --git a/benchmarking/IdentifiableFunctions/run_benchmarks.jl b/benchmarking/IdentifiableFunctions/run_benchmarks.jl index b70d03718..5c840f88e 100644 --- a/benchmarking/IdentifiableFunctions/run_benchmarks.jl +++ b/benchmarking/IdentifiableFunctions/run_benchmarks.jl @@ -1,44 +1,23 @@ #= -This script runs benchmarks and collects useful statistics. +Run identifiable-function benchmarks and compare the returned generators. -To benchmark a particular function `StructuralIdentifiability.XYZ`, execute the -following command in your favorite terminal from this directory: +From this directory, for example: +```sh +julia run_benchmarks.jl find_identifiable_functions \ + --timeout=3600 \ + --memorylimit=8192 \ + --keywords="(with_states=true, cmp=:default); (with_states=true, cmp=:cmp_lie)" ``` -$ julia run_benchmarks.jl XYZ -``` - -`XYZ` can be one of the following: -- `find_identifiable_functions`, -- `reparametrize_global`. - -Going into details, this will: -1. Load all benchmark models from benchmarks.jl and create a directory - `benchmark_results`. -2. Run the function `StructuralIdentifiability.XYZ`, collect useful runtime - statistics, and write them to `benchmark_results`. -3. Read the statistics from `benchmark_results` and produce a table with - results. - -Below in the file there is a description of some potentially useful command line -arguments. -A more involved way of running benchmarks could be the following: -``` -julia run_benchmarks.jl XYZ --target="id_total, id_are_polynomial" --timeout=3600 --keywords="(strategy=(:hybrid,),with_states=true); (strategy=(:normalforms,3),with_states=true)" -``` +Generated model files and results are written below `results/`; subprocess logs +are written below `logs/`. Both directories are ignored by Git. The summarized +comparison is written to `benchmark_result.md`. =# using ArgParse -using CpuId, Logging, Pkg, Printf -using Base.Threads -using Distributed -using Dates -using ProgressMeter - +using CpuId, Dates, Logging, Pkg, ProgressMeter using StructuralIdentifiability -using StructuralIdentifiability: _runtime_logger, ODE -using StructuralIdentifiability.ParamPunPam global_logger(Logging.ConsoleLogger(stdout, Logging.Warn)) include("../benchmarks.jl") @@ -51,7 +30,9 @@ progressbar_enabled() = include("utils.jl") -const BENCHMARK_TABLE = "benchmark_result.md" +const RESULTS_DIR = joinpath(@__DIR__, BENCHMARK_RESULTS) +const LOGS_DIR = joinpath(@__DIR__, BENCHMARK_LOGS) +const BENCHMARK_TABLE = joinpath(@__DIR__, "benchmark_result.md") function parse_commandline() s = ArgParseSettings() @@ -67,19 +48,14 @@ function parse_commandline() Semicolon-separated list of named tuples.""" arg_type = String default = "" - # default = "(strategy=(:gb, ),); (strategy=(:gb, ),with_states=true); (strategy=(:normalforms, 2),); (strategy=(:normalforms, 2),with_states=true); (strategy=(:normalforms, 3),); (strategy=(:normalforms, 3),with_states=true); (strategy=(:hybrid, ),); (strategy=(:hybrid, ),with_states=true)" - "--target" - help = """ - Statistics to be displayed in the table. - Comma-separated list of entities. - NOTE: each of these must be present in the list of - `ALL_POSSIBLE_CATEGORIES` in `utils.jl`""" - arg_type = String - default = join(ALL_CATEGORIES, ",") "--timeout" help = "Timeout, s." arg_type = Int default = 3600 + "--memorylimit" + help = "Optional memory limit in MiB for each benchmark subprocess. Use 0 to disable." + arg_type = Int + default = 0 "--workers" help = "The number of available worker processes." arg_type = Int @@ -87,25 +63,23 @@ function parse_commandline() "--skip" help = "Skip specified benchmark models." arg_type = Vector{Symbol} - default = [:NFkB,:Covid2, :Lincomp1, :Lincomp2] - # "JAK-STAT 1", - # "LeukaemiaLeon2021", - # "MAPK model (5 outputs bis)", - # "QWWC", - # "TumorHu2019", - # "TumorPillis2007", - # "cLV1 (1o)"] + default = [ + :NFkB, + :Covid2, + :Lincomp1, + :Lincomp2, + Symbol("cLV1 (1o)"), + Symbol("cLV1 (2o)"), + :TumorHu2019, + :TumorPillis2007, + ] "--models" help = """ Run specified benchmark models. - A comma-separated list of names. + A comma-separated list of benchmark IDs. Leave empty for selecting all models.""" arg_type = String default = "" - "--augment" - help = "Augment the benchmark dataset with similar models." - arg_type = Bool - default = false "--regen" help = "Re-generate the folder with benchmarks from scratch." arg_type = Bool @@ -122,71 +96,91 @@ function parse_commandline() return parse_args(s) end -function populate_benchmarks(args, kwargs) - regen = args["regen"] - dir_present = isdir((@__DIR__) * "/$BENCHMARK_RESULTS/") - dir_present && !regen && return false - @debug "Re-generating the benchmarks folder" - try - if isdir((@__DIR__) * "/$BENCHMARK_RESULTS/") - rm((@__DIR__) * "/$BENCHMARK_RESULTS/", recursive = true, force = true) - end - catch err - @info "Something went wrong when deleting the benchmarks folder" - showerror(stdout, err) - println(stdout) +function selected_benchmark_names(args) + requested_ids = filter(!isempty, strip.(split(args["models"], ","))) + names = if isempty(requested_ids) + filter(name -> isdir(joinpath(RESULTS_DIR, name)), readdir(RESULTS_DIR)) + else + [benchmarks[Symbol(id)][:name] for id in requested_ids] end - prog = Progress( + skipped_names = [ + benchmarks[id][:name] for id in args["skip"] if haskey(benchmarks, id) + ] + return sort(setdiff(names, skipped_names)) +end + +function populate_benchmarks(args) + results_present = isdir(RESULTS_DIR) + results_present && !args["regen"] && return + results_present && rm(RESULTS_DIR; recursive = true, force = true) + mkpath(RESULTS_DIR) + + @debug "Re-generating the benchmarks folder" + progress = Progress( length(benchmarks), - "Generating benchmarks", - # spinner = true, + "Generating benchmarks"; dt = 0.1, enabled = progressbar_enabled(), color = _progressbar_color, ) - for (id, bmark) in benchmarks - next!(prog) # , spinner = "⌜⌝⌟⌞") - name = bmark[:name] - system = bmark[:ode] + for (_, benchmark) in benchmarks + next!(progress) + name = benchmark[:name] + system = benchmark[:ode] @debug "Generating $name" - mkpath((@__DIR__) * "/$BENCHMARK_RESULTS/$name/") - fd = open((@__DIR__) * "/$BENCHMARK_RESULTS/$name/$name.jl", "w") - println(fd, "# $name") - println(fd, "#! format: off") - println(fd, "using StructuralIdentifiability") - println(fd, "") - ode = join(map(s -> "\t" * s, split(repr(system), "\n")[1:(end - 1)]), ",\n") - println(fd, "system = @ODEmodel(\n$ode\n)") - close(fd) + model_dir = joinpath(RESULTS_DIR, name) + mkpath(model_dir) + open(joinpath(model_dir, "$name.jl"), "w") do io + println(io, "# $name") + println(io, "#! format: off") + println(io, "using StructuralIdentifiability\n") + equations = join( + map(line -> "\t" * line, split(repr(system), "\n")[1:(end - 1)]), + ",\n", + ) + println(io, "system = @ODEmodel(\n$equations\n)") + end end - finish!(prog) - return true + finish!(progress) + return +end + +function benchmark_command(function_name, problem_name, function_kwargs, memory_limit) + project_dir = dirname(Base.active_project()) + script = joinpath(@__DIR__, "run_single_benchmark.jl") + command = `$(Base.julia_cmd()) --project=$project_dir $script $function_name $problem_name $(string(function_kwargs))` + command = Cmd(command; ignorestatus = true, detach = false, env = copy(ENV)) + memory_limit == 0 && return command + + Sys.isunix() || throw(ArgumentError("--memorylimit is supported only on Unix")) + shell_command = string( + "ulimit -v ", + memory_limit * 1024, + "; exec ", + join(Base.shell_escape.(command.exec), " "), + ) + shell = Cmd(["bash", "-lc", shell_command]) + return Cmd(shell; ignorestatus = true, detach = false, env = copy(ENV)) end function run_benchmarks(args, kwargs) - to_skip = args["skip"] timeout = args["timeout"] @assert timeout > 0 + memorylimit = args["memorylimit"] + @assert memorylimit >= 0 function_name = args["function"] @assert function_name in ("find_identifiable_functions", "reparametrize_global") nworkers = args["workers"] @assert nworkers > 0 - dirnames = first(walkdir((@__DIR__) * "/$BENCHMARK_RESULTS/"))[2] - models_from_args = filter(s -> !isempty(s), map(strip, split(args["models"], ","))) - to_run_names = if isempty(models_from_args) - dirnames - else - [benchmarks[Symbol(id)][:name] for id in models_from_args] - end - to_run_names = setdiff(to_run_names, to_skip) - to_run_indices = collect(1:length(to_run_names)) + to_run_names = selected_benchmark_names(args) + to_run_indices = eachindex(to_run_names) @info """ Benchmarking `$function_name`.""" @info """ Passing these keyword arguments to `$function_name`: - \t$(join(map(string, kwargs), "\n\t"))""" + $(join(map(string, kwargs), "\n\t"))""" @info """ Number of benchmark systems: $(length(to_run_indices)) Workers: $(nworkers) @@ -223,7 +217,6 @@ function run_benchmarks(args, kwargs) prog = Progress( length(queue), "Running benchmarks", - # spinner = true, dt = 0.3, enabled = progressbar_enabled(), color = _progressbar_color, @@ -236,21 +229,17 @@ function run_benchmarks(args, kwargs) problem_name = to_run_names[problem_id] global_id = keywords_to_global_id(function_kwargs) @debug "Running $problem_name / $global_id. Kwargs:\n$function_kwargs" - logfn = generic_filename("logs", global_id) - logs = open((@__DIR__) * "/$BENCHMARK_RESULTS/$problem_name/$logfn", "w") - # escaping ':' for Windows - function_kwargs_esc = replace(string(function_kwargs), ":" => "\\:") - cmd = Cmd( - [ - "julia", - (@__DIR__) * "/run_single_benchmark.jl", - "$function_name", - "$problem_name", - "$function_kwargs_esc", - ] + log_dir = joinpath(LOGS_DIR, problem_name) + mkpath(log_dir) + log_path = joinpath(log_dir, generic_filename("run", global_id) * ".log") + logs = open(log_path, "w") + command = benchmark_command( + function_name, + problem_name, + function_kwargs, + memorylimit, ) - cmd = Cmd(cmd, ignorestatus = true, detach = false, env = copy(ENV)) - proc = run(pipeline(cmd, stdout = logs, stderr = logs), wait = false) + proc = run(pipeline(command; stdout = logs, stderr = logs); wait = false) push!( processes, ( @@ -262,7 +251,6 @@ function run_benchmarks(args, kwargs) start_time = time_ns(), global_run_id = global_id, logfile = logs, - # errfile = errs, ), ) push!(running, processes[end]) @@ -271,7 +259,6 @@ function run_benchmarks(args, kwargs) showvalues = generate_showvalues(processes), step = 0, valuecolor = _progressbar_value_color, - # spinner = "⌜⌝⌟⌞", ) end @@ -285,7 +272,6 @@ function run_benchmarks(args, kwargs) push!(errored, proc) end close(proc.logfile) - # close(proc.errfile) kw = proc.function_kwargs start_time = proc.start_time next!( @@ -300,8 +286,9 @@ function run_benchmarks(args, kwargs) if seconds_passed(start_time) > timeout push!(to_be_removed, i) kill(proc.julia_process) + wait(proc.julia_process) close(proc.logfile) - # close(proc.errfile) + push!(errored, proc) kw = proc.function_kwargs next!( prog, @@ -321,7 +308,7 @@ function run_benchmarks(args, kwargs) finish!(prog) if !isempty(errored) - printstyled("(!) Maybe errored:\n", color = :red) + printstyled("Failed or timed out:\n", color = :red) for proc in errored println("\t$(proc.problem_name) / $(proc.global_run_id)") end @@ -330,90 +317,70 @@ function run_benchmarks(args, kwargs) return to_run_names end -function collect_timings(args, kwargs, names; content = :compare) +function compact_result_string(content) + content = strip(content) + isempty(content) && return "-" + first_bracket = findfirst(==('['), content) + if first_bracket !== nothing + content = content[first_bracket:end] + end + content = replace(content, "\n" => " ") + content = replace(content, r"\s+" => " ") + content = replace(content, "|" => "\\|") + return strip(content) +end + +function display_keyword_label(keyword) + isempty(keyword) && return "default" + if haskey(keyword, :cmp) + cmp_label = string(keyword.cmp) + if haskey(keyword, :with_states) && keyword.with_states + return "funcs ($cmp_label)" + end + return string(cmp_label) + end + return string(keywords_to_global_id(keyword)) +end + +function collect_results(args, kwargs, names) function_name = args["function"] - targets = map(Symbol, map(strip, split(args["target"], ","))) - @assert all(target -> target in ALL_CATEGORIES, targets) - @assert length(targets) > 0 @info """ Collecting benchmark results for `$function_name`. Keyword arguments of interest: - \t$(join(map(string, kwargs), "\n\t")) - - Statistics of interest: - \t$(join(map(string, targets), "\n\t")) + $(join(map(string, kwargs), "\n\t")) """ cannot_collect = [] names = sort(names) kwids = map(keywords_to_global_id, kwargs) + labels = map(display_keyword_label, kwargs) + memory_limit = args["memorylimit"] + memory_description = memory_limit == 0 ? "disabled" : "$memory_limit MiB" - # Collect timings and data from directory BENCHMARK_RESULTS. - data = Dict() + data = Dict{String, Dict{Symbol, String}}() for name in names - @debug "==== Reading $name" - data[name] = Dict() - for kwid in kwids - timings_file = nothing - ##### - data[name][kwid] = Dict() - timingsfn = timings_filename(kwid) - try - @debug "==== Opening /$BENCHMARK_RESULTS/$name/$timingsfn" - timings_file = - open((@__DIR__) * "/$BENCHMARK_RESULTS/$name/$timingsfn", "r") - catch e - @debug "Cannot collect timings for $name / $kwid" - push!(cannot_collect, (name, kwid)) - continue + model_data = Dict{Symbol, String}() + data[name] = model_data + for id in kwids + result_path = joinpath(RESULTS_DIR, name, result_filename(id)) + if isfile(result_path) + model_data[id] = compact_result_string(read(result_path, String)) + else + @debug "Cannot collect result for $name / $id" + push!(cannot_collect, (name, id)) + model_data[id] = "-" end - lines = readlines(timings_file) - if isempty(lines) - @debug "Cannot collect timings for $name / $kwid" - push!(cannot_collect, (name, kwid)) - continue - end - @assert lines[1] == name - for line in lines[2:end] - k, v = split(line, ", ") - data[name][kwid][Symbol(k)] = parse(Float64, v) - end - close(timings_file) - ##### - datafn = data_filename(kwid) - data_file = nothing - try - @debug "==== Opening /$BENCHMARK_RESULTS/$name/$datafn" - data_file = open((@__DIR__) * "/$BENCHMARK_RESULTS/$name/$datafn", "r") - catch e - @debug "Cannot collect data for $name / $kwid" - push!(cannot_collect, (name, kwid)) - continue - end - lines = readlines(data_file) - if isempty(lines) - @debug "Cannot collect data for $name / $kwid" - push!(cannot_collect, (name, kwid)) - continue - end - @assert lines[1] == name - for line in lines[2:end] - k, v = map(strip, split(line, ",")) - data[name][kwid][Symbol(k)] = v - end - close(data_file) end end if !isempty(cannot_collect) - printstyled("(!) Cannot collect timings for:\n", color = :red) + printstyled("(!) Cannot collect results for:\n", color = :red) for (name, kwid) in cannot_collect println("\t$name / $kwid") end end - # Print the table to BENCHMARK_TABLE. resulting_md = "" resulting_md *= """ ## Benchmark results @@ -423,33 +390,17 @@ function collect_timings(args, kwargs, names; content = :compare) - Benchmarked function: `$(args["function"])` - Workers: $(args["workers"]) - Timeout: $(args["timeout"]) s - - **All timings in seconds.** + - Memory limit: $memory_description """ - makecolname(kw, target) = - kw === Symbol("") ? HUMAN_READABLE_CATEGORIES[target] : - Symbol(Symbol(kw), Symbol(" / "), HUMAN_READABLE_CATEGORIES[target]) - columns = [makecolname(kwid, target) for kwid in kwids for target in targets] - resulting_md *= "|Model|" * join(map(string, columns), "|") * "|\n" - resulting_md *= "|-----|" * join(["---" for _ in columns], "|") * "|\n" + columns = ["Model"; labels] + resulting_md *= "|" * join(columns, "|") * "|\n" + resulting_md *= "|" * join(["---" for _ in columns], "|") * "|\n" for name in names - model_data = data[name] resulting_md *= "|$name|" for kwid in kwids - if !haskey(model_data, kwid) - resulting_md *= (" - " * "|")^length(columns) - continue - end - for target in targets - if !haskey(model_data[kwid], target) - resulting_md *= " - " * "|" - else - formatting_style = CATEGORY_FORMAT[target] - resulting_md *= formatting_style(model_data[kwid][target]) * "|" - end - end + resulting_md *= data[name][kwid] * "|" end resulting_md *= "\n" end @@ -461,14 +412,13 @@ function collect_timings(args, kwargs, names; content = :compare) resulting_md *= "Versions of the dependencies:\n\n" deps = Pkg.dependencies() - stid_info = deps[findfirst(x -> x.name == "StructuralIdentifiability", deps)] - for (s, uid) in stid_info.dependencies - if deps[uid].version !== nothing - resulting_md *= "* $s : $(deps[uid].version)\n" - end + dependencies = sort(collect(Pkg.project().dependencies); by = first) + for (name, uuid) in dependencies + version = deps[uuid].version + version === nothing || (resulting_md *= "* $name : $version\n") end - return open((@__DIR__) * "/$BENCHMARK_TABLE", "w") do io + return open(BENCHMARK_TABLE, "w") do io write(io, resulting_md) end end @@ -484,32 +434,27 @@ function main() @debug "Keywords for `$(args["function"])`" @debug kwargs if args["tableonly"] - dirnames = first(walkdir((@__DIR__) * "/$BENCHMARK_RESULTS/"))[2] - models_from_args = filter(s -> !isempty(s), map(strip, split(args["models"], ","))) - problems = if isempty(models_from_args) - dirnames - else - [benchmarks[Symbol(id)][:name] for id in models_from_args] - end - problems = setdiff(problems, args["skip"]) + problems = selected_benchmark_names(args) else - flag = populate_benchmarks(args, kwargs) + populate_benchmarks(args) problems = run_benchmarks(args, kwargs) printstyled( """ Benchmarking had finished in $(round((time_ns() - timestamp) / 1.0e9, digits = 2)) seconds. - Results are written to /$BENCHMARK_RESULTS + Results are written to $RESULTS_DIR """, color = :light_green, ) end - collect_timings(args, kwargs, problems) + collect_results(args, kwargs, problems) return printstyled( """ - Table with results is written to /$BENCHMARK_TABLE + Table with results is written to $BENCHMARK_TABLE """, color = :light_green, ) end -main() +if abspath(PROGRAM_FILE) == @__FILE__ + main() +end diff --git a/benchmarking/IdentifiableFunctions/run_single_benchmark.jl b/benchmarking/IdentifiableFunctions/run_single_benchmark.jl index 42dd28cf2..15cd2aa03 100644 --- a/benchmarking/IdentifiableFunctions/run_single_benchmark.jl +++ b/benchmarking/IdentifiableFunctions/run_single_benchmark.jl @@ -1,22 +1,16 @@ -using CpuId, Logging, Pkg, Printf -using Statistics - +using Logging using StructuralIdentifiability -using StructuralIdentifiability: _runtime_logger, ODE -using StructuralIdentifiability.ParamPunPam -logger = Logging.ConsoleLogger(stdout, Logging.Info) -global_logger(logger) +global_logger(Logging.ConsoleLogger(stdout, Logging.Info)) include("utils.jl") -const data = Dict() - -const NUM_RUNS = 5 +const NUM_RUNS = 1 const FUNCTION_NAME = ARGS[1] const PROBLEM_NAME = ARGS[2] -const KWARGS = parse_keywords(ARGS[3])[1] +const KWARGS = only(parse_keywords(ARGS[3])) const GLOBAL_ID = keywords_to_global_id(KWARGS) +const RESULT_DIR = joinpath(@__DIR__, BENCHMARK_RESULTS, PROBLEM_NAME) @info "" FUNCTION_NAME @info "" PROBLEM_NAME @@ -25,98 +19,79 @@ const GLOBAL_ID = keywords_to_global_id(KWARGS) flush(stdout) flush(stderr) -# Load the system -path = (@__DIR__) * "/$BENCHMARK_RESULTS/$PROBLEM_NAME/$PROBLEM_NAME.jl" -include(path) +include(joinpath(RESULT_DIR, "$PROBLEM_NAME.jl")) + +const RESOLVED_KWARGS = if haskey(KWARGS, :cmp) && KWARGS.cmp isa Symbol + merge(KWARGS, (cmp = benchmark_comparator(KWARGS.cmp, system),)) +else + KWARGS +end +const BENCHMARK_FUNCTION = getproperty(StructuralIdentifiability, Symbol(FUNCTION_NAME)) -# Compile -const FUNC_SYM = Symbol(FUNCTION_NAME) -ode = @ODEmodel(x1'(t) = a * x1(t) + x2, x2'(t) = b * x2(t) + c * d, y(t) = x1(t)) -eval(:($FUNC_SYM($system; ($KWARGS...)))) +# Compile before collecting timings. +BENCHMARK_FUNCTION(system; RESOLVED_KWARGS...) function process_system() @info "Processing $PROBLEM_NAME" @info """ Averaging over $NUM_RUNS runs. Using keyword arguments: - $(typeof(KWARGS)) - $KWARGS + $RESOLVED_KWARGS ID: $GLOBAL_ID""" - data[PROBLEM_NAME] = Dict{Any, Any}(c => 0.0 for c in ALL_CATEGORIES) + model_data = Dict{Symbol, Any}(category => 0.0 for category in ID_TIME_CATEGORIES) for _ in 1:NUM_RUNS - timing = @timed result = eval(:($FUNC_SYM($system; ($KWARGS...)))) - data[PROBLEM_NAME][:return_value] = result - @info "Result is" result - for cat in ID_TIME_CATEGORIES - if haskey(StructuralIdentifiability._runtime_logger, cat) - data[PROBLEM_NAME][cat] = StructuralIdentifiability._runtime_logger[cat] + timing = @timed BENCHMARK_FUNCTION(system; RESOLVED_KWARGS...) + model_data[:return_value] = timing.value + @info "Result is" timing.value + + for category in ID_TIME_CATEGORIES + category === :id_total && continue + if haskey(StructuralIdentifiability._runtime_logger, category) + model_data[category] += StructuralIdentifiability._runtime_logger[category] end end - for cat in ID_DATA_CATEGORIES - if haskey(StructuralIdentifiability._runtime_logger, cat) - data[PROBLEM_NAME][cat] = - deepcopy(StructuralIdentifiability._runtime_logger[cat]) + for category in ID_DATA_CATEGORIES + if haskey(StructuralIdentifiability._runtime_logger, category) + model_data[category] = + deepcopy(StructuralIdentifiability._runtime_logger[category]) end end - data[PROBLEM_NAME][:id_total] = timing.time + model_data[:id_total] += timing.time end - for cat in ID_TIME_CATEGORIES - if haskey(data[PROBLEM_NAME], cat) - data[PROBLEM_NAME][cat] = data[PROBLEM_NAME][cat] / NUM_RUNS - end + + for category in ID_TIME_CATEGORIES + model_data[category] /= NUM_RUNS end - return + return model_data end -function dump_timings() - timings = "" - timings *= "$PROBLEM_NAME\n" - for (key, model_data) in data - for c in ID_TIME_CATEGORIES - timings *= "$c, " - timings *= string(model_data[c]) * "\n" +function dump_timings(model_data) + filename = joinpath(RESULT_DIR, timings_filename(GLOBAL_ID)) + return open(filename, "w") do io + println(io, PROBLEM_NAME) + for category in ID_TIME_CATEGORIES + println(io, "$category, $(model_data[category])") end end - filename = timings_filename(GLOBAL_ID) - return open((@__DIR__) * "/$BENCHMARK_RESULTS/$PROBLEM_NAME/$filename", "w") do io - write(io, timings) - end end -function dump_results() - filename = result_filename(GLOBAL_ID) - open((@__DIR__) * "/$BENCHMARK_RESULTS/$PROBLEM_NAME/$filename", "w") do io - if haskey(data, PROBLEM_NAME) - println(io, data[PROBLEM_NAME][:return_value]) - end +function dump_results(model_data) + result_path = joinpath(RESULT_DIR, result_filename(GLOBAL_ID)) + open(result_path, "w") do io + println(io, model_data[:return_value]) end - filename = data_filename(GLOBAL_ID) - open((@__DIR__) * "/$BENCHMARK_RESULTS/$PROBLEM_NAME/$filename", "w") do io - write(io, "$PROBLEM_NAME\n") - end - for cat in ID_DATA_CATEGORIES - if !haskey(data[PROBLEM_NAME], cat) - continue - end - if cat === :something_important - # make a separate file for it - filename_cat = generic_filename(cat, GLOBAL_ID) - open((@__DIR__) * "/$BENCHMARK_RESULTS/$PROBLEM_NAME/$filename_cat", "w") do io - # print something - end - continue - end - # otherwise, print in the data file - open((@__DIR__) * "/$BENCHMARK_RESULTS/$PROBLEM_NAME/$filename", "a+") do io - write(io, "$cat, ") - write(io, string(data[PROBLEM_NAME][cat])) - write(io, "\n") + + data_path = joinpath(RESULT_DIR, data_filename(GLOBAL_ID)) + return open(data_path, "w") do io + println(io, PROBLEM_NAME) + for category in ID_DATA_CATEGORIES + haskey(model_data, category) || continue + println(io, "$category, $(model_data[category])") end end - return end -process_system() -dump_timings() -dump_results() +const MODEL_DATA = process_system() +dump_timings(MODEL_DATA) +dump_results(MODEL_DATA) diff --git a/benchmarking/IdentifiableFunctions/utils.jl b/benchmarking/IdentifiableFunctions/utils.jl index b1906b99b..a394b0b84 100644 --- a/benchmarking/IdentifiableFunctions/utils.jl +++ b/benchmarking/IdentifiableFunctions/utils.jl @@ -1,7 +1,7 @@ # shared utilities for benchmarking -using Printf const BENCHMARK_RESULTS = "results" +const BENCHMARK_LOGS = "logs" const ID_TIME_CATEGORIES = [ :id_total, @@ -18,108 +18,63 @@ const ID_TIME_CATEGORIES = [ const ID_DATA_CATEGORIES = [:id_npoints_degree, :id_npoints_interpolation, :are_id_funcs_polynomial] const ALL_CATEGORIES = union(ID_TIME_CATEGORIES, ID_DATA_CATEGORIES) -const ALL_POSSIBLE_CATEGORIES = union(ALL_CATEGORIES, Symbol[]) -const HUMAN_READABLE_CATEGORIES = Dict( - :id_io_time => "io", - :id_primality_evaluate => "io/primality-evaluate", - :id_uncertain_factorization => "io/uncertain-factor", - :id_global_time => "global id.", - :id_ideal_time => "gen. ideal", - :id_inclusion_check => "inclusion", - :id_inclusion_check_mod_p => "inclusion Zp", - :id_groebner_time => "ParamPunPam.jl", - :id_total => "Runtime", - :id_beautifulization => "beautifulization", - :id_normalforms_time => "normal forms", - :id_gbfan_time => "GB fan", - :id_ranking => "Score", - :implicit_relations => "Algebraic relations", - :dim_before => "Dim. before", - :dim_after => "Dim. after", - :are_id_funcs_polynomial => "Polynomial?", - :id_npoints_degree => "# Points, degree", - :id_npoints_interpolation => "# Points, interpolation", -) +function parse_keywords(input) + isempty(strip(input)) && return [NamedTuple()] -const CATEGORY_FORMAT = Dict() -for cat in ALL_POSSIBLE_CATEGORIES - CATEGORY_FORMAT[cat] = (val) -> if val isa Real - @sprintf("%.2f", val) - else - string(val) + return map(split(input, ";")) do keyword_set + expression = Meta.parse(strip(keyword_set)) + keywords = Core.eval(@__MODULE__, expression) + keywords == () && return NamedTuple() + keywords isa NamedTuple || + throw(ArgumentError("Expected a named tuple of keyword arguments, got $keywords")) + return keywords end end -CATEGORY_FORMAT[:are_id_funcs_polynomial] = (val) -> string(val) == "true" ? "yes" : "no" -function parse_keywords(keywords) - keywords = replace(keywords, "\\:" => ":") - keywords = replace(keywords, "\"" => "") - if isempty(keywords) - return [[]] +function benchmark_comparator(label, ode) + if label in (:default, :cmp_default) + return StructuralIdentifiability.default_cmp(ode) + elseif label === :cmp_lie + return StructuralIdentifiability.cmp_lie( + ode, + StructuralIdentifiability.RationalFunctionFields.rational_function_cmp, + ) end - if isempty(eval(Meta.parse(keywords))) - return [[]] - end - sets_of_keywords = map(strip, split(keywords, ";")) - @assert !isempty(sets_of_keywords) - kws_named_tuples = [] - for kwset in sets_of_keywords - @debug "" kwset - nt = eval(Meta.parse(kwset)) - push!(kws_named_tuples, nt) - end - return kws_named_tuples + throw(ArgumentError("Unknown comparator label: $label")) +end + +function benchmark_keyword_label(value) + value isa Symbol && return value + label = value isa Function ? typeof(value) : value + return Symbol(replace(string(label), r"[^A-Za-z0-9_]+" => "_")) end function keywords_to_global_id(keywords) - if isempty(keywords) - return Symbol() - end - id = get(keywords, :strategy, Symbol("")) - if haskey(keywords, :with_states) - if keywords.with_states - id = if id !== Symbol("") - Symbol(id, :_with_states) - else - Symbol(id, :with_states) - end - end + isempty(keywords) && return Symbol() + + id = get(keywords, :strategy, Symbol()) + if get(keywords, :with_states, false) + id = isempty(string(id)) ? :with_states : Symbol(id, :_with_states) end if haskey(keywords, :rational_interpolator) - interpolator = keywords.rational_interpolator - id = Symbol(id, :_, interpolator) + id = Symbol(id, :_, keywords.rational_interpolator) + end + if haskey(keywords, :cmp) + comparator = benchmark_keyword_label(keywords.cmp) + id = isempty(string(id)) ? comparator : Symbol(id, :_, comparator) end - if haskey(keywords, :adjoin_identifiable) - if keywords.adjoin_identifiable - id = if id !== Symbol("") - Symbol(id, :_adjoin_identifiable) - else - Symbol(id, :adjoin_identifiable) - end - end + if get(keywords, :adjoin_identifiable, false) + id = isempty(string(id)) ? :adjoin_identifiable : Symbol(id, :_adjoin_identifiable) end return id end -function timings_filename(kwid) - return generic_filename("timings", kwid) -end - -function result_filename(kwid) - return generic_filename("result", kwid) -end +timings_filename(id) = generic_filename("timings", id) +result_filename(id) = generic_filename("result", id) +data_filename(id) = generic_filename("data", id) -function data_filename(kwid) - return generic_filename("data", kwid) -end - -function generic_filename(name, kwid) - str = if kwid === Symbol("") - "$name" - else - "$(name)_$kwid" - end - str = replace(str, ":" => "") - return str +function generic_filename(name, id) + filename = isempty(string(id)) ? name : "$(name)_$id" + return replace(filename, ":" => "") end diff --git a/benchmarking/benchmark_result.md b/benchmarking/benchmark_result.md index 7eeed7d4a..86518c171 100644 --- a/benchmarking/benchmark_result.md +++ b/benchmarking/benchmark_result.md @@ -1,38 +1,28 @@ -| Model | loc_time | glob_time | ioeq_time | wrnsk_time | rank_time | check_time | total | -|:-------------------------- |:-------- |:--------- |:--------- |:---------- |:--------- |:---------- |:------- | -| Goodwin oscillator | 0.03 | 0.09 | 0.02 | 0.04 | 0.00 | 0.02 | 0.12 | -| Akt pathway | 0.64 | 4.37 | 0.18 | 0.95 | 0.00 | 3.23 | 5.01 | -| SIWR original | 0.08 | 13.08 | 2.97 | 2.86 | 0.14 | 7.10 | 13.16 | -| SEAIJRC Covid model | 0.07 | 76.34 | 27.58 | 16.23 | 0.63 | 31.90 | 76.41 | -| MAPK model (6 outputs) | 4.96 | 31.59 | 4.24 | 10.70 | 0.00 | 16.65 | 36.55 | -| Pharm | 0.05 | 100.86 | 19.30 | 39.82 | 5.64 | 36.10 | 100.90 | -| SIRS forced | 0.04 | 8.97 | 0.89 | 5.71 | 0.22 | 2.14 | 9.01 | -| SIWR with extra output | 0.09 | 0.67 | 0.44 | 0.09 | 0.00 | 0.13 | 0.76 | -| CD8 T cell differentiation | 0.48 | 0.20 | 0.02 | 0.09 | 0.00 | 0.09 | 0.69 | -| HIV | 0.11 | 0.41 | 0.04 | 0.06 | 0.00 | 0.31 | 0.52 | -| MAPK model (5 outputs bis) | 10.33 | 1665.97 | 461.87 | 214.31 | 0.05 | 989.75 | 1676.31 | -| Chemical reaction network | 0.04 | 0.14 | 0.05 | 0.05 | 0.00 | 0.04 | 0.18 | -| Modified LV for testing | 1.17 | 2.46 | 1.50 | 0.63 | 0.01 | 0.24 | 11.59 | -| MAPK model (5 outputs) | 3.91 | 98.06 | 59.62 | 9.97 | 0.00 | 28.46 | 101.97 | +|Model|loc_time|glob_time|wrnsk_time|rank_time|check_time|total| +|-----|---|---|---|---|---|---| +|SEUIR|0.61|5.91|0.17|0.00|4.78|11.21| *Benchmarking environment:* - - Total RAM (Mb): 16384.0 - - Processor: Intel(R) Core(TM) i5-8210Y CPU @ 1.60GHz - - Julia version: 1.7.2 +* Total RAM (Mb): 257647.125 +* Processor: INTEL(R) XEON(R) GOLD 6538Y+ +* Julia version: 1.12.6 Versions of the dependencies: - - Primes : 0.5.2 - - BenchmarkTools : 1.3.1 - - Singular : 0.10.1 - - IterTools : 1.4.0 - - DataStructures : 0.18.11 - - Groebner : 0.2.5 - - ModelingToolkit : 8.9.0 - - AbstractAlgebra : 0.25.3 - - MacroTools : 0.5.9 - - Nemo : 0.30.0 - - Hecke : 0.13.4 - - SpecialFunctions : 2.1.4 - - TestSetExtensions : 2.0.0 +* Random : 1.11.0 +* Primes : 0.5.7 +* TimerOutputs : 0.5.29 +* IterTools : 1.10.0 +* PrecompileTools : 1.3.4 +* LinearAlgebra : 1.12.0 +* Combinatorics : 1.1.0 +* DataStructures : 0.18.22 +* Groebner : 0.10.6 +* ParamPunPam : 0.5.8 +* RationalFunctionFields : 0.3.2 +* AbstractAlgebra : 0.46.5 +* MacroTools : 0.5.16 +* Dates : 1.11.0 +* Nemo : 0.51.1 +* Logging : 1.11.0 From fb10a485379b8008b0571a79648f332b71f40a6c Mon Sep 17 00:00:00 2001 From: Alexander Demin <60229118+sumiya11@users.noreply.github.com> Date: Mon, 3 Aug 2026 03:40:10 +0200 Subject: [PATCH 11/17] Add default_cmp to exported functions --- src/StructuralIdentifiability.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index e40797159..3d1e351d1 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -44,7 +44,7 @@ export assess_local_identifiability, assess_identifiability export set_parameter_values, rename_variables # extra functionality -export find_ioequations, find_identifiable_functions +export find_ioequations, find_identifiable_functions, default_cmp # exporting to other formats export print_for_maple, print_for_DAISY, print_for_COMBOS, print_for_GenSSI From bb73df485983f2a5178431079acb886acaa8710d Mon Sep 17 00:00:00 2001 From: Sasha Date: Wed, 19 Aug 2026 17:24:24 +0200 Subject: [PATCH 12/17] add docstring to default_cmp --- src/util.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/util.jl b/src/util.jl index 1dbd80a59..3eee2d55e 100644 --- a/src/util.jl +++ b/src/util.jl @@ -303,6 +303,31 @@ end # ----------------------------------------------------------------------------- +""" + default_cmp(ode) + +Returns a comparator for rational functions. +Used as a default measure of simplicity in StructuralIdentifiability.jl. + +## Example + +```jldoctest +using StructuralIdentifiability + +ode = @ODEmodel( + x0'(t) = -(a01 + a21) * x0(t) + a12 * x1(t), + x1'(t) = a21 * x0(t) - a12 * x1(t), + y(t) = x0(t) +) + +f1, f2 = (a01 + a21) // a12, x0 // a12 + +cmp = default_cmp(ode) + +cmp(f1, f2) # returns true, so f1 is considered simpler than f2 +``` + +""" default_cmp(ode) = cmp_prefer_params(ode, rational_function_cmp) """ From 19a0a8e1190626f8f79f1a9c39167080734fc3f8 Mon Sep 17 00:00:00 2001 From: Sasha Date: Wed, 19 Aug 2026 18:36:16 +0200 Subject: [PATCH 13/17] fix jldoctests --- src/identifiable_functions.jl | 3 ++- src/util.jl | 4 ++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index ed21486a9..bbbff8af4 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -41,7 +41,8 @@ ode = @ODEmodel( find_identifiable_functions(ode) -# prints +# Output + 2-element Vector{AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}: a12 + a01 + a21 a12*a01 diff --git a/src/util.jl b/src/util.jl index 3eee2d55e..f55eeef93 100644 --- a/src/util.jl +++ b/src/util.jl @@ -325,6 +325,10 @@ f1, f2 = (a01 + a21) // a12, x0 // a12 cmp = default_cmp(ode) cmp(f1, f2) # returns true, so f1 is considered simpler than f2 + +# Output + +true ``` """ From 8d8a9c0fa1efdec7ad22a8af71d47d55b3693a40 Mon Sep 17 00:00:00 2001 From: Sasha Date: Wed, 19 Aug 2026 18:39:43 +0200 Subject: [PATCH 14/17] do not import generators --- src/StructuralIdentifiability.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index b69530b56..6aa9fedbf 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -67,7 +67,6 @@ using RationalFunctionFields: RationalFunctionField, check_algebraicity, field_contains, - generators, simplified_generating_set using RationalFunctionFields: eval_at_dict, From 0b323a40f7605846e1e1eb8f50924db2a7f6ee75 Mon Sep 17 00:00:00 2001 From: Alexander Demin <60229118+sumiya11@users.noreply.github.com> Date: Wed, 19 Aug 2026 18:56:16 +0200 Subject: [PATCH 15/17] Correct comment formatting in util.jl Fix comment formatting for clarity. --- src/util.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/util.jl b/src/util.jl index f55eeef93..578c5a437 100644 --- a/src/util.jl +++ b/src/util.jl @@ -326,7 +326,7 @@ cmp = default_cmp(ode) cmp(f1, f2) # returns true, so f1 is considered simpler than f2 -# Output +# output true ``` From 4f8297e3270d6b9327b554e9153dbccf26f7e0fc Mon Sep 17 00:00:00 2001 From: Alexander Demin <60229118+sumiya11@users.noreply.github.com> Date: Wed, 19 Aug 2026 18:56:38 +0200 Subject: [PATCH 16/17] Fix output comment in identifiable_functions.jl Updated comment to correct output formatting. --- src/identifiable_functions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index 081b7002f..81eb9d6db 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -41,7 +41,7 @@ ode = @ODEmodel( find_identifiable_functions(ode) -# Output +# output 2-element Vector{AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}: a12 + a01 + a21 From cca1416021d2b5042648188368de6acb82cd6633 Mon Sep 17 00:00:00 2001 From: Sasha Date: Wed, 19 Aug 2026 19:54:18 +0200 Subject: [PATCH 17/17] fix doctest of default_cmp --- src/util.jl | 23 ++++++++++------------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/src/util.jl b/src/util.jl index 578c5a437..46cdf095d 100644 --- a/src/util.jl +++ b/src/util.jl @@ -311,23 +311,20 @@ Used as a default measure of simplicity in StructuralIdentifiability.jl. ## Example -```jldoctest -using StructuralIdentifiability +```jldoctest; setup = :(using Logging; Logging.disable_logging(Logging.Info);) +julia> using StructuralIdentifiability -ode = @ODEmodel( - x0'(t) = -(a01 + a21) * x0(t) + a12 * x1(t), - x1'(t) = a21 * x0(t) - a12 * x1(t), - y(t) = x0(t) -) +julia> ode = @ODEmodel( + x0'(t) = -(a01 + a21) * x0(t) + a12 * x1(t), + x1'(t) = a21 * x0(t) - a12 * x1(t), + y(t) = x0(t) + ); -f1, f2 = (a01 + a21) // a12, x0 // a12 +julia> f1, f2 = (a01 + a21) // a12, x0 // a12; -cmp = default_cmp(ode) - -cmp(f1, f2) # returns true, so f1 is considered simpler than f2 - -# output +julia> cmp = default_cmp(ode); +julia> cmp(f1, f2) # returns true, so f1 is considered simpler than f2 true ```