diff --git a/docs/src/tutorials/identifiability.md b/docs/src/tutorials/identifiability.md index c88ba3d9..dc77cc59 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 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: + +```@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 c7e830dc..c24cd0d8 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 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: + +```@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/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index fbfea666..f2583d9c 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 dennums_to_fractions(F.dennums) +end + # ------------------------------------------------------------------------------ """ @@ -169,6 +173,55 @@ 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 + 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 + push!(result, LinearAlgebra.rank(J) == rank) + end + return result +end + +# ------------------------------------------------------------------------------ + function issubfield( F::RationalFunctionField{T}, E::RationalFunctionField{T}, diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index 66713cca..1b709a60 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -70,6 +70,7 @@ include("lincomp.jl") include("pb_representation.jl") include("submodels.jl") include("discrete.jl") +include("known_ic.jl") include("input_macro.jl") function __init__() @@ -87,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) @@ -98,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 ec116843..743c7fba 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 new file mode 100644 index 00000000..1fee878c --- /dev/null +++ b/src/known_ic.jl @@ -0,0 +1,116 @@ +""" + _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. +- `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) + +**This is experimental functionality** + +``` + +""" +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, +) where {T <: MPolyElem{fmpq}} + Random.seed!(seed) + @assert simplify in (:standard, :weak, :strong, :absent) + half_p = 0.5 + prob_threshold / 2 + runtime_start = time_ns() + id_funcs_general = _find_identifiable_functions( + ode, + prob_threshold = half_p, + with_states = true, + simplify = simplify, + rational_interpolator = rational_interpolator, + seed = seed, + ) + + id_funcs = simplified_generating_set( + RationalFunctionField( + vcat(id_funcs_general, [f // one(parent(ode)) for f in known_ic]), + ), + prob_threshold = half_p, + seed = seed, + simplify = simplify, + rational_interpolator = rational_interpolator, + ) + + @info "The search for identifiable functions with known initial conditions concluded in $((time_ns() - runtime_start) / 1e9) seconds" + + return replace_with_ic(ode, id_funcs) +end + +""" + _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 +- `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 `prob_threshold`. +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, +) 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 + 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 + 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], + 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 diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index a42a2409..8b9c8998 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -180,6 +180,7 @@ function _assess_local_identifiability( prob_threshold::Float64 = 0.99, type = :SE, trbasis = nothing, + known_ic::Array{<:Any, 1} = Array{Any, 1}(), ) where {P <: MPolyRingElem{Nemo.QQFieldElem}} if isempty(funcs_to_check) funcs_to_check = ode.parameters diff --git a/src/util.jl b/src/util.jl index 7acc3bdf..d68a5962 100644 --- a/src/util.jl +++ b/src/util.jl @@ -517,3 +517,27 @@ function difforder(diffpoly::MPolyRingElem, prefix::String) end return max(orders...) 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/check_algebraicity.jl b/test/check_algebraicity.jl new file mode 100644 index 00000000..ecb742f1 --- /dev/null +++ b/test/check_algebraicity.jl @@ -0,0 +1,26 @@ +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/known_ic.jl b/test/known_ic.jl new file mode 100644 index 00000000..ae098ae1 --- /dev/null +++ b/test/known_ic.jl @@ -0,0 +1,111 @@ +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], + :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)) + +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), + ), +) + +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), + ), +) + +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] + known = case[: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(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 + end +end diff --git a/test/runtests.jl b/test/runtests.jl index e5ed9fe7..50936512 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, @@ -59,6 +60,7 @@ using StructuralIdentifiability: lie_derivative, states_generators, RationalFunctionField, + replace_with_ic, x_vars, y_vars, x_equations,