update to MTKv11

CHANGELOG.md

@@ -3,6 +3,10 @@
 ProcessBasedModelling.jl follows semver 2.0.
 Changelog is kept with respect to v1 release.
 
+## 1.9
+
+Updated to ModelingToolkit.jl v11, which also comes with a dependency to AGPL license.
+
 ## 1.8
 
 - Updated to ModelingToolkit.jl v10. `type` keyword in `processes_to_mtkmodel` is now no longer used.

Project.toml

@@ -1,13 +1,15 @@
 name = "ProcessBasedModelling"
 uuid = "ca969041-2cf3-4b10-bc21-86f4417093eb"
 authors = ["Datseris <datseris.george@gmail.com>"]
-version = "1.8.0"
+version = "1.9.0"
 
 [deps]
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [compat]
-ModelingToolkit = "10.0"
+ModelingToolkit = "11"
+Symbolics = "7"
 Reexport = "1.2"
 julia = "1.9.0"

README.md

@@ -6,7 +6,7 @@
 [![codecov](https://codecov.io/gh/JuliaDynamics/ProcessBasedModelling.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/JuliaDynamics/ProcessBasedModelling.jl)
 [![Package Downloads](https://shields.io/endpoint?url=https://pkgs.genieframework.com/api/v1/badge/ProcessBasedModelling)](https://pkgs.genieframework.com?packages=ProcessBasedModelling)
 
-ProcessBasedModelling.jl is an extension to [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) (MTK) for building a model of equations using symbolic expressions.
+ProcessBasedModelling.jl is an extension to [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) (and the wider MTK ecosystem) for building a model of equations using symbolic expressions.
 It is an alternative framework to MTK's [native component-based modelling](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/acausal_components/), but, instead of components, there are "processes".
 This approach is useful in the modelling of physical/biological/whatever systems, where each variable corresponds to a particular physical concept or observable and there are few (or none) duplicate variables to make the definition of MTK "factories" worthwhile.
 On the other hand, there plenty of different physical representations, or _processes_ to represent a given physical concept in equation form.

docs/Project.toml

@@ -4,4 +4,3 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
 DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634"
-ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"

docs/make.jl

@@ -15,4 +15,5 @@ pages =  [
 
 build_docs_with_style(pages, ProcessBasedModelling;
     authors = "George Datseris <datseris.george@gmail.com>",
+    warnonly = [:doctest, :linkcheck],
 )

docs/src/index.md

@@ -31,10 +31,10 @@ Let's say we want to build the system of equations
 y = z - x
 ```
 
-symbolically using ModelingToolkit.jl (**MTK**). We define
+symbolically using ModelingToolkit.jl (**MTKBase**). We define
 
 ```@example MAIN
-using ModelingToolkit
+using ProcessBasedModelling # re-exports ModelingToolkit
 
 @variables t # independent variable _without_ units
 @variables z(t) = 0.0
@@ -69,7 +69,7 @@ model = ODESystem(eqs[1:2], t; name = :example)
 
 ```julia
 # here is the error
-model = structural_simplify(model)
+model = mtkcompile(model)
 ```
 ```
 ERROR: ExtraVariablesSystemException: The system is unbalanced.

src/API.jl

@@ -86,7 +86,7 @@ end
 rhs(e::Equation) = e.rhs
 lhs(e::Equation) = e.lhs
 lhs_variable(e::Equation) = lhs_variable(lhs(e))
-lhs_variable(x::Num) = Num(lhs_variable(x.val))
+lhs_variable(x::Num) = Num(lhs_variable(Symbolics.unwrap(x)))
 function lhs_variable(x) # basically x is SymbolicUtils.BasicSymbolic{Real}
     # check whether `x` is a single variable already
     if is_variable(x)
@@ -95,15 +95,15 @@ function lhs_variable(x) # basically x is SymbolicUtils.BasicSymbolic{Real}
     # check Differential(t)(x)
     if hasproperty(x, :f)
         if x.f isa Differential
-            return x.arguments[1]
+            return x.args[1]
         end
     end
     # check Differential(t)(x)*parameter
-    if hasproperty(x, :arguments)
-        args = x.arguments
+    if hasproperty(x, :args)
+        args = x.args
         di = findfirst(a -> hasproperty(a, :f) && a.f isa Differential, args)
         if !isnothing(di)
-            return args[di].arguments[1]
+            return args[di].args[1]
         end
     end
     # error if all failed

src/ProcessBasedModelling.jl

@@ -11,6 +11,8 @@ using Reexport
 using ModelingToolkit: t_nounits as t, D_nounits as D
 @reexport using ModelingToolkit
 
+getname = ModelingToolkit.SymbolicIndexingInterface.getname
+
 include("API.jl")
 include("utils.jl")
 include("default.jl")

src/make.jl

@@ -18,7 +18,7 @@ During construction, the following automations improve user experience:
 
 `processes` is a `Vector` whose elements can be:
 
-1. Any instance of a subtype of [`Process`](@ref). `Process` is like a
+1. Any instance of a subtype of [`Process`](@ref). `Process` is a
    wrapper around `Equation` that provides some conveniences, e.g., handling of timescales
    or not having limitations on the left-hand-side (LHS) form.
 1. An `Equation`. The LHS format of the equation is limited.
@@ -28,7 +28,7 @@ During construction, the following automations improve user experience:
 2. A `Vector` of the above two, which is then expanded. This allows the convenience of
    functions representing a physical process that may require many equations to be defined
    (because e.g., they may introduce more variables).
-3. A ModelingToolkit.jl `XDESystem`, in which case the `equations` of the system are expanded
+3. A ModelingToolkit.jl `System`, in which case the `equations` of the system are expanded
    as if they were given as a vector of equations like above. This allows the convenience
    of straightforwardly coupling with already existing `XDESystem`s.
 
@@ -59,7 +59,7 @@ These registered default processes are used when `default` is a `Module`.
   (has happened to me many times!).
 """
 function processes_to_mtkmodel(args...;
-        type = System, name = :model, independent = t, kw...,
+        name = :model, independent = t, kw...,
     )
     eqs = processes_to_mtkeqs(args...; kw...)
     sys = System(eqs, independent; name)
@@ -122,7 +122,7 @@ function processes_to_mtkeqs(_processes::Vector, default::Dict{Num, Any};
                     However, a process for $(added_var) was not provided,
                     and there is no default process for it either.
                     Since it has a default value, we make it a parameter by adding a process:
-                    `ParameterProcess($(ModelingToolkit.getname(added_var)))`.
+                    `ParameterProcess($(getname(added_var)))`.
                     """)
                 end
                 parproc = ParameterProcess(added_var)

src/processes_basic.jl

@@ -111,7 +111,7 @@ struct AdditionProcess <: Process
     function AdditionProcess(process, added::Vector)
         for add in added
             if typeof(add) <: Union{Process, Equation}
-                v1, v2 = ModelingToolkit.getname(lhs_variable(process)), ModelingToolkit.getname(lhs_variable(add))
+                v1, v2 = getname(lhs_variable(process)), getname(lhs_variable(add))
                 if v1 ≠ v2
                     throw(ArgumentError(
                         "Added processes do not have the same lhs variable. Got: $(v1), $(v2)"

src/utils.jl

@@ -1,3 +1,5 @@
+import Symbolics
+
 """
     LiteralParameter(p)
 
@@ -17,8 +19,9 @@ _literalvalue(p::LiteralParameter) = p.p
     has_symbolic_var(eqs, var)
 
 Return `true` if symbolic variable `var` exists in the equation(s) `eq`, `false` otherwise.
-This works for either `@parameters` or `@variables`.
-If `var` is a `Symbol` isntead of a `Num`, all variables are converted to their names
+This works for either `@parameters` or `@variables` and currently only compares versus these
+objects (i.e., it ignores `Differentials` or other complex structures).
+If `var` is a `Symbol` instead of a `Num`, all variables are converted to their names
 and equality is checked on the basis of the name only.
 
     has_symbolic_var(model, var)
@@ -27,7 +30,8 @@ When given a MTK model (such as `ODESystem`) search in _all_ the equations of th
 including observed variables.
 """
 function has_symbolic_var(eq::Equation, var)
-    vars = get_variables(eq)
+    vars = get_variables(eq) # this includes differentials
+    vars = filter(x -> is_parameter(x) || is_variable(x), vars)
     return _has_thing(var, vars)
 end
 has_symbolic_var(eqs::Vector{Equation}, var) = any(eq -> has_symbolic_var(eq, var), eqs)
@@ -37,8 +41,8 @@ function _has_thing(var::Num, vars)
     return any(isequal(var), vars)
 end
 function _has_thing(var::Symbol, vars)
-    vars = ModelingToolkit.getname.(vars)
-    var = ModelingToolkit.getname(var)
+    vars = getname.(vars)
+    var = getname(var)
     return any(isequal(var), vars)
 end
 
@@ -55,34 +59,30 @@ all_equations(model) = vcat(equations(model), observed(model))
 Return the default value of a symbolic variable `x` or `nothing`
 if it doesn't have any. Return `x` if `x` is not a symbolic variable.
 The difference with `ModelingToolkit.getdefault` is that this function will
-not error on the absence of a default value.
+not error in the absence of a default value.
 """
 default_value(x) = x
-default_value(x::Num) = default_value(x.val)
-function default_value(x::ModelingToolkit.SymbolicUtils.Symbolic)
-    if haskey(x.metadata, ModelingToolkit.Symbolics.VariableDefaultValue)
-        return x.metadata[ModelingToolkit.Symbolics.VariableDefaultValue]
-    else
-        @warn("No default value assigned to variable/parameter $(x).")
-        return nothing
-    end
+default_value(x::Num) = default_value(Symbolics.unwrap(x))
+function default_value(x::Symbolics.SymbolicT)
+    val = Symbolics.getmetadata(x, Symbolics.VariableDefaultValue, nothing)
+    isnothing(val) && @warn("No default value assigned to variable/parameter $(x).")
+    return val
 end
+
 # trick to get default values for state variables:
 # Base.Fix1(getindex, ModelingToolkit.defaults(ssys)).(states(ssys))
 # while `defaults` returns all default assignments.
 
-is_variable(x::Num) = is_variable(x.val)
-function is_variable(x)
-    if x isa ModelingToolkit.SymbolicUtils.Symbolic
-        if isnothing(x.metadata)
-            return false
-        end
-        if haskey(x.metadata, ModelingToolkit.Symbolics.VariableSource)
-            src = x.metadata[ModelingToolkit.Symbolics.VariableSource]
-            return first(src) == :variables
-        end
-    end
-    return false
+is_variable(x::Num) = is_variable(Symbolics.unwrap(x))
+function is_variable(x::Symbolics.SymbolicT)
+    value = getmetadata(x, Symbolics.VariableSource, nothing)
+    return value isa Tuple{Symbol, Symbol} && value[1] == :variables
+end
+
+is_parameter(x::Num) = is_parameter(Symbolics.unwrap(x))
+function is_parameter(x::Symbolics.SymbolicT)
+    value = getmetadata(x, Symbolics.VariableSource, nothing)
+    return value isa Tuple{Symbol, Symbol} && value[1] == :parameters
 end
 
 """
@@ -111,7 +111,7 @@ Now `p` will be a parameter with name `:τ_x` and default value `0.5`.
 new_derived_named_parameter(v, value::Num, extra::String; kw...) = value
 new_derived_named_parameter(v, value::LiteralParameter, extra::String; kw...) = value.p
 function new_derived_named_parameter(v, value::Real, extra; connector = "_", prefix = true)
-    n = string(ModelingToolkit.getname(v))
+    n = string(Symbolics.getname(v))
     newstring = if prefix
         extra*connector*n
     else
@@ -134,7 +134,7 @@ Convert all variables `vars` into `@parameters` with name the same as `vars`
 and default value the same as the value of `vars`. The macro leaves unaltered
 inputs that are of type `Num`, assumming they are already parameters.
 It also replaces [`LiteralParameter`](@ref) inputs with its literal values.
-This macro is extremely useful to convert e.g., keyword arguments into named parameters,
+This macro is useful to convert e.g., keyword arguments into named parameters,
 while also allowing the user to give custom parameter names,
 or to leave some keywords as numeric literals.
 
@@ -174,7 +174,7 @@ macro convert_to_parameters(vars...)
                 $binding isa Num, $binding,
                 # Else, convert to modeling toolkit param.
                 # This syntax was obtained by doing @macroexpand @parameters A = 0.5
-                (ModelingToolkit.toparam)((Symbolics.wrap)((SymbolicUtils.setmetadata)((Symbolics.setdefaultval)((Sym){Real}($varname), $binding), Symbolics.VariableSource, (:parameters, $varname))))
+                (ModelingToolkit.toparam_validate)((Symbolics.wrap)((SymbolicUtils.setmetadata)((Symbolics.setdefaultval)((SymbolicUtils.Sym){SymbolicUtils.SymReal}($varname; type = Real), $binding), Symbolics.VariableSource, (:parameters, $varname))))
                 ))
             )
         )

test/runtests.jl

@@ -19,7 +19,7 @@ end
     # If that's the case, we are sure model construction was valid
 
     # First, make some default processes
-    @variables T(t) = 300.0       # temperature, in Kelvin
+    @variables T(t) = 300.0 # temperature, in Kelvin
     @variables α(t)         # albedo of ice, unitless
     @variables ε(t)         # effective emissivity, unitless
     solar_constant = 340.25 # W/m^2, already divided by 4
@@ -66,7 +66,7 @@ end
     @test length(unknowns(sys)) == 1
     @test has_symbolic_var(equations(sys), T)
 
-    u0s = [[300.0], [100.0]]
+    u0s = [[T => 300.0], [T => 100.0]]
     ufs = []
     for u0 in u0s
         p = ODEProblem(sys, u0, (0.0, 1000.0*365*24*60*60.0))
@@ -128,7 +128,7 @@ end
         @test length(unknowns(sys)) == 3
         sys = processes_to_mtkmodel(procs[1:3])
         @test length(unknowns(sys)) == 3
-        @test length(unknowns(structural_simplify(sys))) == 2
+        @test length(unknowns(mtkcompile(sys))) == 2
     end
 end
 
@@ -139,10 +139,10 @@ end
     @testset "derived" begin
         @variables x(t) = 0.5
         p = new_derived_named_parameter(x, 0.2, "t")
-        @test ModelingToolkit.getname(p) == :t_x
+        @test ModelingToolkit.SymbolicIndexingInterface.getname(p) == :t_x
         @test default_value(p) == 0.2
         p = new_derived_named_parameter(x, 0.2, "t"; prefix = false, connector = "")
-        @test ModelingToolkit.getname(p) == :xt
+        @test ModelingToolkit.SymbolicIndexingInterface.getname(p) == :xt
     end
 
     @testset "convert" begin
@@ -151,7 +151,7 @@ end
         @convert_to_parameters A B C
         @test A isa Num
         @test default_value(A) == 0.5
-        @test ModelingToolkit.getname(C) == :X
+        @test ModelingToolkit.SymbolicIndexingInterface.getname(C) == :X
     end
 
     @testset "literal in derived" begin
@@ -219,10 +219,10 @@ end
     sys = processes_to_mtkmodel(procs)
     sys2 = processes_to_mtkmodel([sys, w ~ x*y])
     @test length(equations(sys2)) == 4
-    @test sort(ModelingToolkit.getname.(unknowns(sys2))) == [:w, :x, :y, :z]
+    @test sort(ModelingToolkit.SymbolicIndexingInterface.getname.(unknowns(sys2))) == [:w, :x, :y, :z]
 end
 
-@testset "equation in RHS" begin
+@testset "duplcate equation in RHS" begin
     @variables z(t) = 0.0
     @variables x(t) = 0.0
     @variables y(t) = 0.0
@@ -231,7 +231,7 @@ end
         y ~ z-x,                    # is an equation, not a process!
         z ~ (z ~ x^2),
     ]
-    @test_throws ["an `<: Equation` type"] processes_to_mtkeqs(procs)
+    @test_throws ["more than one"] processes_to_mtkeqs(procs)
 end
 
 @testset "not actual process" begin