From 3e0840e2119e7b66d38d26e4b6d4f21c00292866 Mon Sep 17 00:00:00 2001 From: gpogudin Date: Thu, 25 Jan 2024 15:11:20 +0100 Subject: [PATCH] Updating docs, especially the DDS part --- docs/Project.toml | 2 +- docs/src/assets/Project.toml | 4 +- docs/src/input/input.md | 8 ++- docs/src/tutorials/creating_ode.md | 4 +- docs/src/tutorials/discrete_time.md | 97 +++++++++++++++++++++-------- src/discrete.jl | 3 + 6 files changed, 86 insertions(+), 32 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 21e7b1811..e449e2a58 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -7,5 +7,5 @@ StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544" [compat] BenchmarkTools = "1.3" Documenter = "0.27, 1" -ModelingToolkit = "8.34" +ModelingToolkit = "8.74" StructuralIdentifiability = "0.5" diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml index e4c519b6b..e449e2a58 100644 --- a/docs/src/assets/Project.toml +++ b/docs/src/assets/Project.toml @@ -7,5 +7,5 @@ StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544" [compat] BenchmarkTools = "1.3" Documenter = "0.27, 1" -ModelingToolkit = "8.34" -StructuralIdentifiability = "0.4" +ModelingToolkit = "8.74" +StructuralIdentifiability = "0.5" diff --git a/docs/src/input/input.md b/docs/src/input/input.md index d41bf7374..5fa6ad68c 100644 --- a/docs/src/input/input.md +++ b/docs/src/input/input.md @@ -1,4 +1,4 @@ -# Parsing input ODE system +# Parsing input system ```@docs @ODEmodel(ex::Expr...) @@ -11,3 +11,9 @@ set_parameter_values ```@docs linear_compartment_model ``` + +## Discrete-time systems + +```@docs +@DDSmodel +``` diff --git a/docs/src/tutorials/creating_ode.md b/docs/src/tutorials/creating_ode.md index ec7d039ff..00d362846 100644 --- a/docs/src/tutorials/creating_ode.md +++ b/docs/src/tutorials/creating_ode.md @@ -54,7 +54,9 @@ assess_identifiability(ode) ## Defining using `ModelingToolkit` -Alternatively, one can use `ModelingToolkit`: encode the equations for the states as `ODESystem` and specify the outputs separately. +`StructuralIdentifiability` has an extension `ModelingToolkitExt` which allows to use `ODESystem` from `ModelingToolkit` to descibe +a model. The extension is loaded automatically once `ModelingToolkit` is loaded via `using ModelingToolkit`. +In this case, one should encode the equations for the states as `ODESystem` and specify the outputs separately. In order to do this, we first introduce all functions and scalars: ```@example 2; continued = true diff --git a/docs/src/tutorials/discrete_time.md b/docs/src/tutorials/discrete_time.md index 212f59db8..c8f0571f8 100644 --- a/docs/src/tutorials/discrete_time.md +++ b/docs/src/tutorials/discrete_time.md @@ -1,13 +1,20 @@ # Identifiability of Discrete-Time Models (Local) -Now we consider a discrete-time model in the state-space form +Now we consider a discrete-time model in the state-space form. Such a model is typically written either in terms of **shift**: + +$\begin{cases} +\mathbf{x}(t + 1) = \mathbf{f}(\mathbf{x}(t), \mathbf{p}, \mathbf{u}(t)),\\ +\mathbf{y}(t) = \mathbf{g}(\mathbf{x}(t), \mathbf{p}, \mathbf{u(t)}), +\end{cases}$ + +or in terms of **difference** $\begin{cases} \Delta\mathbf{x}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{p}, \mathbf{u}(t)),\\ \mathbf{y}(t) = \mathbf{g}(\mathbf{x}(t), \mathbf{p}, \mathbf{u(t)}), -\end{cases} \quad \text{where} \quad \Delta \mathbf{x}(t) := \mathbf{x}(t + 1) - \mathbf{x}(t)$ +\end{cases} \quad \text{where} \quad \Delta \mathbf{x}(t) := \mathbf{x}(t + 1) - \mathbf{x}(t).$ -and $\mathbf{x}(t), \mathbf{y}(t)$, and $\mathbf{u}(t)$ are time-dependent states, outputs, and inputs, respectively, +In both cases,$\mathbf{x}(t), \mathbf{y}(t)$, and $\mathbf{u}(t)$ are time-dependent states, outputs, and inputs, respectively, and $\mathbf{p}$ are scalar parameters. As in the ODE case, we will call that a parameter or a states (or a function of them) is **identifiable** if its value can be recovered from time series for inputs and outputs (in the generic case, see Definition 3 in [^1] for details). @@ -21,50 +28,53 @@ Currently, `StructuralIdentifiability.jl` allows to assess only local identifiab and below we will describe how this can be done. As a running example, we will use the following discrete version of the [SIR](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model) model: -$\begin{cases} -\Delta S(t) = S(t) - \beta S(t) I(t),\\ -\Delta I(t) = I(t) + \beta S(t) I(t) - \alpha I(t),\\ -\Delta R(t) = R(t) + \alpha I(t),\\ +$ +\begin{cases} +S(t + 1) = S(t) - \beta S(t) I(t),\\ +I(t + 1) = I(t) + \beta S(t) I(t) - \alpha I(t),\\ +R(t + 1) = R(t) + \alpha I(t),\\ +y(t) = I(t), +\end{cases} +\quad \text{or} +\quad +\begin{cases} +\Delta S(t) = -\beta S(t) I(t),\\ +\Delta I(t) = \beta S(t) I(t) - \alpha I(t),\\ +\Delta R(t) = \alpha I(t),\\ y(t) = I(t), \end{cases}$ where the observable is `I`, the number of infected people. -We start with creating a system as a `DiscreteSystem` from `ModelingToolkit`: +The native way to define such a model in `StructuralIdentifiability` is to use `@DDSmodel` macro which +uses the shift notation: -```@example discrete -using ModelingToolkit +```@example discrete_dds using StructuralIdentifiability -@parameters α β -@variables t S(t) I(t) R(t) y(t) -D = Difference(t; dt = 1.0) - -eqs = [D(S) ~ S - β * S * I, D(I) ~ I + β * S * I - α * I, D(R) ~ R + α * I] -@named sir = DiscreteSystem(eqs) +dds = @DDSmodel( + S(t + 1) = S(t) - β * S(t) * I(t), + I(t + 1) = I(t) + β * S(t) * I(t) - α * I(t), + R(t + 1) = R(t) + α * I(t), + y(t) = I(t) +) ``` -Once the model is defined, we can assess identifiability by providing the formula for the observable: +Then local identifiability can be assessed using `assess_local_identifiability` function: -```@example discrete -assess_local_identifiability(sir; measured_quantities = [y ~ I]) +```@example discrete_dds +assess_local_identifiability(dds) ``` For each parameter or state, the value in the returned dictionary is `true` (`1`) if the parameter is locally identifiable and `false` (`0`) otherwise. We see that `R(t)` is not identifiable, which makes sense: it does not affect the dynamics of the observable in any way. -In principle, it is not required to give a name to the observable, so one can write this shorter - -```@example discrete -assess_local_identifiability(sir; measured_quantities = [I]) -``` - The `assess_local_identifiability` function has three important keyword arguments: - `funcs_to_check` is a list of functions for which one want to assess identifiability, for example, the following code will check if `β * S` is locally identifiable. -```@example discrete -assess_local_identifiability(sir; measured_quantities = [I], funcs_to_check = [β * S]) +```@example discrete_dds +assess_local_identifiability(dds; funcs_to_check = [β * S]) ``` - `prob_threshold` is the probability of correctness (default value `0.99`, i.e., 99%). The underlying algorithm is a Monte-Carlo algorithm, so in @@ -77,6 +87,39 @@ assess_local_identifiability(sir; measured_quantities = [I], funcs_to_check = [ As other main functions in the package, `assess_local_identifiability` accepts an optional parameter `loglevel` (default: `Logging.Info`) to adjust the verbosity of logging. +If one loads `ModelingToolkit` (and thus the `ModelingToolkitExt` extension), one can use `DiscreteSystem` from `ModelingToolkit` to +describe the input model (now in terms of difference!): + +```@example discrete_mtk +using ModelingToolkit +using StructuralIdentifiability + +@parameters α β +@variables t S(t) I(t) R(t) y(t) +D = Difference(t; dt = 1.0) + +eqs = [D(S) ~ S - β * S * I, D(I) ~ I + β * S * I - α * I, D(R) ~ R + α * I] +@named sir = DiscreteSystem(eqs) +``` + +Then the same computation can be carried out with the models defined this way: + +```@example discrete_mtk +assess_local_identifiability(sir; measured_quantities = [y ~ I]) +``` + +In principle, it is not required to give a name to the observable, so one can write this shorter + +```@example discrete_mtk +assess_local_identifiability(sir; measured_quantities = [I]) +``` + +The same example but with specified functions to check + +```@example discrete_mtk +assess_local_identifiability(sir; measured_quantities = [I], funcs_to_check = [β * S]) +``` + The implementation is based on a version of the observability rank criterion and will be described in a forthcoming paper. [^1]: > S. Nõmm, C. Moog, [*Identifiability of discrete-time nonlinear systems*](https://doi.org/10.1016/S1474-6670(17)31245-4), IFAC Proceedings Volumes, 2004. diff --git a/src/discrete.jl b/src/discrete.jl index 4b96e6123..a25004994 100644 --- a/src/discrete.jl +++ b/src/discrete.jl @@ -391,6 +391,9 @@ function assess_local_identifiability( restart_logging(loglevel = loglevel) reset_timings() with_logger(_si_logger[]) do + if isempty(funcs_to_check) + funcs_to_check = vcat(parameters(dds), x_vars(dds)) + end return _assess_local_identifiability_discrete_aux( dds, funcs_to_check,