Merge pull request #249 from SciML/known_ic_generic

WIP: Identifiability with known generic initial conditions

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.

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;

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},

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
 

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
 

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

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

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

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