From d6f7927cc5154e22eec72a2d016a7e86c5268ff6 Mon Sep 17 00:00:00 2001 From: gpogudin Date: Tue, 21 Nov 2023 15:02:44 +0100 Subject: [PATCH 01/12] first version of find_identifiable_functions_kic --- src/StructuralIdentifiability.jl | 1 + src/known_ic.jl | 85 ++++++++++++++++++++++++++++++++ test/known_ic.jl | 38 ++++++++++++++ test/runtests.jl | 3 +- 4 files changed, 126 insertions(+), 1 deletion(-) create mode 100644 src/known_ic.jl create mode 100644 test/known_ic.jl diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index 165570f3a..853845aaf 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -72,6 +72,7 @@ include("lincomp.jl") include("pb_representation.jl") include("submodels.jl") include("discrete.jl") +include("known_ic.jl") function __init__() _si_logger[] = @static if VERSION >= v"1.7.0" diff --git a/src/known_ic.jl b/src/known_ic.jl new file mode 100644 index 000000000..6ba2d3572 --- /dev/null +++ b/src/known_ic.jl @@ -0,0 +1,85 @@ +""" + find_identifiable_functions_kic(ode::ODE, known_ic; options...) + +Finds all functions of parameters/initial conditions that are identifiable in the given ODE +system under assumptions that the initial conditions for functions in the list +`known_ic` are known and generic. + +## Options + +This functions takes the following optional arguments: +- `simplify`: The extent to which the output functions are simplified. Stronger + simplification may require more time. Possible options are: + - `:standard`: Default simplification. + - `:weak`: Weak simplification. This option is the fastest, but the output + functions can be quite complex. + - `:strong`: Strong simplification. This option is the slowest, but the output + functions are nice and simple. + - `:absent`: No simplification. +- `p`: A float in the range from 0 to 1, the probability of correctness. Default + is `0.99`. +- `seed`: The rng seed. Default value is `42`. +- `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) + +**This is experimental functionality** + +``` + +""" +function find_identifiable_functions_kic( + ode::ODE{T}, + known_ic::Vector{<: Union{T, Generic.Frac{T}}}; + p::Float64 = 0.99, + seed = 42, + simplify = :standard, + rational_interpolator = :VanDerHoevenLecerf, + loglevel = Logging.Info, +) where {T <: MPolyElem{fmpq}} + restart_logging(loglevel = loglevel) + reset_timings() + with_logger(_si_logger[]) do + return _find_identifiable_functions_kic( + ode, + known_ic, + p = p, + seed = seed, + simplify = simplify, + rational_interpolator = rational_interpolator, + ) + end +end + +function _find_identifiable_functions_kic( + ode::ODE{T}, + known_ic::Vector{<: Union{T, Generic.Frac{T}}}; + p::Float64 = 0.99, + seed = 42, + simplify = :standard, + rational_interpolator = :VanDerHoevenLecerf, +) where {T <: MPolyElem{fmpq}} + Random.seed!(seed) + @assert simplify in (:standard, :weak, :strong, :absent) + half_p = 0.5 + p / 2 + runtime_start = time_ns() + id_funcs_general = find_identifiable_functions( + ode, + p = half_p, + with_states = true, + simplify = simplify, + rational_interpolator = rational_interpolator, + seed = seed, + ) + @info "The search for identifiable functions with known initial conditions concluded in $((time_ns() - runtime_start) / 1e9) seconds" + + id_funcs = simplified_generating_set( + RationalFunctionField(vcat(id_funcs_general, [f // one(parent(ode)) for f in known_ic])), + p = half_p, + seed = seed, + simplify = simplify, + rational_interpolator = rational_interpolator, + ) + + return id_funcs +end + + diff --git a/test/known_ic.jl b/test/known_ic.jl new file mode 100644 index 000000000..2d6db1f58 --- /dev/null +++ b/test/known_ic.jl @@ -0,0 +1,38 @@ +cases = [] + +ode = @ODEmodel( + x1'(t) = a * x1(t) - b * x1(t) * x2(t), + x2'(t) = -c * x2(t) + d * x1(t) * x2(t), + y(t) = x1(t) +) + +push!(cases, Dict( + :ode => ode, + :known => [x2], + :correct => [a, b, c, d, x1, x2], +)) + +ode = @ODEmodel( + x1'(t) = a + x2(t) + x3(t), + x2'(t) = b, + x3'(t) = c, + y(t) = x1(t) +) + +push!(cases, Dict( + :ode => ode, + :known => [x2, x3], + :correct => [a, b + c, x1, x2, x3], +)) + + +@testset "Identifiable functions with known generic initial conditions" begin + p = 0.99 + for case in test_cases + ode = case[:ode] + known = case[:known] + result = find_identifiable_functions_kic(ode, known) + correct = [f // one(parent(ode)) for f in case[:correct]] + @test Set(result) == Set(correct) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 82e6b0dd7..9661d5d5a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -58,7 +58,8 @@ using StructuralIdentifiability: extract_coefficients_ratfunc, lie_derivative, states_generators, - RationalFunctionField + RationalFunctionField, + find_identifiable_functions_kic, function random_ps(ps_ring, range = 1000) result = zero(ps_ring) From 32022c2461e10276fff2e572c2a78f6592ed2ce5 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Tue, 21 Nov 2023 21:43:25 +0100 Subject: [PATCH 02/12] tests fixed --- test/known_ic.jl | 3 +-- test/runtests.jl | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/test/known_ic.jl b/test/known_ic.jl index 2d6db1f58..c64bf5dda 100644 --- a/test/known_ic.jl +++ b/test/known_ic.jl @@ -27,8 +27,7 @@ push!(cases, Dict( @testset "Identifiable functions with known generic initial conditions" begin - p = 0.99 - for case in test_cases + for case in cases ode = case[:ode] known = case[:known] result = find_identifiable_functions_kic(ode, known) diff --git a/test/runtests.jl b/test/runtests.jl index 9661d5d5a..c98edd821 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -59,7 +59,7 @@ using StructuralIdentifiability: lie_derivative, states_generators, RationalFunctionField, - find_identifiable_functions_kic, + find_identifiable_functions_kic function random_ps(ps_ring, range = 1000) result = zero(ps_ring) From 5572d001666177cc688cff22163547a8b1948d17 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Wed, 22 Nov 2023 00:23:46 +0100 Subject: [PATCH 03/12] Algebraicity check implemented --- .../RationalFunctionField.jl | 49 +++++++++++++++++++ src/local_identifiability.jl | 1 + test/check_algebraicity.jl | 20 ++++++++ test/runtests.jl | 1 + 4 files changed, 71 insertions(+) create mode 100644 test/check_algebraicity.jl diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index 9c8e68fbf..f34729e2c 100644 --- a/src/RationalFunctionFields/RationalFunctionField.jl +++ b/src/RationalFunctionFields/RationalFunctionField.jl @@ -41,6 +41,10 @@ function poly_ring(F::RationalFunctionField) return parent(first(first(F.dennums))) end +function generators(F::RationalFunctionField) + return vcat([[f[i] // f[1] for i in 2:length(f)] for f in F.dennums]...) +end + # ------------------------------------------------------------------------------ """ @@ -165,6 +169,51 @@ end # ------------------------------------------------------------------------------ +""" + check_algebraicity(field, ratfuncs, p) + +Checks whether given rational function `ratfuncs` are algebraic over the field `field` +The result is correct with probability at least `p` + +Inputs: +- `field` - a rational function field +- `ratfuncs` - a list of lists of rational functions. +- `p` real number from (0, 1) + +Output: +- a list `L[i]` of bools of length `length(rat_funcs)` such that `L[i]` is true iff + the i-th function is algebraic over the `field` +""" +function check_algebraicity(field, ratfuncs, p) + fgens = generators(field) + base_vars = gens(poly_ring(field)) + maxdeg = maximum([max(total_degree(numerator(f)), total_degree(denominator(f))) for f in vcat(ratfuncs, fgens)]) + # degree of the polynomial whose nonvanishing will be needed for correct result + D = Int(ceil(2 * maxdeg * (length(fgens) + 1)^3 * length(ratfuncs) / (1 - p))) + eval_point = [Nemo.QQ(rand(1:D)) for x in base_vars] + + # Filling the jacobain for generators + S = MatrixSpace(Nemo.QQ, length(base_vars), length(fgens) + 1) + J = zero(S) + for (i, f) in enumerate(fgens) + for (j, x) in enumerate(base_vars) + J[j, i] = evaluate(derivative(f, x), eval_point) + end + end + rank = LinearAlgebra.rank(J) + + result = Bool[] + for f in ratfuncs + for (j, x) in enumerate(base_vars) + J[j, end] = evaluate(derivative(f, x), eval_point) + end + push!(result, LinearAlgebra.rank(J) == rank) + end + return result +end + +# ------------------------------------------------------------------------------ + function issubfield(F::RationalFunctionField{T}, E::RationalFunctionField{T}, p) where {T} return all(field_contains(E, F.dennums, p)) end diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index 3d680c336..e9bc99f8f 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -278,6 +278,7 @@ function _assess_local_identifiability( p::Float64 = 0.99, type = :SE, trbasis = nothing, + known_ic::Array{<:Any, 1} = Array{Any, 1}(), ) where {P <: MPolyElem{Nemo.fmpq}} if isempty(funcs_to_check) funcs_to_check = ode.parameters diff --git a/test/check_algebraicity.jl b/test/check_algebraicity.jl new file mode 100644 index 000000000..842789d09 --- /dev/null +++ b/test/check_algebraicity.jl @@ -0,0 +1,20 @@ +cases = [] + +R, (x, y, z) = PolynomialRing(QQ, ["x", "y", "z"]) +push!(cases, Dict( + :F => RationalFunctionField([[one(R), x + y, x * y]]), + :funcs => [x // one(R), z // one(R), x^3 - y^3 // one(R), x + z // one(R)], + :correct => [true, false, true, false], +)) + +push!(cases, Dict( + :F => RationalFunctionField([[x, y], [y, z]]), + :funcs => [x // z, (x + y) // z, x // one(R), y // one(R), z // one(R)], + :correct => [true, true, false, false, false], +)) + +@testset "Algebraicity over a field" begin + for case in cases + @test check_algebraicity(case[:F], case[:funcs], 0.99) == case[:correct] + end +end diff --git a/test/runtests.jl b/test/runtests.jl index c98edd821..446fcbd07 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,6 +11,7 @@ using StructuralIdentifiability: check_identifiability, check_primality_zerodim, check_primality, + check_algebraicity, det_minor_expansion, ExpVectTrie, get_max_below, From e720f37731eca21667a4746077524ac77552f78f Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Wed, 22 Nov 2023 00:53:00 +0100 Subject: [PATCH 04/12] assess_identifiability_kic but no tests yet --- src/known_ic.jl | 74 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 73 insertions(+), 1 deletion(-) diff --git a/src/known_ic.jl b/src/known_ic.jl index 6ba2d3572..4043478d1 100644 --- a/src/known_ic.jl +++ b/src/known_ic.jl @@ -69,7 +69,6 @@ function _find_identifiable_functions_kic( rational_interpolator = rational_interpolator, seed = seed, ) - @info "The search for identifiable functions with known initial conditions concluded in $((time_ns() - runtime_start) / 1e9) seconds" id_funcs = simplified_generating_set( RationalFunctionField(vcat(id_funcs_general, [f // one(parent(ode)) for f in known_ic])), @@ -79,7 +78,80 @@ function _find_identifiable_functions_kic( rational_interpolator = rational_interpolator, ) + @info "The search for identifiable functions with known initial conditions concluded in $((time_ns() - runtime_start) / 1e9) seconds" + return id_funcs end +""" + assess_identifiability_kic(ode, known_ic; funcs_to_check = [], p=0.99, loglevel=Logging.Info) + +Input: +- `ode` - the ODE model +- `known_ic` - a list of functions for which initial conditions are assumed to be known and generic +- `funcs_to_check` - list of functions to check identifiability for; if empty, all parameters + and states are taken +- `p` - probability of correctness. +- `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) + +Assesses identifiability of parameters and initial conditions of a given ODE model. +The result is guaranteed to be correct with the probability at least `p`. +The function returns an (ordered) dictionary from the functions to check to their identifiability properties +(one of `:nonidentifiable`, `:locally`, `:globally`). +""" +function assess_identifiability_kic( + ode::ODE{P}, + known_ic::Vector{<: Union{P, Generic.Frac{P}}}; + funcs_to_check = Vector(), + p::Float64 = 0.99, + loglevel = Logging.Info, +) where {P <: MPolyElem{fmpq}} + restart_logging(loglevel = loglevel) + reset_timings() + with_logger(_si_logger[]) do + return _assess_identifiability_kic( + ode, + known_ic, + p = p, + funcs_to_check = funcs_to_check, + ) + end +end + +function _assess_identifiability_kic( + ode::ODE{P}, + known_ic::Vector{<: Union{P, Generic.Frac{P}}}; + funcs_to_check = Vector(), + p::Float64 = 0.99, +) where {P <: MPolyElem{fmpq}} + runtime_start = time_ns() + if length(funcs_to_check) == 0 + funcs_to_check = vcat(ode.x_vars, ode.parameters) + end + half_p = 0.5 + p / 2 + id_funcs = _find_identifiable_functions_kic(ode, known_ic, p = half_p) + result = OrderedDict(f => :globally for f in funcs_to_check) + + half_p = 0.5 + p / 2 + local_result = check_algebraicity( + RationalFunctionField([[denominator(f), numerator(f)] for f in id_funcs]), + [f // one(parent(f)) for f in funcs_to_check], + half_p, + ) + global_result = field_contains( + RationalFunctionField([[denominator(f), numerator(f)] for f in id_funcs]), + [f // one(parent(f)) for f in funcs_to_check], + half_p, + ) + for (i, f) in enumerate(funcs_to_check) + if !local_result[i] + result[f] = :nonidentifiable + elseif !global_result[i] + result[f] = :locally + end + end + @info "Assessing identifiability with known initial conditions concluded in $((time_ns() - runtime_start) / 1e9) seconds" + return result +end + From d7868589e7f46650791f21f961a9b90d3963bc77 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Wed, 22 Nov 2023 14:28:47 +0100 Subject: [PATCH 05/12] tests --- test/known_ic.jl | 29 ++++++++++++++++++++++------- test/runtests.jl | 3 ++- 2 files changed, 24 insertions(+), 8 deletions(-) diff --git a/test/known_ic.jl b/test/known_ic.jl index c64bf5dda..73363f13c 100644 --- a/test/known_ic.jl +++ b/test/known_ic.jl @@ -9,29 +9,44 @@ ode = @ODEmodel( push!(cases, Dict( :ode => ode, :known => [x2], - :correct => [a, b, c, d, x1, x2], + :to_check => [], + :correct_funcs => [a, b, c, d, x1, x2], + :correct_ident => OrderedDict(x => :globally for x in [x1 , x2, a, b, c, d]), )) ode = @ODEmodel( x1'(t) = a + x2(t) + x3(t), - x2'(t) = b, - x3'(t) = c, + x2'(t) = b^2 + c, + x3'(t) = -c, y(t) = x1(t) ) push!(cases, Dict( :ode => ode, :known => [x2, x3], - :correct => [a, b + c, x1, x2, x3], + :to_check => [], + :correct_funcs => [a, b^2, x1, x2, x3], + :correct_ident => OrderedDict(x1 => :globally, x2 => :globally, x3 => :globally, a => :globally, b => :locally, c => :nonidentifiable), )) +push!(cases, Dict( + :ode => ode, + :known => [x2, x3], + :to_check => [b^2, x2 * c], + :correct_funcs => [a, b^2, x1, x2, x3], + :correct_ident => OrderedDict(b^2 => :globally, x2 * c => :nonidentifiable), +)) @testset "Identifiable functions with known generic initial conditions" begin for case in cases ode = case[:ode] known = case[:known] - result = find_identifiable_functions_kic(ode, known) - correct = [f // one(parent(ode)) for f in case[:correct]] - @test Set(result) == Set(correct) + + result_funcs = find_identifiable_functions_kic(ode, known) + correct_funcs = [f // one(parent(ode)) for f in case[:correct_funcs]] + @test Set(result_funcs) == Set(correct_funcs) + + result_ident = assess_identifiability_kic(ode, known, funcs_to_check = case[:to_check]) + @test case[:correct_ident] == result_ident end end diff --git a/test/runtests.jl b/test/runtests.jl index 446fcbd07..a62851fdc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -60,7 +60,8 @@ using StructuralIdentifiability: lie_derivative, states_generators, RationalFunctionField, - find_identifiable_functions_kic + find_identifiable_functions_kic, + assess_identifiability_kic function random_ps(ps_ring, range = 1000) result = zero(ps_ring) From 53aaf298e48f80365152fc04a806cca53ab72e87 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Wed, 22 Nov 2023 16:07:24 +0100 Subject: [PATCH 06/12] returning in the form x(0) --- .../RationalFunctionField.jl | 1 + src/known_ic.jl | 5 ++-- src/util.jl | 23 +++++++++++++++++++ test/known_ic.jl | 4 ++-- test/runtests.jl | 3 ++- 5 files changed, 31 insertions(+), 5 deletions(-) diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index f34729e2c..bd9de26d1 100644 --- a/src/RationalFunctionFields/RationalFunctionField.jl +++ b/src/RationalFunctionFields/RationalFunctionField.jl @@ -204,6 +204,7 @@ function check_algebraicity(field, ratfuncs, p) result = Bool[] for f in ratfuncs + f = parent_ring_change(f, poly_ring(field)) for (j, x) in enumerate(base_vars) J[j, end] = evaluate(derivative(f, x), eval_point) end diff --git a/src/known_ic.jl b/src/known_ic.jl index 4043478d1..a69192486 100644 --- a/src/known_ic.jl +++ b/src/known_ic.jl @@ -80,7 +80,7 @@ function _find_identifiable_functions_kic( @info "The search for identifiable functions with known initial conditions concluded in $((time_ns() - runtime_start) / 1e9) seconds" - return id_funcs + return replace_with_ic(ode, id_funcs) end """ @@ -130,6 +130,7 @@ function _assess_identifiability_kic( end half_p = 0.5 + p / 2 id_funcs = _find_identifiable_functions_kic(ode, known_ic, p = half_p) + funcs_to_check = replace_with_ic(ode, funcs_to_check) result = OrderedDict(f => :globally for f in funcs_to_check) half_p = 0.5 + p / 2 @@ -151,7 +152,7 @@ function _assess_identifiability_kic( end end @info "Assessing identifiability with known initial conditions concluded in $((time_ns() - runtime_start) / 1e9) seconds" - return result + return result end diff --git a/src/util.jl b/src/util.jl index a55d5571c..cbea62bc4 100644 --- a/src/util.jl +++ b/src/util.jl @@ -579,3 +579,26 @@ function get_measured_quantities(ode::ModelingToolkit.ODESystem) ) end end + +# ----------------------------------------------------------------------------- + +""" + replace_with_ic(ode::ODE, funcs) + +Takes an ode and a list of functions in the states and parameters and makes a change of variable +names `x(t) -> x(0)`. Function is used to prepare the output for the case of known initial conditions +""" +function replace_with_ic(ode, funcs) + varnames = [(var_to_str(p), var_to_str(p)) for p in ode.parameters] + for x in ode.x_vars + s = var_to_str(x) + if endswith(s, "(t)") + push!(varnames, (s, s[1:(end - 3)] * "(0)")) + else + push!(varnames, (s, s * "(0)")) + end + end + R0, vars0 = PolynomialRing(base_ring(ode.poly_ring), [p[2] for p in varnames]) + eval_dict = Dict(str_to_var(p[1], ode.poly_ring) => str_to_var(p[2], R0) for p in varnames) + return [eval_at_dict(f, eval_dict) for f in funcs] +end diff --git a/test/known_ic.jl b/test/known_ic.jl index 73363f13c..344c32286 100644 --- a/test/known_ic.jl +++ b/test/known_ic.jl @@ -43,10 +43,10 @@ push!(cases, Dict( known = case[:known] result_funcs = find_identifiable_functions_kic(ode, known) - correct_funcs = [f // one(parent(ode)) for f in case[:correct_funcs]] + correct_funcs = replace_with_ic(ode, [f // one(parent(ode)) for f in case[:correct_funcs]]) @test Set(result_funcs) == Set(correct_funcs) result_ident = assess_identifiability_kic(ode, known, funcs_to_check = case[:to_check]) - @test case[:correct_ident] == result_ident + @test OrderedDict(replace_with_ic(ode, [k])[1] => v for (k, v) in case[:correct_ident]) == result_ident end end diff --git a/test/runtests.jl b/test/runtests.jl index a62851fdc..ee88fa79e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -61,7 +61,8 @@ using StructuralIdentifiability: states_generators, RationalFunctionField, find_identifiable_functions_kic, - assess_identifiability_kic + assess_identifiability_kic, + replace_with_ic function random_ps(ps_ring, range = 1000) result = zero(ps_ring) From 37d72f47aa3eeb538005a1e839ced7c531d25e2b Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Wed, 22 Nov 2023 16:10:39 +0100 Subject: [PATCH 07/12] formatting --- .../RationalFunctionField.jl | 5 +- src/known_ic.jl | 14 ++-- src/util.jl | 3 +- test/check_algebraicity.jl | 26 ++++--- test/known_ic.jl | 75 +++++++++++-------- 5 files changed, 74 insertions(+), 49 deletions(-) diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index bd9de26d1..e5af1f040 100644 --- a/src/RationalFunctionFields/RationalFunctionField.jl +++ b/src/RationalFunctionFields/RationalFunctionField.jl @@ -187,7 +187,10 @@ Output: function check_algebraicity(field, ratfuncs, p) fgens = generators(field) base_vars = gens(poly_ring(field)) - maxdeg = maximum([max(total_degree(numerator(f)), total_degree(denominator(f))) for f in vcat(ratfuncs, fgens)]) + maxdeg = maximum([ + max(total_degree(numerator(f)), total_degree(denominator(f))) for + f in vcat(ratfuncs, fgens) + ]) # degree of the polynomial whose nonvanishing will be needed for correct result D = Int(ceil(2 * maxdeg * (length(fgens) + 1)^3 * length(ratfuncs) / (1 - p))) eval_point = [Nemo.QQ(rand(1:D)) for x in base_vars] diff --git a/src/known_ic.jl b/src/known_ic.jl index a69192486..7dbaca06f 100644 --- a/src/known_ic.jl +++ b/src/known_ic.jl @@ -28,7 +28,7 @@ This functions takes the following optional arguments: """ function find_identifiable_functions_kic( ode::ODE{T}, - known_ic::Vector{<: Union{T, Generic.Frac{T}}}; + known_ic::Vector{<:Union{T, Generic.Frac{T}}}; p::Float64 = 0.99, seed = 42, simplify = :standard, @@ -51,7 +51,7 @@ end function _find_identifiable_functions_kic( ode::ODE{T}, - known_ic::Vector{<: Union{T, Generic.Frac{T}}}; + known_ic::Vector{<:Union{T, Generic.Frac{T}}}; p::Float64 = 0.99, seed = 42, simplify = :standard, @@ -71,7 +71,9 @@ function _find_identifiable_functions_kic( ) id_funcs = simplified_generating_set( - RationalFunctionField(vcat(id_funcs_general, [f // one(parent(ode)) for f in known_ic])), + RationalFunctionField( + vcat(id_funcs_general, [f // one(parent(ode)) for f in known_ic]), + ), p = half_p, seed = seed, simplify = simplify, @@ -101,7 +103,7 @@ The function returns an (ordered) dictionary from the functions to check to thei """ function assess_identifiability_kic( ode::ODE{P}, - known_ic::Vector{<: Union{P, Generic.Frac{P}}}; + known_ic::Vector{<:Union{P, Generic.Frac{P}}}; funcs_to_check = Vector(), p::Float64 = 0.99, loglevel = Logging.Info, @@ -120,7 +122,7 @@ end function _assess_identifiability_kic( ode::ODE{P}, - known_ic::Vector{<: Union{P, Generic.Frac{P}}}; + known_ic::Vector{<:Union{P, Generic.Frac{P}}}; funcs_to_check = Vector(), p::Float64 = 0.99, ) where {P <: MPolyElem{fmpq}} @@ -154,5 +156,3 @@ function _assess_identifiability_kic( @info "Assessing identifiability with known initial conditions concluded in $((time_ns() - runtime_start) / 1e9) seconds" return result end - - diff --git a/src/util.jl b/src/util.jl index cbea62bc4..38eb9d9e2 100644 --- a/src/util.jl +++ b/src/util.jl @@ -599,6 +599,7 @@ function replace_with_ic(ode, funcs) end end R0, vars0 = PolynomialRing(base_ring(ode.poly_ring), [p[2] for p in varnames]) - eval_dict = Dict(str_to_var(p[1], ode.poly_ring) => str_to_var(p[2], R0) for p in varnames) + eval_dict = + Dict(str_to_var(p[1], ode.poly_ring) => str_to_var(p[2], R0) for p in varnames) return [eval_at_dict(f, eval_dict) for f in funcs] end diff --git a/test/check_algebraicity.jl b/test/check_algebraicity.jl index 842789d09..ecb742f1a 100644 --- a/test/check_algebraicity.jl +++ b/test/check_algebraicity.jl @@ -1,17 +1,23 @@ cases = [] R, (x, y, z) = PolynomialRing(QQ, ["x", "y", "z"]) -push!(cases, Dict( - :F => RationalFunctionField([[one(R), x + y, x * y]]), - :funcs => [x // one(R), z // one(R), x^3 - y^3 // one(R), x + z // one(R)], - :correct => [true, false, true, false], -)) +push!( + cases, + Dict( + :F => RationalFunctionField([[one(R), x + y, x * y]]), + :funcs => [x // one(R), z // one(R), x^3 - y^3 // one(R), x + z // one(R)], + :correct => [true, false, true, false], + ), +) -push!(cases, Dict( - :F => RationalFunctionField([[x, y], [y, z]]), - :funcs => [x // z, (x + y) // z, x // one(R), y // one(R), z // one(R)], - :correct => [true, true, false, false, false], -)) +push!( + cases, + Dict( + :F => RationalFunctionField([[x, y], [y, z]]), + :funcs => [x // z, (x + y) // z, x // one(R), y // one(R), z // one(R)], + :correct => [true, true, false, false, false], + ), +) @testset "Algebraicity over a field" begin for case in cases diff --git a/test/known_ic.jl b/test/known_ic.jl index 344c32286..4c59df835 100644 --- a/test/known_ic.jl +++ b/test/known_ic.jl @@ -6,36 +6,47 @@ ode = @ODEmodel( y(t) = x1(t) ) -push!(cases, Dict( - :ode => ode, - :known => [x2], - :to_check => [], - :correct_funcs => [a, b, c, d, x1, x2], - :correct_ident => OrderedDict(x => :globally for x in [x1 , x2, a, b, c, d]), -)) +push!( + cases, + Dict( + :ode => ode, + :known => [x2], + :to_check => [], + :correct_funcs => [a, b, c, d, x1, x2], + :correct_ident => OrderedDict(x => :globally for x in [x1, x2, a, b, c, d]), + ), +) -ode = @ODEmodel( - x1'(t) = a + x2(t) + x3(t), - x2'(t) = b^2 + c, - x3'(t) = -c, - y(t) = x1(t) +ode = @ODEmodel(x1'(t) = a + x2(t) + x3(t), x2'(t) = b^2 + c, x3'(t) = -c, y(t) = x1(t)) + +push!( + cases, + Dict( + :ode => ode, + :known => [x2, x3], + :to_check => [], + :correct_funcs => [a, b^2, x1, x2, x3], + :correct_ident => OrderedDict( + x1 => :globally, + x2 => :globally, + x3 => :globally, + a => :globally, + b => :locally, + c => :nonidentifiable, + ), + ), ) -push!(cases, Dict( - :ode => ode, - :known => [x2, x3], - :to_check => [], - :correct_funcs => [a, b^2, x1, x2, x3], - :correct_ident => OrderedDict(x1 => :globally, x2 => :globally, x3 => :globally, a => :globally, b => :locally, c => :nonidentifiable), -)) - -push!(cases, Dict( - :ode => ode, - :known => [x2, x3], - :to_check => [b^2, x2 * c], - :correct_funcs => [a, b^2, x1, x2, x3], - :correct_ident => OrderedDict(b^2 => :globally, x2 * c => :nonidentifiable), -)) +push!( + cases, + Dict( + :ode => ode, + :known => [x2, x3], + :to_check => [b^2, x2 * c], + :correct_funcs => [a, b^2, x1, x2, x3], + :correct_ident => OrderedDict(b^2 => :globally, x2 * c => :nonidentifiable), + ), +) @testset "Identifiable functions with known generic initial conditions" begin for case in cases @@ -43,10 +54,14 @@ push!(cases, Dict( known = case[:known] result_funcs = find_identifiable_functions_kic(ode, known) - correct_funcs = replace_with_ic(ode, [f // one(parent(ode)) for f in case[:correct_funcs]]) + correct_funcs = + replace_with_ic(ode, [f // one(parent(ode)) for f in case[:correct_funcs]]) @test Set(result_funcs) == Set(correct_funcs) - result_ident = assess_identifiability_kic(ode, known, funcs_to_check = case[:to_check]) - @test OrderedDict(replace_with_ic(ode, [k])[1] => v for (k, v) in case[:correct_ident]) == result_ident + result_ident = + assess_identifiability_kic(ode, known, funcs_to_check = case[:to_check]) + @test OrderedDict( + replace_with_ic(ode, [k])[1] => v for (k, v) in case[:correct_ident] + ) == result_ident end end From b0fc44369955d1d6c4c302b907ab197a6c900203 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Thu, 23 Nov 2023 22:22:16 +0100 Subject: [PATCH 08/12] using dennums_to_fracs --- src/RationalFunctionFields/RationalFunctionField.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index e5af1f040..365e1f493 100644 --- a/src/RationalFunctionFields/RationalFunctionField.jl +++ b/src/RationalFunctionFields/RationalFunctionField.jl @@ -42,7 +42,7 @@ function poly_ring(F::RationalFunctionField) end function generators(F::RationalFunctionField) - return vcat([[f[i] // f[1] for i in 2:length(f)] for f in F.dennums]...) + return dennums_to_fractions(F.dennums) end # ------------------------------------------------------------------------------ From dab25d5639d4ecfbf87a8bd1df201167a975dbb4 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Sun, 4 Feb 2024 23:51:38 +0100 Subject: [PATCH 09/12] fixing the tests --- src/known_ic.jl | 32 ++++++++++++++++---------------- src/util.jl | 17 ----------------- 2 files changed, 16 insertions(+), 33 deletions(-) diff --git a/src/known_ic.jl b/src/known_ic.jl index 7dbaca06f..c5874cd37 100644 --- a/src/known_ic.jl +++ b/src/known_ic.jl @@ -16,7 +16,7 @@ This functions takes the following optional arguments: - `:strong`: Strong simplification. This option is the slowest, but the output functions are nice and simple. - `:absent`: No simplification. -- `p`: A float in the range from 0 to 1, the probability of correctness. Default +- `prob_threshold`: A float in the range from 0 to 1, the probability of correctness. Default is `0.99`. - `seed`: The rng seed. Default value is `42`. - `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) @@ -29,7 +29,7 @@ This functions takes the following optional arguments: function find_identifiable_functions_kic( ode::ODE{T}, known_ic::Vector{<:Union{T, Generic.Frac{T}}}; - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, seed = 42, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, @@ -41,7 +41,7 @@ function find_identifiable_functions_kic( return _find_identifiable_functions_kic( ode, known_ic, - p = p, + prob_threshold = prob_threshold, seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, @@ -52,18 +52,18 @@ end function _find_identifiable_functions_kic( ode::ODE{T}, known_ic::Vector{<:Union{T, Generic.Frac{T}}}; - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, seed = 42, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, ) where {T <: MPolyElem{fmpq}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) - half_p = 0.5 + p / 2 + half_p = 0.5 + prob_threshold / 2 runtime_start = time_ns() id_funcs_general = find_identifiable_functions( ode, - p = half_p, + prob_threshold = half_p, with_states = true, simplify = simplify, rational_interpolator = rational_interpolator, @@ -74,7 +74,7 @@ function _find_identifiable_functions_kic( RationalFunctionField( vcat(id_funcs_general, [f // one(parent(ode)) for f in known_ic]), ), - p = half_p, + prob_threshold = half_p, seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, @@ -86,18 +86,18 @@ function _find_identifiable_functions_kic( end """ - assess_identifiability_kic(ode, known_ic; funcs_to_check = [], p=0.99, loglevel=Logging.Info) + assess_identifiability_kic(ode, known_ic; funcs_to_check = [], prob_threshold=0.99, loglevel=Logging.Info) Input: - `ode` - the ODE model - `known_ic` - a list of functions for which initial conditions are assumed to be known and generic - `funcs_to_check` - list of functions to check identifiability for; if empty, all parameters and states are taken -- `p` - probability of correctness. +- `prob_threshold` - probability of correctness. - `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) Assesses identifiability of parameters and initial conditions of a given ODE model. -The result is guaranteed to be correct with the probability at least `p`. +The result is guaranteed to be correct with the probability at least `prob_threshold`. The function returns an (ordered) dictionary from the functions to check to their identifiability properties (one of `:nonidentifiable`, `:locally`, `:globally`). """ @@ -105,7 +105,7 @@ function assess_identifiability_kic( ode::ODE{P}, known_ic::Vector{<:Union{P, Generic.Frac{P}}}; funcs_to_check = Vector(), - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, loglevel = Logging.Info, ) where {P <: MPolyElem{fmpq}} restart_logging(loglevel = loglevel) @@ -114,7 +114,7 @@ function assess_identifiability_kic( return _assess_identifiability_kic( ode, known_ic, - p = p, + prob_threshold = prob_threshold, funcs_to_check = funcs_to_check, ) end @@ -124,18 +124,18 @@ function _assess_identifiability_kic( ode::ODE{P}, known_ic::Vector{<:Union{P, Generic.Frac{P}}}; funcs_to_check = Vector(), - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, ) where {P <: MPolyElem{fmpq}} runtime_start = time_ns() if length(funcs_to_check) == 0 funcs_to_check = vcat(ode.x_vars, ode.parameters) end - half_p = 0.5 + p / 2 - id_funcs = _find_identifiable_functions_kic(ode, known_ic, p = half_p) + half_p = 0.5 + prob_threshold / 2 + id_funcs = _find_identifiable_functions_kic(ode, known_ic, prob_threshold = half_p) funcs_to_check = replace_with_ic(ode, funcs_to_check) result = OrderedDict(f => :globally for f in funcs_to_check) - half_p = 0.5 + p / 2 + half_p = 0.5 + half_p / 2 local_result = check_algebraicity( RationalFunctionField([[denominator(f), numerator(f)] for f in id_funcs]), [f // one(parent(f)) for f in funcs_to_check], diff --git a/src/util.jl b/src/util.jl index c85533365..d68a59628 100644 --- a/src/util.jl +++ b/src/util.jl @@ -518,22 +518,6 @@ function difforder(diffpoly::MPolyRingElem, prefix::String) return max(orders...) end -function get_measured_quantities(ode::ModelingToolkit.ODESystem) - if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(ode)) - @info "Measured quantities are not provided, trying to find the outputs in input ODE." - return filter( - eq -> (ModelingToolkit.isoutput(eq.lhs)), - ModelingToolkit.equations(ode), - ) - else - throw( - error( - "Measured quantities (output functions) were not provided and no outputs were found.", - ), - ) - end -end - # ----------------------------------------------------------------------------- """ @@ -557,4 +541,3 @@ function replace_with_ic(ode, funcs) Dict(str_to_var(p[1], ode.poly_ring) => str_to_var(p[2], R0) for p in varnames) return [eval_at_dict(f, eval_dict) for f in funcs] end - From 762c568dc20b2032f99c4f544e0320bc4c78bb94 Mon Sep 17 00:00:00 2001 From: gpogudin Date: Mon, 5 Feb 2024 13:21:48 +0100 Subject: [PATCH 10/12] refactoring the dispatch --- src/StructuralIdentifiability.jl | 23 +++++++++++---- src/identifiable_functions.jl | 31 +++++++++++++++------ src/known_ic.jl | 48 ++------------------------------ test/known_ic.jl | 17 +++++++++-- test/runtests.jl | 2 -- 5 files changed, 59 insertions(+), 62 deletions(-) diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index a0d15e336..1b709a60b 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -88,6 +88,9 @@ Input: - `ode` - the ODE model - `funcs_to_check` - list of functions to check identifiability for; if empty, all parameters and states are taken +- `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). - `prob_threshold` - probability of correctness. - `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) @@ -99,17 +102,27 @@ The function returns an (ordered) dictionary from the functions to check to thei function assess_identifiability( ode::ODE{P}; funcs_to_check = Vector(), + known_ic::Vector{<:Union{P, Generic.Frac{P}}} = Vector{Union{P, Generic.Frac{P}}}(), prob_threshold::Float64 = 0.99, loglevel = Logging.Info, ) where {P <: MPolyRingElem{QQFieldElem}} restart_logging(loglevel = loglevel) reset_timings() with_logger(_si_logger[]) do - return _assess_identifiability( - ode, - funcs_to_check = funcs_to_check, - prob_threshold = prob_threshold, - ) + if isempty(known_ic) + return _assess_identifiability( + ode, + funcs_to_check = funcs_to_check, + prob_threshold = prob_threshold, + ) + else + return _assess_identifiability_kic( + ode, + known_ic, + funcs_to_check = funcs_to_check, + prob_threshold = prob_threshold, + ) + end end end diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index ec116843e..743c7fba5 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -18,6 +18,9 @@ This functions takes the following optional arguments: - `:strong`: Strong simplification. This option is the slowest, but the output functions are nice and simple. - `:absent`: No simplification. +- `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). - `prob_threshold`: A float in the range from 0 to 1, the probability of correctness. Default is `0.99`. - `seed`: The rng seed. Default value is `42`. @@ -45,6 +48,7 @@ find_identifiable_functions(ode) """ function find_identifiable_functions( ode::ODE{T}; + known_ic::Vector{<:Union{T, Generic.Frac{T}}} = Vector{Union{T, Generic.Frac{T}}}(), prob_threshold::Float64 = 0.99, seed = 42, with_states = false, @@ -55,14 +59,25 @@ function find_identifiable_functions( restart_logging(loglevel = loglevel) reset_timings() with_logger(_si_logger[]) do - return _find_identifiable_functions( - ode, - prob_threshold = prob_threshold, - seed = seed, - with_states = with_states, - simplify = simplify, - rational_interpolator = rational_interpolator, - ) + if isempty(known_ic) + return _find_identifiable_functions( + ode, + prob_threshold = prob_threshold, + seed = seed, + with_states = with_states, + simplify = simplify, + rational_interpolator = rational_interpolator, + ) + else + return _find_identifiable_functions_kic( + ode, + known_ic, + prob_threshold = prob_threshold, + seed = seed, + simplify = simplify, + rational_interpolator = rational_interpolator, + ) + end end end diff --git a/src/known_ic.jl b/src/known_ic.jl index c5874cd37..1fee878cb 100644 --- a/src/known_ic.jl +++ b/src/known_ic.jl @@ -1,5 +1,5 @@ """ - find_identifiable_functions_kic(ode::ODE, known_ic; options...) + _find_identifiable_functions_kic(ode::ODE, known_ic; options...) Finds all functions of parameters/initial conditions that are identifiable in the given ODE system under assumptions that the initial conditions for functions in the list @@ -26,29 +26,6 @@ This functions takes the following optional arguments: ``` """ -function find_identifiable_functions_kic( - ode::ODE{T}, - known_ic::Vector{<:Union{T, Generic.Frac{T}}}; - prob_threshold::Float64 = 0.99, - seed = 42, - simplify = :standard, - rational_interpolator = :VanDerHoevenLecerf, - loglevel = Logging.Info, -) where {T <: MPolyElem{fmpq}} - restart_logging(loglevel = loglevel) - reset_timings() - with_logger(_si_logger[]) do - return _find_identifiable_functions_kic( - ode, - known_ic, - prob_threshold = prob_threshold, - seed = seed, - simplify = simplify, - rational_interpolator = rational_interpolator, - ) - end -end - function _find_identifiable_functions_kic( ode::ODE{T}, known_ic::Vector{<:Union{T, Generic.Frac{T}}}; @@ -61,7 +38,7 @@ function _find_identifiable_functions_kic( @assert simplify in (:standard, :weak, :strong, :absent) half_p = 0.5 + prob_threshold / 2 runtime_start = time_ns() - id_funcs_general = find_identifiable_functions( + id_funcs_general = _find_identifiable_functions( ode, prob_threshold = half_p, with_states = true, @@ -86,7 +63,7 @@ function _find_identifiable_functions_kic( end """ - assess_identifiability_kic(ode, known_ic; funcs_to_check = [], prob_threshold=0.99, loglevel=Logging.Info) + _assess_identifiability_kic(ode; known_ic, funcs_to_check = [], prob_threshold=0.99, loglevel=Logging.Info) Input: - `ode` - the ODE model @@ -101,25 +78,6 @@ The result is guaranteed to be correct with the probability at least `prob_thres The function returns an (ordered) dictionary from the functions to check to their identifiability properties (one of `:nonidentifiable`, `:locally`, `:globally`). """ -function assess_identifiability_kic( - ode::ODE{P}, - known_ic::Vector{<:Union{P, Generic.Frac{P}}}; - funcs_to_check = Vector(), - prob_threshold::Float64 = 0.99, - loglevel = Logging.Info, -) where {P <: MPolyElem{fmpq}} - restart_logging(loglevel = loglevel) - reset_timings() - with_logger(_si_logger[]) do - return _assess_identifiability_kic( - ode, - known_ic, - prob_threshold = prob_threshold, - funcs_to_check = funcs_to_check, - ) - end -end - function _assess_identifiability_kic( ode::ODE{P}, known_ic::Vector{<:Union{P, Generic.Frac{P}}}; diff --git a/test/known_ic.jl b/test/known_ic.jl index 4c59df835..bf26d1555 100644 --- a/test/known_ic.jl +++ b/test/known_ic.jl @@ -48,18 +48,31 @@ push!( ), ) +ode = @ODEmodel(x1'(t) = a * x1(t), x2'(t) = x1(t) + 1 / x1(t), y(t) = x2(t)) + +push!( + cases, + Dict( + :ode => ode, + :known => [x1], + :to_check => [], + :correct_funcs => [a, x1, x2], + :correct_ident => OrderedDict(x1 => :globally, x2 => :globally, a => :globally), + ), +) + @testset "Identifiable functions with known generic initial conditions" begin for case in cases ode = case[:ode] known = case[:known] - result_funcs = find_identifiable_functions_kic(ode, known) + result_funcs = find_identifiable_functions(ode, known_ic = known) correct_funcs = replace_with_ic(ode, [f // one(parent(ode)) for f in case[:correct_funcs]]) @test Set(result_funcs) == Set(correct_funcs) result_ident = - assess_identifiability_kic(ode, known, funcs_to_check = case[:to_check]) + assess_identifiability(ode, known_ic = known, funcs_to_check = case[:to_check]) @test OrderedDict( replace_with_ic(ode, [k])[1] => v for (k, v) in case[:correct_ident] ) == result_ident diff --git a/test/runtests.jl b/test/runtests.jl index 1fdaa5a2b..509365125 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -60,8 +60,6 @@ using StructuralIdentifiability: lie_derivative, states_generators, RationalFunctionField, - find_identifiable_functions_kic, - assess_identifiability_kic, replace_with_ic, x_vars, y_vars, From bf930bb9c7f1d3384ae0aa1e1382ccd916ff4ac6 Mon Sep 17 00:00:00 2001 From: gpogudin Date: Mon, 5 Feb 2024 13:52:26 +0100 Subject: [PATCH 11/12] adding docs for known_ic --- docs/src/tutorials/identifiability.md | 13 ++++++++ docs/src/tutorials/identifiable_functions.md | 11 +++++++ test/known_ic.jl | 31 ++++++++++++++++++++ 3 files changed, 55 insertions(+) diff --git a/docs/src/tutorials/identifiability.md b/docs/src/tutorials/identifiability.md index c88ba3d94..a5e0e4888 100644 --- a/docs/src/tutorials/identifiability.md +++ b/docs/src/tutorials/identifiability.md @@ -123,6 +123,19 @@ And we see that they indeed are. This means, in particular, that the reason why can be exchanged. One may wonder how could we guess these functions `beta + delta, beta * delta`. In fact, they can be just computed using `find_identifiable_functions` function as we will explain in the next tutorial. Stay tuned! +## Assuming known initial conditions + +An experimental feature allows to provide an additional keyword argument `known_ic` to inidcate functions of states and parameters for which the +initial conditions are assumed to be known (while the initial conditions of the system are still assumed to be genric). In this case, +the identifiability will be assessed for parameters and all the initial conditions or for the initial conditions of `funcs_to_check`. +Let us add an assumption that the initial conditions `x2(0)` and `x3(0)` are known: + +```@example global +assess_identifiability(ode, known_ic = [x2, x3]) +``` + +And we see that now `alpha` and `gama` become locally identifiable. + [^1]: > A. Sedoglavic, [*A probabilistic algorithm to test local algebraic observability in polynomial time*](https://doi.org/10.1006/jsco.2002.0532), Journal of Symbolic Computation, 2002. [^2]: > A. Ovchinnikov, A. Pillay, G. Pogudin, T. Scanlon, [*Multi-experiment Parameter Identifiability of ODEs and Model Theory*](https://doi.org/10.1137/21M1389845), SIAM Journal on Applied Algebra and Geometry, 2022. [^3]: > D. Gonze, P. Ruoff, [*The Goodwin Oscillator and its Legacy*](https://doi.org/10.1007/s10441-020-09379-8), Acta Biotheoretica, 2020. diff --git a/docs/src/tutorials/identifiable_functions.md b/docs/src/tutorials/identifiable_functions.md index c7e830dca..1e1986385 100644 --- a/docs/src/tutorials/identifiable_functions.md +++ b/docs/src/tutorials/identifiable_functions.md @@ -43,4 +43,15 @@ the degree of simplification. The default value is `:standard` but one could use As `assess_identifiability` and `assess_local_identifiability`, `find_identifiable_functions` accepts an optional parameter `loglevel` (default: `Logging.Info`) to adjust the verbosity of logging. +Finally, as for `assess_identifiability`, an experimental feature allows to provide an additional keyword argument `known_ic` to inidcate functions of states and parameters for which the +initial conditions are assumed to be known (while the initial conditions of the system are still assumed to be genric). +In this case, the function will find identifiable functions of parameters and initial conditions rather than of parameters and states. +Let us add an assumption that the initial condition `x1(0)` is known: + +```@example funcs +find_identifiable_functions(LLW1987, known_ic = [x1]) +``` + +We see that `x2(0)` becomes an identifiable function as well, which is natural since `x1(t) * x2(t)` was an identifiable function before. + [^1]: > Y. Lecourtier, F. Lamnabhi-Lagarrigue, and E. Walter, [*A method to prove that nonlinear models can be unidentifiable*](https://doi.org/10.1109/CDC.1987.272467), Proceedings of the 26th Conference on Decision and Control, December 1987, 2144-2145; diff --git a/test/known_ic.jl b/test/known_ic.jl index bf26d1555..ae098ae1d 100644 --- a/test/known_ic.jl +++ b/test/known_ic.jl @@ -61,6 +61,37 @@ push!( ), ) +ode = @ODEmodel( + x1'(t) = -b * x1(t) + 1 / (c + x4(t)), + x2'(t) = alpha * x1(t) - beta * x2(t), + x3'(t) = gama * x2(t) - delta * x3(t), + x4'(t) = sigma * x4(t) * (gama * x2(t) - delta * x3(t)) / x3(t), + y(t) = x1(t) +) + +push!( + cases, + Dict( + :ode => ode, + :known => [x2, x3], + :to_check => [alpha, alpha * gama], + :correct_funcs => [ + sigma, + c, + b, + x4, + x3, + x2, + x1, + beta * delta, + alpha * gama, + beta + delta, + -delta * x3 + gama * x2, + ], + :correct_ident => OrderedDict(alpha => :locally, alpha * gama => :globally), + ), +) + @testset "Identifiable functions with known generic initial conditions" begin for case in cases ode = case[:ode] From 3b4d971bc02a442a6509fcf240f242d9439078f7 Mon Sep 17 00:00:00 2001 From: gpogudin Date: Mon, 5 Feb 2024 15:34:19 +0100 Subject: [PATCH 12/12] fixing typos --- docs/src/tutorials/identifiability.md | 2 +- docs/src/tutorials/identifiable_functions.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/tutorials/identifiability.md b/docs/src/tutorials/identifiability.md index a5e0e4888..dc77cc59a 100644 --- a/docs/src/tutorials/identifiability.md +++ b/docs/src/tutorials/identifiability.md @@ -126,7 +126,7 @@ can be exchanged. One may wonder how could we guess these functions `beta + delt ## Assuming known initial conditions An experimental feature allows to provide an additional keyword argument `known_ic` to inidcate functions of states and parameters for which the -initial conditions are assumed to be known (while the initial conditions of the system are still assumed to be genric). In this case, +initial conditions are assumed to be known (while the initial conditions of the system are still assumed to be generic). In this case, the identifiability will be assessed for parameters and all the initial conditions or for the initial conditions of `funcs_to_check`. Let us add an assumption that the initial conditions `x2(0)` and `x3(0)` are known: diff --git a/docs/src/tutorials/identifiable_functions.md b/docs/src/tutorials/identifiable_functions.md index 1e1986385..c24cd0d85 100644 --- a/docs/src/tutorials/identifiable_functions.md +++ b/docs/src/tutorials/identifiable_functions.md @@ -44,7 +44,7 @@ As `assess_identifiability` and `assess_local_identifiability`, `find_identifiab to adjust the verbosity of logging. Finally, as for `assess_identifiability`, an experimental feature allows to provide an additional keyword argument `known_ic` to inidcate functions of states and parameters for which the -initial conditions are assumed to be known (while the initial conditions of the system are still assumed to be genric). +initial conditions are assumed to be known (while the initial conditions of the system are still assumed to be generic). In this case, the function will find identifiable functions of parameters and initial conditions rather than of parameters and states. Let us add an assumption that the initial condition `x1(0)` is known: