Merge pull request #266 from SciML/mtk_extension

Render ModelingToolkit as an extension

.github/workflows/CI.yml

@@ -13,6 +13,7 @@ jobs:
       matrix:
         group:
           - Core
+          - ModelingToolkitExt
         version:
           - '1'
           - '1.6'

Project.toml

@@ -14,45 +14,55 @@ IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 ParamPunPam = "3e851597-e36f-45a9-af0a-b7781937992f"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+
+[weakdeps]
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+
+[extensions]
+ModelingToolkitExt = ["ModelingToolkit", "SymbolicUtils", "Symbolics"]
 
 [compat]
-AbstractAlgebra = "0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34, 0.35"
+AbstractAlgebra = "0.34.5, 0.35"
 BenchmarkTools = "1"
 Combinatorics = "1"
 DataStructures = "0.18"
 Dates = "1.6, 1.7"
-Groebner = "0.4, 0.5, 0.6"
+Groebner = "0.6.3"
 IterTools = "1"
 LinearAlgebra = "1.6, 1.7"
 Logging = "1.6, 1.7"
 MacroTools = "0.5"
-ModelingToolkit = "8.51"
-Nemo = "0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39"
-ParamPunPam = "0.2"
-PrecompileTools = "1.1, 1.2"
+ModelingToolkit = "8.74"
+Nemo = "0.38.3, 0.39"
+ParamPunPam = "0.3.1"
+PrecompileTools = "1.2"
 Primes = "0.5"
 Random = "1.6, 1.7"
-SpecialFunctions = "1, 2"
-SymbolicUtils = "1"
-Symbolics = "5.5"
+SpecialFunctions = "2"
+SymbolicUtils = "1.4, 1.5"
+Symbolics = "5.16"
 TestSetExtensions = "2"
 TimerOutputs = "0.5"
 julia = "1.6, 1.7"
 
 [extras]
 CPUSummary = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
 
 [targets]
-test = ["CPUSummary", "Test", "TestSetExtensions"]
+test = ["CPUSummary", "Pkg", "Test", "TestSetExtensions"]

benchmarking/IdentifiableFunctions/param.jl

@@ -31,7 +31,7 @@ function check_constructive_field_membership(generators::AbstractVector, to_be_r
     $(join(map(x -> string(x[1]) * " -> " * string(x[2]),  zip(fracs_gen, tag_strings)), "\t\n"))
     """
     var_strings = vcat(sat_string, map(string, gens(ring)), tag_strings)
-    ring_tag, xs_tag = PolynomialRing(K, var_strings, ordering = Nemo.ordering(ring))
+    ring_tag, xs_tag = polynomial_ring(K, var_strings, ordering = Nemo.ordering(ring))
     orig_vars = xs_tag[2:(nvars(ring) + 1)]
     tag_vars = xs_tag[(nvars(ring) + 2):end]
     sat_var = xs_tag[1]
@@ -121,12 +121,12 @@ rem_tags, tag_to_gen = check_constructive_field_membership(fracs_generators, to_
 #=
 ┌ Info: 
 │   rem_tags =
-│    3-element Vector{AbstractAlgebra.Generic.Frac{fmpq_mpoly}}:
+│    3-element Vector{AbstractAlgebra.Generic.Frac{QQMPolyRingElem}}:
 │     T1^2
 │     -5*T1 + T2
 │     T2//T1^10
 │   tag_to_gen =
-│    Dict{fmpq_mpoly, AbstractAlgebra.Generic.Frac{fmpq_mpoly}} with 2 entries:
+│    Dict{QQMPolyRingElem, AbstractAlgebra.Generic.Frac{QQMPolyRingElem}} with 2 entries:
 │      T1 => a^2
 └      T2 => (a + b)//b
 =#
@@ -205,10 +205,10 @@ new_vector_field, new_outputs, new_vars =
 │    Dict{Any, Any} with 1 entry:
 │      T1 => 2*T1 + T2
 │   new_outputs =
-│    Dict{fmpq_mpoly, AbstractAlgebra.Generic.Frac{fmpq_mpoly}} with 1 entry:
+│    Dict{QQMPolyRingElem, AbstractAlgebra.Generic.Frac{QQMPolyRingElem}} with 1 entry:
 │      y => T1
 │   new_vars =
-│    Dict{fmpq_mpoly, AbstractAlgebra.Generic.Frac{fmpq_mpoly}} with 2 entries:
+│    Dict{QQMPolyRingElem, AbstractAlgebra.Generic.Frac{QQMPolyRingElem}} with 2 entries:
 │      T2 => a + b
 └      T1 => x2 + x1
 =#
@@ -230,7 +230,7 @@ id_funcs = StructuralIdentifiability.find_identifiable_functions(
 #=
 ┌ Info: 
 │   id_funcs =
-│    2-element Vector{AbstractAlgebra.Generic.Frac{fmpq_mpoly}}:
+│    2-element Vector{AbstractAlgebra.Generic.Frac{QQMPolyRingElem}}:
 │     x2*x1
 └     a + b
 =#
@@ -246,10 +246,10 @@ new_vector_field, new_outputs, new_vars = reparametrize_with_respect_to(ode, new
 │    Dict{Any, Any} with 1 entry:
 │      T1 => T1*T2
 │   new_outputs =
-│    Dict{fmpq_mpoly, AbstractAlgebra.Generic.Frac{fmpq_mpoly}} with 1 entry:
+│    Dict{QQMPolyRingElem, AbstractAlgebra.Generic.Frac{QQMPolyRingElem}} with 1 entry:
 │      y => T1
 │   new_vars =
-│    Dict{fmpq_mpoly, AbstractAlgebra.Generic.Frac{fmpq_mpoly}} with 2 entries:
+│    Dict{QQMPolyRingElem, AbstractAlgebra.Generic.Frac{QQMPolyRingElem}} with 2 entries:
 │      T2 => a + b
 └      T1 => x2*x1
 =#
@@ -271,7 +271,7 @@ id_funcs = StructuralIdentifiability.find_identifiable_functions(
 #=
 ┌ Info: 
 │   id_funcs =
-│    5-element Vector{AbstractAlgebra.Generic.Frac{fmpq_mpoly}}:
+│    5-element Vector{AbstractAlgebra.Generic.Frac{QQMPolyRingElem}}:
 │     x2*x1
 │     a*b
 │     x2 + x1
@@ -292,10 +292,10 @@ new_vector_field, new_iutputs, new_vars = reparametrize_with_respect_to(ode, new
 │      T2 => T2*T4
 │      T3 => -1//2*T1*T4^2 + 2*T1*T5 + 1//2*T3*T4
 │   new_outputs =
-│    Dict{fmpq_mpoly, AbstractAlgebra.Generic.Frac{fmpq_mpoly}} with 1 entry:
+│    Dict{QQMPolyRingElem, AbstractAlgebra.Generic.Frac{QQMPolyRingElem}} with 1 entry:
 │      y => T1
 │   new_vars =
-│    Dict{fmpq_mpoly, AbstractAlgebra.Generic.Frac{fmpq_mpoly}} with 5 entries:
+│    Dict{QQMPolyRingElem, AbstractAlgebra.Generic.Frac{QQMPolyRingElem}} with 5 entries:
 │      T1 => x2 + x1
 │      T2 => x2*x1
 │      T3 => a*x2 - a*x1 - b*x2 + b*x1

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"

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"

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
+```

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

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.