diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 13e7d89b9..997b5426e 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -13,6 +13,7 @@ jobs: matrix: group: - Core + - ModelingToolkitExt version: - '1' - '1.6' diff --git a/Project.toml b/Project.toml index 32bc21b06..8513a0819 100644 --- a/Project.toml +++ b/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"] diff --git a/benchmarking/IdentifiableFunctions/param.jl b/benchmarking/IdentifiableFunctions/param.jl index 44b2f54fa..1a4fd93f5 100644 --- a/benchmarking/IdentifiableFunctions/param.jl +++ b/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 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/ext/ModelingToolkitExt.jl b/ext/ModelingToolkitExt.jl new file mode 100644 index 000000000..4ca8de85f --- /dev/null +++ b/ext/ModelingToolkitExt.jl @@ -0,0 +1,595 @@ +module ModelingToolkitExt + +using DataStructures +using Logging +using Nemo +using Random +using StructuralIdentifiability +using StructuralIdentifiability: str_to_var, parent_ring_change, eval_at_dict +using StructuralIdentifiability: restart_logging, _si_logger, reset_timings, _to +using TimerOutputs + +if isdefined(Base, :get_extension) + using ModelingToolkit +else + using ..ModelingToolkit +end + +export assess_local_identifiability, assess_identifiability, find_identifiable_functions + +# ------------------------------------------------------------------------------ + +function eval_at_nemo(e::Num, vals::Dict) + e = Symbolics.value(e) + return eval_at_nemo(e, vals) +end + +function eval_at_nemo(e::SymbolicUtils.BasicSymbolic, vals::Dict) + if Symbolics.istree(e) + # checking if it is a function of the form x(t), a bit dirty + if length(Symbolics.arguments(e)) == 1 && "$(first(Symbolics.arguments(e)))" == "t" + return vals[e] + end + # checking if this is a vector entry like x(t)[1] + if Symbolics.operation(e) == getindex + return vals[e] + end + # otherwise, this is a term + args = map(a -> eval_at_nemo(a, vals), Symbolics.arguments(e)) + if Symbolics.operation(e) in (+, -, *) + return Symbolics.operation(e)(args...) + elseif isequal(Symbolics.operation(e), /) + return //(args...) + elseif isequal(Symbolics.operation(e), ^) + if args[2] >= 0 + return args[1]^args[2] + end + return 1 // args[1]^(-args[2]) + end + throw(Base.ArgumentError("Function $(Symbolics.operation(e)) is not supported")) + elseif e isa Symbolics.Symbolic + return get(vals, e, e) + end +end + +function eval_at_nemo(e::Union{Integer, Rational}, vals::Dict) + return e +end + +function eval_at_nemo(e::Union{Float16, Float32, Float64}, vals::Dict) + if isequal(e % 1, 0) + out = Int(e) + else + out = rationalize(e) + end + @warn "Floating point value $e will be converted to $(out)." + return out +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 + +""" + function mtk_to_si(de::ModelingToolkit.AbstractTimeDependentSystem, measured_quantities::Array{ModelingToolkit.Equation}) + function mtk_to_si(de::ModelingToolkit.AbstractTimeDependentSystem, measured_quantities::Array{SymbolicUtils.BasicSymbolic}) + +Input: +- `de` - ModelingToolkit.AbstractTimeDependentSystem, a system for identifiability query +- `measured_quantities` - array of output functions (as equations of just functions) + +Output: +- `ODE` object containing required data for identifiability assessment +- `conversion` dictionary from the symbols in the input MTK model to the variable + involved in the produced `ODE` object +""" +function mtk_to_si( + de::ModelingToolkit.AbstractTimeDependentSystem, + measured_quantities::Array{ModelingToolkit.Equation}, +) + return __mtk_to_si( + de, + [(replace(string(e.lhs), "(t)" => ""), e.rhs) for e in measured_quantities], + ) +end + +function mtk_to_si( + de::ModelingToolkit.AbstractTimeDependentSystem, + measured_quantities::Array{<:Symbolics.Num}, +) + return __mtk_to_si( + de, + [("y$i", Symbolics.value(e)) for (i, e) in enumerate(measured_quantities)], + ) +end + +function mtk_to_si( + de::ModelingToolkit.AbstractTimeDependentSystem, + measured_quantities::Array{<:SymbolicUtils.BasicSymbolic}, +) + return __mtk_to_si(de, [("y$i", e) for (i, e) in enumerate(measured_quantities)]) +end + +#------------------------------------------------------------------------------ +# old name kept for compatibility purposes + +function preprocess_ode( + de::ModelingToolkit.AbstractTimeDependentSystem, + measured_quantities::Array{ModelingToolkit.Equation}, +) + @warn "Function `preprocess_ode` has been renamed to `mtk_to_si`. The old name can be still used but will disappear in the future releases." + return mtk_to_si(de, measured_quantities) +end + +function preprocess_ode( + de::ModelingToolkit.AbstractTimeDependentSystem, + measured_quantities::Array{<:Symbolics.Num}, +) + @warn "Function `preprocess_ode` has been renamed to `mtk_to_si`. The old name can be still used but will disappear in the future releases." + return mtk_to_si(de, measured_quantities) +end + +function preprocess_ode( + de::ModelingToolkit.AbstractTimeDependentSystem, + measured_quantities::Array{<:SymbolicUtils.BasicSymbolic}, +) + @warn "Function `preprocess_ode` has been renamed to `mtk_to_si`. The old name can be still used but will disappear in the future releases." + return mtk_to_si(de, measured_quantities) +end + +#------------------------------------------------------------------------------ +""" + function __mtk_to_si(de::ModelingToolkit.AbstractTimeDependentSystem, measured_quantities::Array{Tuple{String, SymbolicUtils.BasicSymbolic}}) + +Input: +- `de` - ModelingToolkit.AbstractTimeDependentSystem, a system for identifiability query +- `measured_quantities` - array of input function in the form (name, expression) + +Output: +- `ODE` object containing required data for identifiability assessment +- `conversion` dictionary from the symbols in the input MTK model to the variable + involved in the produced `ODE` object +""" +function __mtk_to_si( + de::ModelingToolkit.AbstractTimeDependentSystem, + measured_quantities::Array{<:Tuple{String, <:SymbolicUtils.BasicSymbolic}}, +) + polytype = StructuralIdentifiability.Nemo.QQMPolyRingElem + fractype = StructuralIdentifiability.Nemo.Generic.Frac{polytype} + diff_eqs = + filter(eq -> !(ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(de)) + # performing full structural simplification + if length(observed(de)) > 0 + rules = Dict(s.lhs => s.rhs for s in observed(de)) + while any([ + length(intersect(get_variables(r), keys(rules))) > 0 for r in values(rules) + ]) + rules = Dict(k => SymbolicUtils.substitute(v, rules) for (k, v) in rules) + end + diff_eqs = [SymbolicUtils.substitute(eq, rules) for eq in diff_eqs] + end + + y_functions = [each[2] for each in measured_quantities] + inputs = filter(v -> ModelingToolkit.isinput(v), ModelingToolkit.states(de)) + state_vars = filter( + s -> !(ModelingToolkit.isinput(s) || ModelingToolkit.isoutput(s)), + ModelingToolkit.states(de), + ) + params = ModelingToolkit.parameters(de) + t = ModelingToolkit.arguments(diff_eqs[1].lhs)[1] + params_from_measured_quantities = union( + [filter(s -> !istree(s), get_variables(y[2])) for y in measured_quantities]..., + ) + params = union(params, params_from_measured_quantities) + + input_symbols = vcat(state_vars, inputs, params) + generators = vcat(string.(input_symbols), [e[1] for e in measured_quantities]) + generators = map(g -> replace(g, "(t)" => ""), generators) + R, gens_ = Nemo.polynomial_ring(Nemo.QQ, generators) + y_vars = Vector{polytype}([str_to_var(e[1], R) for e in measured_quantities]) + symb2gens = Dict(input_symbols .=> gens_[1:length(input_symbols)]) + + x_vars = Vector{polytype}() + + state_eqn_dict = Dict{polytype, Union{polytype, fractype}}() + out_eqn_dict = Dict{polytype, Union{polytype, fractype}}() + + for i in 1:length(diff_eqs) + x = substitute(state_vars[i], symb2gens) + push!(x_vars, x) + if !(diff_eqs[i].rhs isa Number) + state_eqn_dict[x] = eval_at_nemo(diff_eqs[i].rhs, symb2gens) + else + state_eqn_dict[x] = R(diff_eqs[i].rhs) + end + end + for i in 1:length(measured_quantities) + out_eqn_dict[y_vars[i]] = eval_at_nemo(measured_quantities[i][2], symb2gens) + end + + inputs_ = [substitute(each, symb2gens) for each in inputs] + if isequal(length(inputs_), 0) + inputs_ = Vector{polytype}() + end + return ( + StructuralIdentifiability.ODE{polytype}( + x_vars, + y_vars, + state_eqn_dict, + out_eqn_dict, + inputs_, + ), + symb2gens, + ) +end +# ----------------------------------------------------------------------------- +""" + function assess_local_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=Array{}[], prob_threshold::Float64=0.99, type=:SE, loglevel=Logging.Info) + +Input: +- `ode` - the ODESystem object from ModelingToolkit +- `measured_quantities` - the measurable outputs of the model +- `funcs_to_check` - functions of parameters for which to check identifiability +- `prob_threshold` - probability of correctness +- `type` - identifiability type (`:SE` for single-experiment, `:ME` for multi-experiment) +- `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) + +Output: +- for `type=:SE`, the result is an (ordered) dictionary from each parameter to boolean; +- for `type=:ME`, the result is a tuple with the dictionary as in `:SE` case and array of number of experiments. + +The function determines local identifiability of parameters in `funcs_to_check` or all possible parameters if `funcs_to_check` is empty + +The result is correct with probability at least `prob_threshold`. + +`type` can be either `:SE` (single-experiment identifiability) or `:ME` (multi-experiment identifiability). +The return value is a tuple consisting of the array of bools and the number of experiments to be performed. +""" +function StructuralIdentifiability.assess_local_identifiability( + ode::ModelingToolkit.ODESystem; + measured_quantities = Array{ModelingToolkit.Equation}[], + funcs_to_check = Array{}[], + prob_threshold::Float64 = 0.99, + type = :SE, + loglevel = Logging.Info, +) + restart_logging(loglevel = loglevel) + with_logger(_si_logger[]) do + return _assess_local_identifiability( + ode, + measured_quantities = measured_quantities, + funcs_to_check = funcs_to_check, + prob_threshold = prob_threshold, + type = type, + ) + end +end + +@timeit _to function _assess_local_identifiability( + ode::ModelingToolkit.ODESystem; + measured_quantities = Array{ModelingToolkit.Equation}[], + funcs_to_check = Array{}[], + prob_threshold::Float64 = 0.99, + type = :SE, +) + if length(measured_quantities) == 0 + 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." + measured_quantities = 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 + if length(funcs_to_check) == 0 + funcs_to_check = vcat( + [e for e in ModelingToolkit.states(ode) if !ModelingToolkit.isoutput(e)], + ModelingToolkit.parameters(ode), + ) + end + ode, conversion = mtk_to_si(ode, measured_quantities) + funcs_to_check_ = [eval_at_nemo(x, conversion) for x in funcs_to_check] + + if isequal(type, :SE) + result = StructuralIdentifiability._assess_local_identifiability( + ode, + funcs_to_check = funcs_to_check_, + prob_threshold = prob_threshold, + type = type, + ) + nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) + out_dict = + OrderedDict(nemo2mtk[param] => result[param] for param in funcs_to_check_) + return out_dict + elseif isequal(type, :ME) + result, bd = StructuralIdentifiability._assess_local_identifiability( + ode, + funcs_to_check = funcs_to_check_, + prob_threshold = prob_threshold, + type = type, + ) + nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) + out_dict = + OrderedDict(nemo2mtk[param] => result[param] for param in funcs_to_check_) + return (out_dict, bd) + end +end + +# ------------------------------------------------------------------------------ + +""" + assess_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=[], prob_threshold = 0.99, loglevel=Logging.Info) + +Input: +- `ode` - the ModelingToolkit.ODESystem object that defines the model +- `measured_quantities` - the output functions of the model +- `funcs_to_check` - functions of parameters for which to check the identifiability +- `prob_threshold` - probability of correctness. +- `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) + +Assesses identifiability (both local and global) of a given ODE model (parameters detected automatically). The result is guaranteed to be correct with the probability +at least `prob_threshold`. +""" +function StructuralIdentifiability.assess_identifiability( + ode::ModelingToolkit.ODESystem; + measured_quantities = Array{ModelingToolkit.Equation}[], + funcs_to_check = [], + prob_threshold = 0.99, + loglevel = Logging.Info, +) + restart_logging(loglevel = loglevel) + with_logger(_si_logger[]) do + return _assess_identifiability( + ode, + measured_quantities = measured_quantities, + funcs_to_check = funcs_to_check, + prob_threshold = prob_threshold, + ) + end +end + +function _assess_identifiability( + ode::ModelingToolkit.ODESystem; + measured_quantities = Array{ModelingToolkit.Equation}[], + funcs_to_check = [], + prob_threshold = 0.99, +) + if isempty(measured_quantities) + measured_quantities = get_measured_quantities(ode) + end + + ode, conversion = mtk_to_si(ode, measured_quantities) + conversion_back = Dict(v => k for (k, v) in conversion) + if isempty(funcs_to_check) + funcs_to_check = [conversion_back[x] for x in [ode.x_vars..., ode.parameters...]] + end + funcs_to_check_ = [eval_at_nemo(each, conversion) for each in funcs_to_check] + + result = StructuralIdentifiability._assess_identifiability( + ode, + funcs_to_check = funcs_to_check_, + prob_threshold = prob_threshold, + ) + nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) + out_dict = OrderedDict(nemo2mtk[param] => result[param] for param in funcs_to_check_) + return out_dict +end + +# ------------------------------------------------------------------------------ + +""" + function assess_local_identifiability( + dds::ModelingToolkit.DiscreteSystem; + measured_quantities=Array{ModelingToolkit.Equation}[], + funcs_to_check=Array{}[], + known_ic=Array{}[], + prob_threshold::Float64=0.99) + +Input: +- `dds` - the DiscreteSystem object from ModelingToolkit (with **difference** operator in the left-hand side) +- `measured_quantities` - the measurable outputs of the model +- `funcs_to_check` - functions of parameters for which to check identifiability (all parameters and states if not specified) +- `known_ic` - functions (of states and parameter) whose initial conditions are assumed to be known +- `prob_threshold` - probability of correctness + +Output: +- the result is an (ordered) dictionary from each function to to boolean; + +The result is correct with probability at least `prob_threshold`. +""" +function StructuralIdentifiability.assess_local_identifiability( + dds::ModelingToolkit.DiscreteSystem; + measured_quantities = Array{ModelingToolkit.Equation}[], + funcs_to_check = Array{}[], + known_ic = Array{}[], + prob_threshold::Float64 = 0.99, + loglevel = Logging.Info, +) + restart_logging(loglevel = loglevel) + with_logger(_si_logger[]) do + return _assess_local_identifiability( + dds, + measured_quantities = measured_quantities, + funcs_to_check = funcs_to_check, + known_ic = known_ic, + prob_threshold = prob_threshold, + ) + end +end + +function _assess_local_identifiability( + dds::ModelingToolkit.DiscreteSystem; + measured_quantities = Array{ModelingToolkit.Equation}[], + funcs_to_check = Array{}[], + known_ic = Array{}[], + prob_threshold::Float64 = 0.99, +) + if length(measured_quantities) == 0 + if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(dds)) + @info "Measured quantities are not provided, trying to find the outputs in input dynamical system." + measured_quantities = filter( + eq -> (ModelingToolkit.isoutput(eq.lhs)), + ModelingToolkit.equations(dds), + ) + else + throw( + error( + "Measured quantities (output functions) were not provided and no outputs were found.", + ), + ) + end + end + + # Converting the finite difference operator in the right-hand side to + # the corresponding shift operator + eqs = filter(eq -> !(ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(dds)) + deltas = [Symbolics.operation(e.lhs).dt for e in eqs] + @assert length(Set(deltas)) == 1 + eqs_shift = [e.lhs ~ e.rhs + first(Symbolics.arguments(e.lhs)) for e in eqs] + dds_shift = DiscreteSystem(eqs_shift, name = gensym()) + @debug "System transformed from difference to shift: $dds_shift" + + dds_aux_ode, conversion = mtk_to_si(dds_shift, measured_quantities) + dds_aux = StructuralIdentifiability.DDS{QQMPolyRingElem}(dds_aux_ode) + if length(funcs_to_check) == 0 + params = parameters(dds) + params_from_measured_quantities = union( + [filter(s -> !istree(s), get_variables(y)) for y in measured_quantities]..., + ) + funcs_to_check = vcat( + [ + x for x in states(dds) if + conversion[x] in StructuralIdentifiability.x_vars(dds_aux) + ], + union(params, params_from_measured_quantities), + ) + end + funcs_to_check_ = [eval_at_nemo(x, conversion) for x in funcs_to_check] + known_ic_ = [eval_at_nemo(x, conversion) for x in known_ic] + + result = StructuralIdentifiability._assess_local_identifiability_discrete_aux( + dds_aux, + funcs_to_check_, + known_ic_, + prob_threshold, + ) + nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) + out_dict = OrderedDict(nemo2mtk[param] => result[param] for param in funcs_to_check_) + if length(known_ic) > 0 + @warn "Since known initial conditions were provided, identifiability of states (e.g., `x(t)`) is at t = 0 only !" + t = SymbolicUtils.Sym{Real}(:t) + out_dict = OrderedDict(substitute(k, Dict(t => 0)) => v for (k, v) in out_dict) + end + return out_dict +end +# ------------------------------------------------------------------------------ + +""" + find_identifiable_functions(ode::ModelingToolkit.ODESystem; measured_quantities=[], options...) + +Finds all functions of parameters/states that are identifiable in the given ODE +system. + +## Options + +This functions takes the following optional arguments: +- `measured_quantities` - the output functions of the model. +- `loglevel` - the verbosity of the logging + (can be Logging.Error, Logging.Warn, Logging.Info, Logging.Debug) + +## Example + +```jldoctest +using StructuralIdentifiability +using ModelingToolkit + +@parameters a01 a21 a12 +@variables t x0(t) x1(t) y1(t) [output = true] +D = Differential(t) + +eqs = [ + D(x0) ~ -(a01 + a21) * x0 + a12 * x1, + D(x1) ~ a21 * x0 - a12 * x1, y1 ~ x0 +] +de = ODESystem(eqs, t, name = :Test) + +find_identifiable_functions(de, measured_quantities = [y1 ~ x0]) + +# prints +2-element Vector{Num}: + a01*a12 + a01 + a12 + a21 +``` +""" +function StructuralIdentifiability.find_identifiable_functions( + ode::ModelingToolkit.ODESystem; + measured_quantities = Array{ModelingToolkit.Equation}[], + prob_threshold::Float64 = 0.99, + seed = 42, + with_states = false, + simplify = :standard, + rational_interpolator = :VanDerHoevenLecerf, + loglevel = Logging.Info, +) + restart_logging(loglevel = loglevel) + reset_timings() + with_logger(_si_logger[]) do + return _find_identifiable_functions( + ode, + measured_quantities = measured_quantities, + prob_threshold = prob_threshold, + seed = seed, + with_states = with_states, + simplify = simplify, + rational_interpolator = rational_interpolator, + ) + end +end + +function _find_identifiable_functions( + ode::ModelingToolkit.ODESystem; + measured_quantities = Array{ModelingToolkit.Equation}[], + prob_threshold::Float64 = 0.99, + seed = 42, + with_states = false, + simplify = :standard, + rational_interpolator = :VanDerHoevenLecerf, +) + Random.seed!(seed) + if isempty(measured_quantities) + measured_quantities = get_measured_quantities(ode) + end + ode, conversion = mtk_to_si(ode, measured_quantities) + result = StructuralIdentifiability._find_identifiable_functions( + ode, + simplify = simplify, + prob_threshold = prob_threshold, + seed = seed, + with_states = with_states, + rational_interpolator = rational_interpolator, + ) + result = [parent_ring_change(f, ode.poly_ring) for f in result] + nemo2mtk = Dict(v => Num(k) for (k, v) in conversion) + out_funcs = [eval_at_dict(func, nemo2mtk) for func in result] + return out_funcs +end + +end diff --git a/src/ODE.jl b/src/ODE.jl index 625c0d21e..addf9f933 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -20,7 +20,7 @@ struct ODE{P} # P is the type of polynomials in the rhs of the ODE system x_eqs::Dict{P, <:Union{P, Generic.Frac{P}}}, y_eqs::Dict{P, <:Union{P, Generic.Frac{P}}}, inputs::Array{P, 1}, - ) where {P <: MPolyElem{<:FieldElem}} + ) where {P <: MPolyRingElem{<:FieldElem}} # Initialize ODE # x_eqs is a dictionary x_i => f_i(x, u, params) # y_eqs is a dictionary y_i => g_i(x, u, params) @@ -40,7 +40,7 @@ struct ODE{P} # P is the type of polynomials in the rhs of the ODE system x_eqs::Dict{P, <:Union{P, Generic.Frac{P}}}, y_eqs::Dict{P, <:Union{P, Generic.Frac{P}}}, inputs::Array{P, 1}, - ) where {P <: MPolyElem{<:FieldElem}} + ) where {P <: MPolyRingElem{<:FieldElem}} x_vars = collect(keys(x_eqs)) y_vars = collect(keys(y_eqs)) return ODE{P}(x_vars, y_vars, x_eqs, y_eqs, inputs) @@ -53,10 +53,13 @@ end #------------------------------------------------------------------------------ -function add_outputs(ode::ODE{P}, extra_y::Dict{String, <:RingElem}) where {P <: MPolyElem} +function add_outputs( + ode::ODE{P}, + extra_y::Dict{String, <:RingElem}, +) where {P <: MPolyRingElem} new_var_names = vcat(collect(map(var_to_str, gens(ode.poly_ring))), collect(keys(extra_y))) - new_ring, new_vars = Nemo.PolynomialRing(base_ring(ode.poly_ring), new_var_names) + new_ring, new_vars = Nemo.polynomial_ring(base_ring(ode.poly_ring), new_var_names) new_x = Array{P, 1}([parent_ring_change(x, new_ring) for x in ode.x_vars]) new_x_eqs = Dict{P, Union{P, Generic.Frac{P}}}( @@ -94,10 +97,10 @@ Output: function set_parameter_values( ode::ODE{P}, param_values::Dict{P, T}, -) where {T <: FieldElem, P <: MPolyElem{T}} +) where {T <: FieldElem, P <: MPolyRingElem{T}} new_vars = map(var_to_str, [v for v in gens(ode.poly_ring) if !(v in keys(param_values))]) - small_ring, small_vars = Nemo.PolynomialRing(base_ring(ode.poly_ring), new_vars) + small_ring, small_vars = Nemo.polynomial_ring(base_ring(ode.poly_ring), new_vars) eval_dict = Dict(str_to_var(v, ode.poly_ring) => str_to_var(v, small_ring) for v in new_vars) merge!(eval_dict, Dict(p => small_ring(val) for (p, val) in param_values)) @@ -130,8 +133,8 @@ end function set_parameter_values( ode::ODE{P}, param_values::Dict{P, <:Number}, -) where {P <: MPolyElem} - new_values = Dict{P, fmpq}(x => _to_rational(v) for (x, v) in param_values) +) where {P <: MPolyRingElem} + new_values = Dict{P, QQFieldElem}(x => _to_rational(v) for (x, v) in param_values) return set_parameter_values(ode, new_values) end @@ -156,11 +159,11 @@ function power_series_solution( initial_conditions::Dict{P, T}, input_values::Dict{P, Array{T, 1}}, prec::Int, -) where {T <: FieldElem, P <: MPolyElem{T}} +) where {T <: FieldElem, P <: MPolyRingElem{T}} new_varnames = map(var_to_str, vcat(ode.x_vars, ode.u_vars)) append!(new_varnames, map(v -> var_to_str(v) * "_dot", ode.x_vars)) - new_ring, new_vars = Nemo.PolynomialRing(base_ring(ode.poly_ring), new_varnames) + new_ring, new_vars = Nemo.polynomial_ring(base_ring(ode.poly_ring), new_varnames) equations = Array{P, 1}() evaluation = Dict(k => new_ring(v) for (k, v) in param_values) for v in vcat(ode.x_vars, ode.u_vars) @@ -196,7 +199,7 @@ function power_series_solution( initial_conditions::Dict{P, Int}, input_values::Dict{P, Array{Int, 1}}, prec::Int, -) where {P <: MPolyElem{<:FieldElem}} +) where {P <: MPolyRingElem{<:FieldElem}} bring = base_ring(ode.poly_ring) return power_series_solution( ode, @@ -213,16 +216,16 @@ end Reduces a polynomial/rational function over Q modulo p """ -function _reduce_mod_p(poly::fmpq_mpoly, p::Int) +function _reduce_mod_p(poly::QQMPolyRingElem, p::Int) den = denominator(poly) num = change_base_ring(Nemo.ZZ, den * poly) - if Nemo.GF(p)(den) == 0 + if Nemo.Native.GF(p)(den) == 0 throw(Base.ArgumentError("Prime $p divides the denominator of $poly")) end - return change_base_ring(Nemo.GF(p), num) * (1 // Nemo.GF(p)(den)) + return change_base_ring(Nemo.Native.GF(p), num) * (1 // Nemo.Native.GF(p)(den)) end -function _reduce_mod_p(rat::Generic.Frac{fmpq_mpoly}, p::Int) +function _reduce_mod_p(rat::Generic.Frac{QQMPolyRingElem}, p::Int) num, den = map(poly -> _reduce_mod_p(poly, p), [numerator(rat), denominator(rat)]) if den == 0 throw(Base.ArgumentError("Prime $p divides the denominator of $rat")) @@ -238,9 +241,9 @@ end Input: ode is an ODE over QQ, p is a prime number Output: the reduction mod p, throws an exception if p divides one of the denominators """ -function reduce_ode_mod_p(ode::ODE{<:MPolyElem{Nemo.fmpq}}, p::Int) +function reduce_ode_mod_p(ode::ODE{<:MPolyRingElem{Nemo.QQFieldElem}}, p::Int) new_ring, new_vars = - Nemo.PolynomialRing(Nemo.GF(p), map(var_to_str, gens(ode.poly_ring))) + Nemo.polynomial_ring(Nemo.Native.GF(p), map(var_to_str, gens(ode.poly_ring))) new_type = typeof(new_vars[1]) new_inputs = map(u -> switch_ring(u, new_ring), ode.u_vars) new_x = map(x -> switch_ring(x, new_ring), ode.x_vars) @@ -258,245 +261,6 @@ end #------------------------------------------------------------------------------ -function _extract_aux!(funcs, all_symb, eq, ders_ok = false) - aux_symb = Set([:(+), :(-), :(=), :(*), :(^), :t, :(/), :(//)]) - MacroTools.postwalk( - x -> begin - if @capture(x, f_'(t)) - if !ders_ok - throw( - Base.ArgumentError( - "Derivative are not allowed in the right-hand side", - ), - ) - end - push!(all_symb, f) - elseif @capture(x, f_(t)) - push!(funcs, f) - elseif (x isa Symbol) && !(x in aux_symb) - push!(all_symb, x) - end - return x - end, - eq, - ) -end - -""" - For an expression of the form f'(t) or f(t) returns (f, true) and (f, false), resp -""" -function _get_var(expr) - if @capture(expr, f_'(t)) - return (f, true) - end - if @capture(expr, f_(t)) - return (f, false) - end - error("cannot extract the single function name from $expr") -end - -function macrohelper_extract_vars(equations::Array{Expr, 1}) - funcs, all_symb = Set(), Set() - x_vars, y_vars = Vector(), Vector() - aux_symb = Set([:(+), :(-), :(=), :(*), :(^), :t, :(/), :(//)]) - for eq in equations - if eq.head != :(=) - _extract_aux!(funcs, all_symb, eq) - else - lhs, rhs = eq.args[1:2] - _extract_aux!(funcs, all_symb, lhs, true) - _extract_aux!(funcs, all_symb, rhs) - (v, is_state) = _get_var(lhs) - if is_state - push!(x_vars, v) - else - push!(y_vars, v) - end - end - end - u_vars = setdiff(funcs, vcat(x_vars, y_vars)) - all_symb = collect(all_symb) - return x_vars, y_vars, collect(u_vars), collect(all_symb) -end - -function macrohelper_extract_vars(equations::Array{Symbol, 1}) - return macrohelper_extract_vars(map(Expr, equations)) -end - -#------------------------------------------------------------------------------ - -function macrohelper_clean(ex::Expr) - ex = MacroTools.postwalk(x -> @capture(x, f_'(t)) ? f : x, ex) - ex = MacroTools.postwalk(x -> @capture(x, f_(t)) ? f : x, ex) - ex = MacroTools.postwalk(x -> x == :(/) ? :(//) : x, ex) - ex = MacroTools.postwalk(x -> x isa Float64 ? rationalize(x) : x, ex) - return ex -end - -#------------------------------------------------------------------------------ - -""" - macro ODEmodel - -Macro for creating an ODE from a list of equations. -It also injects all variables into the global scope. - -## Example - -Creating a simple `ODE`: - -```jldoctest -using StructuralIdentifiability - -ode = @ODEmodel( - x1'(t) = a * x1(t) + u(t), - x2'(t) = b * x2(t) + c*x1(t)*x2(t), - y(t) = x1(t) -) -``` - -Here, -- `x1`, `x2` are state variables -- `y` is an output variable -- `u` is an input variable -- `a`, `b`, `c` are time-indepdendent parameters - -""" -macro ODEmodel(ex::Expr...) - equations = [ex...] - x_vars, y_vars, u_vars, all_symb = macrohelper_extract_vars(equations) - time_dependent = vcat(x_vars, y_vars, u_vars) - params = sort([s for s in all_symb if !(s in time_dependent)]) - all_symb_no_t = vcat(time_dependent, params) - all_symb_with_t = vcat([:($s(t)) for s in time_dependent], params) - - # creating the polynomial ring - vars_list = :([$(all_symb_with_t...)]) - R = gensym() - vars_aux = gensym() - exp_ring = :( - ($R, $vars_aux) = StructuralIdentifiability.Nemo.PolynomialRing( - StructuralIdentifiability.Nemo.QQ, - map(string, $all_symb_with_t), - ) - ) - assignments = [:($(all_symb_no_t[i]) = $vars_aux[$i]) for i in 1:length(all_symb_no_t)] - - # setting x_vars and y_vars in the right order - vx = gensym() - vy = gensym() - x_var_expr = :($vx = Vector{StructuralIdentifiability.Nemo.fmpq_mpoly}([$(x_vars...)])) - y_var_expr = :($vy = Vector{StructuralIdentifiability.Nemo.fmpq_mpoly}([$(y_vars...)])) - - # preparing equations - equations = map(macrohelper_clean, equations) - x_dict = gensym() - y_dict = gensym() - x_dict_create_expr = :( - $x_dict = Dict{ - StructuralIdentifiability.Nemo.fmpq_mpoly, - Union{ - StructuralIdentifiability.Nemo.fmpq_mpoly, - StructuralIdentifiability.AbstractAlgebra.Generic.Frac{ - StructuralIdentifiability.Nemo.fmpq_mpoly, - }, - }, - }() - ) - y_dict_create_expr = :( - $y_dict = Dict{ - StructuralIdentifiability.Nemo.fmpq_mpoly, - Union{ - StructuralIdentifiability.Nemo.fmpq_mpoly, - StructuralIdentifiability.AbstractAlgebra.Generic.Frac{ - StructuralIdentifiability.Nemo.fmpq_mpoly, - }, - }, - }() - ) - eqs_expr = [] - for eq in equations - if eq.head != :(=) - throw("Problem with parsing at $eq") - end - lhs, rhs = eq.args[1:2] - loc_all_symb = macrohelper_extract_vars([rhs])[4] - to_insert = undef - if lhs in x_vars - to_insert = x_dict - elseif lhs in y_vars - to_insert = y_dict - else - throw("Unknown left-hand side $lhs") - end - - uniqueness_check_expr = quote - if haskey($to_insert, $lhs) - throw( - DomainError( - $lhs, - "The variable occurs twice in the left-hand-side of the ODE system", - ), - ) - end - end - push!(eqs_expr, uniqueness_check_expr) - if isempty(loc_all_symb) - push!(eqs_expr, :($to_insert[$lhs] = $R($rhs))) - else - push!(eqs_expr, :($to_insert[$lhs] = ($rhs))) - end - end - - for n in all_symb_no_t - if !Base.isidentifier(n) - throw( - ArgumentError( - "The names of the variables will be injected into the global scope, so their name must be allowed Julia names, $n is not", - ), - ) - end - end - - logging_exprs = [ - :( - StructuralIdentifiability.Logging.with_logger( - StructuralIdentifiability._si_logger[], - ) do - @info "Summary of the model:" - @info "State variables: " * $(join(map(string, collect(x_vars)), ", ")) - @info "Parameters: " * $(join(map(string, collect(params)), ", ")) - @info "Inputs: " * $(join(map(string, collect(u_vars)), ", ")) - @info "Outputs: " * $(join(map(string, collect(y_vars)), ", ")) - end - ), - ] - # creating the ode object - ode_expr = :(StructuralIdentifiability.ODE{StructuralIdentifiability.Nemo.fmpq_mpoly}( - $vx, - $vy, - $x_dict, - $y_dict, - Array{StructuralIdentifiability.Nemo.fmpq_mpoly}([$(u_vars...)]), - )) - - result = Expr( - :block, - logging_exprs..., - exp_ring, - assignments..., - x_var_expr, - y_var_expr, - x_dict_create_expr, - y_dict_create_expr, - eqs_expr..., - ode_expr, - ) - return esc(result) -end - -#------------------------------------------------------------------------------ - function Base.show(io::IO, ode::ODE) for x in ode.x_vars if endswith(var_to_str(x), "(t)") @@ -515,155 +279,3 @@ function Base.show(io::IO, ode::ODE) end #------------------------------------------------------------------------------ -""" - function mtk_to_si(de::ModelingToolkit.AbstractTimeDependentSystem, measured_quantities::Array{ModelingToolkit.Equation}) - function mtk_to_si(de::ModelingToolkit.AbstractTimeDependentSystem, measured_quantities::Array{SymbolicUtils.BasicSymbolic}) - -Input: -- `de` - ModelingToolkit.AbstractTimeDependentSystem, a system for identifiability query -- `measured_quantities` - array of output functions (as equations of just functions) - -Output: -- `ODE` object containing required data for identifiability assessment -- `conversion` dictionary from the symbols in the input MTK model to the variable - involved in the produced `ODE` object -""" -function mtk_to_si( - de::ModelingToolkit.AbstractTimeDependentSystem, - measured_quantities::Array{ModelingToolkit.Equation}, -) - return __mtk_to_si( - de, - [(replace(string(e.lhs), "(t)" => ""), e.rhs) for e in measured_quantities], - ) -end - -function mtk_to_si( - de::ModelingToolkit.AbstractTimeDependentSystem, - measured_quantities::Array{<:Symbolics.Num}, -) - return __mtk_to_si( - de, - [("y$i", Symbolics.value(e)) for (i, e) in enumerate(measured_quantities)], - ) -end - -function mtk_to_si( - de::ModelingToolkit.AbstractTimeDependentSystem, - measured_quantities::Array{<:SymbolicUtils.BasicSymbolic}, -) - return __mtk_to_si(de, [("y$i", e) for (i, e) in enumerate(measured_quantities)]) -end - -#------------------------------------------------------------------------------ -# old name kept for compatibility purposes - -function preprocess_ode( - de::ModelingToolkit.AbstractTimeDependentSystem, - measured_quantities::Array{ModelingToolkit.Equation}, -) - @warn "Function `preprocess_ode` has been renamed to `mtk_to_si`. The old name can be still used but will disappear in the future releases." - return mtk_to_si(de, measured_quantities) -end - -function preprocess_ode( - de::ModelingToolkit.AbstractTimeDependentSystem, - measured_quantities::Array{<:Symbolics.Num}, -) - @warn "Function `preprocess_ode` has been renamed to `mtk_to_si`. The old name can be still used but will disappear in the future releases." - return mtk_to_si(de, measured_quantities) -end - -function preprocess_ode( - de::ModelingToolkit.AbstractTimeDependentSystem, - measured_quantities::Array{<:SymbolicUtils.BasicSymbolic}, -) - @warn "Function `preprocess_ode` has been renamed to `mtk_to_si`. The old name can be still used but will disappear in the future releases." - return mtk_to_si(de, measured_quantities) -end - -#------------------------------------------------------------------------------ -""" - function __mtk_to_si(de::ModelingToolkit.AbstractTimeDependentSystem, measured_quantities::Array{Tuple{String, SymbolicUtils.BasicSymbolic}}) - -Input: -- `de` - ModelingToolkit.AbstractTimeDependentSystem, a system for identifiability query -- `measured_quantities` - array of input function in the form (name, expression) - -Output: -- `ODE` object containing required data for identifiability assessment -- `conversion` dictionary from the symbols in the input MTK model to the variable - involved in the produced `ODE` object -""" -function __mtk_to_si( - de::ModelingToolkit.AbstractTimeDependentSystem, - measured_quantities::Array{<:Tuple{String, <:SymbolicUtils.BasicSymbolic}}, -) - polytype = StructuralIdentifiability.Nemo.fmpq_mpoly - fractype = StructuralIdentifiability.Nemo.Generic.Frac{polytype} - diff_eqs = - filter(eq -> !(ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(de)) - # performing full structural simplification - if length(observed(de)) > 0 - rules = Dict(s.lhs => s.rhs for s in observed(de)) - while any([ - length(intersect(get_variables(r), keys(rules))) > 0 for r in values(rules) - ]) - rules = Dict(k => SymbolicUtils.substitute(v, rules) for (k, v) in rules) - end - diff_eqs = [SymbolicUtils.substitute(eq, rules) for eq in diff_eqs] - end - - y_functions = [each[2] for each in measured_quantities] - inputs = filter(v -> ModelingToolkit.isinput(v), ModelingToolkit.states(de)) - state_vars = filter( - s -> !(ModelingToolkit.isinput(s) || ModelingToolkit.isoutput(s)), - ModelingToolkit.states(de), - ) - params = ModelingToolkit.parameters(de) - t = ModelingToolkit.arguments(diff_eqs[1].lhs)[1] - params_from_measured_quantities = union( - [filter(s -> !istree(s), get_variables(y[2])) for y in measured_quantities]..., - ) - params = union(params, params_from_measured_quantities) - - input_symbols = vcat(state_vars, inputs, params) - generators = vcat(string.(input_symbols), [e[1] for e in measured_quantities]) - generators = map(g -> replace(g, "(t)" => ""), generators) - R, gens_ = Nemo.PolynomialRing(Nemo.QQ, generators) - y_vars = Vector{polytype}([str_to_var(e[1], R) for e in measured_quantities]) - symb2gens = Dict(input_symbols .=> gens_[1:length(input_symbols)]) - - x_vars = Vector{polytype}() - - state_eqn_dict = Dict{polytype, Union{polytype, fractype}}() - out_eqn_dict = Dict{polytype, Union{polytype, fractype}}() - - for i in 1:length(diff_eqs) - x = substitute(state_vars[i], symb2gens) - push!(x_vars, x) - if !(diff_eqs[i].rhs isa Number) - state_eqn_dict[x] = eval_at_nemo(diff_eqs[i].rhs, symb2gens) - else - state_eqn_dict[x] = R(diff_eqs[i].rhs) - end - end - for i in 1:length(measured_quantities) - out_eqn_dict[y_vars[i]] = eval_at_nemo(measured_quantities[i][2], symb2gens) - end - - inputs_ = [substitute(each, symb2gens) for each in inputs] - if isequal(length(inputs_), 0) - inputs_ = Vector{polytype}() - end - return ( - StructuralIdentifiability.ODE{polytype}( - x_vars, - y_vars, - state_eqn_dict, - out_eqn_dict, - inputs_, - ), - symb2gens, - ) -end diff --git a/src/ODEexport.jl b/src/ODEexport.jl index aa4f8d63d..c74898476 100644 --- a/src/ODEexport.jl +++ b/src/ODEexport.jl @@ -10,7 +10,7 @@ function print_for_maple(ode::ODE, package = :SIAN) varstr = Dict(x => var_to_str(x) * "(t)" for x in vcat(ode.x_vars, ode.u_vars, ode.y_vars)) merge!(varstr, Dict(p => var_to_str(p) for p in ode.parameters)) - R_print, vars_print = Nemo.PolynomialRing( + R_print, vars_print = Nemo.polynomial_ring( base_ring(ode.poly_ring), [varstr[v] for v in gens(ode.poly_ring)], ) @@ -189,7 +189,7 @@ function print_for_COMBOS(ode::ODE) merge!(varstr, Dict(u => "u" * string(ind) for (ind, u) in enumerate(ode.u_vars))) merge!(varstr, Dict(y => "y" * string(ind) for (ind, y) in enumerate(ode.y_vars))) merge!(varstr, Dict(p => var_to_str(p) for p in ode.parameters)) - R_print, vars_print = Nemo.PolynomialRing( + R_print, vars_print = Nemo.polynomial_ring( base_ring(ode.poly_ring), [varstr[v] for v in gens(ode.poly_ring)], ) diff --git a/src/RationalFunctionFields/IdealMQS.jl b/src/RationalFunctionFields/IdealMQS.jl index 6d8ddb81c..95b9fa060 100644 --- a/src/RationalFunctionFields/IdealMQS.jl +++ b/src/RationalFunctionFields/IdealMQS.jl @@ -100,7 +100,7 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal end varnames = append_at_index(ystrs, sat_var_index, sat_varname) @debug "Saturating variable is $sat_varname, index is $sat_var_index" - R_sat, v_sat = Nemo.PolynomialRing(K, varnames, ordering = ordering) + R_sat, v_sat = Nemo.polynomial_ring(K, varnames, ordering = ordering) # Saturation t_sat = v_sat[sat_var_index] den_lcm_orig = den_lcm @@ -144,7 +144,7 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal push!(nums_qq, num) end end - parent_ring_param, _ = PolynomialRing(ring, varnames, ordering = ordering) + parent_ring_param, _ = polynomial_ring(ring, varnames, ordering = ordering) @debug "Constructed MQS ideal in $R_sat with $(length(nums_qq) + 1) elements" @assert length(pivots_indices) == length(dens_indices) == length(dens_qq) @assert length(pivots_indices) == length(funcs_den_nums) @@ -202,7 +202,7 @@ function fractionfree_generators_raw(mqs::IdealMQS) end # NOTE: new variables go first! big_ring, big_vars = - PolynomialRing(K, vcat(new_varnames, old_varnames), ordering = :lex) + polynomial_ring(K, vcat(new_varnames, old_varnames), ordering = :lex) @info "$(mqs.sat_var_index) $(varnames) $ring_params $(parent(mqs.sat_qq))" nums_qq, dens_qq, sat_qq = mqs.nums_qq, mqs.dens_qq, mqs.sat_qq nums_y = map(num -> parent_ring_change(num, big_ring, matching = :byindex), nums_qq) @@ -229,7 +229,7 @@ end function ParamPunPam.reduce_mod_p!( mqs::IdealMQS, ff::Field, -) where {Field <: Union{Nemo.GaloisField, Nemo.GaloisFmpzField}} +) where {Field <: Union{Nemo.fpField, Nemo.FpField}} @debug "Reducing MQS ideal modulo $(ff)" # If there is a reduction modulo this field already, if haskey(mqs.cached_nums_gf, ff) @@ -251,7 +251,7 @@ function ParamPunPam.specialize_mod_p( mqs::IdealMQS, point::Vector{T}; saturated = true, -) where {T <: Union{gfp_elem, gfp_fmpz_elem}} +) where {T <: Union{fpFieldElem, FpFieldElem}} K_1 = parent(first(point)) @debug "Evaluating MQS ideal over $K_1 at $point" @assert haskey(mqs.cached_nums_gf, K_1) @@ -283,7 +283,7 @@ function ParamPunPam.specialize_mod_p( return polys end -function specialize(mqs::IdealMQS, point::Vector{Nemo.fmpq}; saturated = true) +function specialize(mqs::IdealMQS, point::Vector{Nemo.QQFieldElem}; saturated = true) @debug "Evaluating MQS ideal over QQ at $point" nums_qq, dens_qq, sat_qq = mqs.nums_qq, mqs.dens_qq, mqs.sat_qq dens_indices = mqs.dens_indices diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index cf38a9da4..100252e17 100644 --- a/src/RationalFunctionFields/RationalFunctionField.jl +++ b/src/RationalFunctionFields/RationalFunctionField.jl @@ -194,7 +194,7 @@ end ) where {T} mqs_generators = generators.mqs mqs_tobereduced = tobereduced.mqs - ff = Nemo.GF(2^31 - 1) + ff = Nemo.Native.GF(2^31 - 1) reduce_mod_p!(mqs_generators, ff) reduce_mod_p!(mqs_tobereduced, ff) param_ring = ParamPunPam.parent_params(mqs_generators) diff --git a/src/RationalFunctionFields/normalforms.jl b/src/RationalFunctionFields/normalforms.jl index 80c6329c0..9477fe7f0 100644 --- a/src/RationalFunctionFields/normalforms.jl +++ b/src/RationalFunctionFields/normalforms.jl @@ -301,7 +301,7 @@ is not specified but is assumed to be close to 1. ring_param = ParamPunPam.parent_params(mqs) xs_param = gens(ring_param) nparams = nvars(ring_param) - finite_field = Nemo.GF(2^30 + 3) + finite_field = Nemo.Native.GF(2^30 + 3) ParamPunPam.reduce_mod_p!(mqs, finite_field) @info "Computing normal forms of degree $up_to_degree in $nparams variables" @debug """Variables ($nparams in total): $xs_param @@ -319,7 +319,7 @@ is not specified but is assumed to be close to 1. push!(relations_ff_1, monoms_ff_1[i]) push!(tref, exponent_vector(monoms_ff_1[i], 1)) end - complete_intersection_relations_ff = Vector{Nemo.gfp_mpoly}(undef, 0) + complete_intersection_relations_ff = Vector{Nemo.fpMPolyRingElem}(undef, 0) iters = 0 # Compute relations at several random points until a consensus is reached while true diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index 86c83993e..a18c98285 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -20,10 +20,8 @@ using ParamPunPam using ParamPunPam: reduce_mod_p!, specialize_mod_p, AbstractBlackboxIdeal ParamPunPam.enable_progressbar(false) -using ModelingToolkit - # defining a model -export ODE, @ODEmodel, mtk_to_si +export ODE, @ODEmodel, @DDSmodel, mtk_to_si # assessing identifiability export assess_local_identifiability, assess_identifiability @@ -72,6 +70,7 @@ include("lincomp.jl") include("pb_representation.jl") include("submodels.jl") include("discrete.jl") +include("input_macro.jl") function __init__() _si_logger[] = @static if VERSION >= v"1.7.0" @@ -101,7 +100,7 @@ function assess_identifiability( funcs_to_check = Vector(), prob_threshold::Float64 = 0.99, loglevel = Logging.Info, -) where {P <: MPolyElem{fmpq}} +) where {P <: MPolyRingElem{QQFieldElem}} restart_logging(loglevel = loglevel) reset_timings() with_logger(_si_logger[]) do @@ -117,7 +116,7 @@ function _assess_identifiability( ode::ODE{P}; funcs_to_check = Vector(), prob_threshold::Float64 = 0.99, -) where {P <: MPolyElem{fmpq}} +) where {P <: MPolyRingElem{QQFieldElem}} p_glob = 1 - (1 - prob_threshold) * 0.9 p_loc = 1 - (1 - prob_threshold) * 0.1 @@ -126,7 +125,7 @@ function _assess_identifiability( end @info "Assessing local identifiability" - trbasis = Array{fmpq_mpoly, 1}() + trbasis = Array{QQMPolyRingElem, 1}() runtime = @elapsed local_result = _assess_local_identifiability( ode, funcs_to_check = funcs_to_check, @@ -170,64 +169,6 @@ function _assess_identifiability( return result end -""" - assess_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=[], prob_threshold = 0.99, loglevel=Logging.Info) - -Input: -- `ode` - the ModelingToolkit.ODESystem object that defines the model -- `measured_quantities` - the output functions of the model -- `funcs_to_check` - functions of parameters for which to check the identifiability -- `prob_threshold` - probability of correctness. -- `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) - -Assesses identifiability (both local and global) of a given ODE model (parameters detected automatically). The result is guaranteed to be correct with the probability -at least `prob_threshold`. -""" -function assess_identifiability( - ode::ModelingToolkit.ODESystem; - measured_quantities = Array{ModelingToolkit.Equation}[], - funcs_to_check = [], - prob_threshold = 0.99, - loglevel = Logging.Info, -) - restart_logging(loglevel = loglevel) - with_logger(_si_logger[]) do - return _assess_identifiability( - ode, - measured_quantities = measured_quantities, - funcs_to_check = funcs_to_check, - prob_threshold = prob_threshold, - ) - end -end - -function _assess_identifiability( - ode::ModelingToolkit.ODESystem; - measured_quantities = Array{ModelingToolkit.Equation}[], - funcs_to_check = [], - prob_threshold = 0.99, -) - if isempty(measured_quantities) - measured_quantities = get_measured_quantities(ode) - end - - ode, conversion = mtk_to_si(ode, measured_quantities) - conversion_back = Dict(v => k for (k, v) in conversion) - if isempty(funcs_to_check) - funcs_to_check = [conversion_back[x] for x in [ode.x_vars..., ode.parameters...]] - end - funcs_to_check_ = [eval_at_nemo(each, conversion) for each in funcs_to_check] - - result = _assess_identifiability( - ode, - funcs_to_check = funcs_to_check_, - prob_threshold = prob_threshold, - ) - nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) - out_dict = OrderedDict(nemo2mtk[param] => result[param] for param in funcs_to_check_) - return out_dict -end - using PrecompileTools include("precompile.jl") diff --git a/src/discrete.jl b/src/discrete.jl index 193d4dcb3..3134fb5e2 100644 --- a/src/discrete.jl +++ b/src/discrete.jl @@ -1,8 +1,95 @@ +""" +The structure to represent a discrete dynamical system +with respect to *shift*. Internally just stores an ODE structure + +Can be constructed with @DDSmodel macro +""" +struct DDS{P} # P is the type of polynomials in the rhs of the DDS system + ode::ODE{P} + + function DDS{P}( + x_vars::Array{P, 1}, + y_vars::Array{P, 1}, + x_eqs::Dict{P, <:Union{P, Generic.Frac{P}}}, + y_eqs::Dict{P, <:Union{P, Generic.Frac{P}}}, + u_vars::Array{P, 1}, + ) where {P <: MPolyRingElem{<:FieldElem}} + new{P}(ODE{P}(x_vars, y_vars, x_eqs, y_eqs, u_vars)) + end + + function DDS{P}(ode::ODE{P}) where {P <: MPolyRingElem{<:FieldElem}} + new{P}(ode) + end +end + +#------------------------------------------------------------------------------ + +# getters + +function x_vars(dds::DDS) + return dds.ode.x_vars +end + +function y_vars(dds::DDS) + return dds.ode.y_vars +end + +function parameters(dds::DDS) + return dds.ode.parameters +end + +function inputs(dds::DDS) + return dds.ode.u_vars +end + +function x_equations(dds::DDS) + return dds.ode.x_equations +end + +function y_equations(dds::DDS) + return dds.ode.y_equations +end + +function Base.parent(dds::DDS) + return parent(dds.ode) +end + +#------------------------------------------------------------------------------ +# Some functions to transform DDS's + +function add_outputs( + dds::DDS{P}, + extra_y::Dict{String, <:RingElem}, +) where {P <: MPolyRingElem} + return DDS{P}(add_outputs(dds.ode, extra_y)) +end + +#------------------------------------------------------------------------------ + +function Base.show(io::IO, dds::DDS) + for x in x_vars(dds) + if endswith(var_to_str(x), "(t)") + print(io, var_to_str(x)[1:(end - 3)] * "(t + 1) = ") + else + print(io, var_to_str(x) * "(t + 1) = ") + end + print(io, x_equations(dds)[x]) + print(io, "\n") + end + for y in y_vars(dds) + print(io, var_to_str(y) * " = ") + print(io, y_equations(dds)[y]) + print(io, "\n") + end +end + +#------------------------------------------------------------------------------ + """ sequence_solution(dds, param_values, initial_conditions, input_values, num_terms) Input: -- `dds` - a discrete dynamical system to solve (represented as an ODE struct) +- `dds` - a discrete dynamical system to solve - `param_values` - parameter values, must be a dictionary mapping parameter to a value - `initial_conditions` - initial conditions of `ode`, must be a dictionary mapping state variable to a value - `input_values` - input sequences in the form input => list of terms; length of the lists must be at least @@ -13,21 +100,21 @@ Output: - computes a sequence solution with teh required number of terms prec presented as a dictionary state_variable => corresponding sequence """ function sequence_solution( - dds::ODE{P}, + dds::DDS{P}, param_values::Dict{P, T}, initial_conditions::Dict{P, T}, input_values::Dict{P, Array{T, 1}}, num_terms::Int, -) where {T <: FieldElem, P <: MPolyElem{T}} - result = Dict(x => [initial_conditions[x]] for x in dds.x_vars) +) where {T <: FieldElem, P <: MPolyRingElem{T}} + result = Dict(x => [initial_conditions[x]] for x in x_vars(dds)) for i in 2:num_terms eval_dict = merge( param_values, Dict(k => v[end] for (k, v) in result), Dict(u => val[i - 1] for (u, val) in input_values), ) - for x in dds.x_vars - push!(result[x], eval_at_dict(dds.x_equations[x], eval_dict)) + for x in x_vars(dds) + push!(result[x], eval_at_dict(x_equations(dds)[x], eval_dict)) end end return result @@ -36,13 +123,13 @@ end #------------------------------------------------------------------------------ function sequence_solution( - dds::ODE{P}, + dds::DDS{P}, param_values::Dict{P, Int}, initial_conditions::Dict{P, Int}, input_values::Dict{P, Array{Int, 1}}, num_terms::Int, -) where {P <: MPolyElem{<:FieldElem}} - bring = base_ring(dds.poly_ring) +) where {P <: MPolyRingElem{<:FieldElem}} + bring = base_ring(parent(dds)) return sequence_solution( dds, Dict(p => bring(v) for (p, v) in param_values), @@ -55,7 +142,7 @@ end #----------------------------------------------------------------------------- """ - differentiate_sequence_solution(dds, params, ic, inputs, num_terms) + differentiate_sequence_solution(dds, params, ic, input_values, num_terms) Input: - the same as for `sequence_solutions` @@ -66,40 +153,40 @@ Output: the function `u` w.r.t. `v` evaluated at the solution """ function differentiate_sequence_solution( - dds::ODE{P}, + dds::DDS{P}, params::Dict{P, T}, ic::Dict{P, T}, - inputs::Dict{P, Array{T, 1}}, + input_values::Dict{P, Array{T, 1}}, num_terms::Int, -) where {T <: Generic.FieldElem, P <: MPolyElem{T}} +) where {T <: Generic.FieldElem, P <: MPolyRingElem{T}} @debug "Computing the power series solution of the system" - seq_sol = sequence_solution(dds, params, ic, inputs, num_terms) - generalized_params = vcat(dds.x_vars, dds.parameters) - bring = base_ring(dds.poly_ring) + seq_sol = sequence_solution(dds, params, ic, input_values, num_terms) + generalized_params = vcat(x_vars(dds), parameters(dds)) + bring = base_ring(parent(dds)) @debug "Solving the variational system at the solution" part_diffs = Dict( - (x, p) => derivative(dds.x_equations[x], p) for x in dds.x_vars for + (x, p) => derivative(x_equations(dds)[x], p) for x in x_vars(dds) for p in generalized_params ) result = Dict( - (x, p) => [x == p ? one(bring) : zero(bring)] for x in dds.x_vars for + (x, p) => [x == p ? one(bring) : zero(bring)] for x in x_vars(dds) for p in generalized_params ) for i in 2:num_terms eval_dict = merge( params, Dict(k => v[i - 1] for (k, v) in seq_sol), - Dict(u => val[i - 1] for (u, val) in inputs), + Dict(u => val[i - 1] for (u, val) in input_values), ) for p in generalized_params - local_eval = Dict(x => result[(x, p)][end] for x in dds.x_vars) - for x in dds.x_vars + local_eval = Dict(x => result[(x, p)][end] for x in x_vars(dds)) + for x in x_vars(dds) res = sum([ eval_at_dict(part_diffs[(x, x2)], eval_dict) * local_eval[x2] for - x2 in dds.x_vars + x2 in x_vars(dds) ]) - if p in dds.parameters + if p in parameters(dds) res += eval_at_dict(part_diffs[(x, p)], eval_dict) end push!(result[(x, p)], res) @@ -113,45 +200,46 @@ end # ------------------------------------------------------------------------------ """ - differentiate_sequence_output(dds, params, ic, inputs, num_terms) + differentiate_sequence_output(dds, params, ic, input_values, num_terms) Similar to `differentiate_sequence_solution` but computes partial derivatives of prescribed outputs returns a dictionary of the form `y_function => Dict(var => dy/dvar)` where `dy/dvar` is the derivative of `y_function` with respect to `var`. """ function differentiate_sequence_output( - dds::ODE{P}, + dds::DDS{P}, params::Dict{P, T}, ic::Dict{P, T}, - inputs::Dict{P, Array{T, 1}}, + input_values::Dict{P, Array{T, 1}}, num_terms::Int, -) where {T <: Generic.FieldElem, P <: MPolyElem{T}} +) where {T <: Generic.FieldElem, P <: MPolyRingElem{T}} @debug "Computing partial derivatives of the solution" - seq_sol, sol_diff = differentiate_sequence_solution(dds, params, ic, inputs, num_terms) + seq_sol, sol_diff = + differentiate_sequence_solution(dds, params, ic, input_values, num_terms) @debug "Evaluating the partial derivatives of the outputs" - generalized_params = vcat(dds.x_vars, dds.parameters) + generalized_params = vcat(x_vars(dds), parameters(dds)) part_diffs = Dict( - (y, p) => derivative(dds.y_equations[y], p) for y in dds.y_vars for + (y, p) => derivative(y_equations(dds)[y], p) for y in y_vars(dds) for p in generalized_params ) - result = Dict((y, p) => [] for y in dds.y_vars for p in generalized_params) + result = Dict((y, p) => [] for y in y_vars(dds) for p in generalized_params) for i in 1:num_terms eval_dict = merge( params, Dict(k => v[i] for (k, v) in seq_sol), - Dict(u => val[i] for (u, val) in inputs), + Dict(u => val[i] for (u, val) in input_values), ) - for p in vcat(dds.x_vars, dds.parameters) - local_eval = Dict(x => sol_diff[(x, p)][i] for x in dds.x_vars) - for (y, y_eq) in dds.y_equations + for p in generalized_params + local_eval = Dict(x => sol_diff[(x, p)][i] for x in x_vars(dds)) + for (y, y_eq) in y_equations(dds) res = sum([ eval_at_dict(part_diffs[(y, x)], eval_dict) * local_eval[x] for - x in dds.x_vars + x in x_vars(dds) ]) - if p in dds.parameters + if p in parameters(dds) res += eval_at_dict(part_diffs[(y, p)], eval_dict) end push!(result[(y, p)], res) @@ -178,10 +266,9 @@ function _degree_with_common_denom(polys) end """ - _assess_local_identifiability_discrete_aux(dds::ODE{P}, funcs_to_check::Array{<: Any, 1}, known_ic, prob_threshold::Float64=0.99) where P <: MPolyElem{Nemo.fmpq} + _assess_local_identifiability_discrete_aux(dds::DDS{P}, funcs_to_check::Array{<: Any, 1}, known_ic, prob_threshold::Float64=0.99) where P <: MPolyRingElem{Nemo.QQFieldElem} -Checks the local identifiability/observability of the functions in `funcs_to_check` treating `dds` as a discrete-time system with **shift** -instead of derivative in the right-hand side. +Checks the local identifiability/observability of the functions in `funcs_to_check`. The result is correct with probability at least `prob_threshold`. `known_ic` can take one of the following * `:none` - no initial conditions are assumed to be known @@ -189,12 +276,12 @@ The result is correct with probability at least `prob_threshold`. * a list of rational functions in states and parameters assumed to be known at t = 0 """ function _assess_local_identifiability_discrete_aux( - dds::ODE{P}, + dds::DDS{P}, funcs_to_check::Array{<:Any, 1}, known_ic = :none, prob_threshold::Float64 = 0.99, -) where {P <: MPolyElem{Nemo.fmpq}} - bring = base_ring(dds.poly_ring) +) where {P <: MPolyRingElem{Nemo.QQFieldElem}} + bring = base_ring(parent(dds)) @debug "Extending the model" dds_ext = @@ -204,16 +291,16 @@ function _assess_local_identifiability_discrete_aux( known_ic = [] end if known_ic == :all - known_ic = dds_ext.x_vars + known_ic = x_vars(dds_ext) end @debug "Computing the observability matrix" - prec = length(dds.x_vars) + length(dds.parameters) + prec = length(x_vars(dds)) + length(parameters(dds)) @debug "The truncation order is $prec" # Computing the bound from the Schwartz-Zippel-DeMilo-Lipton lemma - deg_x = _degree_with_common_denom(values(dds.x_equations)) - deg_y = _degree_with_common_denom(values(dds.y_equations)) + deg_x = _degree_with_common_denom(values(x_equations(dds))) + deg_y = _degree_with_common_denom(values(y_equations(dds))) deg_known = reduce(+, map(total_degree, known_ic), init = 0) deg_to_check = max(map(total_degree, funcs_to_check)...) Jac_degree = deg_to_check + deg_known @@ -226,43 +313,44 @@ function _assess_local_identifiability_discrete_aux( @debug "Sampling range $D" # Parameter values are the same across all the replicas - params_vals = Dict(p => bring(rand(1:D)) for p in dds_ext.parameters) - ic = Dict(x => bring(rand(1:D)) for x in dds_ext.x_vars) - # TODO: parametric type instead of fmpq - inputs = Dict{P, Array{fmpq, 1}}( - u => [bring(rand(1:D)) for i in 1:prec] for u in dds_ext.u_vars + params_vals = Dict(p => bring(rand(1:D)) for p in parameters(dds_ext)) + ic = Dict(x => bring(rand(1:D)) for x in x_vars(dds_ext)) + # TODO: parametric type instead of QQFieldElem + input_values = Dict{P, Array{QQFieldElem, 1}}( + u => [bring(rand(1:D)) for i in 1:prec] for + u in StructuralIdentifiability.inputs(dds_ext) ) @debug "Computing the output derivatives" output_derivatives = - differentiate_sequence_output(dds_ext, params_vals, ic, inputs, prec) + differentiate_sequence_output(dds_ext, params_vals, ic, input_values, prec) @debug "Building the matrices" Jac = zero( - Nemo.MatrixSpace( + Nemo.matrix_space( bring, - length(dds.x_vars) + length(dds.parameters), - 1 + prec * length(dds.y_vars) + length(known_ic), + length(x_vars(dds)) + length(parameters(dds)), + 1 + prec * length(y_vars(dds)) + length(known_ic), ), ) - xs_params = vcat(dds_ext.x_vars, dds_ext.parameters) - for (i, y) in enumerate(dds.y_vars) - y = switch_ring(y, dds_ext.poly_ring) + xs_params = vcat(x_vars(dds_ext), parameters(dds_ext)) + for (i, y) in enumerate(y_vars(dds)) + y = switch_ring(y, parent(dds_ext)) for j in 1:prec - for (k, p) in enumerate(dds_ext.parameters) + for (k, p) in enumerate(parameters(dds_ext)) Jac[k, 1 + (i - 1) * prec + j] = output_derivatives[(y, p)][j] end - for (k, x) in enumerate(dds_ext.x_vars) + for (k, x) in enumerate(x_vars(dds_ext)) Jac[end - k + 1, 1 + (i - 1) * prec + j] = output_derivatives[(y, x)][j] end end end eval_point = merge(params_vals, ic) for (i, v) in enumerate(known_ic) - for (k, p) in enumerate(dds_ext.parameters) + for (k, p) in enumerate(parameters(dds_ext)) Jac[k, end - i + 1] = eval_at_dict(derivative(v, p), eval_point) end - for (k, x) in enumerate(dds_ext.x_vars) + for (k, x) in enumerate(x_vars(dds_ext)) Jac[end - k + 1, end - i + 1] = eval_at_dict(derivative(v, x), eval_point) end end @@ -271,13 +359,13 @@ function _assess_local_identifiability_discrete_aux( base_rank = LinearAlgebra.rank(Jac) result = OrderedDict{Any, Bool}() for i in 1:length(funcs_to_check) - for (k, p) in enumerate(dds_ext.parameters) + for (k, p) in enumerate(parameters(dds_ext)) Jac[k, 1] = - output_derivatives[(str_to_var("loc_aux_$i", dds_ext.poly_ring), p)][1] + output_derivatives[(str_to_var("loc_aux_$i", parent(dds_ext)), p)][1] end - for (k, x) in enumerate(dds_ext.x_vars) + for (k, x) in enumerate(x_vars(dds_ext)) Jac[end - k + 1, 1] = - output_derivatives[(str_to_var("loc_aux_$i", dds_ext.poly_ring), x)][1] + output_derivatives[(str_to_var("loc_aux_$i", parent(dds_ext)), x)][1] end result[funcs_to_check[i]] = LinearAlgebra.rank(Jac) == base_rank end @@ -288,103 +376,31 @@ end # ------------------------------------------------------------------------------ """ - function assess_local_identifiability( - dds::ModelingToolkit.DiscreteSystem; - measured_quantities=Array{ModelingToolkit.Equation}[], - funcs_to_check=Array{}[], - known_ic=Array{}[], - prob_threshold::Float64=0.99) - -Input: -- `dds` - the DiscreteSystem object from ModelingToolkit (with **difference** operator in the right-hand side) -- `measured_quantities` - the measurable outputs of the model -- `funcs_to_check` - functions of parameters for which to check identifiability (all parameters and states if not specified) -- `known_ic` - functions (of states and parameter) whose initial conditions are assumed to be known -- `prob_threshold` - probability of correctness - -Output: -- the result is an (ordered) dictionary from each function to to boolean; + assess_local_identifiability(dds::DDS{P}; funcs_to_check::Array{<: Any, 1}, known_ic, prob_threshold::Float64=0.99, loglevel=Logging.Info) where P <: MPolyRingElem{Nemo.QQFieldElem} -The result is correct with probability at least `prob_threshold`. +Checks the local identifiability/observability of the functions in `funcs_to_check`. The result is correct with probability at least `prob_threshold`. +A list of quantities can be provided as `known_ic` for which the initial conditions can be assumed to be known and generic. """ function assess_local_identifiability( - dds::ModelingToolkit.DiscreteSystem; - measured_quantities = Array{ModelingToolkit.Equation}[], - funcs_to_check = Array{}[], - known_ic = Array{}[], + dds::DDS{P}; + funcs_to_check::Array{<:Any, 1} = Array{Any, 1}(), + known_ic = :none, prob_threshold::Float64 = 0.99, loglevel = Logging.Info, -) +) where {P <: MPolyRingElem{Nemo.QQFieldElem}} restart_logging(loglevel = loglevel) + reset_timings() with_logger(_si_logger[]) do - return _assess_local_identifiability( + if isempty(funcs_to_check) + funcs_to_check = vcat(parameters(dds), x_vars(dds)) + end + return _assess_local_identifiability_discrete_aux( dds, - measured_quantities = measured_quantities, - funcs_to_check = funcs_to_check, - known_ic = known_ic, - prob_threshold = prob_threshold, + funcs_to_check, + known_ic, + prob_threshold, ) end end -function _assess_local_identifiability( - dds::ModelingToolkit.DiscreteSystem; - measured_quantities = Array{ModelingToolkit.Equation}[], - funcs_to_check = Array{}[], - known_ic = Array{}[], - prob_threshold::Float64 = 0.99, -) - if length(measured_quantities) == 0 - if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(dds)) - @info "Measured quantities are not provided, trying to find the outputs in input dynamical system." - measured_quantities = filter( - eq -> (ModelingToolkit.isoutput(eq.lhs)), - ModelingToolkit.equations(dds), - ) - else - throw( - error( - "Measured quantities (output functions) were not provided and no outputs were found.", - ), - ) - end - end - - # Converting the finite difference operator in the right-hand side to - # the corresponding shift operator - eqs = filter(eq -> !(ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(dds)) - deltas = [Symbolics.operation(e.lhs).dt for e in eqs] - @assert length(Set(deltas)) == 1 - eqs_shift = [e.lhs ~ e.rhs + first(Symbolics.arguments(e.lhs)) for e in eqs] - dds_shift = DiscreteSystem(eqs_shift, name = gensym()) - @debug "System transformed from difference to shift: $dds_shift" - - dds_aux, conversion = mtk_to_si(dds_shift, measured_quantities) - if length(funcs_to_check) == 0 - params = parameters(dds) - params_from_measured_quantities = union( - [filter(s -> !istree(s), get_variables(y)) for y in measured_quantities]..., - ) - funcs_to_check = vcat( - [x for x in states(dds) if conversion[x] in dds_aux.x_vars], - union(params, params_from_measured_quantities), - ) - end - funcs_to_check_ = [eval_at_nemo(x, conversion) for x in funcs_to_check] - known_ic_ = [eval_at_nemo(x, conversion) for x in known_ic] - - result = _assess_local_identifiability_discrete_aux( - dds_aux, - funcs_to_check_, - known_ic_, - prob_threshold, - ) - nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) - out_dict = OrderedDict(nemo2mtk[param] => result[param] for param in funcs_to_check_) - if length(known_ic) > 0 - @warn "Since known initial conditions were provided, identifiability of states (e.g., `x(t)`) is at t = 0 only !" - out_dict = OrderedDict(substitute(k, Dict(t => 0)) => v for (k, v) in out_dict) - end - return out_dict -end # ------------------------------------------------------------------------------ diff --git a/src/elimination.jl b/src/elimination.jl index 6798889e9..426bd924a 100644 --- a/src/elimination.jl +++ b/src/elimination.jl @@ -58,14 +58,14 @@ Inputs: Output: - `M::MatrixElem` - The Bezout matrix """ -function Bezout_matrix(f::P, g::P, var_elim::P) where {P <: MPolyElem} +function Bezout_matrix(f::P, g::P, var_elim::P) where {P <: MPolyRingElem} parent_ring = parent(f) deg_f = Nemo.degree(f, var_elim) deg_g = Nemo.degree(g, var_elim) n = max(deg_f, deg_g) coeffs_f = [coeff(f, [var_elim], [i]) for i in 0:n] coeffs_g = [coeff(g, [var_elim], [i]) for i in 0:n] - GL = AbstractAlgebra.MatrixSpace(parent_ring, n, n) + GL = AbstractAlgebra.matrix_space(parent_ring, n, n) M = zero(GL) for i in 1:n for j in 1:n @@ -91,12 +91,12 @@ Inputs: Output: - `M::MatrixElem` - The Sylvester matrix """ -function Sylvester_matrix(f::P, g::P, var_elim::P) where {P <: MPolyElem} +function Sylvester_matrix(f::P, g::P, var_elim::P) where {P <: MPolyRingElem} parent_ring = parent(f) deg_f = Nemo.degree(f, var_elim) deg_g = Nemo.degree(g, var_elim) n = deg_f + deg_g - GL = AbstractAlgebra.MatrixSpace(parent_ring, n, n) + GL = AbstractAlgebra.matrix_space(parent_ring, n, n) M = zero(GL) for i in 1:deg_f for j in 0:deg_g @@ -124,9 +124,9 @@ Input: Output: - `M::MatrixElem` - Simplified matrix -- `extra_factors::Vector{AbstractAlgebra.MPolyElem}` - array of GCDs eliminated from `M`. +- `extra_factors::Vector{AbstractAlgebra.MPolyRingElem}` - array of GCDs eliminated from `M`. """ -function simplify_matrix(M::MatElem{P}) where {P <: MPolyElem} +function simplify_matrix(M::MatElem{P}) where {P <: MPolyRingElem} """ An auxiliary function taking a list of coordinates of cells and dividing them by their gcd. @@ -188,7 +188,10 @@ mutable struct ODEPointGenerator{P} <: PointGenerator{P} cached_points::Array{Dict{P, <:FieldElem}, 1} number_type::Type - function ODEPointGenerator{P}(ode::ODE{P}, big_ring::MPolyRing) where {P <: MPolyElem} + function ODEPointGenerator{P}( + ode::ODE{P}, + big_ring::MPolyRing, + ) where {P <: MPolyRingElem} prec = length(ode.x_vars) + 1 number_type = typeof(one(base_ring(big_ring))) return new(ode, big_ring, prec, Array{Dict{P, number_type}}[], number_type) @@ -200,7 +203,7 @@ end function Base.iterate( gpg::ODEPointGenerator{P}, i::Int = 1, -) where {P <: MPolyElem{<:FieldElem}} +) where {P <: MPolyRingElem{<:FieldElem}} if i > length(gpg.cached_points) @debug "Generating new point on the variety" sample_max = i * 50 @@ -273,7 +276,7 @@ Output: function choose( polys::Array{P, 1}, generic_point_generator, -) where {P <: MPolyElem{<:FieldElem}} +) where {P <: MPolyRingElem{<:FieldElem}} vars = gens(parent(polys[1])) for p in generic_point_generator # get accounts for the fact that the big ring may contain some auxiliary variables, e.g. rand_proj_var @@ -307,7 +310,7 @@ Output: g::P, var_elim::P, generic_point_generator, -) where {P <: MPolyElem{<:FieldElem}} +) where {P <: MPolyRingElem{<:FieldElem}} # Linear comb while f != 0 && g != 0 if Nemo.degree(f, var_elim) > Nemo.degree(g, var_elim) @@ -381,7 +384,7 @@ Output: flush(_si_logger[].stream) M_simp, matrix_factors = simplify_matrix(M) @debug "Removed factors $(map(length, matrix_factors))" - M_size = zero(Nemo.MatrixSpace(Nemo.ZZ, ncols(M_simp), ncols(M_simp))) + M_size = zero(Nemo.matrix_space(Nemo.ZZ, ncols(M_simp), ncols(M_simp))) for i in 1:ncols(M_simp) for j in 1:ncols(M_simp) M_size[i, j] = length(M_simp[i, j]) diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl index 984747e3b..da93e6d74 100644 --- a/src/global_identifiability.jl +++ b/src/global_identifiability.jl @@ -17,14 +17,14 @@ are identifiable functions containing or not the state variables ode::ODE{P}, known::Array{P, 1}, with_states::Bool, -) where {P <: MPolyElem{fmpq}} +) where {P <: MPolyRingElem{QQFieldElem}} coeff_lists = Dict(:with_states => Array{Array{P, 1}, 1}(), :no_states => Array{Array{P, 1}, 1}()) varnames = [var_to_str(p) for p in ode.parameters] if with_states append!(varnames, map(var_to_str, ode.x_vars)) end - bring, _ = Nemo.PolynomialRing(base_ring(ode.poly_ring), varnames) + bring, _ = Nemo.polynomial_ring(base_ring(ode.poly_ring), varnames) if with_states @debug "Computing Lie derivatives" @@ -67,7 +67,7 @@ are identifiable functions containing or not the state variables if with_states new_vars = vcat(new_vars, ode.x_vars) end - new_ring, _ = PolynomialRing(Nemo.QQ, map(Symbol, new_vars)) + new_ring, _ = polynomial_ring(Nemo.QQ, map(Symbol, new_vars)) new_coeff_lists = empty(coeff_lists) for (key, coeff_list) in coeff_lists new_coeff_lists[key] = @@ -145,7 +145,7 @@ The function returns a tuple containing the following: if with_states && !isempty(ode.parameters) @debug "Generators of identifiable functions involve states, the parameter-only part is getting simplified" # NOTE: switching to a ring without states for a moment - param_ring, _ = PolynomialRing( + param_ring, _ = polynomial_ring( base_ring(bring), map(string, ode.parameters), ordering = Nemo.ordering(bring), @@ -195,7 +195,7 @@ Output: a list L of booleans with L[i] being the identifiability status of the i known::Array{P, 1} = Array{P, 1}(), prob_threshold::Float64 = 0.99, var_change_policy = :default, -) where {P <: MPolyElem{fmpq}} +) where {P <: MPolyRingElem{QQFieldElem}} states_needed = false for f in funcs_to_check num, den = unpack_fraction(f) @@ -236,7 +236,7 @@ function check_identifiability( known::Array{P, 1} = Array{P, 1}(), prob_threshold::Float64 = 0.99, var_change_policy = :default, -) where {P <: MPolyElem{fmpq}} +) where {P <: MPolyRingElem{QQFieldElem}} return check_identifiability( ode, ode.parameters, @@ -248,7 +248,7 @@ end #------------------------------------------------------------------------------ """ - assess_global_identifiability(ode::ODE{P}, prob_threshold::Float64=0.99; var_change=:default) where P <: MPolyElem{fmpq} + assess_global_identifiability(ode::ODE{P}, prob_threshold::Float64=0.99; var_change=:default) where P <: MPolyRingElem{QQFieldElem} Input: - `ode` - the ODE model @@ -266,7 +266,7 @@ function assess_global_identifiability( known::Array{P, 1} = Array{P, 1}(), prob_threshold::Float64 = 0.99; var_change = :default, -) where {P <: MPolyElem{fmpq}} +) where {P <: MPolyRingElem{QQFieldElem}} result_list = assess_global_identifiability( ode, ode.parameters, @@ -303,7 +303,7 @@ Checks global identifiability of functions of parameters specified in `funcs_to_ known::Array{P, 1} = Array{P, 1}(), prob_threshold::Float64 = 0.99; var_change = :default, -) where {P <: MPolyElem{fmpq}} +) where {P <: MPolyRingElem{QQFieldElem}} submodels = find_submodels(ode) if length(submodels) > 0 @info "Note: the input model has nontrivial submodels. If the computation for the full model will be too heavy, you may want to try to first analyze one of the submodels. They can be produced using function `find_submodels`" diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index 79db0ffbd..ec116843e 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -37,7 +37,7 @@ ode = @ODEmodel( find_identifiable_functions(ode) # prints -3-element Vector{AbstractAlgebra.Generic.Frac{Nemo.fmpq_mpoly}}: +3-element Vector{AbstractAlgebra.Generic.Frac{Nemo.QQMPolyRingElem}}: a12 + a01 + a21 a12*a01 ``` @@ -51,7 +51,7 @@ function find_identifiable_functions( simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, loglevel = Logging.Info, -) where {T <: MPolyElem{fmpq}} +) where {T <: MPolyRingElem{QQFieldElem}} restart_logging(loglevel = loglevel) reset_timings() with_logger(_si_logger[]) do @@ -73,7 +73,7 @@ function _find_identifiable_functions( with_states = false, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, -) where {T <: MPolyElem{fmpq}} +) where {T <: MPolyRingElem{QQFieldElem}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) runtime_start = time_ns() @@ -117,93 +117,3 @@ function _find_identifiable_functions( return id_funcs_fracs end - -""" - find_identifiable_functions(ode::ModelingToolkit.ODESystem; measured_quantities=[], options...) - -Finds all functions of parameters/states that are identifiable in the given ODE -system. - -## Options - -This functions takes the following optional arguments: -- `measured_quantities` - the output functions of the model. -- `loglevel` - the verbosity of the logging - (can be Logging.Error, Logging.Warn, Logging.Info, Logging.Debug) - -## Example - -```jldoctest -using StructuralIdentifiability -using ModelingToolkit - -@parameters a01 a21 a12 -@variables t x0(t) x1(t) y1(t) [output = true] -D = Differential(t) - -eqs = [ - D(x0) ~ -(a01 + a21) * x0 + a12 * x1, - D(x1) ~ a21 * x0 - a12 * x1, y1 ~ x0 -] -de = ODESystem(eqs, t, name = :Test) - -find_identifiable_functions(de, measured_quantities = [y1 ~ x0]) - -# prints -2-element Vector{Num}: - a01*a12 - a01 + a12 + a21 -``` -""" -function find_identifiable_functions( - ode::ModelingToolkit.ODESystem; - measured_quantities = Array{ModelingToolkit.Equation}[], - prob_threshold::Float64 = 0.99, - seed = 42, - with_states = false, - simplify = :standard, - rational_interpolator = :VanDerHoevenLecerf, - loglevel = Logging.Info, -) - restart_logging(loglevel = loglevel) - reset_timings() - with_logger(_si_logger[]) do - return _find_identifiable_functions( - ode, - measured_quantities = measured_quantities, - prob_threshold = prob_threshold, - seed = seed, - with_states = with_states, - simplify = simplify, - rational_interpolator = rational_interpolator, - ) - end -end - -function _find_identifiable_functions( - ode::ModelingToolkit.ODESystem; - measured_quantities = Array{ModelingToolkit.Equation}[], - prob_threshold::Float64 = 0.99, - seed = 42, - with_states = false, - simplify = :standard, - rational_interpolator = :VanDerHoevenLecerf, -) - Random.seed!(seed) - if isempty(measured_quantities) - measured_quantities = get_measured_quantities(ode) - end - ode, conversion = mtk_to_si(ode, measured_quantities) - result = _find_identifiable_functions( - ode, - simplify = simplify, - prob_threshold = prob_threshold, - seed = seed, - with_states = with_states, - rational_interpolator = rational_interpolator, - ) - result = [parent_ring_change(f, ode.poly_ring) for f in result] - nemo2mtk = Dict(v => Num(k) for (k, v) in conversion) - out_funcs = [eval_at_dict(func, nemo2mtk) for func in result] - return out_funcs -end diff --git a/src/input_macro.jl b/src/input_macro.jl new file mode 100644 index 000000000..3743117c7 --- /dev/null +++ b/src/input_macro.jl @@ -0,0 +1,307 @@ +function _extract_aux!(funcs, all_symb, eq; ders_ok = false, type = :ode) + aux_symb = Set([:(+), :(-), :(=), :(*), :(^), :t, :(/), :(//)]) + MacroTools.postwalk( + x -> begin + if @capture(x, f_'(t)) + if !ders_ok + throw( + Base.ArgumentError( + "Derivative are not allowed in the right-hand side", + ), + ) + end + if type != :ode + throw( + Base.ArgumentError( + "Derivative are not expected in the discrete case", + ), + ) + end + push!(all_symb, f) + elseif @capture(x, f_(t + 1)) + if !ders_ok + throw(Base.ArgumentError("Shifts are not allowed in the right-hand side")) + end + if type != :dds + throw( + Base.ArgumentError( + "Shifts are not expected in the differential case", + ), + ) + end + push!(all_symb, f) + elseif @capture(x, f_(t)) + push!(funcs, f) + elseif (x isa Symbol) && !(x in aux_symb) + push!(all_symb, x) + end + return x + end, + eq, + ) +end + +""" + For an expression of the form f'(t)/f(t + 1) or f(t) returns (f, true) and (f, false), resp +""" +function _get_var(expr, type = :ode) + if @capture(expr, f_'(t)) + @assert type == :ode + return (f, true) + end + if @capture(expr, f_(t + 1)) + @assert type == :dds + return (f, true) + end + if @capture(expr, f_(t)) + return (f, false) + end + error("cannot extract the single function name from $expr") +end + +function macrohelper_extract_vars(equations::Array{Expr, 1}, type = :ode) + funcs, all_symb = Set(), Set() + x_vars, y_vars = Vector(), Vector() + for eq in equations + if eq.head != :(=) + _extract_aux!(funcs, all_symb, eq, type = type) + else + lhs, rhs = eq.args[1:2] + _extract_aux!(funcs, all_symb, lhs, ders_ok = true, type = type) + _extract_aux!(funcs, all_symb, rhs, type = type) + (v, is_state) = _get_var(lhs, type) + if is_state + push!(x_vars, v) + else + push!(y_vars, v) + end + end + end + u_vars = setdiff(funcs, vcat(x_vars, y_vars)) + all_symb = collect(all_symb) + return x_vars, y_vars, collect(u_vars), collect(all_symb) +end + +function macrohelper_extract_vars(equations::Array{Symbol, 1}, type = :ode) + return macrohelper_extract_vars(map(Expr, equations), type) +end + +#------------------------------------------------------------------------------ + +function macrohelper_clean(ex::Expr) + ex = MacroTools.postwalk(x -> @capture(x, f_'(t)) ? f : x, ex) + ex = MacroTools.postwalk(x -> @capture(x, f_(t + 1)) ? f : x, ex) + ex = MacroTools.postwalk(x -> @capture(x, f_(t)) ? f : x, ex) + ex = MacroTools.postwalk(x -> x == :(/) ? :(//) : x, ex) + ex = MacroTools.postwalk(x -> x isa Float64 ? rationalize(x) : x, ex) + return ex +end + +#------------------------------------------------------------------------------ + +function generate_model_code(type, ex::Expr...) + @assert type in (:ode, :dds) + equations = [ex...] + x_vars, y_vars, u_vars, all_symb = macrohelper_extract_vars(equations, type) + time_dependent = vcat(x_vars, y_vars, u_vars) + params = sort([s for s in all_symb if !(s in time_dependent)]) + all_symb_no_t = vcat(time_dependent, params) + all_symb_with_t = vcat([:($s(t)) for s in time_dependent], params) + + # creating the polynomial ring + R = gensym() + vars_aux = gensym() + exp_ring = :( + ($R, $vars_aux) = StructuralIdentifiability.Nemo.polynomial_ring( + StructuralIdentifiability.Nemo.QQ, + map(string, $all_symb_with_t), + ) + ) + assignments = [:($(all_symb_no_t[i]) = $vars_aux[$i]) for i in 1:length(all_symb_no_t)] + + # setting x_vars and y_vars in the right order + vx = gensym() + vy = gensym() + x_var_expr = + :($vx = Vector{StructuralIdentifiability.Nemo.QQMPolyRingElem}([$(x_vars...)])) + y_var_expr = + :($vy = Vector{StructuralIdentifiability.Nemo.QQMPolyRingElem}([$(y_vars...)])) + + # preparing equations + equations = map(macrohelper_clean, equations) + x_dict = gensym() + y_dict = gensym() + x_dict_create_expr = :( + $x_dict = Dict{ + StructuralIdentifiability.Nemo.QQMPolyRingElem, + Union{ + StructuralIdentifiability.Nemo.QQMPolyRingElem, + StructuralIdentifiability.AbstractAlgebra.Generic.Frac{ + StructuralIdentifiability.Nemo.QQMPolyRingElem, + }, + }, + }() + ) + y_dict_create_expr = :( + $y_dict = Dict{ + StructuralIdentifiability.Nemo.QQMPolyRingElem, + Union{ + StructuralIdentifiability.Nemo.QQMPolyRingElem, + StructuralIdentifiability.AbstractAlgebra.Generic.Frac{ + StructuralIdentifiability.Nemo.QQMPolyRingElem, + }, + }, + }() + ) + eqs_expr = [] + for eq in equations + if eq.head != :(=) + throw("Problem with parsing at $eq") + end + lhs, rhs = eq.args[1:2] + loc_all_symb = macrohelper_extract_vars([rhs], type)[4] + to_insert = undef + if lhs in x_vars + to_insert = x_dict + elseif lhs in y_vars + to_insert = y_dict + else + throw("Unknown left-hand side $lhs") + end + + uniqueness_check_expr = quote + if haskey($to_insert, $lhs) + throw( + DomainError( + $lhs, + "The variable occurs twice in the left-hand-side of the ODE system", + ), + ) + end + end + push!(eqs_expr, uniqueness_check_expr) + if isempty(loc_all_symb) + push!(eqs_expr, :($to_insert[$lhs] = $R($rhs))) + else + push!(eqs_expr, :($to_insert[$lhs] = ($rhs))) + end + end + + for n in all_symb_no_t + if !Base.isidentifier(n) + throw( + ArgumentError( + "The names of the variables will be injected into the global scope, so their name must be allowed Julia names, $n is not", + ), + ) + end + end + + logging_exprs = [ + :( + StructuralIdentifiability.Logging.with_logger( + StructuralIdentifiability._si_logger[], + ) do + @info "Summary of the model:" + @info "State variables: " * $(join(map(string, collect(x_vars)), ", ")) + @info "Parameters: " * $(join(map(string, collect(params)), ", ")) + @info "Inputs: " * $(join(map(string, collect(u_vars)), ", ")) + @info "Outputs: " * $(join(map(string, collect(y_vars)), ", ")) + end + ), + ] + # creating the ode/dds object + obj_type = Dict(:ode => :ODE, :dds => :DDS) + ds_expr = :(StructuralIdentifiability.$(obj_type[type]){ + StructuralIdentifiability.Nemo.QQMPolyRingElem, + }( + $vx, + $vy, + $x_dict, + $y_dict, + Array{StructuralIdentifiability.Nemo.QQMPolyRingElem}([$(u_vars...)]), + )) + + result = Expr( + :block, + logging_exprs..., + exp_ring, + assignments..., + x_var_expr, + y_var_expr, + x_dict_create_expr, + y_dict_create_expr, + eqs_expr..., + ds_expr, + ) + return result +end + +#------------------------------------------------------------------------------ + +""" + macro ODEmodel + +Macro for creating an ODE from a list of equations. +It also injects all variables into the global scope. + +## Example + +Creating a simple `ODE`: + +```jldoctest +using StructuralIdentifiability + +ode = @ODEmodel( + x1'(t) = a * x1(t) + u(t), + x2'(t) = b * x2(t) + c*x1(t)*x2(t), + y(t) = x1(t) +) +``` + +Here, +- `x1`, `x2` are state variables +- `y` is an output variable +- `u` is an input variable +- `a`, `b`, `c` are time-indepdendent parameters + +""" +macro ODEmodel(ex::Expr...) + return esc(generate_model_code(:ode, ex...)) +end + +#------------------------------------------------------------------------------ + +""" + macro DDSmodel + +Macro for creating a DDS (discrete dynamical system) +from a list of equations. +It also injects all variables into the global scope. + +## Example + +Creating a simple `DDS`: + +```jldoctest +using StructuralIdentifiability + +dds = @DDSmodel( + x1(t + 1) = a * x1(t) + u(t), + x2(t + 1) = b * x2(t) + c*x1(t)*x2(t), + y(t) = x1(t) +) +``` + +Here, +- `x1`, `x2` are state variables +- `y` is an output variable +- `u` is an input variable +- `a`, `b`, `c` are time-indepdendent parameters + +""" +macro DDSmodel(ex::Expr...) + return esc(generate_model_code(:dds, ex...)) +end + +#------------------------------------------------------------------------------ diff --git a/src/io_equation.jl b/src/io_equation.jl index ed1b1b1c4..92d90ff4e 100644 --- a/src/io_equation.jl +++ b/src/io_equation.jl @@ -2,7 +2,12 @@ const PROJECTION_VARNAME = "rand_proj_var" # ------------------------------------------------------------------------------ -function generator_var_change(generator, var::P, numer::P, denom::P) where {P <: MPolyElem} +function generator_var_change( + generator, + var::P, + numer::P, + denom::P, +) where {P <: MPolyRingElem} return IterTools.imap( point -> begin result = copy(point) @@ -15,7 +20,7 @@ end # ------------------------------------------------------------------------------ -function diff_poly(poly::P, derivation::Dict{P, T}) where {P <: MPolyElem, T} +function diff_poly(poly::P, derivation::Dict{P, T}) where {P <: MPolyRingElem, T} return sum(derivative(poly, x) * xd for (x, xd) in derivation) end @@ -28,7 +33,7 @@ end # ------------------------------------------------------------------------------ -function generate_io_equation_problem(ode::ODE{P}) where {P <: MPolyElem{<:FieldElem}} +function generate_io_equation_problem(ode::ODE{P}) where {P <: MPolyRingElem{<:FieldElem}} dim_x = length(ode.x_vars) # Creating a ring @@ -40,7 +45,7 @@ function generate_io_equation_problem(ode::ODE{P}) where {P <: MPolyElem{<:Field [var_to_str(u) * "_$i" for i in 0:dim_x for u in ode.u_vars], [PROJECTION_VARNAME], ) - ring, ring_vars = Nemo.PolynomialRing(base_ring(ode.poly_ring), var_names) + ring, ring_vars = Nemo.polynomial_ring(base_ring(ode.poly_ring), var_names) # Definiting a (partial) derivation on it derivation = Dict{P, P}() @@ -105,7 +110,7 @@ Output: ode::ODE{P}, auto_var_change::Bool, extra_projection = nothing, -) where {P <: MPolyElem{<:FieldElem}} +) where {P <: MPolyRingElem{<:FieldElem}} # Initialization ring, derivation, x_equations, y_equations, point_generator = generate_io_equation_problem(ode) @@ -341,7 +346,7 @@ Output: ode::ODE{P}; var_change_policy = :default, loglevel = Logging.Info, -) where {P <: MPolyElem{<:FieldElem}} +) where {P <: MPolyRingElem{<:FieldElem}} # Setting the var_change policy if (var_change_policy == :yes) || (var_change_policy == :default && length(ode.y_vars) < 3) diff --git a/src/lincomp.jl b/src/lincomp.jl index 72dd6a327..0603e39a9 100644 --- a/src/lincomp.jl +++ b/src/lincomp.jl @@ -45,14 +45,15 @@ function linear_compartment_model( push!(edges_vars_names, "a_0_$(s)") end - R, vars = StructuralIdentifiability.Nemo.PolynomialRing( + R, vars = StructuralIdentifiability.Nemo.polynomial_ring( StructuralIdentifiability.Nemo.QQ, vcat(x_vars_names, y_vars_names, u_vars_names, edges_vars_names), ) x_vars = @view vars[1:n] - x_equations = Dict{fmpq_mpoly, Union{fmpq_mpoly, Generic.Frac{fmpq_mpoly}}}( - x => R(0) for x in x_vars - ) + x_equations = + Dict{QQMPolyRingElem, Union{QQMPolyRingElem, Generic.Frac{QQMPolyRingElem}}}( + x => R(0) for x in x_vars + ) for i in 1:n for j in graph[i] rate = str_to_var("a_$(j)_$(i)", R) @@ -68,14 +69,15 @@ function linear_compartment_model( end end - y_equations = Dict{fmpq_mpoly, Union{fmpq_mpoly, Generic.Frac{fmpq_mpoly}}}( - str_to_var("y$i", R) => str_to_var("x$i", R) for i in outputs - ) + y_equations = + Dict{QQMPolyRingElem, Union{QQMPolyRingElem, Generic.Frac{QQMPolyRingElem}}}( + str_to_var("y$i", R) => str_to_var("x$i", R) for i in outputs + ) - return ODE{fmpq_mpoly}( + return ODE{QQMPolyRingElem}( x_equations, y_equations, - Array{fmpq_mpoly}([str_to_var("u$i", R) for i in inputs]), + Array{QQMPolyRingElem}([str_to_var("u$i", R) for i in inputs]), ) end diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index ae2a3cbba..a42a24098 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -30,7 +30,7 @@ function differentiate_solution( ic::Dict{P, T}, inputs::Dict{P, Array{T, 1}}, prec::Int, -) where {T <: Generic.FieldElem, P <: MPolyElem{T}} +) where {T <: Generic.FieldElem, P <: MPolyRingElem{T}} @debug "Computing the power series solution of the system" ps_sol = power_series_solution(ode, params, ic, inputs, prec) ps_ring = parent(first(values(ps_sol))) @@ -41,12 +41,12 @@ function differentiate_solution( @debug "Building the variational system at the solution" # Y' = AY + B vars = vcat(ode.x_vars, ode.parameters) - SA = AbstractAlgebra.MatrixSpace(ps_ring, length(ode.x_vars), length(ode.x_vars)) + SA = AbstractAlgebra.matrix_space(ps_ring, length(ode.x_vars), length(ode.x_vars)) A = SA([ eval_at_dict(derivative(ode.x_equations[vars[i]], vars[j]), ps_sol) for i in 1:length(ode.x_vars), j in 1:length(ode.x_vars) ]) - SB = AbstractAlgebra.MatrixSpace(ps_ring, length(ode.x_vars), length(vars)) + SB = AbstractAlgebra.matrix_space(ps_ring, length(ode.x_vars), length(vars)) B = zero(SB) for i in 1:length(ode.x_vars) for j in (length(ode.x_vars) + 1):length(vars) @@ -55,7 +55,7 @@ function differentiate_solution( end # TODO: make use of one() function (problems modulo prime) initial_condition = - zero(Nemo.MatrixSpace(base_ring(ode.poly_ring), length(ode.x_vars), length(vars))) + zero(Nemo.matrix_space(base_ring(ode.poly_ring), length(ode.x_vars), length(vars))) for i in 1:length(ode.x_vars) initial_condition[i, i] = 1 end @@ -86,7 +86,7 @@ function differentiate_output( ic::Dict{P, T}, inputs::Dict{P, Array{T, 1}}, prec::Int, -) where {T <: Generic.FieldElem, P <: MPolyElem{T}} +) where {T <: Generic.FieldElem, P <: MPolyRingElem{T}} @debug "Computing partial derivatives of the solution" ps_sol, sol_diff = differentiate_solution(ode, params, ic, inputs, prec) ps_ring = parent(first(values(ps_sol))) @@ -124,7 +124,7 @@ end for `f` being a polynomial/rational function over rationals (`QQ`) returns a tuple `(degree, max_coef_size)` """ -function get_degree_and_coeffsize(f::MPolyElem{Nemo.fmpq}) +function get_degree_and_coeffsize(f::MPolyRingElem{Nemo.QQFieldElem}) if length(f) == 0 return (0, 1) end @@ -135,114 +135,16 @@ function get_degree_and_coeffsize(f::MPolyElem{Nemo.fmpq}) return (total_degree(f), max_coef) end -function get_degree_and_coeffsize(f::Generic.Frac{<:MPolyElem{Nemo.fmpq}}) +function get_degree_and_coeffsize(f::Generic.Frac{<:MPolyRingElem{Nemo.QQFieldElem}}) num_deg, num_coef = get_degree_and_coeffsize(numerator(f)) den_deg, den_coef = get_degree_and_coeffsize(denominator(f)) return (max(num_deg, den_deg), max(num_coef, den_coef)) end -# ------------------------------------------------------------------------------ -""" - function assess_local_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=Array{}[], prob_threshold::Float64=0.99, type=:SE, loglevel=Logging.Info) - -Input: -- `ode` - the ODESystem object from ModelingToolkit -- `measured_quantities` - the measurable outputs of the model -- `funcs_to_check` - functions of parameters for which to check identifiability -- `prob_threshold` - probability of correctness -- `type` - identifiability type (`:SE` for single-experiment, `:ME` for multi-experiment) -- `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) - -Output: -- for `type=:SE`, the result is an (ordered) dictionary from each parameter to boolean; -- for `type=:ME`, the result is a tuple with the dictionary as in `:SE` case and array of number of experiments. - -The function determines local identifiability of parameters in `funcs_to_check` or all possible parameters if `funcs_to_check` is empty - -The result is correct with probability at least `prob_threshold`. - -`type` can be either `:SE` (single-experiment identifiability) or `:ME` (multi-experiment identifiability). -The return value is a tuple consisting of the array of bools and the number of experiments to be performed. -""" -function assess_local_identifiability( - ode::ModelingToolkit.ODESystem; - measured_quantities = Array{ModelingToolkit.Equation}[], - funcs_to_check = Array{}[], - prob_threshold::Float64 = 0.99, - type = :SE, - loglevel = Logging.Info, -) - restart_logging(loglevel = loglevel) - with_logger(_si_logger[]) do - return _assess_local_identifiability( - ode, - measured_quantities = measured_quantities, - funcs_to_check = funcs_to_check, - prob_threshold = prob_threshold, - type = type, - ) - end -end - -@timeit _to function _assess_local_identifiability( - ode::ModelingToolkit.ODESystem; - measured_quantities = Array{ModelingToolkit.Equation}[], - funcs_to_check = Array{}[], - prob_threshold::Float64 = 0.99, - type = :SE, -) - if length(measured_quantities) == 0 - 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." - measured_quantities = 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 - if length(funcs_to_check) == 0 - funcs_to_check = vcat( - [e for e in ModelingToolkit.states(ode) if !ModelingToolkit.isoutput(e)], - ModelingToolkit.parameters(ode), - ) - end - ode, conversion = mtk_to_si(ode, measured_quantities) - funcs_to_check_ = [eval_at_nemo(x, conversion) for x in funcs_to_check] - - if isequal(type, :SE) - result = _assess_local_identifiability( - ode, - funcs_to_check = funcs_to_check_, - prob_threshold = prob_threshold, - type = type, - ) - nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) - out_dict = - OrderedDict(nemo2mtk[param] => result[param] for param in funcs_to_check_) - return out_dict - elseif isequal(type, :ME) - result, bd = _assess_local_identifiability( - ode, - funcs_to_check = funcs_to_check_, - prob_threshold = prob_threshold, - type = type, - ) - nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) - out_dict = - OrderedDict(nemo2mtk[param] => result[param] for param in funcs_to_check_) - return (out_dict, bd) - end -end # ------------------------------------------------------------------------------ """ - assess_local_identifiability(ode::ODE{P}; funcs_to_check::Array{<: Any, 1}, prob_threshold::Float64=0.99, type=:SE, loglevel=Logging.Info) where P <: MPolyElem{Nemo.fmpq} + assess_local_identifiability(ode::ODE{P}; funcs_to_check::Array{<: Any, 1}, prob_threshold::Float64=0.99, type=:SE, loglevel=Logging.Info) where P <: MPolyRingElem{Nemo.QQFieldElem} Checks the local identifiability/observability of the functions in `funcs_to_check`. The result is correct with probability at least `prob_threshold`. @@ -258,7 +160,7 @@ function assess_local_identifiability( type = :SE, trbasis = nothing, loglevel = Logging.Info, -) where {P <: MPolyElem{Nemo.fmpq}} +) where {P <: MPolyRingElem{Nemo.QQFieldElem}} restart_logging(loglevel = loglevel) reset_timings() with_logger(_si_logger[]) do @@ -278,7 +180,7 @@ function _assess_local_identifiability( prob_threshold::Float64 = 0.99, type = :SE, trbasis = nothing, -) where {P <: MPolyElem{Nemo.fmpq}} +) where {P <: MPolyRingElem{Nemo.QQFieldElem}} if isempty(funcs_to_check) funcs_to_check = ode.parameters if type == :SE @@ -326,7 +228,7 @@ function _assess_local_identifiability( Dprime = max(Dprime, 1.0) prime = Primes.nextprime(Int(ceil(2 * mu * Dprime))) @debug "The prime is $prime" - F = Nemo.GF(prime) + F = Nemo.Native.GF(prime) @debug "Extending the model" ode_ext = @@ -345,13 +247,13 @@ function _assess_local_identifiability( num_exp = 0 # rows are the "parameters": parameters and initial conditions # columns are "observations": derivatives of the outputs - Jac = zero(Nemo.MatrixSpace(F, length(ode.parameters), 1)) + Jac = zero(Nemo.matrix_space(F, length(ode.parameters), 1)) output_derivatives = undef # the while loop is primarily for ME-deintifiability, it is adding replicas until the rank stabilizes # in the SE case, it will exit right away while true ic = Dict(x => F(rand(1:prime)) for x in ode_red.x_vars) - inputs = Dict{Nemo.gfp_mpoly, Array{Nemo.gfp_elem, 1}}( + inputs = Dict{Nemo.fpMPolyRingElem, Array{Nemo.fpFieldElem, 1}}( u => [F(rand(1:prime)) for i in 1:prec] for u in ode_red.u_vars ) @@ -359,10 +261,10 @@ function _assess_local_identifiability( output_derivatives = differentiate_output(ode_red, params_vals, ic, inputs, prec) @debug "Building the matrices" - newJac = vcat(Jac, zero(Nemo.MatrixSpace(F, length(ode.x_vars), ncols(Jac)))) + newJac = vcat(Jac, zero(Nemo.matrix_space(F, length(ode.x_vars), ncols(Jac)))) newJac = hcat( newJac, - zero(Nemo.MatrixSpace(F, nrows(newJac), prec * length(ode.y_vars))), + zero(Nemo.matrix_space(F, nrows(newJac), prec * length(ode.y_vars))), ) xs_params = vcat(ode_red.x_vars, ode_red.parameters) for (i, y) in enumerate(ode.y_vars) @@ -398,7 +300,7 @@ function _assess_local_identifiability( if !isnothing(trbasis) @debug "Transcendence basis computation requested" - reverted_Jac = zero(Nemo.MatrixSpace(F, size(Jac)[2], size(Jac)[1])) + reverted_Jac = zero(Nemo.matrix_space(F, size(Jac)[2], size(Jac)[1])) for i in 1:size(Jac)[1] for j in 1:size(Jac)[2] reverted_Jac[j, i] = Jac[size(Jac)[1] - i + 1, j] diff --git a/src/parametrizations.jl b/src/parametrizations.jl index a8b6e6ec4..b0d7f6fe3 100644 --- a/src/parametrizations.jl +++ b/src/parametrizations.jl @@ -67,7 +67,7 @@ function check_constructive_field_membership( # If norm_form(Num) // norm_form(Den) does not belongs to K(T), then # the fraction does not belong to the field if iszero(den_rem) - @warn """ + @debug """ The element $tagged_num // $tagged_den is not in the sub-field Normal form, numerator: $num_rem Normal form, denominator: $den_rem @@ -76,7 +76,7 @@ function check_constructive_field_membership( return false, zero(ring_of_tags) // one(ring_of_tags) end if total_degree(num_rem) > 0 || total_degree(den_rem) > 0 - @warn """ + @debug """ The element $tagged_num // $tagged_den is not in the sub-field Normal form, numerator: $num_rem Normal form, denominator: $den_rem @@ -174,7 +174,8 @@ $(join(map(x -> string(x[1]) * " -> " * string(x[2]), zip(fracs_gen, tag_string Saturation tag: $sat_string """ - poly_ring_tag, vars_tag = PolynomialRing(K, vcat(sat_string, orig_strings, tag_strings)) + poly_ring_tag, vars_tag = + polynomial_ring(K, vcat(sat_string, orig_strings, tag_strings)) sat_var = vars_tag[1] orig_vars = vars_tag[2:(nvars(poly_ring) + 1)] tag_vars = vars_tag[(nvars(poly_ring) + 2):end] @@ -227,7 +228,7 @@ $sat_string # # NOTE: reduction actually happens in K(T)[x]. So we map polynomials to the # parametric ring K(T)[x]. - ring_of_tags, tags = PolynomialRing(K, tag_strings) + ring_of_tags, tags = polynomial_ring(K, tag_strings) tag_to_gen = Dict(tags[i] => fracs_gen[i] for i in 1:length(fracs_gen)) if !isempty(intersect(tag_strings, orig_strings)) @warn """ @@ -236,7 +237,7 @@ $sat_string Original vars: $orig_strings""" end parametric_ring, _ = - PolynomialRing(FractionField(ring_of_tags), orig_strings, ordering = :degrevlex) + polynomial_ring(fraction_field(ring_of_tags), orig_strings, ordering = :degrevlex) relations_between_tags = map(poly -> parent_ring_change(poly, ring_of_tags), relations_between_tags) param_var_mapping = merge( @@ -380,7 +381,7 @@ function reparametrize_with_respect_to(ode, new_states, new_params) input_variable_names, output_variable_names, ) - ring_output, _ = PolynomialRing( + ring_output, _ = polynomial_ring( base_ring(ring_of_tags), all_variable_names, ordering = Nemo.ordering(ring_of_tags), @@ -492,7 +493,7 @@ function _reparametrize_global(ode::ODE{P}; prob_threshold = 0.99, seed = 42) wh ode_ring = parent(ode) @assert base_ring(parent(first(id_funcs))) == ode_ring @info "Constructing a new parametrization" - contains_states(poly::MPolyElem) = any(x -> degree(poly, x) > 0, ode.x_vars) + contains_states(poly::MPolyRingElem) = any(x -> degree(poly, x) > 0, ode.x_vars) contains_states(func) = contains_states(numerator(func)) || contains_states(denominator(func)) id_funcs_contains_states = filter(contains_states, id_funcs) diff --git a/src/pb_representation.jl b/src/pb_representation.jl index 3594b143a..34a618bd4 100644 --- a/src/pb_representation.jl +++ b/src/pb_representation.jl @@ -13,7 +13,7 @@ struct PBRepresentation u_names::Array{String} # variables with infinite orders in the profile param_names::Array{String} # scalar parameters profile::Dict{String, Int} # the profile restricted on the y-variables - projections::Dict{String, <:MPolyElem} + projections::Dict{String, <:MPolyRingElem} function PBRepresentation(ode::ODE, io_equations) if length(keys(io_equations)) > length(ode.y_vars) @@ -33,10 +33,10 @@ struct PBRepresentation decompose_derivative(v, vcat(y_names, u_names)) != nothing, map(var_to_str, gens(old_ring)), ) - newring, _ = Nemo.PolynomialRing(base_ring(old_ring), new_varnames) + newring, _ = Nemo.polynomial_ring(base_ring(old_ring), new_varnames) profile = Dict{String, Int}() - projections = Dict{String, MPolyElem}() + projections = Dict{String, MPolyRingElem}() for (y, eq) in io_equations (name, ord) = decompose_derivative(var_to_str(y), y_names) profile[name] = ord @@ -60,7 +60,7 @@ Among the variables `vars`, determines the leading derivative if the y-variable (if exists) with respect to the ordering defined by the PB-representation (see Remark 2.20 in https://arxiv.org/abs/2111.00991) """ -function find_leader(vars::Array{<:MPolyElem}, pbr::PBRepresentation) +function find_leader(vars::Array{<:MPolyRingElem}, pbr::PBRepresentation) y_ders = filter(v -> decompose_derivative(var_to_str(v), pbr.y_names) != nothing, vars) if length(y_ders) == 0 return nothing @@ -83,7 +83,7 @@ For a polynomial `poly` in the same differential variables as `pbr`, finds a polynomial ring sufficient for carrying out the reduction and the corresponding differentiation mapping on the variables """ -function common_ring(poly::MPolyElem, pbr::PBRepresentation) +function common_ring(poly::MPolyRingElem, pbr::PBRepresentation) max_ords = Dict{String, Int}(v => 0 for v in vcat(pbr.y_names, pbr.u_names)) new_params = Array{String, 1}() for v in vars(poly) @@ -118,8 +118,8 @@ function common_ring(poly::MPolyElem, pbr::PBRepresentation) append!(varnames, new_params) newring, _ = - StructuralIdentifiability.Nemo.PolynomialRing(base_ring(parent(poly)), varnames) - derivation = Dict{MPolyElem, MPolyElem}() + StructuralIdentifiability.Nemo.polynomial_ring(base_ring(parent(poly)), varnames) + derivation = Dict{MPolyRingElem, MPolyRingElem}() for v in varnames d = decompose_derivative(v, vcat(pbr.y_names, pbr.u_names)) if d == nothing @@ -143,7 +143,7 @@ end Computes the leading coefficient of `f` viewed as a univariate polynomiall in variable `x` """ -function lc_univariate(f::MPolyElem, x::MPolyElem) +function lc_univariate(f::MPolyRingElem, x::MPolyRingElem) FieldType = typeof(one(base_ring(parent(f)))) dict_result = Dict{Array{Int, 1}, FieldType}() x_ind = findfirst(v -> v == x, gens(parent(f))) @@ -174,7 +174,7 @@ Input: Output: the pseudoremainder of `f` divided by `g` w.r.t. `x` """ -function pseudodivision(f::MPolyElem, g::MPolyElem, x::MPolyElem) +function pseudodivision(f::MPolyRingElem, g::MPolyRingElem, x::MPolyRingElem) result = f lcg = lc_univariate(g, x) while Nemo.degree(result, x) >= Nemo.degree(g, x) @@ -190,7 +190,7 @@ end # ----------------------------------------------------------------------------- -function diff(p::MPolyElem, derivation::Dict{<:MPolyElem, <:MPolyElem}, i::Int) +function diff(p::MPolyRingElem, derivation::Dict{<:MPolyRingElem, <:MPolyRingElem}, i::Int) if i == 0 return p end @@ -213,7 +213,7 @@ Input: Output: the result of differential reduction of `diffpoly` by `pbr` considered as a characteristic set (see Remark 2.20 in the paper) """ -function diffreduce(diffpoly::MPolyElem, pbr::PBRepresentation) +function diffreduce(diffpoly::MPolyRingElem, pbr::PBRepresentation) (ring, der) = common_ring(diffpoly, pbr) result = parent_ring_change(diffpoly, ring) diff --git a/src/power_series_utils.jl b/src/power_series_utils.jl index 3124d47c0..ac4429243 100644 --- a/src/power_series_utils.jl +++ b/src/power_series_utils.jl @@ -1,18 +1,18 @@ #------------------------------------------------------------------------------ -function truncate_matrix(M::MatElem{<:Generic.AbsSeriesElem}, prec::Int) +function truncate_matrix(M::MatElem{<:AbsPowerSeriesRingElem}, prec::Int) return map(e -> truncate(e, prec), M) end #------------------------------------------------------------------------------ -function matrix_set_precision!(M::MatElem{<:Generic.AbsSeriesElem}, prec::Int) +function matrix_set_precision!(M::MatElem{<:AbsPowerSeriesRingElem}, prec::Int) map(e -> set_precision!(e, prec), M) end #------------------------------------------------------------------------------ -function ps_matrix_const_term(M::MatElem{<:Generic.AbsSeriesElem}) +function ps_matrix_const_term(M::MatElem{<:AbsPowerSeriesRingElem}) return map(e -> coeff(e, 0), M) end @@ -26,7 +26,7 @@ Performs a single step of Newton iteration for inverting `M` with `Minv` being a function _matrix_inv_newton_iteration( M::MatElem{T}, Minv::MatElem{T}, -) where {T <: Generic.AbsSeriesElem{<:Generic.FieldElem}} +) where {T <: AbsPowerSeriesRingElem{<:Generic.FieldElem}} return 2 * Minv - Minv * M * Minv end @@ -44,7 +44,7 @@ Output: - the inverse of `M` computed up to `prec` """ function ps_matrix_inv( - M::MatElem{<:Generic.AbsSeriesElem{<:Generic.FieldElem}}, + M::MatElem{<:AbsPowerSeriesRingElem{<:Generic.FieldElem}}, prec::Int = -1, ) const_term = ps_matrix_const_term(M) @@ -70,7 +70,7 @@ Input: Output: - the derivative of `ps` """ -function ps_diff(ps::Generic.AbsSeriesElem{<:Generic.RingElem}) +function ps_diff(ps::AbsPowerSeriesRingElem{<:Generic.RingElem}) result = zero(parent(ps)) set_precision!(result, precision(ps)) for exp in 1:(precision(ps) - 1) @@ -89,7 +89,7 @@ Input: Output: - the integral of `ps` without constant term """ -function ps_integrate(ps::Generic.AbsSeriesElem{<:Generic.FieldElem}) +function ps_integrate(ps::AbsPowerSeriesRingElem{<:Generic.FieldElem}) result = zero(parent(ps)) set_precision!(result, precision(ps) + 1) for exp in 0:(precision(ps) - 1) @@ -109,7 +109,7 @@ Input: Output: - the natural log of `M` """ -function ps_matrix_log(M::MatElem{<:Generic.AbsSeriesElem{<:Generic.FieldElem}}) +function ps_matrix_log(M::MatElem{<:AbsPowerSeriesRingElem{<:Generic.FieldElem}}) const_term = ps_matrix_const_term(M) if const_term != one(parent(const_term)) throw(Base.DomainError("Constant term must be the identity matrix")) @@ -123,9 +123,9 @@ end #------------------------------------------------------------------------------ function _matrix_homlinear_de_newton_iteration( - A::MatElem{<:Generic.AbsSeriesElem{T}}, - Y::MatElem{<:Generic.AbsSeriesElem{T}}, - Z::MatElem{<:Generic.AbsSeriesElem{T}}, + A::MatElem{<:AbsPowerSeriesRingElem{T}}, + Y::MatElem{<:AbsPowerSeriesRingElem{T}}, + Z::MatElem{<:AbsPowerSeriesRingElem{T}}, cur_prec::Int, ) where {T <: Generic.FieldElem} Yprime = map(ps_diff, Y) @@ -150,7 +150,7 @@ Output: - matrix `Y` such that `Y' = AY` up to precision of `A - 1` and `Y(0) = Y0` """ function ps_matrix_homlinear_de( - A::MatElem{<:Generic.AbsSeriesElem{T}}, + A::MatElem{<:AbsPowerSeriesRingElem{T}}, Y0::MatElem{<:T}, prec::Int = -1, ) where {T <: Generic.FieldElem} @@ -173,10 +173,10 @@ end #------------------------------------------------------------------------------ function _variation_of_constants( - A::MatElem{<:Generic.AbsSeriesElem{T}}, - B::MatElem{<:Generic.AbsSeriesElem{T}}, - Yh::MatElem{<:Generic.AbsSeriesElem{T}}, - Zh::MatElem{<:Generic.AbsSeriesElem{T}}, + A::MatElem{<:AbsPowerSeriesRingElem{T}}, + B::MatElem{<:AbsPowerSeriesRingElem{T}}, + Yh::MatElem{<:AbsPowerSeriesRingElem{T}}, + Zh::MatElem{<:AbsPowerSeriesRingElem{T}}, Y0::MatElem{<:T}, prec::Int, ) where {T <: Generic.FieldElem} @@ -198,14 +198,14 @@ Output: - matrix `Y` such that `Y' = AY + B` up to precision of `A - 1` and `Y(0) = Y0` """ function ps_matrix_linear_de( - A::MatElem{<:Generic.AbsSeriesElem{T}}, - B::MatElem{<:Generic.AbsSeriesElem{T}}, + A::MatElem{<:AbsPowerSeriesRingElem{T}}, + B::MatElem{<:AbsPowerSeriesRingElem{T}}, Y0::MatElem{<:T}, prec::Int = -1, ) where {T <: Generic.FieldElem} prec = (prec == -1) ? precision(A[1, 1]) : prec n = nrows(A) - identity = one(AbstractAlgebra.MatrixSpace(base_ring(parent(Y0)), n, n)) + identity = one(AbstractAlgebra.matrix_space(base_ring(parent(Y0)), n, n)) Yh, Zh = ps_matrix_homlinear_de(A, identity, prec) matrix_set_precision!(Zh, prec) return _variation_of_constants(A, B, Yh, Zh, Y0, prec) @@ -231,12 +231,12 @@ function ps_ode_solution( ic::Dict{P, T}, inputs::Dict{P, Array{T, 1}}, prec::Int, -) where {T <: Generic.FieldElem, P <: MPolyElem{T}} +) where {T <: Generic.FieldElem, P <: MPolyRingElem{T}} n = length(equations) ring = parent(equations[1]) - S = AbstractAlgebra.MatrixSpace(ring, n, n) - Sv = AbstractAlgebra.MatrixSpace(ring, n, 1) - Svconst = AbstractAlgebra.MatrixSpace(base_ring(ring), n, 1) + S = AbstractAlgebra.matrix_space(ring, n, n) + Sv = AbstractAlgebra.matrix_space(ring, n, 1) + Svconst = AbstractAlgebra.matrix_space(base_ring(ring), n, 1) eqs = Sv(equations) x_vars = filter(v -> ("$(v)_dot" in map(string, gens(ring))), gens(ring)) @@ -246,7 +246,7 @@ function ps_ode_solution( Jac_dots = S([derivative(p, xd) for p in equations, xd in x_dot_vars]) Jac_xs = S([derivative(p, x) for p in equations, x in x_vars]) - ps_ring, t = PowerSeriesRing(base_ring(ring), prec, "t"; model = :capped_absolute) + ps_ring, t = power_series_ring(base_ring(ring), prec, "t"; model = :capped_absolute) solution = Dict() for (u, coeffs) in inputs solution[u] = sum(coeffs[i] * t^(i - 1) for i in 1:length(coeffs)) @@ -292,7 +292,7 @@ function ps_ode_solution( ic::Dict{P, Int}, inputs::Dict{P, Array{Int, 1}}, prec::Int, -) where {P <: MPolyElem{<:Generic.FieldElem}} +) where {P <: MPolyRingElem{<:Generic.FieldElem}} bring = base_ring(parent(equations[1])) ps_ode_solution( equations, diff --git a/src/precompile.jl b/src/precompile.jl index fb2c2094f..8930129e9 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -3,12 +3,6 @@ # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the # precompile file and potentially make loading faster. using Logging - using ModelingToolkit - @parameters a01 a21 a12 - @variables t x0(t) x1(t) y(t) - D = Differential(t) - eqs = [D(x0) ~ -(a01 + a21) * x0 + a12 * x1, D(x1) ~ a21 * x0 - a12 * x1] - de = ODESystem(eqs, t, name = :Test) @compile_workload begin restart_logging(loglevel = Logging.Warn) with_logger(_si_logger[]) do @@ -21,12 +15,6 @@ y(t) = x2(t) ) assess_identifiability(ode, loglevel = Logging.Warn) - assess_identifiability(de; measured_quantities = [x0], loglevel = Logging.Warn) - assess_identifiability( - de; - measured_quantities = [y ~ x0], - loglevel = Logging.Warn, - ) find_identifiable_functions(ode, with_states = true, loglevel = Logging.Warn) end restart_logging(loglevel = Logging.Info) diff --git a/src/primality_check.jl b/src/primality_check.jl index ae98f006b..acccf0ebe 100644 --- a/src/primality_check.jl +++ b/src/primality_check.jl @@ -1,10 +1,10 @@ # ------------------------------------------------------------------------------ -function check_primality_zerodim(J::Array{fmpq_mpoly, 1}) +function check_primality_zerodim(J::Array{QQMPolyRingElem, 1}) J = Groebner.groebner(J, loglevel = _groebner_loglevel[]) basis = Groebner.kbase(J, loglevel = _groebner_loglevel[]) dim = length(basis) - S = Nemo.MatrixSpace(Nemo.QQ, dim, dim) + S = Nemo.matrix_space(Nemo.QQ, dim, dim) matrices = [] @debug "$J $basis" @debug "Dim is $dim" @@ -22,15 +22,15 @@ function check_primality_zerodim(J::Array{fmpq_mpoly, 1}) generic_multiplication = sum(Nemo.QQ(rand(1:100)) * M for M in matrices) @debug generic_multiplication - R, t = Nemo.PolynomialRing(Nemo.QQ, "t") + R, t = Nemo.polynomial_ring(Nemo.QQ, "t") @debug "$(Nemo.charpoly(R, generic_multiplication))" - return Nemo.isirreducible(Nemo.charpoly(R, generic_multiplication)) + return Nemo.is_irreducible(Nemo.charpoly(R, generic_multiplication)) end #------------------------------------------------------------------------------ """ - check_primality(polys::Dict{fmpq_mpoly, fmpq_mpoly}, extra_relations::Array{fmpq_mpoly, 1}) + check_primality(polys::Dict{QQMPolyRingElem, QQMPolyRingElem}, extra_relations::Array{QQMPolyRingElem, 1}) The function checks if the ideal generated by the polynomials and saturated at the leading coefficient with respect to the corresponding variables is prime @@ -39,13 +39,13 @@ over rationals. The `extra_relations` allows adding more polynomials to the generators (not affecting the saturation). """ function check_primality( - polys::Dict{fmpq_mpoly, fmpq_mpoly}, - extra_relations::Array{fmpq_mpoly, 1}, + polys::Dict{QQMPolyRingElem, QQMPolyRingElem}, + extra_relations::Array{QQMPolyRingElem, 1}, ) leaders = collect(keys(polys)) ring = parent(leaders[1]) - Rspec, vspec = Nemo.PolynomialRing(Nemo.QQ, [var_to_str(l) for l in leaders]) + Rspec, vspec = Nemo.polynomial_ring(Nemo.QQ, [var_to_str(l) for l in leaders]) eval_point = [v in keys(polys) ? v : ring(rand(1:100)) for v in gens(ring)] all_polys = vcat(collect(values(polys)), extra_relations) zerodim_ideal = @@ -56,14 +56,14 @@ end #------------------------------------------------------------------------------ """ - check_primality(polys::Dict{fmpq_mpoly, fmpq_mpoly}) + check_primality(polys::Dict{QQMPolyRingElem, QQMPolyRingElem}) The function checks if the ideal generated by the polynomials and saturated at the leading coefficient with respect to the corresponding variables is prime over rationals. """ -function check_primality(polys::Dict{fmpq_mpoly, fmpq_mpoly}) - return check_primality(polys, Array{fmpq_mpoly, 1}()) +function check_primality(polys::Dict{QQMPolyRingElem, QQMPolyRingElem}) + return check_primality(polys, Array{QQMPolyRingElem, 1}()) end # ------------------------------------------------------------------------------ diff --git a/src/states.jl b/src/states.jl index c706a9ab5..86e0269d8 100644 --- a/src/states.jl +++ b/src/states.jl @@ -15,7 +15,7 @@ Output: function extract_coefficients_ratfunc( f::AbstractAlgebra.Generic.Frac{<:P}, vars::Vector{<:P}, -) where {P <: MPolyElem{<:FieldElem}} +) where {P <: MPolyRingElem{<:FieldElem}} num, denom = unpack_fraction(f) total_coeffs = Vector{P}() for p in (num, denom) @@ -32,7 +32,7 @@ end function extract_coefficients_ratfunc( f::P, vars::Vector{<:P}, -) where {P <: MPolyElem{<:FieldElem}} +) where {P <: MPolyRingElem{<:FieldElem}} return extract_coefficients_ratfunc(f // 1, vars) end @@ -50,7 +50,7 @@ Output: function lie_derivative( f::Generic.Frac{<:P}, ode::ODE{<:P}, -) where {P <: MPolyElem{<:FieldElem}} +) where {P <: MPolyRingElem{<:FieldElem}} @assert all([(x in ode.parameters) || (x in ode.x_vars) for x in vars(f)]) res = zero(parent(ode)) // 1 for (x, eq) in ode.x_equations @@ -59,7 +59,7 @@ function lie_derivative( return res end -function lie_derivative(f::P, ode::ODE{<:P}) where {P <: MPolyElem{<:FieldElem}} +function lie_derivative(f::P, ode::ODE{<:P}) where {P <: MPolyRingElem{<:FieldElem}} return lie_derivative(f // 1, ode) end @@ -79,7 +79,7 @@ identifiable functions of parameters only @timeit _to function states_generators( ode::ODE{P}, io_equations::Dict{P, P}, -) where {P <: MPolyElem{<:FieldElem}} +) where {P <: MPolyRingElem{<:FieldElem}} y_to_ord = Dict{P, Int}() ynames = [var_to_str(y) for y in ode.y_vars] for (leader, ioeq) in io_equations diff --git a/src/submodels.jl b/src/submodels.jl index 107ceac4e..793e77ba4 100644 --- a/src/submodels.jl +++ b/src/submodels.jl @@ -12,15 +12,15 @@ Output: - Dictionary where each key is a variable and each value is a list of variables on which the key depends """ -function construct_graph(ode::ODE{P}) where {P <: MPolyElem} - graph = Dict{fmpq_mpoly, Set{fmpq_mpoly}}() +function construct_graph(ode::ODE{P}) where {P <: MPolyRingElem} + graph = Dict{QQMPolyRingElem, Set{QQMPolyRingElem}}() for (x, f) in ode.x_equations temp = unpack_fraction(f) - graph[x] = Set{fmpq_mpoly}(union(vars(temp[1]), vars(temp[2]))) + graph[x] = Set{QQMPolyRingElem}(union(vars(temp[1]), vars(temp[2]))) end for (y, f) in ode.y_equations temp = unpack_fraction(f) - graph[y] = Set{fmpq_mpoly}(union(vars(temp[1]), vars(temp[2]))) + graph[y] = Set{QQMPolyRingElem}(union(vars(temp[1]), vars(temp[2]))) end return graph @@ -29,9 +29,9 @@ end # ------------------------------------------------------------------------------ function dfs( - graph::Dict{fmpq_mpoly, Set{fmpq_mpoly}}, - start::fmpq_mpoly, - visited::Set{fmpq_mpoly}, + graph::Dict{QQMPolyRingElem, Set{QQMPolyRingElem}}, + start::QQMPolyRingElem, + visited::Set{QQMPolyRingElem}, ) push!(visited, start) if start in keys(graph) @@ -47,12 +47,12 @@ end # ------------------------------------------------------------------------------ function traverse_outputs( - graph::Dict{fmpq_mpoly, Set{fmpq_mpoly}}, - ys::Array{fmpq_mpoly, 1}, + graph::Dict{QQMPolyRingElem, Set{QQMPolyRingElem}}, + ys::Array{QQMPolyRingElem, 1}, ) - raw_models = Dict{fmpq_mpoly, Set{fmpq_mpoly}}() + raw_models = Dict{QQMPolyRingElem, Set{QQMPolyRingElem}}() for y in ys - model = dfs(graph, y, Set{fmpq_mpoly}()) + model = dfs(graph, y, Set{QQMPolyRingElem}()) raw_models[y] = model end return raw_models @@ -61,10 +61,10 @@ end # ------------------------------------------------------------------------------ function saturate_ys( - unions::Array{Set{fmpq_mpoly}, 1}, - Y::Array{fmpq_mpoly, 1}, - graph::Dict{fmpq_mpoly, Set{fmpq_mpoly}}, - X::Array{fmpq_mpoly, 1}, + unions::Array{Set{QQMPolyRingElem}, 1}, + Y::Array{QQMPolyRingElem, 1}, + graph::Dict{QQMPolyRingElem, Set{QQMPolyRingElem}}, + X::Array{QQMPolyRingElem, 1}, ) for element in unions for y in Y @@ -78,8 +78,8 @@ end # ------------------------------------------------------------------------------ -function search_add_unions(submodels::Array{Set{fmpq_mpoly}, 1}) - result = Array{Set{fmpq_mpoly}, 1}([Set{fmpq_mpoly}()]) +function search_add_unions(submodels::Array{Set{QQMPolyRingElem}, 1}) + result = Array{Set{QQMPolyRingElem}, 1}([Set{QQMPolyRingElem}()]) for model in submodels for index in 1:length(result) push!(result, union(result[index], model)) @@ -93,10 +93,10 @@ end # filters the models containin all states or no states function filter_max_empty( ode::ODE{P}, - submodels::Array{Set{fmpq_mpoly}, 1}, -) where {P <: MPolyElem} + submodels::Array{Set{QQMPolyRingElem}, 1}, +) where {P <: MPolyRingElem} n = length(ode.x_vars) - new_sub = Array{Set{fmpq_mpoly}, 1}([]) + new_sub = Array{Set{QQMPolyRingElem}, 1}([]) for submodel in submodels list_x = [x for x in submodel if x in ode.x_vars] if !(length(list_x) == n) && (length(list_x) > 0) @@ -108,10 +108,10 @@ end # ------------------------------------------------------------------------------ -function ode_aux(ode::ODE{P}, submodel::Set{fmpq_mpoly}) where {P <: MPolyElem} +function ode_aux(ode::ODE{P}, submodel::Set{QQMPolyRingElem}) where {P <: MPolyRingElem} new_y = copy(ode.y_equations) new_x = copy(ode.x_equations) - new_u = Array{fmpq_mpoly, 1}([u for u in ode.u_vars if u in submodel]) + new_u = Array{QQMPolyRingElem, 1}([u for u in ode.u_vars if u in submodel]) for (x, f) in ode.x_equations if !(issubset(vars(x), submodel) && issubset(vars(f), submodel)) delete!(new_x, x) @@ -125,8 +125,8 @@ function ode_aux(ode::ODE{P}, submodel::Set{fmpq_mpoly}) where {P <: MPolyElem} end sub_str = map(var_to_str, collect(submodel)) - S, _ = Nemo.PolynomialRing(Nemo.QQ, sub_str) - dict_type = Dict{fmpq_mpoly, Union{fmpq_mpoly, Generic.Frac{fmpq_mpoly}}} + S, _ = Nemo.polynomial_ring(Nemo.QQ, sub_str) + dict_type = Dict{QQMPolyRingElem, Union{QQMPolyRingElem, Generic.Frac{QQMPolyRingElem}}} fin_x = dict_type(parent_ring_change(x, S) => parent_ring_change(f, S) for (x, f) in new_x) fin_y = @@ -142,7 +142,7 @@ function ode_aux(ode::ODE{P}, submodel::Set{fmpq_mpoly}) where {P <: MPolyElem} y in ode.y_vars if var_to_str(y) in map(var_to_str, collect(keys(fin_y))) ] - return ODE{fmpq_mpoly}(new_x_vars, new_y_vars, fin_x, fin_y, fin_u) + return ODE{QQMPolyRingElem}(new_x_vars, new_y_vars, fin_x, fin_y, fin_u) end # ------------------------------------------------------------------------------ @@ -163,8 +163,8 @@ Output: function submodel2ode( ode::ODE{P}, - submodels::Array{Set{fmpq_mpoly}, 1}, -) where {P <: MPolyElem} + submodels::Array{Set{QQMPolyRingElem}, 1}, +) where {P <: MPolyRingElem} return [ode_aux(ode, submodel) for submodel in submodels] end @@ -187,14 +187,14 @@ Example: y1(t) = x1(t), y2(t) = x2(t)) >find_submodels(ode) - ODE{fmpq_mpoly}[ + ODE{QQMPolyRingElem}[ x1'(t) = a(t)*x2(t)^2 + x1(t) y1(t) = x1(t) ] ``` """ -function find_submodels(ode::ODE{P}) where {P <: MPolyElem} +function find_submodels(ode::ODE{P}) where {P <: MPolyRingElem} graph = construct_graph(ode) ys = ode.y_vars xs = ode.x_vars diff --git a/src/util.jl b/src/util.jl index a55d5571c..919c2bb69 100644 --- a/src/util.jl +++ b/src/util.jl @@ -1,10 +1,10 @@ # ------------------------------------------------------------------------------ -function Nemo.vars(f::Generic.Frac{<:MPolyElem}) +function Nemo.vars(f::Generic.Frac{<:MPolyRingElem}) return collect(union(Set(vars(numerator(f))), Set(vars(denominator(f))))) end -function Nemo.total_degree(f::Generic.Frac{<:MPolyElem}) +function Nemo.total_degree(f::Generic.Frac{<:MPolyRingElem}) return sum(map(total_degree, unpack_fraction(f))) end @@ -36,13 +36,13 @@ end Evaluates a polynomial/rational function on a dictionary of type `var => val` and missing values are replaced with zeroes """ -function eval_at_dict(poly::P, d::Dict{P, <:RingElem}) where {P <: MPolyElem} +function eval_at_dict(poly::P, d::Dict{P, <:RingElem}) where {P <: MPolyRingElem} R = parent(first(values(d))) point = [get(d, v, zero(R)) for v in gens(parent(poly))] return evaluate(poly, point) end -function eval_at_dict(poly::P, d::Dict{P, S}) where {P <: MPolyElem, S <: Real} +function eval_at_dict(poly::P, d::Dict{P, S}) where {P <: MPolyRingElem, S <: Real} R = parent(poly) @assert R == parent(first(keys(d))) xs = gens(parent(first(keys(d)))) @@ -59,7 +59,10 @@ function eval_at_dict(poly::P, d::Dict{P, S}) where {P <: MPolyElem, S <: Real} return accum end -function eval_at_dict(rational::Generic.Frac{T}, d::Dict{T, V}) where {T <: MPolyElem, V} +function eval_at_dict( + rational::Generic.Frac{T}, + d::Dict{T, V}, +) where {T <: MPolyRingElem, V} f, g = unpack_fraction(rational) return eval_at_dict(f, d) / eval_at_dict(g, d) end @@ -67,7 +70,7 @@ end function eval_at_dict( rational::Generic.Frac{<:T}, d::Dict{T, <:RingElem}, -) where {T <: MPolyElem} +) where {T <: MPolyRingElem} f, g = unpack_fraction(rational) return eval_at_dict(f, d) * inv(eval_at_dict(g, d)) end @@ -75,24 +78,24 @@ end function eval_at_dict( rational::Generic.Frac{<:P}, d::Dict{<:P, <:Union{<:Generic.Frac, <:P}}, -) where {P <: MPolyElem} +) where {P <: MPolyRingElem} f, g = unpack_fraction(rational) return eval_at_dict(f, d) // eval_at_dict(g, d) end # ------------------------------------------------------------------------------ -function unpack_fraction(f::MPolyElem) +function unpack_fraction(f::MPolyRingElem) return (f, one(parent(f))) end -function unpack_fraction(f::Generic.Frac{<:MPolyElem}) +function unpack_fraction(f::Generic.Frac{<:MPolyRingElem}) return (numerator(f), denominator(f)) end # ------------------------------------------------------------------------------ -function simplify_frac(numer::P, denom::P) where {P <: MPolyElem} +function simplify_frac(numer::P, denom::P) where {P <: MPolyRingElem} gcd_sub = gcd(numer, denom) sub_numer = divexact(numer, gcd_sub) sub_denom = divexact(denom, gcd_sub) @@ -163,7 +166,7 @@ function make_substitution( var_sub::P, val_numer::P, val_denom::P, -) where {P <: MPolyElem} +) where {P <: MPolyRingElem} d = Nemo.degree(f, var_sub) result = 0 @@ -180,7 +183,7 @@ end function homogenize(fs) ring = parent(fs[1]) - newring, hom_vars = PolynomialRing( + newring, hom_vars = polynomial_ring( base_ring(ring), vcat("X0", map(string, gens(ring))), ordering = ordering(ring), @@ -202,7 +205,7 @@ end function dehomogenize(Fs) ring = parent(Fs[1]) - newring, dehom_vars = PolynomialRing( + newring, dehom_vars = polynomial_ring( base_ring(ring), map(string, gens(ring)[2:end]), ordering = ordering(ring), @@ -229,7 +232,7 @@ Output: - a polynomial in `new_ring` “equal” to `poly` """ function parent_ring_change( - poly::MPolyElem, + poly::MPolyRingElem, new_ring::MPolyRing; matching = :byname, shift = 0, @@ -279,7 +282,7 @@ function parent_ring_change( end function parent_ring_change( - f::Generic.Frac{<:MPolyElem}, + f::Generic.Frac{<:MPolyRingElem}, new_ring::MPolyRing; matching = :byname, ) @@ -301,7 +304,7 @@ Output: - `div`'s are divisors of `f` such that `f` is their product with certain powers - if `certainty` is true, `div` is ``Q``-irreducible """ -function uncertain_factorization(f::MPolyElem{fmpq}) +function uncertain_factorization(f::MPolyRingElem{QQFieldElem}) vars_f = vars(f) if isempty(vars_f) return Array{Tuple{typeof(f), Bool}, 1}() @@ -319,7 +322,7 @@ function uncertain_factorization(f::MPolyElem{fmpq}) while true plugin = rand(-3:3, length(gens(parent(f)))) if evaluate(mainvar_coeffs[end], plugin) != 0 - uni_ring, T = Nemo.PolynomialRing(base_ring(f), "T") + uni_ring, T = Nemo.polynomial_ring(base_ring(f), "T") f_uni = sum([evaluate(mainvar_coeffs[i + 1], plugin) * T^i for i in 0:d]) if !isone(gcd(f_uni, derivative(f_uni))) factor_out = gcd(f, derivative(f, main_var)) @@ -330,7 +333,7 @@ function uncertain_factorization(f::MPolyElem{fmpq}) f = divexact(f, factor_out) f_uni = divexact(f_uni, gcd(f_uni, derivative(f_uni))) end - is_irr = Nemo.isirreducible(f_uni) + is_irr = Nemo.is_irreducible(f_uni) break end end @@ -342,7 +345,7 @@ end # ------------------------------------------------------------------------------ -function fast_factor(poly::MPolyElem{fmpq}) +function fast_factor(poly::MPolyRingElem{QQFieldElem}) prelim_factors = uncertain_factorization(poly) cert_factors = map(pair -> pair[1], filter(f -> f[2], prelim_factors)) uncert_factors = map(pair -> pair[1], filter(f -> !f[2], prelim_factors)) @@ -375,7 +378,7 @@ Input: Output: - dictionary with keys being tuples of length `lenght(variables)` and values being polynomials in the variables other than those which are the coefficients at the corresponding monomials (in a smaller polynomial ring) """ -function extract_coefficients(poly::P, variables::Array{P, 1}) where {P <: MPolyElem} +function extract_coefficients(poly::P, variables::Array{P, 1}) where {P <: MPolyRingElem} xs = gens(parent(poly)) @assert all(in(xs), variables) cut_indices = map(v -> findfirst(x -> x == v, xs), variables) @@ -383,7 +386,7 @@ function extract_coefficients(poly::P, variables::Array{P, 1}) where {P <: MPoly coeff_vars = xs[coeff_indices] K = base_ring(parent(poly)) - new_ring, _ = Nemo.PolynomialRing(K, map(vv -> var_to_str(vv, xs = xs), coeff_vars)) + new_ring, _ = Nemo.polynomial_ring(K, map(vv -> var_to_str(vv, xs = xs), coeff_vars)) FieldType = elem_type(K) result = Dict{Vector{Int}, Tuple{Vector{Vector{Int}}, Vector{FieldType}}}() @@ -422,7 +425,7 @@ end # ------------------------------------------------------------------------------ -function var_to_str(v::MPolyElem; xs = gens(parent(v))) +function var_to_str(v::MPolyRingElem; xs = gens(parent(v))) ind = findfirst(vv -> vv == v, xs) return string(symbols(parent(v))[ind]) end @@ -470,7 +473,7 @@ end For a variable `v`, returns a variable in `ring` with the same name """ -function switch_ring(v::MPolyElem, ring::MPolyRing) +function switch_ring(v::MPolyRingElem, ring::MPolyRing) ind = findfirst(vv -> vv == v, gens(parent(v))) return str_to_var(string(symbols(parent(v))[ind]), ring) @@ -478,55 +481,6 @@ end # ------------------------------------------------------------------------------ -function eval_at_nemo(e::Num, vals::Dict) - e = Symbolics.value(e) - return eval_at_nemo(e, vals) -end - -function eval_at_nemo(e::SymbolicUtils.BasicSymbolic, vals::Dict) - if Symbolics.istree(e) - # checking if it is a function of the form x(t), a bit dirty - if length(Symbolics.arguments(e)) == 1 && "$(first(Symbolics.arguments(e)))" == "t" - return vals[e] - end - # checking if this is a vector entry like x(t)[1] - if Symbolics.operation(e) == getindex - return vals[e] - end - # otherwise, this is a term - args = map(a -> eval_at_nemo(a, vals), Symbolics.arguments(e)) - if Symbolics.operation(e) in (+, -, *) - return Symbolics.operation(e)(args...) - elseif isequal(Symbolics.operation(e), /) - return //(args...) - elseif isequal(Symbolics.operation(e), ^) - if args[2] >= 0 - return args[1]^args[2] - end - return 1 // args[1]^(-args[2]) - end - throw(Base.ArgumentError("Function $(Symbolics.operation(e)) is not supported")) - elseif e isa Symbolics.Symbolic - return get(vals, e, e) - end -end - -function eval_at_nemo(e::Union{Integer, Rational}, vals::Dict) - return e -end - -function eval_at_nemo(e::Union{Float16, Float32, Float64}, vals::Dict) - if isequal(e % 1, 0) - out = Int(e) - else - out = rationalize(e) - end - @warn "Floating point value $e will be converted to $(out)." - return out -end - -# ----------------------------------------------------------------------------- - """ decompose_derivative(varname, prefixes) @@ -553,7 +507,7 @@ Finds the order of a differential polynomial `diffpoly` in a variable `prefix`, returns -1 is the variable does not appear """ -function difforder(diffpoly::MPolyElem, prefix::String) +function difforder(diffpoly::MPolyRingElem, prefix::String) orders = [-1] for v in vars(diffpoly) d = decompose_derivative(var_to_str(v), [prefix]) @@ -563,19 +517,3 @@ function difforder(diffpoly::MPolyElem, prefix::String) end 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 diff --git a/src/wronskian.jl b/src/wronskian.jl index ac5750a12..cca62bc5a 100644 --- a/src/wronskian.jl +++ b/src/wronskian.jl @@ -20,7 +20,7 @@ function monomial_compress(io_equation, ode::ODE) return monomial_compress(io_equation, ode.parameters) end -function monomial_compress(io_equation, params::Array{<:MPolyElem, 1}) +function monomial_compress(io_equation, params::Array{<:MPolyRingElem, 1}) params_xs = isempty(params) ? empty(params) : gens(parent(first(params))) other_vars = [ v for v in gens(parent(io_equation)) if @@ -197,7 +197,10 @@ Output: Computes the Wronskians of io_equations """ -@timeit _to function wronskian(io_equations::Dict{P, P}, ode::ODE{P}) where {P <: MPolyElem} +@timeit _to function wronskian( + io_equations::Dict{P, P}, + ode::ODE{P}, +) where {P <: MPolyRingElem} @debug "Compressing monomials" termlists = [monomial_compress(ioeq, ode)[2] for ioeq in values(io_equations)] @debug "Matrix sizes $(map(length, termlists))" @@ -207,9 +210,9 @@ Computes the Wronskians of io_equations # reducing everything modulo prime PRIME = 2^31 - 1 - F = Nemo.GF(PRIME) + F = Nemo.Native.GF(PRIME) polyring_red, gens_red = - Nemo.PolynomialRing(F, map(var_to_str, gens(parent(termlists[1][1])))) + Nemo.polynomial_ring(F, map(var_to_str, gens(parent(termlists[1][1])))) termlists = [map(p -> parent_ring_change(p, polyring_red), tlist) for tlist in termlists] ode_red = reduce_ode_mod_p(ode, PRIME) @@ -223,7 +226,7 @@ Computes the Wronskians of io_equations ord, ) @debug "Computing the derivatives of the solution" - ps_ext = Dict{MPolyElem, Nemo.gfp_abs_series}()# Generic.AbsSeries}() + ps_ext = Dict{MPolyRingElem, Nemo.gfp_abs_series}()# Generic.AbsSeries}() for v in vcat(ode_red.y_vars, ode_red.u_vars) cur_ps = ps[v] for i in 0:length(ode_red.x_vars) @@ -237,7 +240,7 @@ Computes the Wronskians of io_equations for (i, tlist) in enumerate(termlists) n = length(tlist) evaled = massive_eval(tlist, ps_ext) - S = Nemo.MatrixSpace(F, n, n) + S = Nemo.matrix_space(F, n, n) W = zero(S) for (i, ps) in enumerate(evaled) for j in 1:n diff --git a/test/RationalFunctionFields/normalforms.jl b/test/RationalFunctionFields/normalforms.jl index abce0d95f..a6635b669 100644 --- a/test/RationalFunctionFields/normalforms.jl +++ b/test/RationalFunctionFields/normalforms.jl @@ -62,7 +62,7 @@ eq_up_to_the_order(a, b) = issubset(a, b) && issubset(b, a) ### # Some arbitrary generators for the SLIQR model R, (b, e, In, S, Ninv, s, Q, g, u, a, y, L) = - PolynomialRing(QQ, [:b, :e, :In, :S, :Ninv, :s, :Q, :g, :u, :a, :y, :L]) + polynomial_ring(QQ, [:b, :e, :In, :S, :Ninv, :s, :Q, :g, :u, :a, :y, :L]) f = [ In // one(R), s // one(R), diff --git a/test/check_field_membership.jl b/test/check_field_membership.jl index 6521db4d7..a72ed5d23 100644 --- a/test/check_field_membership.jl +++ b/test/check_field_membership.jl @@ -1,5 +1,5 @@ @testset "Check field membership" begin - R, (x, y, z) = Nemo.PolynomialRing(Nemo.QQ, ["x", "y", "z"]) + R, (x, y, z) = Nemo.polynomial_ring(Nemo.QQ, ["x", "y", "z"]) @test field_contains( RationalFunctionField([[R(1), x + y], [R(1), x * y], [z, (x + y)^2]]), diff --git a/test/check_primality_zerodim.jl b/test/check_primality_zerodim.jl index 6178efa89..b9e880f1f 100644 --- a/test/check_primality_zerodim.jl +++ b/test/check_primality_zerodim.jl @@ -1,11 +1,13 @@ -@testset "Primality check (zerodim subroutine)" begin - R, (x, y) = Nemo.PolynomialRing(Nemo.QQ, ["x", "y"]) +if GROUP == "All" || GROUP == "Core" + @testset "Primality check (zerodim subroutine)" begin + R, (x, y) = Nemo.polynomial_ring(Nemo.QQ, ["x", "y"]) - @test check_primality_zerodim([x^2 - 1, y^2 - 4]) == false + @test check_primality_zerodim([x^2 - 1, y^2 - 4]) == false - @test check_primality_zerodim([(x + 5) * (x^3 - 7), y - 3]) == false + @test check_primality_zerodim([(x + 5) * (x^3 - 7), y - 3]) == false - @test check_primality_zerodim([x^3 - 5, y - 1]) == true + @test check_primality_zerodim([x^3 - 5, y - 1]) == true - @test check_primality_zerodim([x^2 + 1, y^3 - 3 * x + x + 5]) == true + @test check_primality_zerodim([x^2 + 1, y^3 - 3 * x + x + 5]) == true + end end diff --git a/test/common_ring.jl b/test/common_ring.jl index 91512b439..acec21ecf 100644 --- a/test/common_ring.jl +++ b/test/common_ring.jl @@ -1,39 +1,42 @@ -@testset "Computing common ring for the PB-reduction" begin - ode = @ODEmodel(x1'(t) = x2(t), x2'(t) = a * x1(t), y(t) = x1(t)) - ioeqs = find_ioequations(ode) - pbr = PBRepresentation(ode, ioeqs) - R, (y_2, y_5, c) = Nemo.PolynomialRing(Nemo.QQ, ["y(t)_2", "y(t)_5", "c"]) - p = y_2^2 + c * y_5 - (r, der) = common_ring(p, pbr) - @test Set(map(var_to_str, gens(r))) == - Set(["y(t)_0", "y(t)_1", "y(t)_2", "y(t)_3", "y(t)_4", "y(t)_5", "c", "a"]) +if GROUP == "All" || GROUP == "Core" + @testset "Computing common ring for the PB-reduction" begin + ode = @ODEmodel(x1'(t) = x2(t), x2'(t) = a * x1(t), y(t) = x1(t)) + ioeqs = find_ioequations(ode) + pbr = PBRepresentation(ode, ioeqs) + R, (y_2, y_5, c) = Nemo.polynomial_ring(Nemo.QQ, ["y(t)_2", "y(t)_5", "c"]) + p = y_2^2 + c * y_5 + (r, der) = common_ring(p, pbr) + @test Set(map(var_to_str, gens(r))) == + Set(["y(t)_0", "y(t)_1", "y(t)_2", "y(t)_3", "y(t)_4", "y(t)_5", "c", "a"]) - ode = @ODEmodel( - x1'(t) = x3(t), - x2'(t) = a * x2(t), - x3'(t) = x1(t), - y1(t) = x1(t), - y2(t) = x2(t) + u(t) - ) - ioeqs = find_ioequations(ode) - pbr = PBRepresentation(ode, ioeqs) - R, (y1_0, y2_3, u_3) = Nemo.PolynomialRing(Nemo.QQ, ["y1(t)_0", "y2(t)_3", "u(t)_3"]) - p = y1_0 + y2_3 + u_3 - (r, der) = common_ring(p, pbr) - @test Set([var_to_str(v) for v in gens(r)]) == Set([ - "y1(t)_0", - "y1(t)_1", - "y1(t)_2", - "y1(t)_3", - "y1(t)_4", - "y2(t)_0", - "y2(t)_1", - "y2(t)_2", - "y2(t)_3", - "u(t)_0", - "u(t)_1", - "u(t)_2", - "u(t)_3", - "a", - ]) + ode = @ODEmodel( + x1'(t) = x3(t), + x2'(t) = a * x2(t), + x3'(t) = x1(t), + y1(t) = x1(t), + y2(t) = x2(t) + u(t) + ) + ioeqs = find_ioequations(ode) + pbr = PBRepresentation(ode, ioeqs) + R, (y1_0, y2_3, u_3) = + Nemo.polynomial_ring(Nemo.QQ, ["y1(t)_0", "y2(t)_3", "u(t)_3"]) + p = y1_0 + y2_3 + u_3 + (r, der) = common_ring(p, pbr) + @test Set([var_to_str(v) for v in gens(r)]) == Set([ + "y1(t)_0", + "y1(t)_1", + "y1(t)_2", + "y1(t)_3", + "y1(t)_4", + "y2(t)_0", + "y2(t)_1", + "y2(t)_2", + "y2(t)_3", + "u(t)_0", + "u(t)_1", + "u(t)_2", + "u(t)_3", + "a", + ]) + end end diff --git a/test/constructive_membership.jl b/test/constructive_membership.jl index d80aca25f..56bfe0636 100644 --- a/test/constructive_membership.jl +++ b/test/constructive_membership.jl @@ -1,125 +1,131 @@ -@testset "Constructive field membership" begin - R, (x,) = PolynomialRing(Nemo.QQ, ["x"]) +if GROUP == "All" || GROUP == "Core" + @testset "Constructive field membership" begin + R, (x,) = polynomial_ring(Nemo.QQ, ["x"]) - generators = [x^2, x^3] - to_be_reduced = [x^2, x, 3one(R), zero(R)] + generators = [x^2, x^3] + to_be_reduced = [x^2, x, 3one(R), zero(R)] - memberships, remainders, relations_between_tags, tag_to_gen = - StructuralIdentifiability.check_constructive_field_membership( - StructuralIdentifiability.RationalFunctionField(generators), - map(f -> f // one(f), to_be_reduced), - tag_names = ["T1", "T2"], - ) - tags = gens(base_ring(parent(first(remainders)))) - - @test length(tags) == 2 - @test all(memberships) - @test map(string, remainders) == ["T1", "T1^2//T2", "3", "0"] - @test tag_to_gen == Dict(tags[1] => x^2, tags[2] => x^3) - @test length(relations_between_tags) == 1 - @test string(relations_between_tags[1]) == "T1^3 - T2^2" + memberships, remainders, relations_between_tags, tag_to_gen = + StructuralIdentifiability.check_constructive_field_membership( + StructuralIdentifiability.RationalFunctionField(generators), + map(f -> f // one(f), to_be_reduced), + tag_names = ["T1", "T2"], + ) + tags = gens(base_ring(parent(first(remainders)))) - cases = [] + @test length(tags) == 2 + @test all(memberships) + @test map(string, remainders) == ["T1", "T1^2//T2", "3", "0"] + @test tag_to_gen == Dict(tags[1] => x^2, tags[2] => x^3) + @test length(relations_between_tags) == 1 + @test string(relations_between_tags[1]) == "T1^3 - T2^2" - R, (T1,) = PolynomialRing(Nemo.QQ, ["T1"]) - append!( - cases, - [(generators = [T1^2], to_be_reduced = [T1, T1^2], memberships = Bool[0, 1])], - ) + cases = [] - R, (T1, t, _t) = PolynomialRing(Nemo.QQ, ["T1", "t", "_t"]) - append!( - cases, - [( - generators = [T1, t, _t], - to_be_reduced = [_t, t, T1 * t * _t], - memberships = Bool[1, 1, 1], - )], - ) + R, (T1,) = polynomial_ring(Nemo.QQ, ["T1"]) + append!( + cases, + [(generators = [T1^2], to_be_reduced = [T1, T1^2], memberships = Bool[0, 1])], + ) - R, (x,) = PolynomialRing(Nemo.QQ, ["x"]) - append!( - cases, - [ - ( - generators = [(x - 1) // R(1), R(1) // (x^5 - 1), x // R(1)], - to_be_reduced = [ - (x^4 + x^3 + x^2 + x + 1) // one(R), - x // R(1), - R(33) // x^2, - ], + R, (T1, t, _t) = polynomial_ring(Nemo.QQ, ["T1", "t", "_t"]) + append!( + cases, + [( + generators = [T1, t, _t], + to_be_reduced = [_t, t, T1 * t * _t], memberships = Bool[1, 1, 1], - ), - ( - generators = [(x^10 + x^9 + x^2 + 1) // (x^7 - x^6 - x^3 + 1)], - to_be_reduced = [x // one(R), 2x // one(R), -3x // one(R)], - memberships = Bool[0, 0, 0], - ), - (generators = [x^2], to_be_reduced = [x, x^88], memberships = Bool[0, 1]), - ], - ) + )], + ) - R, (x, y, z) = PolynomialRing(Nemo.QQ, ["x", "y", "z"]) - append!( - cases, - [ - (generators = [x, y], to_be_reduced = [x^2 + y^2, z], memberships = Bool[1, 0]), - ( - generators = [x^2 + y^2, x^3 + y^3, x^4 + y^4], - to_be_reduced = [x * y, x + y, x + y + 1, x + y + z], - memberships = Bool[1, 1, 1, 0], - ), - ( - generators = [(x + y + z)^2, (x + y + z)^3, (x + y + z)^4], - to_be_reduced = [(x + y + z)^18, x + 1, y + 2, z + 3], - memberships = Bool[1, 0, 0, 0], - ), - ], - ) + R, (x,) = polynomial_ring(Nemo.QQ, ["x"]) + append!( + cases, + [ + ( + generators = [(x - 1) // R(1), R(1) // (x^5 - 1), x // R(1)], + to_be_reduced = [ + (x^4 + x^3 + x^2 + x + 1) // one(R), + x // R(1), + R(33) // x^2, + ], + memberships = Bool[1, 1, 1], + ), + ( + generators = [(x^10 + x^9 + x^2 + 1) // (x^7 - x^6 - x^3 + 1)], + to_be_reduced = [x // one(R), 2x // one(R), -3x // one(R)], + memberships = Bool[0, 0, 0], + ), + (generators = [x^2], to_be_reduced = [x, x^88], memberships = Bool[0, 1]), + ], + ) - # NOTE: in this case it actually matter to cancel out the gcd after - # computing the normal forms - R, (a, b, y, x2, c, x1) = PolynomialRing(Nemo.QQ, ["a", "b", "y", "x2", "c", "x1"]) - append!( - cases, - [ - ( - generators = [ - x1 // one(R), - a // one(R), - (a * c + c^2) // one(R), - c // x2, - x2 // (a + b), - ], - to_be_reduced = [ - (a * c + c^2 + x1) // (a * c + c^2), - (a * c + c^2 + x1) // (a^2 + a * b + a * c + b * c), - (a * x2 + a * x1 + b * x1) // x2, - ], - memberships = Bool[1, 1, 1], - ), - ], - ) + R, (x, y, z) = polynomial_ring(Nemo.QQ, ["x", "y", "z"]) + append!( + cases, + [ + ( + generators = [x, y], + to_be_reduced = [x^2 + y^2, z], + memberships = Bool[1, 0], + ), + ( + generators = [x^2 + y^2, x^3 + y^3, x^4 + y^4], + to_be_reduced = [x * y, x + y, x + y + 1, x + y + z], + memberships = Bool[1, 1, 1, 0], + ), + ( + generators = [(x + y + z)^2, (x + y + z)^3, (x + y + z)^4], + to_be_reduced = [(x + y + z)^18, x + 1, y + 2, z + 3], + memberships = Bool[1, 0, 0, 0], + ), + ], + ) - for case in cases - generators = case.generators - to_be_reduced = case.to_be_reduced - memberships, remainders, relations_between_tags, tag_to_gen = - StructuralIdentifiability.check_constructive_field_membership( - StructuralIdentifiability.RationalFunctionField(generators), - map(f -> f // one(f), to_be_reduced), - ) - @test memberships == case.memberships - tags = gens(base_ring(parent(first(remainders)))) - evaluate_tags = poly -> evaluate(poly, [tag_to_gen[tag] for tag in tags]) - for i in 1:length(relations_between_tags) - @test iszero(evaluate_tags(relations_between_tags[i])) - end - for i in 1:length(remainders) - if !memberships[i] - continue + # NOTE: in this case it actually matter to cancel out the gcd after + # computing the normal forms + R, (a, b, y, x2, c, x1) = polynomial_ring(Nemo.QQ, ["a", "b", "y", "x2", "c", "x1"]) + append!( + cases, + [ + ( + generators = [ + x1 // one(R), + a // one(R), + (a * c + c^2) // one(R), + c // x2, + x2 // (a + b), + ], + to_be_reduced = [ + (a * c + c^2 + x1) // (a * c + c^2), + (a * c + c^2 + x1) // (a^2 + a * b + a * c + b * c), + (a * x2 + a * x1 + b * x1) // x2, + ], + memberships = Bool[1, 1, 1], + ), + ], + ) + + for case in cases + generators = case.generators + to_be_reduced = case.to_be_reduced + memberships, remainders, relations_between_tags, tag_to_gen = + StructuralIdentifiability.check_constructive_field_membership( + StructuralIdentifiability.RationalFunctionField(generators), + map(f -> f // one(f), to_be_reduced), + ) + @test memberships == case.memberships + tags = gens(base_ring(parent(first(remainders)))) + evaluate_tags = poly -> evaluate(poly, [tag_to_gen[tag] for tag in tags]) + for i in 1:length(relations_between_tags) + @test iszero(evaluate_tags(relations_between_tags[i])) + end + for i in 1:length(remainders) + if !memberships[i] + continue + end + @test iszero(evaluate_tags(remainders[i]) - to_be_reduced[i]) end - @test iszero(evaluate_tags(remainders[i]) - to_be_reduced[i]) end end end diff --git a/test/decompose_derivative.jl b/test/decompose_derivative.jl index d44311707..7f2c8c060 100644 --- a/test/decompose_derivative.jl +++ b/test/decompose_derivative.jl @@ -1,12 +1,14 @@ -@testset "Decomposing derivative" begin - cases = [ - ["yy_11", ["y", "yy", "yy_"], ("yy", 11)], - ["xx_xx_xx_0", ["xx", "x", "xx_xx_xx"], ("xx_xx_xx", 0)], - ["abc154f", ["ab", "abc"], nothing], - ["c_1542673", ["a", "b", "c"], ("c", 1542673)], - ["a", ["a"], nothing], - ] - for c in cases - @test decompose_derivative(c[1], c[2]) == c[3] +if GROUP == "All" || GROUP == "Core" + @testset "Decomposing derivative" begin + cases = [ + ["yy_11", ["y", "yy", "yy_"], ("yy", 11)], + ["xx_xx_xx_0", ["xx", "x", "xx_xx_xx"], ("xx_xx_xx", 0)], + ["abc154f", ["ab", "abc"], nothing], + ["c_1542673", ["a", "b", "c"], ("c", 1542673)], + ["a", ["a"], nothing], + ] + for c in cases + @test decompose_derivative(c[1], c[2]) == c[3] + end end end diff --git a/test/det_minor_expansion.jl b/test/det_minor_expansion.jl index 4240feae7..42a84069a 100644 --- a/test/det_minor_expansion.jl +++ b/test/det_minor_expansion.jl @@ -1,9 +1,11 @@ -@testset "Determinant by minor expansion" begin - for d in 1:5 - for testcase in 1:10 - mat_space = Nemo.MatrixSpace(Nemo.QQ, d, d) - rnd_matrix = mat_space([mod(rand(Int), 1000) for i in 1:d, j in 1:d]) - @test det(rnd_matrix) == det_minor_expansion(rnd_matrix) +if GROUP == "All" || GROUP == "Core" + @testset "Determinant by minor expansion" begin + for d in 1:5 + for testcase in 1:10 + mat_space = Nemo.matrix_space(Nemo.QQ, d, d) + rnd_matrix = mat_space([mod(rand(Int), 1000) for i in 1:d, j in 1:d]) + @test det(rnd_matrix) == det_minor_expansion(rnd_matrix) + end end end end diff --git a/test/diff_sequence_solution.jl b/test/diff_sequence_solution.jl index 1df86d8f7..912a632cb 100644 --- a/test/diff_sequence_solution.jl +++ b/test/diff_sequence_solution.jl @@ -1,103 +1,117 @@ -@testset "Computing variations around a sequence solution" begin - # Computing sensitivities directly be explicitly writing down Lie derivatives - function use_lie_derivatives( - dds::ODE{P}, - params::Dict{P, T}, - ic::Dict{P, T}, - inputs::Dict{P, Array{T, 1}}, - num_terms::Int, - ) where {T <: Generic.FieldElem, P <: MPolyElem{T}} - newvars = [var_to_str(v) for v in gens(dds.poly_ring)] - append!(newvars, [var_to_str(v) * "$i" for v in dds.u_vars for i in 1:num_terms]) - R, _ = - StructuralIdentifiability.Nemo.PolynomialRing(base_ring(dds.poly_ring), newvars) - explicit_sol = merge( - Dict( - parent_ring_change(x, R) => Vector{Any}([parent_ring_change(x, R)]) for - (x, eq) in dds.x_equations - ), - Dict( - parent_ring_change(y, R) => Vector{Any}([parent_ring_change(eq, R)]) for - (y, eq) in dds.y_equations - ), - ) - time_step = merge( - Dict( - parent_ring_change(x, R) => parent_ring_change(eq, R) for - (x, eq) in dds.x_equations - ), - Dict( - parent_ring_change(p, R) => parent_ring_change(p, R) for - p in dds.parameters - ), - Dict( - parent_ring_change(u, R) => str_to_var(var_to_str(u) * "1", R) for - u in dds.u_vars - ), - Dict( - str_to_var(s * "$i", R) => str_to_var(s * "$(i + 1)", R) for - s in map(var_to_str, dds.u_vars) for i in 1:(num_terms - 1) - ), - ) - eval_dict = merge( - Dict(parent_ring_change(p, R) => v for (p, v) in params), - Dict(parent_ring_change(x, R) => v for (x, v) in ic), - Dict(parent_ring_change(u, R) => val[1] for (u, val) in inputs), - Dict( - str_to_var(var_to_str(u) * "$i", R) => inputs[u][i + 1] for - u in dds.u_vars for i in 1:(num_terms - 1) - ), - ) - generalized_parameters = - [parent_ring_change(p, R) for p in vcat(dds.x_vars, dds.parameters)] - for i in 2:num_terms - for (k, v) in explicit_sol - push!(explicit_sol[k], eval_at_dict(v[end], time_step)) +if GROUP == "All" || GROUP == "Core" + @testset "Computing variations around a sequence solution" begin + # Computing sensitivities directly be explicitly writing down Lie derivatives + function use_lie_derivatives( + dds::StructuralIdentifiability.DDS{P}, + params::Dict{P, T}, + ic::Dict{P, T}, + input_values::Dict{P, Array{T, 1}}, + num_terms::Int, + ) where {T <: Generic.FieldElem, P <: MPolyRingElem{T}} + newvars = [var_to_str(v) for v in gens(parent(dds))] + append!( + newvars, + [var_to_str(v) * "$i" for v in inputs(dds) for i in 1:num_terms], + ) + R, _ = StructuralIdentifiability.Nemo.polynomial_ring( + base_ring(parent(dds)), + newvars, + ) + explicit_sol = merge( + Dict( + parent_ring_change(x, R) => Vector{Any}([parent_ring_change(x, R)]) + for (x, eq) in x_equations(dds) + ), + Dict( + parent_ring_change(y, R) => + Vector{Any}([parent_ring_change(eq, R)]) for + (y, eq) in y_equations(dds) + ), + ) + time_step = merge( + Dict( + parent_ring_change(x, R) => parent_ring_change(eq, R) for + (x, eq) in x_equations(dds) + ), + Dict( + parent_ring_change(p, R) => parent_ring_change(p, R) for + p in StructuralIdentifiability.parameters(dds) + ), + Dict( + parent_ring_change(u, R) => str_to_var(var_to_str(u) * "1", R) for + u in inputs(dds) + ), + Dict( + str_to_var(s * "$i", R) => str_to_var(s * "$(i + 1)", R) for + s in map(var_to_str, inputs(dds)) for i in 1:(num_terms - 1) + ), + ) + eval_dict = merge( + Dict(parent_ring_change(p, R) => v for (p, v) in params), + Dict(parent_ring_change(x, R) => v for (x, v) in ic), + Dict(parent_ring_change(u, R) => val[1] for (u, val) in input_values), + Dict( + str_to_var(var_to_str(u) * "$i", R) => input_values[u][i + 1] for + u in inputs(dds) for i in 1:(num_terms - 1) + ), + ) + generalized_parameters = [ + parent_ring_change(p, R) for + p in vcat(x_vars(dds), StructuralIdentifiability.parameters(dds)) + ] + for i in 2:num_terms + for (k, v) in explicit_sol + push!(explicit_sol[k], eval_at_dict(v[end], time_step)) + end end - end - part_diffs = - Dict((f, p) => [] for f in keys(explicit_sol) for p in generalized_parameters) - for i in 1:num_terms - for (f, ders) in explicit_sol - for p in generalized_parameters - push!( - part_diffs[(f, p)], - eval_at_dict(derivative(ders[i], p), eval_dict), - ) + part_diffs = Dict( + (f, p) => [] for f in keys(explicit_sol) for p in generalized_parameters + ) + for i in 1:num_terms + for (f, ders) in explicit_sol + for p in generalized_parameters + push!( + part_diffs[(f, p)], + eval_at_dict(derivative(ders[i], p), eval_dict), + ) + end end end + return Dict( + ( + parent_ring_change(k[1], parent(dds)), + parent_ring_change(k[2], parent(dds)), + ) => res for (k, res) in part_diffs + ) end - return Dict( - ( - parent_ring_change(k[1], dds.poly_ring), - parent_ring_change(k[2], dds.poly_ring), - ) => res for (k, res) in part_diffs - ) - end - locQQ = StructuralIdentifiability.Nemo.QQ + locQQ = StructuralIdentifiability.Nemo.QQ - ode = @ODEmodel(a'(t) = a(t)^2 + b, y(t) = 1 / (a(t) * c(t))) - params = Dict(b => locQQ(1)) - ic = Dict(a => locQQ(2)) - inputs = Dict(c => [locQQ(1), locQQ(-2), locQQ(3), locQQ(-4), locQQ(5)]) - seq_sol, diff_sol = differentiate_sequence_solution(ode, params, ic, inputs, 4) - diff_y = differentiate_sequence_output(ode, params, ic, inputs, 4) - lie_ders_sol = use_lie_derivatives(ode, params, ic, inputs, 4) - merged = merge(diff_sol, diff_y) - @test merged == lie_ders_sol + dds = @DDSmodel(a(t + 1) = a(t)^2 + b, y(t) = 1 / (a(t) * c(t))) + params = Dict(b => locQQ(1)) + ic = Dict(a => locQQ(2)) + input_values = Dict(c => [locQQ(1), locQQ(-2), locQQ(3), locQQ(-4), locQQ(5)]) + seq_sol, diff_sol = + differentiate_sequence_solution(dds, params, ic, input_values, 4) + diff_y = differentiate_sequence_output(dds, params, ic, input_values, 4) + lie_ders_sol = use_lie_derivatives(dds, params, ic, input_values, 4) + merged = merge(diff_sol, diff_y) + @test merged == lie_ders_sol - ode = @ODEmodel( - a'(t) = (23 * k1 * a(t) - 7 * b(t)^3) // (a(t)^2 + b(t)^2) - c(t)^3 * k1 * b(t), - b'(t) = a(t) + 17 * (b(t) - c(t))^2 + 1 // (a(t) + b(t) - k2), - y(t) = (a(t) + b(t) - c(t)) // (a(t)^2 + k2) - ) - params = Dict(k1 => locQQ(1), k2 => locQQ(2)) - ic = Dict(a => locQQ(3), b => locQQ(-4)) - inputs = Dict(c => [locQQ(5), locQQ(-6), locQQ(7), locQQ(-8)]) - seq_sol, diff_sol = differentiate_sequence_solution(ode, params, ic, inputs, 2) - diff_y = differentiate_sequence_output(ode, params, ic, inputs, 2) - lie_ders_sol = use_lie_derivatives(ode, params, ic, inputs, 2) - merged = merge(diff_sol, diff_y) - @test merged == lie_ders_sol + dds = @DDSmodel( + a(t + 1) = + (23 * k1 * a(t) - 7 * b(t)^3) // (a(t)^2 + b(t)^2) - c(t)^3 * k1 * b(t), + b(t + 1) = a(t) + 17 * (b(t) - c(t))^2 + 1 // (a(t) + b(t) - k2), + y(t) = (a(t) + b(t) - c(t)) // (a(t)^2 + k2) + ) + params = Dict(k1 => locQQ(1), k2 => locQQ(2)) + ic = Dict(a => locQQ(3), b => locQQ(-4)) + input_values = Dict(c => [locQQ(5), locQQ(-6), locQQ(7), locQQ(-8)]) + seq_sol, diff_sol = + differentiate_sequence_solution(dds, params, ic, input_values, 2) + diff_y = differentiate_sequence_output(dds, params, ic, input_values, 2) + lie_ders_sol = use_lie_derivatives(dds, params, ic, input_values, 2) + merged = merge(diff_sol, diff_y) + @test merged == lie_ders_sol + end end diff --git a/test/differentiate_output.jl b/test/differentiate_output.jl index 75d45dc0a..3c8ef8662 100644 --- a/test/differentiate_output.jl +++ b/test/differentiate_output.jl @@ -9,7 +9,7 @@ function diff_sol_Lie_derivatives(ode::ODE, params, ic, inputs, prec::Int) [var_to_str(u) * "_$i" for u in ode.u_vars for i in 0:(prec - 1)], ) end - new_ring, vars = Nemo.PolynomialRing(base_ring(ode.poly_ring), new_varnames) + new_ring, vars = Nemo.polynomial_ring(base_ring(ode.poly_ring), new_varnames) # mapping everything to the new ring eval_point = Dict(v => switch_ring(v, new_ring) for v in gens(ode.poly_ring)) @@ -77,35 +77,9 @@ end #------------------------------------------------------------------------------ -function rand_poly(deg, vars) - if deg == 0 - return parent(vars[1])(1) - end - result = 0 - indices = collect(1:length(vars)) - monomials = [] - for d in 0:deg - for subs in StructuralIdentifiability.IterTools.subsets(indices, d) - push!(monomials, subs) - end - end - - for subs in monomials - monom = rand(-50:50) - for v_ind in subs - monom *= vars[v_ind] - end - result += monom - end - - return result -end - -#------------------------------------------------------------------------------ - @testset "Partial derivatives of an output w.r.t. to initial conditions and parameters" begin test_cases = [] - P = fmpq_mpoly + P = QQMPolyRingElem DType = Union{P, Generic.Frac{P}} ode = @ODEmodel(x'(t) = x(t) + a, y(t) = x(t)^2) @@ -115,7 +89,7 @@ end :ODE => ode, :ic => Dict(x => Nemo.QQ(rand(1:10))), :param_vals => Dict(a => Nemo.QQ(rand(1:10))), - :inputs => Dict{P, Array{fmpq, 1}}(), + :inputs => Dict{P, Array{QQFieldElem, 1}}(), :prec => 20, ), ) @@ -127,7 +101,7 @@ end :ODE => ode, :ic => Dict(x => Nemo.QQ(rand(1:10))), :param_vals => Dict(a => Nemo.QQ(rand(1:10))), - :inputs => Dict{P, Array{fmpq, 1}}(), + :inputs => Dict{P, Array{QQFieldElem, 1}}(), :prec => 20, ), ) @@ -144,7 +118,7 @@ end :ODE => ode, :ic => Dict(x => Nemo.QQ(rand(1:10)), y => Nemo.QQ(rand(1:10))), :param_vals => Dict(a => Nemo.QQ(rand(1:10)), b => Nemo.QQ(rand(1:10))), - :inputs => Dict{P, Array{fmpq, 1}}(), + :inputs => Dict{P, Array{QQFieldElem, 1}}(), :prec => 8, ), ) @@ -161,8 +135,8 @@ end ), ) - F = Nemo.GF(2^31 - 1) - P = gfp_mpoly + F = Nemo.Native.GF(2^31 - 1) + P = fpMPolyRingElem DType = Union{P, Generic.Frac{P}} varnames = vcat( @@ -171,7 +145,7 @@ end ["u_$i" for i in 1:2], ["y_$i" for i in 1:3], ) - R, vars = Nemo.PolynomialRing(F, varnames) + R, vars = Nemo.polynomial_ring(F, varnames) push!( test_cases, Dict( @@ -193,7 +167,7 @@ end ["u_$i" for i in 1:2], ["y_$i" for i in 1:3], ) - R, vars = Nemo.PolynomialRing(F, varnames) + R, vars = Nemo.polynomial_ring(F, varnames) push!( test_cases, Dict( @@ -210,7 +184,7 @@ end ) varnames = vcat(["x_$i" for i in 1:2], ["p_$i" for i in 1:2], "u", ["y_1", "y_2"]) - R, vars = Nemo.PolynomialRing(F, varnames) + R, vars = Nemo.polynomial_ring(F, varnames) push!( test_cases, Dict( diff --git a/test/diffreduction.jl b/test/diffreduction.jl index b624438b2..9b8d11101 100644 --- a/test/diffreduction.jl +++ b/test/diffreduction.jl @@ -1,7 +1,7 @@ @testset "Differential reduction" begin ode = @ODEmodel(x'(t) = a * x(t), y(t) = x(t)) pbr = PBRepresentation(ode, find_ioequations(ode)) - R, (y_5, y_4) = Nemo.PolynomialRing(Nemo.QQ, ["y(t)_5", "y(t)_4"]) + R, (y_5, y_4) = Nemo.polynomial_ring(Nemo.QQ, ["y(t)_5", "y(t)_4"]) res = diffreduce(y_5, pbr) @test res == str_to_var("a", parent(res))^5 * str_to_var("y(t)_0", parent(res)) res = diffreduce(y_4, pbr) @@ -11,7 +11,7 @@ ode = @ODEmodel(x1'(t) = x2(t), x2'(t) = -x1(t), y(t) = x1(t) + u(t)) pbr = PBRepresentation(ode, find_ioequations(ode)) - R, (y_5, y_4, u_10) = Nemo.PolynomialRing(Nemo.QQ, ["y(t)_5", "y(t)_4", "u(t)_10"]) + R, (y_5, y_4, u_10) = Nemo.polynomial_ring(Nemo.QQ, ["y(t)_5", "y(t)_4", "u(t)_10"]) res = diffreduce(y_4 + u_10, pbr) @test res == str_to_var("u(t)_10", parent(res)) + str_to_var("y(t)_0", parent(res)) - @@ -20,7 +20,7 @@ # Next two verified with Maple # Mizuka's example R, (y_0, y_1, y_2, y_3, u_0, u_1, u_2) = - Nemo.PolynomialRing(Nemo.QQ, ["y_0", "y_1", "y_2", "y_3", "u_0", "u_1", "u_2"]) + Nemo.polynomial_ring(Nemo.QQ, ["y_0", "y_1", "y_2", "y_3", "u_0", "u_1", "u_2"]) pbr = PBRepresentation( ["y"], ["u"], @@ -127,7 +127,7 @@ y(t) = x1(t) ) pbr = PBRepresentation(ode, find_ioequations(ode)) - R, (y_0, y_1, y_2, y_3, y_4, u_0, u_3) = Nemo.PolynomialRing( + R, (y_0, y_1, y_2, y_3, y_4, u_0, u_3) = Nemo.polynomial_ring( Nemo.QQ, ["y(t)_0", "y(t)_1", "y(t)_2", "y(t)_3", "y(t)_4", "u(t)_0", "u(t)_3"], ) diff --git a/test/extensions/modelingtoolkit.jl b/test/extensions/modelingtoolkit.jl new file mode 100644 index 000000000..7753e5651 --- /dev/null +++ b/test/extensions/modelingtoolkit.jl @@ -0,0 +1,659 @@ +if GROUP == "All" || GROUP == "ModelingToolkitExt" + @testset "Check identifiability of `ODESystem` object" begin + using ModelingToolkit + using ModelingToolkit: parameters + using Symbolics + + @parameters a01 a21 a12 + @variables t x0(t) x1(t) y1(t) [output = true] + D = Differential(t) + + eqs = [D(x0) ~ -(a01 + a21) * x0 + a12 * x1, D(x1) ~ a21 * x0 - a12 * x1, y1 ~ x0] + de = ODESystem(eqs, t, name = :Test) + + correct = OrderedDict( + a01 => :nonidentifiable, + a21 => :nonidentifiable, + a12 => :nonidentifiable, + x0 => :globally, + x1 => :nonidentifiable, + ) + + @test isequal(correct, assess_identifiability(de; measured_quantities = [y1 ~ x0])) + @test isequal(correct, assess_identifiability(de; measured_quantities = [x0])) + @test isequal( + correct, + assess_identifiability(de; measured_quantities = [(y1 ~ x0).rhs]), + ) + + # check identifiabile functions + correct = [a01 * a12, a01 + a12 + a21] + result = find_identifiable_functions(de, measured_quantities = [y1 ~ x0]) + @test isequal(Set(correct), Set(result)) + + # -------------------------------------------------------------------------- + + # check identifiabile functions + @parameters V_m k_m k01 c + @variables t x(t) y1(t) [output = true] + D = Differential(t) + + eqs = [D(x) ~ (-V_m * x) / (k_m + x) + k01 * x, y1 ~ c * x] + de = ODESystem(eqs, t, name = :Test) + + correct = [k01, c * k_m, V_m * c] + result = find_identifiable_functions(de) + @test isequal(Set(correct), Set(result)) + + correct = [k01, c * x, k_m * c, V_m * c] + result = find_identifiable_functions(de, with_states = true) + @test isequal(Set(correct), Set(result)) + + # -------------------------------------------------------------------------- + @parameters a01 a21 a12 + @variables t x0(t) x1(t) y1(t) [output = true] + D = Differential(t) + + eqs = [D(x0) ~ -(a01 + a21) * x0 + a12 * x1, D(x1) ~ a21 * x0 - a12 * x1, y1 ~ x0] + de = ODESystem(eqs, t, name = :Test) + + correct = OrderedDict( + a01 => :nonidentifiable, + a21 => :nonidentifiable, + a12 => :nonidentifiable, + x0 => :globally, + x1 => :nonidentifiable, + ) + + @test isequal(correct, assess_identifiability(de)) + + # -------------------------------------------------------------------------- + + @parameters a01 a21 a12 + @variables t x0(t) x1(t) y1(t) [output = true] + D = Differential(t) + + eqs = [D(x0) ~ -(a01 + a21) * x0 + a12 * x1, D(x1) ~ a21 * x0 - a12 * x1, y1 ~ x0] + de = ODESystem(eqs, t, name = :Test) + funcs_to_check = [a01, a21, a12, a01 * a12, a01 + a12 + a21] + correct = OrderedDict( + a01 => :nonidentifiable, + a21 => :nonidentifiable, + a12 => :nonidentifiable, + a01 * a12 => :globally, + a01 + a12 + a21 => :globally, + ) + @test isequal(correct, assess_identifiability(de; funcs_to_check = funcs_to_check)) + + # -------------------------------------------------------------------------- + + @parameters a01 a21 a12 + @variables t x0(t) x1(t) y1(t) + D = Differential(t) + + eqs = [D(x0) ~ -(a01 + a21) * x0 + a12 * x1, D(x1) ~ a21 * x0 - a12 * x1] + measured_quantities = [y1 ~ x0] + de = ODESystem(eqs, t, name = :Test) + funcs_to_check = [a01, a21, a12, a01 * a12, a01 + a12 + a21] + correct = OrderedDict( + a01 => :nonidentifiable, + a21 => :nonidentifiable, + a12 => :nonidentifiable, + a01 * a12 => :globally, + a01 + a12 + a21 => :globally, + ) + @test isequal( + correct, + assess_identifiability( + de; + measured_quantities = measured_quantities, + funcs_to_check = funcs_to_check, + ), + ) + + # -------------------------------------------------------------------------- + @parameters μ bi bw a χ γ k + @variables t S(t) I(t) W(t) R(t) y(t) + + eqs = [ + D(S) ~ μ - bi * S * I - bw * S * W - μ * S + a * R, + D(I) ~ bw * S * W + bi * S * I - (γ + μ) * I, + D(W) ~ χ * (I - W), + D(R) ~ γ * I - (μ + a) * R, + ] + de = ODESystem(eqs, t, name = :TestSIWR) + measured_quantities = [y ~ k * I] + # check all parameters (default) + @test isequal( + true, + all( + values( + assess_local_identifiability( + de; + measured_quantities = measured_quantities, + ), + ), + ), + ) + + # check specific parameters + funcs_to_check = [μ, bi, bw, a, χ, γ, γ + μ, k, S, I, W, R] + correct = OrderedDict(f => true for f in funcs_to_check) + @test isequal( + correct, + assess_local_identifiability( + de; + measured_quantities = measured_quantities, + funcs_to_check = funcs_to_check, + ), + ) + + # checking ME identifiability + funcs_to_check = [μ, bi, bw, a, χ, γ, γ + μ, k] + correct = OrderedDict(f => true for f in funcs_to_check) + @test isequal( + (correct, 1), + assess_local_identifiability( + de; + measured_quantities = measured_quantities, + funcs_to_check = funcs_to_check, + prob_threshold = 0.99, + type = :ME, + ), + ) + + # checking identifiabile functions + correct = [a, bw, χ, bi, k, γ, μ] + result = find_identifiable_functions(de, measured_quantities = measured_quantities) + @test isequal(Set(correct), Set(result)) + + # -------------------------------------------------------------------------- + @parameters mu bi bw a xi gm k + @variables t S(t) I(t) W(t) R(t) y(t) [output = true] + + eqs = [ + D(S) ~ mu - bi * S * I - bw * S * W - mu * S + a * R, + D(I) ~ bw * S * W + bi * S * I - (gm + mu) * I, + D(W) ~ xi * (I - W), + D(R) ~ gm * I - (mu + a) * R, + y ~ k * I, + ] + de = ODESystem(eqs, t, name = :TestSIWR) + # check all parameters (default) + @test isequal(true, all(values(assess_local_identifiability(de)))) + + @test isequal( + true, + all( + values(assess_local_identifiability(de; measured_quantities = [y ~ k * I])), + ), + ) + + # check specific parameters + funcs_to_check = [mu, bi, bw, a, xi, gm, gm + mu, k, S, I, W, R] + correct = OrderedDict(f => true for f in funcs_to_check) + @test isequal( + correct, + assess_local_identifiability(de; funcs_to_check = funcs_to_check), + ) + + # checking ME identifiability + funcs_to_check = [mu, bi, bw, a, xi, gm, gm + mu, k] + correct = OrderedDict(f => true for f in funcs_to_check) + @test isequal( + (correct, 1), + assess_local_identifiability( + de; + funcs_to_check = funcs_to_check, + prob_threshold = 0.99, + type = :ME, + ), + ) + + # -------------------------------------------------------------------------- + @parameters mu bi bw a xi gm k + @variables t S(t) I(t) W(t) R(t) y(t) + + eqs = [ + D(S) ~ 2.0 * mu - bi * S * I - bw * S * W - mu * S + a * R, + D(I) ~ bw * S * W + bi * S * I - (gm + mu) * I, + D(W) ~ xi * (I - 0.6 * W), + D(R) ~ gm * I - (mu + a) * R, + ] + de = ODESystem(eqs, t, name = :TestSIWR) + measured_quantities = [y ~ 1.57 * I * k] + funcs_to_check = [mu, bi, bw, a, xi, gm, mu, gm + mu, k, S, I, W, R] + correct = OrderedDict(f => true for f in funcs_to_check) + @test isequal( + correct, + assess_local_identifiability( + de; + measured_quantities = measured_quantities, + funcs_to_check = funcs_to_check, + ), + ) + + # checking ME identifiability + funcs_to_check = [bi, bw, a, xi, gm, mu, gm + mu, k] + correct = OrderedDict(f => true for f in funcs_to_check) + @test isequal( + (correct, 1), + assess_local_identifiability( + de; + measured_quantities = measured_quantities, + funcs_to_check = funcs_to_check, + prob_threshold = 0.99, + type = :ME, + ), + ) + + # ---------- + + @parameters a01 a21 a12 + @variables t x0(t) x1(t) y1(t) + D = Differential(t) + using SpecialFunctions + + eqs = [ + D(x0) ~ -(a01 + a21) * SpecialFunctions.erfc(x0) + a12 * x1, + D(x1) ~ a21 * x0 - a12 * x1, + ] + + de = ODESystem(eqs, t, name = :Test) + measured_quantities = [y1 ~ x0] + funcs_to_check = [a01, a21, a12, a01 * a12, a01 + a12 + a21] + correct = Dict( + a01 => :nonidentifiable, + a21 => :nonidentifiable, + a12 => :nonidentifiable, + a01 * a12 => :globally, + a01 + a12 + a21 => :globally, + ) + @test_throws ArgumentError assess_identifiability( + de; + measured_quantities = measured_quantities, + funcs_to_check = funcs_to_check, + ) + # ---------- + @parameters a b c + @variables t x1(t) x2(t) y(t) + D = Differential(t) + + eqs = [D(x1) ~ -a * x1 + x2 * b / (x1 + b / (c^2 - x2)), D(x2) ~ x2 * c^2 + x1] + de = ODESystem(eqs, t, name = :Test) + measured_quantities = [y ~ x2] + correct = Dict(a => :globally, b => :globally, c => :locally) + to_check = [a, b, c] + @test isequal( + correct, + assess_identifiability( + de; + measured_quantities = measured_quantities, + funcs_to_check = to_check, + ), + ) + + # check identifiabile functions + result = find_identifiable_functions(de, measured_quantities = measured_quantities) + correct = [b, a, c^2] + @test isequal(Set(result), Set(correct)) + + # ---------- + @parameters a b + @variables t c(t) x1(t) x2(t) y1(t) y2(t) + D = Differential(t) + + eqs = [ + D(x1) ~ -a * x1 + x2 * b / (x1 + b / (c^2 - x2)), + D(x2) ~ x2 * c^2 + x1, + D(c) ~ 0, + ] + de = ODESystem(eqs, t, name = :Test) + measured_quantities = [y1 ~ x2, y2 ~ c] + correct = OrderedDict(a => :globally, b => :globally) + to_check = [a, b] + @test isequal( + correct, + assess_identifiability( + de; + measured_quantities = measured_quantities, + funcs_to_check = to_check, + ), + ) + + #---------------------------------- + # Composable models test (from https://github.com/SciML/StructuralIdentifiability.jl/issues/162) + @variables t + function rabbits_creator(; name) + ps = @parameters α = 1.5 + vars = @variables x(t) = 1.0 z(t) = 0.0 [input = true] + D = Differential(t) + equs = [D(x) ~ α^2 * x + z] + + ODESystem(equs, t, vars, ps; name = name) + end + + function wolves_creator(; name) + ps = @parameters δ = 3.0 + vars = @variables y(t) = 1.0 q(t) = 0.0 [input = true] + D = Differential(t) + equs = [D(y) ~ -δ * y + q] + + ODESystem(equs, t, vars, ps; name = name) + end + + function lotka_volterra_creator(; name) + @named wolves = wolves_creator() + @named rabbits = rabbits_creator() + + ps = @parameters β = 1.0 γ = 1.0 + D = Differential(t) + + eqs = + [rabbits.z ~ -β * wolves.y * rabbits.x, wolves.q ~ γ * wolves.y * rabbits.x] + + ModelingToolkit.compose(ODESystem(eqs, t, [], ps; name = name), wolves, rabbits) + end + + function getbyname(sys, name) + println(name) + return first([ + v for v in vcat(states(sys), parameters(sys)) if + replace(string(v), "(t)" => "") == name + ]) + end + + @named ltk_mtk = lotka_volterra_creator() + simp_ltk_mtk = structural_simplify(ltk_mtk) + wolves₊δ = getbyname(simp_ltk_mtk, "wolves₊δ") + rabbits₊α = getbyname(simp_ltk_mtk, "rabbits₊α") + β = getbyname(simp_ltk_mtk, "β") + γ = getbyname(simp_ltk_mtk, "γ") + wolves₊y = getbyname(simp_ltk_mtk, "wolves₊y") + rabbits₊x = getbyname(simp_ltk_mtk, "rabbits₊x") + @variables y(t) + measured_quantities = [y ~ wolves₊y] + result = + assess_identifiability(simp_ltk_mtk, measured_quantities = measured_quantities) + correct = Dict( + rabbits₊α => :locally, + γ => :nonidentifiable, + β => :globally, + wolves₊δ => :globally, + rabbits₊x => :nonidentifiable, + wolves₊y => :globally, + ) + @test Dict(result) == correct + + #---------------------------------- + + @variables t, x(t), y(t), z(t), w(t) + @parameters a + @named sys = ODESystem([D(x) ~ a * y], t, [x], [a]; observed = [y ~ z, z ~ x]) + measured_quantities = [w ~ x] + result = assess_identifiability(sys, measured_quantities = measured_quantities) + @test result[a] == :globally + + result = find_identifiable_functions(sys, measured_quantities = measured_quantities) + @test isequal(result, [a]) + + #---------------------------------- + + # Tensor definition case as reported in + # https://github.com/SciML/StructuralIdentifiability.jl/issues/178 + @variables t, x(t)[1:2], y(t)[1:2] + @parameters k1, k2 + + eqs = [D(x[1]) ~ -k1 * x[2], D(x[2]) ~ -k2 * x[1]] + + sys = ODESystem(eqs, t, name = :example_vector) + correct = OrderedDict(x[1] => true, x[2] => true, k1 => true, k2 => true) + @test assess_local_identifiability(sys, measured_quantities = [x[1], x[2]]) == + correct + end + + @testset "Discrete local identifiability, ModelingToolkit interface" begin + cases = [] + + @parameters α β + @variables t S(t) I(t) R(t) y(t) + D = Difference(t; dt = 1.0) + + eqs = [D(S) ~ -β * S * I, D(I) ~ β * S * I - α * I, D(R) ~ α * I] + @named sir = DiscreteSystem(eqs) + push!( + cases, + Dict( + :dds => sir, + :res => OrderedDict(S => true, I => true, R => false, α => true, β => true), + :y => [y ~ I], + :y2 => [I], + :known_ic => Array{}[], + :to_check => Array{}[], + ), + ) + + @parameters θ + @variables t x(t) y(t) + D = Difference(t; dt = 1.0) + + eqs = [D(x) ~ θ * x^3 - x] + + @named eqs = DiscreteSystem(eqs) + push!( + cases, + Dict( + :dds => eqs, + :res => OrderedDict(x => true, θ => true), + :y => [y ~ x], + :y2 => [x], + :known_ic => Array{}[], + :to_check => Array{}[], + ), + ) + + @parameters θ β + @variables t x1(t) x2(t) y(t) + D = Difference(t; dt = 1.0) + + eqs = [D(x1) ~ x1 + x2, D(x2) ~ θ + β] + + @named eqs = DiscreteSystem(eqs) + push!( + cases, + Dict( + :dds => eqs, + :res => OrderedDict(x1 => true, x2 => true, θ => false, β => false), + :y => [y ~ x1], + :y2 => [x1], + :known_ic => Array{}[], + :to_check => Array{}[], + ), + ) + + @parameters a b c d + @variables t x1(t) x2(t) u(t) y2(t) + D = Difference(t; dt = 1.0) + + eqs = [D(x1) ~ a * x1 - b * x1 * x2 + u, D(x2) ~ -c * x2 + d * x1 * x2] + + @named lv = DiscreteSystem(eqs) + push!( + cases, + Dict( + :dds => lv, + :res => OrderedDict( + x1 => true, + x2 => false, + a => true, + b => false, + c => true, + d => true, + ), + :y => [y ~ x1], + :y2 => [x1], + :known_ic => Array{}[], + :to_check => Array{}[], + ), + ) + + push!( + cases, + Dict( + :dds => lv, + :res => OrderedDict(b * x2 => true), + :y => [y ~ x1], + :y2 => [x1], + :known_ic => Array{}[], + :to_check => [b * x2], + ), + ) + + push!( + cases, + Dict( + :dds => lv, + :res => OrderedDict( + x1 => true, + x2 => true, + a => true, + b => true, + c => true, + d => true, + ), + :y => [y ~ x1, y2 ~ x1 / x2], + :y2 => [x1, x1 / x2], + :known_ic => Array{}[], + :to_check => Array{}[], + ), + ) + + push!( + cases, + Dict( + :dds => lv, + :res => OrderedDict( + substitute(x1, Dict(t => 0)) => true, + substitute(x2, Dict(t => 0)) => true, + a => true, + b => true, + c => true, + d => true, + ), + :y => [y ~ x1], + :y2 => [x1], + :known_ic => [x2], + :to_check => Array{}[], + ), + ) + + # Example 1 from https://doi.org/10.1016/j.automatica.2008.03.019 + @parameters theta1 theta2 + @variables t x1(t) x2(t) u(t) y(t) + D = Difference(t; dt = 1.0) + + eqs = [D(x1) ~ theta1 * x1 + x2, D(x2) ~ (1 - theta2) * x1 + x2^2 + u - x2] + + @named abmd1 = DiscreteSystem(eqs) + push!( + cases, + Dict( + :dds => abmd1, + :res => OrderedDict(x1 => true, x2 => true, theta1 => true, theta2 => true), + :y => [y ~ x1], + :y2 => [x1], + :known_ic => Array{}[], + :to_check => Array{}[], + ), + ) + + # Example 2 from https://doi.org/10.1016/j.automatica.2008.03.019 + @parameters theta1 theta2 theta3 + @variables t x1(t) x2(t) u(t) y(t) y2(t) + D = Difference(t; dt = 1.0) + + eqs = [D(x1) ~ theta1 * x1^2 + theta2 * x2 + u - x1, D(x2) ~ theta3 * x1 - x2] + + @named abmd2 = DiscreteSystem(eqs) + push!( + cases, + Dict( + :dds => abmd2, + :res => OrderedDict( + x1 => true, + x2 => false, + theta1 => true, + theta2 => false, + theta3 => false, + ), + :y => [y ~ x1], + :y2 => [x1], + :known_ic => Array{}[], + :to_check => Array{}[], + ), + ) + push!( + cases, + Dict( + :dds => abmd2, + :res => Dict( + x1 => true, + x2 => true, + theta1 => true, + theta2 => true, + theta3 => true, + ), + :y => [y ~ x1, y2 ~ x2], + :y2 => [x1, x2], + :known_ic => Array{}[], + :to_check => Array{}[], + ), + ) + + @parameters a b + @variables t x1(t) y(t) + D = Difference(t; dt = 1.0) + + eqs = [D(x1) ~ a] + + @named kic = DiscreteSystem(eqs) + push!( + cases, + Dict( + :dds => kic, + :res => OrderedDict(x1 => false, a => true, b => false), + :y => [y ~ x1 + b], + :y2 => [x1 + b], + :known_ic => Array{}[], + :to_check => Array{}[], + ), + ) + push!( + cases, + Dict( + :dds => kic, + :res => + OrderedDict(substitute(x1, Dict(t => 0)) => true, a => true, b => true), + :y => [y ~ x1 + b], + :y2 => [x1 + b], + :known_ic => [x1], + :to_check => Array{}[], + ), + ) + + for c in cases + @test assess_local_identifiability( + c[:dds]; + measured_quantities = c[:y], + known_ic = c[:known_ic], + funcs_to_check = c[:to_check], + ) == c[:res] + @test assess_local_identifiability( + c[:dds]; + measured_quantities = c[:y2], + known_ic = c[:known_ic], + funcs_to_check = c[:to_check], + ) == c[:res] + end + end +end diff --git a/test/extract_coefficients.jl b/test/extract_coefficients.jl index da3bae87f..7a66fd8a7 100644 --- a/test/extract_coefficients.jl +++ b/test/extract_coefficients.jl @@ -1,5 +1,5 @@ @testset "Coefficient extraction for rational fucntions" begin - R, (x, y, z) = PolynomialRing(QQ, ["x", "y", "z"]) + R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]) C = extract_coefficients_ratfunc( (x^2 + y * z - y^2 * z^3 + 3 * x * z^3) // (x + y + z + z^2 * (x^2 + 1)), [z], @@ -14,12 +14,12 @@ (x^2 + 1) // 1, ]) - R, (x, y) = PolynomialRing(QQ, ["x", "y"]) + R, (x, y) = polynomial_ring(QQ, ["x", "y"]) f = (x^2 + y^2) // (1 - x - 3 * y) - @test Set(extract_coefficients_ratfunc(f, Vector{Nemo.fmpq_mpoly}())) == + @test Set(extract_coefficients_ratfunc(f, Vector{Nemo.QQMPolyRingElem}())) == Set([f, one(R) // 1]) - R, (x, y, u, v) = PolynomialRing(QQ, ["x", "y", "u", "v"]) + R, (x, y, u, v) = polynomial_ring(QQ, ["x", "y", "u", "v"]) C = extract_coefficients_ratfunc( (x + (y + 3) * u * v + y^2 * v^3) // (u + 3 * v - (x^2 + y^2) * u^2), [u, v], @@ -35,21 +35,21 @@ end @testset "Coefficient extraction for polynomials" begin - R, (x, y, z) = PolynomialRing(QQ, ["x", "y", "z"]) + R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]) C = extract_coefficients((y + z + 8), [x]) R_coef = parent(first(values(C))) y, z = gens(R_coef) @test symbols(R_coef) == [:y, :z] @test C == Dict([0] => y + z + 8) - R, (x, y, z) = PolynomialRing(QQ, ["x", "y", "z"]) + R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]) C = extract_coefficients((x^2 + y * z - y^2 * z^3 + 3 * x * z^3), [z]) R_coef = parent(first(values(C))) x, y = gens(R_coef) @test symbols(R_coef) == [:x, :y] @test C == Dict([3] => 3x - y^2, [1] => y, [0] => x^2) - R, (x, y, z) = PolynomialRing(QQ, ["x", "y", "z"]) + R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]) C = extract_coefficients((x^2 + y * z - y^2 * z^3 + 3 * x * z^3), [x, z]) R_coef = parent(first(values(C))) y = gens(R_coef)[1] diff --git a/test/find_leader.jl b/test/find_leader.jl index 1994e82d2..5b12feed9 100644 --- a/test/find_leader.jl +++ b/test/find_leader.jl @@ -10,7 +10,7 @@ println("IOEQS: ", ioeqs) pbr = PBRepresentation(ode, ioeqs) - R, (y1_0, y1_1, y1_2, y2_0, y2_1, y2_2) = Nemo.PolynomialRing( + R, (y1_0, y1_1, y1_2, y2_0, y2_1, y2_2) = Nemo.polynomial_ring( Nemo.QQ, ["y1(t)_0", "y1(t)_1", "y1(t)_2", "y2(t)_0", "y2(t)_1", "y2(t)_2"], ) diff --git a/test/io_cases.jl b/test/io_cases.jl index 3242d75c1..9710334cc 100644 --- a/test/io_cases.jl +++ b/test/io_cases.jl @@ -2,7 +2,7 @@ test_cases = [] # 2-compartiment model - R, (y_0, y_1, y_2, u_0, u_1, a01, a12, a21) = PolynomialRing( + R, (y_0, y_1, y_2, u_0, u_1, a01, a12, a21) = polynomial_ring( QQ, ["y(t)_0", "y(t)_1", "y(t)_2", "u(t)_0", "u(t)_1", "a01", "a12", "a21"], ) @@ -22,7 +22,7 @@ #--------------------------------------- # Chen-Lee model R, (y1_0, y1_1, y1_2, y1_3, a, b, c) = - PolynomialRing(QQ, ["y1(t)_0", "y1(t)_1", "y1(t)_2", "y1(t)_3", "a", "b", "c"]) + polynomial_ring(QQ, ["y1(t)_0", "y1(t)_1", "y1(t)_2", "y1(t)_3", "a", "b", "c"]) correct = 9 * y1_3^2 * y1_0^2 - 18 * y1_3 * y1_2 * y1_1 * y1_0 - 18 * y1_3 * y1_2 * y1_0^2 * a - 36 * y1_3 * y1_2 * y1_0^2 * b - @@ -121,7 +121,7 @@ # predator-prey model R, (y1_0, y1_1, y1_2, a, b, c, d) = - PolynomialRing(QQ, ["y1(t)_0", "y1(t)_1", "y1(t)_2", "a", "b", "c", "d"]) + polynomial_ring(QQ, ["y1(t)_0", "y1(t)_1", "y1(t)_2", "a", "b", "c", "d"]) correct = y1_2 * y1_0 - y1_1^2 - y1_1 * y1_0^2 * d + y1_1 * y1_0 * c + y1_0^3 * a * d - y1_0^2 * a * c diff --git a/test/lc_univariate.jl b/test/lc_univariate.jl index 1d0919312..08c2922c8 100644 --- a/test/lc_univariate.jl +++ b/test/lc_univariate.jl @@ -1,5 +1,5 @@ @testset "Univariate leading coefficient" begin - R, (x, y, z) = Nemo.PolynomialRing(Nemo.QQ, ["x", "y", "z"]) + R, (x, y, z) = Nemo.polynomial_ring(Nemo.QQ, ["x", "y", "z"]) p = x^2 * y + x^2 * (z + z^3) + y - 5 @test lc_univariate(p, x) == y + z + z^3 @test lc_univariate(p, z) == x^2 diff --git a/test/lie_derivative.jl b/test/lie_derivative.jl index 64ea774eb..250b22d71 100644 --- a/test/lie_derivative.jl +++ b/test/lie_derivative.jl @@ -1,6 +1,3 @@ -using Test -using TestSetExtensions - @testset "Lie derivative" begin ode = @ODEmodel( x1'(t) = a * x1(t) + b * u(t), diff --git a/test/linear_compartment.jl b/test/linear_compartment.jl index f52fad656..d25ca9f1e 100644 --- a/test/linear_compartment.jl +++ b/test/linear_compartment.jl @@ -157,7 +157,7 @@ case[:leaks], ) bring = ode.poly_ring - correct = Dict{fmpq_mpoly, Symbol}() + correct = Dict{QQMPolyRingElem, Symbol}() for (e, id) in case[:result] correct[str_to_var("a_$(e[2])_$(e[1])", bring)] = id end diff --git a/test/local_identifiability_discrete.jl b/test/local_identifiability_discrete.jl index fd56b77e4..adeeaf8be 100644 --- a/test/local_identifiability_discrete.jl +++ b/test/local_identifiability_discrete.jl @@ -1,243 +1,58 @@ -@testset "Discrete local identifiability, internal function" begin +@testset "Discrete local identifiability, @DDSmodel interface" begin cases = [] - @parameters α β - @variables t S(t) I(t) R(t) y(t) - D = Difference(t; dt = 1.0) + dds = @DDSmodel(a(t + 1) = (b + c) * a(t) + 1, y(t) = a(t)) - eqs = [D(S) ~ S - β * S * I, D(I) ~ I + β * S * I - α * I, D(R) ~ R + α * I] - @named sir = DiscreteSystem(eqs) push!( cases, Dict( - :dds => sir, - :res => OrderedDict(S => true, I => true, R => false, α => true, β => true), - :y => [y ~ I], - :y2 => [I], - :known_ic => Array{}[], - :to_check => Array{}[], - ), - ) - - @parameters θ - @variables t x(t) y(t) - D = Difference(t; dt = 1.0) - - eqs = [D(x) ~ θ * x^3] - - @named eqs = DiscreteSystem(eqs) - push!( - cases, - Dict( - :dds => eqs, - :res => OrderedDict(x => true, θ => true), - :y => [y ~ x], - :y2 => [x], - :known_ic => Array{}[], - :to_check => Array{}[], - ), - ) - - @parameters θ β - @variables t x1(t) x2(t) y(t) - D = Difference(t; dt = 1.0) - - eqs = [D(x1) ~ x1 + x2, D(x2) ~ θ + β] - - @named eqs = DiscreteSystem(eqs) - push!( - cases, - Dict( - :dds => eqs, - :res => OrderedDict(x1 => true, x2 => true, θ => false, β => false), - :y => [y ~ x1], - :y2 => [x1], - :known_ic => Array{}[], - :to_check => Array{}[], - ), - ) - - @parameters a b c d - @variables t x1(t) x2(t) u(t) y2(t) - D = Difference(t; dt = 1.0) - - eqs = [D(x1) ~ a * x1 - b * x1 * x2 + u, D(x2) ~ -c * x2 + d * x1 * x2] - - @named lv = DiscreteSystem(eqs) - push!( - cases, - Dict( - :dds => lv, - :res => OrderedDict( - x1 => true, - x2 => false, - a => true, - b => false, - c => true, - d => true, - ), - :y => [y ~ x1], - :y2 => [x1], - :known_ic => Array{}[], - :to_check => Array{}[], - ), - ) - - push!( - cases, - Dict( - :dds => lv, - :res => OrderedDict(b * x2 => true), - :y => [y ~ x1], - :y2 => [x1], - :known_ic => Array{}[], - :to_check => [b * x2], + :dds => dds, + :res => OrderedDict(a => true, b => false, c => false, b + c => true), + :known => :none, ), ) push!( cases, Dict( - :dds => lv, - :res => OrderedDict( - x1 => true, - x2 => true, - a => true, - b => true, - c => true, - d => true, - ), - :y => [y ~ x1, y2 ~ x1 / x2], - :y2 => [x1, x1 / x2], - :known_ic => Array{}[], - :to_check => Array{}[], + :dds => dds, + :res => OrderedDict(a => true, b => false, c => false, b + c => true), + :known => :all, ), ) - push!( - cases, - Dict( - :dds => lv, - :res => OrderedDict( - substitute(x1, Dict(t => 0)) => true, - substitute(x2, Dict(t => 0)) => true, - a => true, - b => true, - c => true, - d => true, - ), - :y => [y ~ x1], - :y2 => [x1], - :known_ic => [x2], - :to_check => Array{}[], - ), - ) + #--------------------- - # Example 1 from https://doi.org/10.1016/j.automatica.2008.03.019 - @parameters theta1 theta2 - @variables t x1(t) x2(t) u(t) y(t) - D = Difference(t; dt = 1.0) + cases = [] - eqs = [D(x1) ~ theta1 * x1 + x2, D(x2) ~ (1 - theta2) * x1 + x2^2 + u - x2] + dds = @DDSmodel(x1(t + 1) = a * x1(t) * x2(t), x2(t + 1) = b * x1(t), y(t) = x2(t)) - @named abmd1 = DiscreteSystem(eqs) push!( cases, Dict( - :dds => abmd1, - :res => OrderedDict(x1 => true, x2 => true, theta1 => true, theta2 => true), - :y => [y ~ x1], - :y2 => [x1], - :known_ic => Array{}[], - :to_check => Array{}[], + :dds => dds, + :res => + OrderedDict(a => true, b => false, x1 => false, x2 => true, b * x1 => true), + :known => :none, ), ) - # Example 2 from https://doi.org/10.1016/j.automatica.2008.03.019 - @parameters theta1 theta2 theta3 - @variables t x1(t) x2(t) u(t) y(t) y2(t) - D = Difference(t; dt = 1.0) - - eqs = [D(x1) ~ theta1 * x1^2 + theta2 * x2 + u - x1, D(x2) ~ theta3 * x1 - x2] - - @named abmd2 = DiscreteSystem(eqs) push!( cases, Dict( - :dds => abmd2, - :res => OrderedDict( - x1 => true, - x2 => false, - theta1 => true, - theta2 => false, - theta3 => false, - ), - :y => [y ~ x1], - :y2 => [x1], - :known_ic => Array{}[], - :to_check => Array{}[], + :dds => dds, + :res => OrderedDict(a => true, b => true, x1 => true, x2 => true), + :known => :all, ), ) - push!( - cases, - Dict( - :dds => abmd2, - :res => Dict( - x1 => true, - x2 => true, - theta1 => true, - theta2 => true, - theta3 => true, - ), - :y => [y ~ x1, y2 ~ x2], - :y2 => [x1, x2], - :known_ic => Array{}[], - :to_check => Array{}[], - ), - ) - - @parameters a b - @variables t x1(t) y(t) - D = Difference(t; dt = 1.0) - eqs = [D(x1) ~ a] - - @named kic = DiscreteSystem(eqs) - push!( - cases, - Dict( - :dds => kic, - :res => OrderedDict(x1 => false, a => true, b => false), - :y => [y ~ x1 + b], - :y2 => [x1 + b], - :known_ic => Array{}[], - :to_check => Array{}[], - ), - ) - push!( - cases, - Dict( - :dds => kic, - :res => OrderedDict(substitute(x1, Dict(t => 0)) => true, a => true, b => true), - :y => [y ~ x1 + b], - :y2 => [x1 + b], - :known_ic => [x1], - :to_check => Array{}[], - ), - ) + #--------------------- for c in cases @test assess_local_identifiability( c[:dds]; - measured_quantities = c[:y], - known_ic = c[:known_ic], - funcs_to_check = c[:to_check], - ) == c[:res] - @test assess_local_identifiability( - c[:dds]; - measured_quantities = c[:y2], - known_ic = c[:known_ic], - funcs_to_check = c[:to_check], + funcs_to_check = collect(keys(c[:res])), + known_ic = c[:known], ) == c[:res] end end diff --git a/test/local_identifiability_discrete_aux.jl b/test/local_identifiability_discrete_aux.jl index 2a172b1a2..b2566719e 100644 --- a/test/local_identifiability_discrete_aux.jl +++ b/test/local_identifiability_discrete_aux.jl @@ -1,7 +1,7 @@ @testset "Discrete local identifiability, internal function" begin cases = [] - dds = @ODEmodel(a'(t) = (b + c) * a(t) + 1, y(t) = a(t)) + dds = @DDSmodel(a(t + 1) = (b + c) * a(t) + 1, y(t) = a(t)) push!( cases, @@ -23,7 +23,7 @@ #--------------------- - dds = @ODEmodel(a'(t) = b(t) * a(t) + c, b'(t) = d * a(t), y(t) = b(t)) + dds = @DDSmodel(a(t + 1) = b(t) * a(t) + c, b(t + 1) = d * a(t), y(t) = b(t)) push!( cases, @@ -62,7 +62,7 @@ # ------------------- # Example 4 from https://doi.org/10.1016/j.automatica.2016.01.054 - dds = @ODEmodel(x'(t) = theta^3 * x(t), y(t) = x(t)) + dds = @DDSmodel(x(t + 1) = theta^3 * x(t), y(t) = x(t)) push!( cases, diff --git a/test/local_identifiability_me.jl b/test/local_identifiability_me.jl index 01365120c..07065365d 100644 --- a/test/local_identifiability_me.jl +++ b/test/local_identifiability_me.jl @@ -29,12 +29,13 @@ function _linear_compartment_model(graph, sinks) push!(edges_vars_names, "b_$(s)_0") end R, vars = - Nemo.PolynomialRing(Nemo.QQ, vcat(x_vars_names, edges_vars_names, ["y1", "y2"])) + Nemo.polynomial_ring(Nemo.QQ, vcat(x_vars_names, edges_vars_names, ["y1", "y2"])) x_vars = vars[2:(n + 1)] x0 = vars[1] - equations = Dict{fmpq_mpoly, Union{fmpq_mpoly, Generic.Frac{fmpq_mpoly}}}( - x => R(0) for x in x_vars - ) + equations = + Dict{QQMPolyRingElem, Union{QQMPolyRingElem, Generic.Frac{QQMPolyRingElem}}}( + x => R(0) for x in x_vars + ) equations[x0] = R(0) for i in 1:n for j in graph[i] @@ -52,10 +53,10 @@ function _linear_compartment_model(graph, sinks) equations[x_vars[i]] += -x_vars[i] * rate end end - return ODE{fmpq_mpoly}( + return ODE{QQMPolyRingElem}( equations, Dict(vars[end] => x_vars[1], vars[end - 1] => x0), - Array{fmpq_mpoly, 1}(), + Array{QQMPolyRingElem, 1}(), ) end diff --git a/test/monomial_compress.jl b/test/monomial_compress.jl index 74eed4c4d..697d1525c 100644 --- a/test/monomial_compress.jl +++ b/test/monomial_compress.jl @@ -1,5 +1,5 @@ @testset "Monomial compression test" begin - R, (v1, v2, v3, v4) = Nemo.PolynomialRing(Nemo.QQ, ["v1", "v2", "v3", "v4"]) + R, (v1, v2, v3, v4) = Nemo.polynomial_ring(Nemo.QQ, ["v1", "v2", "v3", "v4"]) tests = [ v1 + v2 - 2, v4 + v4 * v1 - 3 * v2 * v4, diff --git a/test/mtk_compat.jl b/test/mtk_compat.jl deleted file mode 100644 index ccbdaf287..000000000 --- a/test/mtk_compat.jl +++ /dev/null @@ -1,397 +0,0 @@ -@testset "Check identifiability of `ODESystem` object" begin - @parameters a01 a21 a12 - @variables t x0(t) x1(t) y1(t) [output = true] - D = Differential(t) - - eqs = [D(x0) ~ -(a01 + a21) * x0 + a12 * x1, D(x1) ~ a21 * x0 - a12 * x1, y1 ~ x0] - de = ODESystem(eqs, t, name = :Test) - - correct = OrderedDict( - a01 => :nonidentifiable, - a21 => :nonidentifiable, - a12 => :nonidentifiable, - x0 => :globally, - x1 => :nonidentifiable, - ) - - @test isequal(correct, assess_identifiability(de; measured_quantities = [y1 ~ x0])) - @test isequal(correct, assess_identifiability(de; measured_quantities = [x0])) - @test isequal( - correct, - assess_identifiability(de; measured_quantities = [(y1 ~ x0).rhs]), - ) - - # check identifiabile functions - correct = [a01 * a12, a01 + a12 + a21] - result = find_identifiable_functions(de, measured_quantities = [y1 ~ x0]) - @test isequal(Set(correct), Set(result)) - - # -------------------------------------------------------------------------- - - # check identifiabile functions - @parameters V_m k_m k01 c - @variables t x(t) y1(t) [output = true] - D = Differential(t) - - eqs = [D(x) ~ (-V_m * x) / (k_m + x) + k01 * x, y1 ~ c * x] - de = ODESystem(eqs, t, name = :Test) - - correct = [k01, c * k_m, V_m * c] - result = find_identifiable_functions(de) - @test isequal(Set(correct), Set(result)) - - correct = [k01, c * x, k_m * c, V_m * c] - result = find_identifiable_functions(de, with_states = true) - @test isequal(Set(correct), Set(result)) - - # -------------------------------------------------------------------------- - @parameters a01 a21 a12 - @variables t x0(t) x1(t) y1(t) [output = true] - D = Differential(t) - - eqs = [D(x0) ~ -(a01 + a21) * x0 + a12 * x1, D(x1) ~ a21 * x0 - a12 * x1, y1 ~ x0] - de = ODESystem(eqs, t, name = :Test) - - correct = OrderedDict( - a01 => :nonidentifiable, - a21 => :nonidentifiable, - a12 => :nonidentifiable, - x0 => :globally, - x1 => :nonidentifiable, - ) - - @test isequal(correct, assess_identifiability(de)) - - # -------------------------------------------------------------------------- - - @parameters a01 a21 a12 - @variables t x0(t) x1(t) y1(t) [output = true] - D = Differential(t) - - eqs = [D(x0) ~ -(a01 + a21) * x0 + a12 * x1, D(x1) ~ a21 * x0 - a12 * x1, y1 ~ x0] - de = ODESystem(eqs, t, name = :Test) - funcs_to_check = [a01, a21, a12, a01 * a12, a01 + a12 + a21] - correct = OrderedDict( - a01 => :nonidentifiable, - a21 => :nonidentifiable, - a12 => :nonidentifiable, - a01 * a12 => :globally, - a01 + a12 + a21 => :globally, - ) - @test isequal(correct, assess_identifiability(de; funcs_to_check = funcs_to_check)) - - # -------------------------------------------------------------------------- - - @parameters a01 a21 a12 - @variables t x0(t) x1(t) y1(t) - D = Differential(t) - - eqs = [D(x0) ~ -(a01 + a21) * x0 + a12 * x1, D(x1) ~ a21 * x0 - a12 * x1] - measured_quantities = [y1 ~ x0] - de = ODESystem(eqs, t, name = :Test) - funcs_to_check = [a01, a21, a12, a01 * a12, a01 + a12 + a21] - correct = OrderedDict( - a01 => :nonidentifiable, - a21 => :nonidentifiable, - a12 => :nonidentifiable, - a01 * a12 => :globally, - a01 + a12 + a21 => :globally, - ) - @test isequal( - correct, - assess_identifiability( - de; - measured_quantities = measured_quantities, - funcs_to_check = funcs_to_check, - ), - ) - - # -------------------------------------------------------------------------- - @parameters μ bi bw a χ γ k - @variables t S(t) I(t) W(t) R(t) y(t) - - eqs = [ - D(S) ~ μ - bi * S * I - bw * S * W - μ * S + a * R, - D(I) ~ bw * S * W + bi * S * I - (γ + μ) * I, - D(W) ~ χ * (I - W), - D(R) ~ γ * I - (μ + a) * R, - ] - de = ODESystem(eqs, t, name = :TestSIWR) - measured_quantities = [y ~ k * I] - # check all parameters (default) - @test isequal( - true, - all( - values( - assess_local_identifiability(de; measured_quantities = measured_quantities), - ), - ), - ) - - # check specific parameters - funcs_to_check = [μ, bi, bw, a, χ, γ, γ + μ, k, S, I, W, R] - correct = OrderedDict(f => true for f in funcs_to_check) - @test isequal( - correct, - assess_local_identifiability( - de; - measured_quantities = measured_quantities, - funcs_to_check = funcs_to_check, - ), - ) - - # checking ME identifiability - funcs_to_check = [μ, bi, bw, a, χ, γ, γ + μ, k] - correct = OrderedDict(f => true for f in funcs_to_check) - @test isequal( - (correct, 1), - assess_local_identifiability( - de; - measured_quantities = measured_quantities, - funcs_to_check = funcs_to_check, - prob_threshold = 0.99, - type = :ME, - ), - ) - - # checking identifiabile functions - correct = [a, bw, χ, bi, k, γ, μ] - result = find_identifiable_functions(de, measured_quantities = measured_quantities) - @test isequal(Set(correct), Set(result)) - - # -------------------------------------------------------------------------- - @parameters mu bi bw a xi gm k - @variables t S(t) I(t) W(t) R(t) y(t) [output = true] - - eqs = [ - D(S) ~ mu - bi * S * I - bw * S * W - mu * S + a * R, - D(I) ~ bw * S * W + bi * S * I - (gm + mu) * I, - D(W) ~ xi * (I - W), - D(R) ~ gm * I - (mu + a) * R, - y ~ k * I, - ] - de = ODESystem(eqs, t, name = :TestSIWR) - # check all parameters (default) - @test isequal(true, all(values(assess_local_identifiability(de)))) - - @test isequal( - true, - all(values(assess_local_identifiability(de; measured_quantities = [y ~ k * I]))), - ) - - # check specific parameters - funcs_to_check = [mu, bi, bw, a, xi, gm, gm + mu, k, S, I, W, R] - correct = OrderedDict(f => true for f in funcs_to_check) - @test isequal( - correct, - assess_local_identifiability(de; funcs_to_check = funcs_to_check), - ) - - # checking ME identifiability - funcs_to_check = [mu, bi, bw, a, xi, gm, gm + mu, k] - correct = OrderedDict(f => true for f in funcs_to_check) - @test isequal( - (correct, 1), - assess_local_identifiability( - de; - funcs_to_check = funcs_to_check, - prob_threshold = 0.99, - type = :ME, - ), - ) - - # -------------------------------------------------------------------------- - @parameters mu bi bw a xi gm k - @variables t S(t) I(t) W(t) R(t) y(t) - - eqs = [ - D(S) ~ 2.0 * mu - bi * S * I - bw * S * W - mu * S + a * R, - D(I) ~ bw * S * W + bi * S * I - (gm + mu) * I, - D(W) ~ xi * (I - 0.6 * W), - D(R) ~ gm * I - (mu + a) * R, - ] - de = ODESystem(eqs, t, name = :TestSIWR) - measured_quantities = [y ~ 1.57 * I * k] - funcs_to_check = [mu, bi, bw, a, xi, gm, mu, gm + mu, k, S, I, W, R] - correct = OrderedDict(f => true for f in funcs_to_check) - @test isequal( - correct, - assess_local_identifiability( - de; - measured_quantities = measured_quantities, - funcs_to_check = funcs_to_check, - ), - ) - - # checking ME identifiability - funcs_to_check = [bi, bw, a, xi, gm, mu, gm + mu, k] - correct = OrderedDict(f => true for f in funcs_to_check) - @test isequal( - (correct, 1), - assess_local_identifiability( - de; - measured_quantities = measured_quantities, - funcs_to_check = funcs_to_check, - prob_threshold = 0.99, - type = :ME, - ), - ) - - # ---------- - - @parameters a01 a21 a12 - @variables t x0(t) x1(t) y1(t) - D = Differential(t) - using SpecialFunctions - - eqs = [ - D(x0) ~ -(a01 + a21) * SpecialFunctions.erfc(x0) + a12 * x1, - D(x1) ~ a21 * x0 - a12 * x1, - ] - - de = ODESystem(eqs, t, name = :Test) - measured_quantities = [y1 ~ x0] - funcs_to_check = [a01, a21, a12, a01 * a12, a01 + a12 + a21] - correct = Dict( - a01 => :nonidentifiable, - a21 => :nonidentifiable, - a12 => :nonidentifiable, - a01 * a12 => :globally, - a01 + a12 + a21 => :globally, - ) - @test_throws ArgumentError assess_identifiability( - de; - measured_quantities = measured_quantities, - funcs_to_check = funcs_to_check, - ) - # ---------- - @parameters a b c - @variables t x1(t) x2(t) y(t) - D = Differential(t) - - eqs = [D(x1) ~ -a * x1 + x2 * b / (x1 + b / (c^2 - x2)), D(x2) ~ x2 * c^2 + x1] - de = ODESystem(eqs, t, name = :Test) - measured_quantities = [y ~ x2] - correct = Dict(a => :globally, b => :globally, c => :locally) - to_check = [a, b, c] - @test isequal( - correct, - assess_identifiability( - de; - measured_quantities = measured_quantities, - funcs_to_check = to_check, - ), - ) - - # check identifiabile functions - result = find_identifiable_functions(de, measured_quantities = measured_quantities) - correct = [b, a, c^2] - @test isequal(Set(result), Set(correct)) - - # ---------- - @parameters a b - @variables t c(t) x1(t) x2(t) y1(t) y2(t) - D = Differential(t) - - eqs = - [D(x1) ~ -a * x1 + x2 * b / (x1 + b / (c^2 - x2)), D(x2) ~ x2 * c^2 + x1, D(c) ~ 0] - de = ODESystem(eqs, t, name = :Test) - measured_quantities = [y1 ~ x2, y2 ~ c] - correct = OrderedDict(a => :globally, b => :globally) - to_check = [a, b] - @test isequal( - correct, - assess_identifiability( - de; - measured_quantities = measured_quantities, - funcs_to_check = to_check, - ), - ) - - #---------------------------------- - # Composable models test (from https://github.com/SciML/StructuralIdentifiability.jl/issues/162) - @variables t - function rabbits_creator(; name) - ps = @parameters α = 1.5 - vars = @variables x(t) = 1.0 z(t) = 0.0 [input = true] - D = Differential(t) - equs = [D(x) ~ α^2 * x + z] - - ODESystem(equs, t, vars, ps; name = name) - end - - function wolves_creator(; name) - ps = @parameters δ = 3.0 - vars = @variables y(t) = 1.0 q(t) = 0.0 [input = true] - D = Differential(t) - equs = [D(y) ~ -δ * y + q] - - ODESystem(equs, t, vars, ps; name = name) - end - - function lotka_volterra_creator(; name) - @named wolves = wolves_creator() - @named rabbits = rabbits_creator() - - ps = @parameters β = 1.0 γ = 1.0 - D = Differential(t) - - eqs = [rabbits.z ~ -β * wolves.y * rabbits.x, wolves.q ~ γ * wolves.y * rabbits.x] - - ModelingToolkit.compose(ODESystem(eqs, t, [], ps; name = name), wolves, rabbits) - end - - function getbyname(sys, name) - println(name) - return first([ - v for v in vcat(states(sys), parameters(sys)) if - replace(string(v), "(t)" => "") == name - ]) - end - - @named ltk_mtk = lotka_volterra_creator() - simp_ltk_mtk = structural_simplify(ltk_mtk) - wolves₊δ = getbyname(simp_ltk_mtk, "wolves₊δ") - rabbits₊α = getbyname(simp_ltk_mtk, "rabbits₊α") - β = getbyname(simp_ltk_mtk, "β") - γ = getbyname(simp_ltk_mtk, "γ") - wolves₊y = getbyname(simp_ltk_mtk, "wolves₊y") - rabbits₊x = getbyname(simp_ltk_mtk, "rabbits₊x") - @variables y(t) - measured_quantities = [y ~ wolves₊y] - result = assess_identifiability(simp_ltk_mtk, measured_quantities = measured_quantities) - correct = Dict( - rabbits₊α => :locally, - γ => :nonidentifiable, - β => :globally, - wolves₊δ => :globally, - rabbits₊x => :nonidentifiable, - wolves₊y => :globally, - ) - @test Dict(result) == correct - - #---------------------------------- - - @variables t, x(t), y(t), z(t), w(t) - @parameters a - @named sys = ODESystem([D(x) ~ a * y], t, [x], [a]; observed = [y ~ z, z ~ x]) - measured_quantities = [w ~ x] - result = assess_identifiability(sys, measured_quantities = measured_quantities) - @test result[a] == :globally - - result = find_identifiable_functions(sys, measured_quantities = measured_quantities) - @test isequal(result, [a]) - - #---------------------------------- - - # Tensor definition case as reported in - # https://github.com/SciML/StructuralIdentifiability.jl/issues/178 - @variables t, x(t)[1:2], y(t)[1:2] - @parameters k1, k2 - - eqs = [D(x[1]) ~ -k1 * x[2], D(x[2]) ~ -k2 * x[1]] - - sys = ODESystem(eqs, t, name = :example_vector) - correct = OrderedDict(x[1] => true, x[2] => true, k1 => true, k2 => true) - @test assess_local_identifiability(sys, measured_quantities = [x[1], x[2]]) == correct -end diff --git a/test/ode_ps_solution.jl b/test/ode_ps_solution.jl index 177bec543..dcaa60739 100644 --- a/test/ode_ps_solution.jl +++ b/test/ode_ps_solution.jl @@ -1,45 +1,22 @@ -# using Profile - -function rand_poly(deg, vars) - result = 0 - indices = vcat(collect(1:length(vars)), collect(1:length(vars))) - monomials = [] - for d in 0:deg - for subs in StructuralIdentifiability.IterTools.subsets(indices, d) - push!(monomials, subs) - end - end - - for subs in monomials - monom = rand(-50:50) - for v_ind in subs - monom *= vars[v_ind] - end - result += monom - end - - return result -end - @testset "Power series solution for an ODE system" begin - R, (x, x_dot) = Nemo.PolynomialRing(Nemo.QQ, ["x", "x_dot"]) + R, (x, x_dot) = Nemo.polynomial_ring(Nemo.QQ, ["x", "x_dot"]) exp_t = ps_ode_solution( [x_dot - x], - Dict{fmpq_mpoly, fmpq}(x => 1), - Dict{fmpq_mpoly, Array{fmpq, 1}}(), + Dict{QQMPolyRingElem, QQFieldElem}(x => 1), + Dict{QQMPolyRingElem, Array{QQFieldElem, 1}}(), 20, )[x] @test valuation(ps_diff(exp_t) - exp_t) == 19 R, (x, y, x_dot, y_dot, u) = - Nemo.PolynomialRing(Nemo.QQ, ["x", "y", "x_dot", "y_dot", "u"]) + Nemo.polynomial_ring(Nemo.QQ, ["x", "y", "x_dot", "y_dot", "u"]) prec = 100 eqs = [x_dot - x + 3 * x * y - u, y_dot + 2 * y - 4 * x * y] u_coeff = [rand(1:5) for i in 1:prec] sol = ps_ode_solution(eqs, Dict(x => 0, y => -2), Dict(u => u_coeff), prec) @test map(e -> valuation(evaluate(e, [sol[v] for v in gens(R)])), eqs) == [prec - 1, prec - 1] - F = Nemo.GF(2^31 - 1) + F = Nemo.Native.GF(2^31 - 1) prec = 100 # Testing the function ps_ode_solution by itself @@ -52,7 +29,7 @@ end ["x_$i" for i in 1:NUMX], ["u_$i" for i in 1:NUMU], ) - R, vars = Nemo.PolynomialRing(F, varnames) + R, vars = Nemo.polynomial_ring(F, varnames) # Generating the initial conditions and inputs ic = Dict(vars[i + NUMX] => F(rand(-5:5)) for i in 1:NUMX) @@ -87,8 +64,8 @@ end ["p_$i" for i in 1:NUMP], ["u_$i" for i in 1:NUMU], ) - R, vars = Nemo.PolynomialRing(F, varnames) - PType = gfp_mpoly + R, vars = Nemo.polynomial_ring(F, varnames) + PType = fpMPolyRingElem TDict = Dict{PType, Union{PType, Generic.Frac{PType}}} # Generating the intial conditions etc diff --git a/test/paradigm_shift.jl b/test/paradigm_shift.jl index f59b0b7da..5c3125c6b 100644 --- a/test/paradigm_shift.jl +++ b/test/paradigm_shift.jl @@ -22,7 +22,7 @@ function test_reparametrization(old_ode, new_ode, var_mapping, implicit_relation ic_point = map(Nemo.QQ, rand(1:bound, length(old_vars))) old_var_ic = Dict(old_vars .=> ic_point) input_ts = map(_ -> [Nemo.QQ(rand(1:bound))], 1:length(old_inputs)) - old_input_ts = Dict{eltype(old_vars), Vector{Nemo.fmpq}}(old_inputs .=> input_ts) + old_input_ts = Dict{eltype(old_vars), Vector{Nemo.QQFieldElem}}(old_inputs .=> input_ts) old_solutions = StructuralIdentifiability.power_series_solution( old_ode, @@ -47,7 +47,7 @@ function test_reparametrization(old_ode, new_ode, var_mapping, implicit_relation merge(old_param_spec, old_var_ic), ) for new_var in new_vars ) - new_input_ts = Dict{eltype(new_vars), Vector{Nemo.fmpq}}( + new_input_ts = Dict{eltype(new_vars), Vector{Nemo.QQFieldElem}}( new_input => old_input_ts[numerator(var_mapping[new_input])] for new_input in new_inputs ) diff --git a/test/parent_ring_change.jl b/test/parent_ring_change.jl index fcb952c75..327d2cae8 100644 --- a/test/parent_ring_change.jl +++ b/test/parent_ring_change.jl @@ -1,9 +1,9 @@ @testset "Parent ring change" begin - for field in [Nemo.QQ, Nemo.GF(2^31 - 1)] - R, (z, x, y) = PolynomialRing(QQ, ["z", "x", "y"], ordering = :degrevlex) - R_, (y_, x_) = PolynomialRing(QQ, ["y", "x"], ordering = :lex) + for field in [Nemo.QQ, Nemo.Native.GF(2^31 - 1)] + R, (z, x, y) = polynomial_ring(QQ, ["z", "x", "y"], ordering = :degrevlex) + R_, (y_, x_) = polynomial_ring(QQ, ["y", "x"], ordering = :lex) R__, (x__, t__, y__, z__) = - PolynomialRing(QQ, ["x", "t", "y", "z"], ordering = :deglex) + polynomial_ring(QQ, ["x", "t", "y", "z"], ordering = :deglex) f = 2x + 3y + x^7 * y @test f == StructuralIdentifiability.parent_ring_change(f, R, matching = :byname) @@ -37,9 +37,9 @@ @test f__ == StructuralIdentifiability.parent_ring_change(f, R__, matching = :byindex) - R, (x,) = PolynomialRing(field, ["x"]) - R2, (x2,) = PolynomialRing(Nemo.FractionField(R), ["x"]) - R3, (x3,) = PolynomialRing(field, ["x"]) + R, (x,) = polynomial_ring(field, ["x"]) + R2, (x2,) = polynomial_ring(Nemo.fraction_field(R), ["x"]) + R3, (x3,) = polynomial_ring(field, ["x"]) f = x3^2 f_ = StructuralIdentifiability.parent_ring_change(f, R2) @test parent(f) == R3 diff --git a/test/ps_diff.jl b/test/ps_diff.jl index 99f56801a..974459614 100644 --- a/test/ps_diff.jl +++ b/test/ps_diff.jl @@ -1,5 +1,5 @@ @testset "Power Series Differentiation" begin - T, t = Nemo.PowerSeriesRing(Nemo.QQ, 300, "t"; model = :capped_absolute) + T, t = Nemo.power_series_ring(Nemo.QQ, 300, "t"; model = :capped_absolute) @test ps_diff(zero(T)) == zero(T) diff --git a/test/ps_integrate.jl b/test/ps_integrate.jl index 53b1968ef..0bfc5c457 100644 --- a/test/ps_integrate.jl +++ b/test/ps_integrate.jl @@ -1,5 +1,5 @@ @testset "Power series integration" begin - T, t = Nemo.PowerSeriesRing(Nemo.QQ, 300, "t"; model = :capped_absolute) + T, t = Nemo.power_series_ring(Nemo.QQ, 300, "t"; model = :capped_absolute) @test ps_integrate(zero(T)) == zero(T) diff --git a/test/ps_inverse.jl b/test/ps_inverse.jl index c5a95d124..0f9aed191 100644 --- a/test/ps_inverse.jl +++ b/test/ps_inverse.jl @@ -1,8 +1,9 @@ @testset "Power series matrix inverse" begin - T, t = Nemo.PowerSeriesRing(Nemo.GF(2^31 - 1), 50, "t"; model = :capped_absolute) + T, t = + Nemo.power_series_ring(Nemo.Native.GF(2^31 - 1), 50, "t"; model = :capped_absolute) for d in 1:5 - S = Nemo.MatrixSpace(T, d, d) + S = Nemo.matrix_space(T, d, d) for case in 1:20 M = S([random_ps(T) for i in 1:d, j in 1:d]) while isequal( diff --git a/test/ps_matrix_homlinear.jl b/test/ps_matrix_homlinear.jl index aa8ab8230..ebdd9fc0b 100644 --- a/test/ps_matrix_homlinear.jl +++ b/test/ps_matrix_homlinear.jl @@ -1,11 +1,12 @@ @testset "Homogeneous linear differential equations" begin - T, t = Nemo.PowerSeriesRing(Nemo.GF(2^31 - 1), 300, "t"; model = :capped_absolute) + T, t = + Nemo.power_series_ring(Nemo.Native.GF(2^31 - 1), 300, "t"; model = :capped_absolute) for d in 1:5 for c in 1:5 - S = Nemo.MatrixSpace(T, d, d) + S = Nemo.matrix_space(T, d, d) A = random_ps_matrix(T, S) - Sconst = Nemo.MatrixSpace(Nemo.GF(2^31 - 1), d, d) + Sconst = Nemo.matrix_space(Nemo.Native.GF(2^31 - 1), d, d) Y0 = Sconst([rand(Int) % 100 for i in 1:d, j in 1:d]) while StructuralIdentifiability.LinearAlgebra.det(Y0) == 0 Y0 = Sconst([rand(Int) % 100 for i in 1:d, j in 1:d]) diff --git a/test/ps_matrix_linear.jl b/test/ps_matrix_linear.jl index e681c54d6..410e2cf76 100644 --- a/test/ps_matrix_linear.jl +++ b/test/ps_matrix_linear.jl @@ -1,12 +1,13 @@ @testset "Linear differential equations" begin - T, t = Nemo.PowerSeriesRing(Nemo.GF(2^31 - 1), 300, "t"; model = :capped_absolute) + T, t = + Nemo.power_series_ring(Nemo.Native.GF(2^31 - 1), 300, "t"; model = :capped_absolute) for d in 1:5 for c in 1:5 - S = Nemo.MatrixSpace(T, d, d) + S = Nemo.matrix_space(T, d, d) A = random_ps_matrix(T, S) B = random_ps_matrix(T, S) - Sconst = Nemo.MatrixSpace(Nemo.GF(2^31 - 1), d, d) + Sconst = Nemo.matrix_space(Nemo.Native.GF(2^31 - 1), d, d) Y0 = Sconst([rand(Int) % 100 for i in 1:d, j in 1:d]) @time sol = ps_matrix_linear_de(A, B, Y0) to_be_zero = map(ps_diff, sol) - A * sol - B diff --git a/test/ps_matrix_log.jl b/test/ps_matrix_log.jl index 5aa089c3f..6a32943ab 100644 --- a/test/ps_matrix_log.jl +++ b/test/ps_matrix_log.jl @@ -1,9 +1,9 @@ @testset "Logarith of power series matrices" begin - T, t = Nemo.PowerSeriesRing(Nemo.QQ, 300, "t"; model = :capped_absolute) + T, t = Nemo.power_series_ring(Nemo.QQ, 300, "t"; model = :capped_absolute) for d in 1:5 diag_elements = [1 + t * random_ps(T) for i in 1:d] - S = Nemo.MatrixSpace(T, d, d) + S = Nemo.matrix_space(T, d, d) M = S([(i == j ? diag_elements[i] : T(0)) for i in 1:d, j in 1:d]) result = ps_matrix_log(M) correct = S([(i == j ? log(diag_elements[i]) : T(0)) for i in 1:d, j in 1:d]) @@ -12,7 +12,7 @@ for c in 1:5 f = random_ps(T) * t - S = Nemo.MatrixSpace(T, 2, 2) + S = Nemo.matrix_space(T, 2, 2) m = S([ T(1) f T(0) T(1) @@ -26,7 +26,7 @@ for c in 1:5 f, g = random_ps(T) * t, random_ps(T) * t - S = Nemo.MatrixSpace(T, 3, 3) + S = Nemo.matrix_space(T, 3, 3) m = S([ T(1) f g T(0) T(1) f diff --git a/test/pseudodivision.jl b/test/pseudodivision.jl index 7c9c5d542..895300ece 100644 --- a/test/pseudodivision.jl +++ b/test/pseudodivision.jl @@ -1,5 +1,5 @@ @testset "Pseudodivision" begin - R, (x, y, z) = Nemo.PolynomialRing(Nemo.QQ, ["x", "y", "z"]) + R, (x, y, z) = Nemo.polynomial_ring(Nemo.QQ, ["x", "y", "z"]) @test iszero(pseudodivision(x^2 - y^2, x - y, x)) @test pseudodivision(y * x^4 + x^3 + x, z * x^2, x) == x diff --git a/test/runtests.jl b/test/runtests.jl index 82e6b0dd7..a2e3dc781 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,7 +5,6 @@ using TestSetExtensions using StructuralIdentifiability.DataStructures using StructuralIdentifiability.Nemo -using StructuralIdentifiability.ModelingToolkit using StructuralIdentifiability: field_contains, check_identifiability, @@ -58,7 +57,24 @@ using StructuralIdentifiability: extract_coefficients_ratfunc, lie_derivative, states_generators, - RationalFunctionField + RationalFunctionField, + x_vars, + y_vars, + x_equations, + y_equations, + inputs + +const GROUP = get(ENV, "GROUP", "All") + +@static if VERSION >= v"1.10.0" + if GROUP == "All" || GROUP == "ModelingToolkitExt" + using Pkg + Pkg.add("ModelingToolkit") + Pkg.add("Symbolics") + using ModelingToolkit + using Symbolics + end +end function random_ps(ps_ring, range = 1000) result = zero(ps_ring) @@ -77,11 +93,59 @@ function random_ps_matrix(ps_ring, matrix_space) return result end +function rand_poly(deg, vars) + result = 0 + indices = vcat(collect(1:length(vars)), collect(1:length(vars))) + monomials = [] + for d in 0:deg + for subs in StructuralIdentifiability.IterTools.subsets(indices, d) + push!(monomials, subs) + end + end + + for subs in monomials + monom = rand(-50:50) + for v_ind in subs + monom *= vars[v_ind] + end + result += monom + end + + return result +end + +function get_test_files(group) + result = Vector{String}() + for (dir, _, files) in walkdir("./") + for fname in files + if fname != "runtests.jl" && endswith(fname, ".jl") + if group == "All" || + (group == "Core" && dir != "./extensions") || + ( + group == "ModelingToolkitExt" && + dir == "./extensions" && + VERSION >= v"1.10.0" + ) + push!(result, dir * "/" * fname) + end + end + end + end + return result +end + @info "Testing started" @test isempty(Test.detect_ambiguities(StructuralIdentifiability)) @test isempty(Test.detect_unbound_args(StructuralIdentifiability)) +all_tests = get_test_files(GROUP) +if !isempty(ARGS) + all_tests = ARGS +end + @time @testset "All the tests" verbose = true begin - @includetests ARGS + for test_file in all_tests + include(test_file) + end end diff --git a/test/sequence_solution.jl b/test/sequence_solution.jl index 5a21b7956..5ebeba3b2 100644 --- a/test/sequence_solution.jl +++ b/test/sequence_solution.jl @@ -1,19 +1,19 @@ @testset "Sequence solutions in the discrete case" begin # Fibonacci example - ode = @ODEmodel(f1'(t) = f1(t) + f0(t), f0'(t) = f1(t), y(t) = f1(t)) + ode = @DDSmodel(f1(t + 1) = f1(t) + f0(t), f0(t + 1) = f1(t), y(t) = f1(t)) sol = sequence_solution( ode, - Dict{fmpq_mpoly, Int}(), + Dict{QQMPolyRingElem, Int}(), Dict(f0 => 1, f1 => 1), - Dict{fmpq_mpoly, Array{Int, 1}}(), + Dict{QQMPolyRingElem, Array{Int, 1}}(), 10, ) @test sol[f1] == [1, 2, 3, 5, 8, 13, 21, 34, 55, 89] @test sol[f0] == [1, 1, 2, 3, 5, 8, 13, 21, 34, 55] # Sum of powers - ode = @ODEmodel(a'(t) = 2 * a(t) + b * u(t), y(t) = a(t)) + ode = @DDSmodel(a(t + 1) = 2 * a(t) + b * u(t), y(t) = a(t)) sol = sequence_solution( ode, Dict(b => 3),