Check for wrong Equation type in RHS. (#11)

* Check for wrong `Equation` type in RHS.

* oops use `isa`

* use new ci file which will hopefully work?

* simlper ode test dependency

* change ODE to Tsit4

* only check RHS when keyword true

* remove obsolete arguments from mtkeqs function

* improve central docstring

* wrap all tests in a testset

* module outside tests

Author: Datseris

.github/workflows/ci.yml

@@ -34,7 +34,7 @@ jobs:
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
+      - uses: actions/cache@v4
         env:
           cache-name: cache-artifacts
         with:

CHANGELOG.md

@@ -3,6 +3,10 @@
 ProcessBasedModelling.jl follows semver 2.0.
 Changelog is kept with respect to v1 release.
 
+## 1.6
+
+- Added an additional step when constructing the raw equation vector to be passed into an MTK model. In this step it is also checked that the RHS for all equations is an `Expression`. Sometimes it is easy to get confused and mess up and make it be an `Equation` (i.e., assigning the LHS-variable twice). This now will give an informative error.
+
 ## 1.5
 
 - Add docstring to `processes_to_mtkeqs` and list it in the documentation.

Project.toml

@@ -1,7 +1,7 @@
 name = "ProcessBasedModelling"
 uuid = "ca969041-2cf3-4b10-bc21-86f4417093eb"
 authors = ["Datseris <datseris.george@gmail.com>"]
-version = "1.5.0"
+version = "1.6.0"
 
 [deps]
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"

src/make.jl

@@ -8,22 +8,22 @@ passed to the MTK model/system like `ODESystem`.
 
 During construction, the following automations improve user experience:
 
-- Variable(s) introduced in `processes` that does not itself have a process obtain
+- Variable(s) introduced in `processes` that does not themselves have a process obtain
   a default process from `default`.
-- If no default exists, but the variable(s) itself has a default numerical value,
-  a [`ParameterProcess`](@ref) is created for said variable and a warning is thrown.
+- If no default exists, but the variable(s) have a default numerical value,
+  a [`ParameterProcess`](@ref) is created for said variable(s) and a warning is thrown.
 - Else, an informative error is thrown.
 - An error is also thrown if any variable has two or more processes assigned to it.
 
 `processes` is a `Vector` whose elements can be:
 
-1. Any instance of a subtype of [`Process`](@ref). `Process` is a
+1. Any instance of a subtype of [`Process`](@ref). `Process` is like 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.
    Let `x` be a `@variable` and `p` be a `@parameter`. Then, the LHS can only be one of:
-   `x`, `Differential(t)(x)`, `Differential(t)(x)*p`, `p*Differential(t)(x)`,
-   however, the versions with `p` may fail unexpectedly. Anything else will error.
+   `x`, `Differential(t)(x)`, `Differential(t)(x)*p`, or `p*Differential(t)(x)`.
+   The versions with `p` may fail unexpectedly. Anything else will error.
 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).
@@ -43,7 +43,7 @@ modelling libraries based on ProcessBasedModelling.jl is to define modules/submo
 that offer a pool of pre-defined variables and processes.
 Modules may register their own default processes via the function
 [`register_default_process!`](@ref).
-These registered processes are used when `default` is a `Module`.
+These registered default processes are used when `default` is a `Module`.
 
 ## Keyword arguments
 
@@ -53,6 +53,10 @@ These registered processes are used when `default` is a `Module`.
   `t` is also exported by ProcessBasedModelling.jl for convenience.
 - `warn_default::Bool = true`: if `true`, throw a warning when a variable does not
   have an assigned process but it has a default value so that it becomes a parameter instead.
+- `check_rhs::Bool = true`: if `true`, check that the RHS of all processes is NOT an
+  `Equation` type. Throw an informative error if there is one. This
+  helps scenarios where the RHS is wrongly an `Equation` assigning the LHS itself
+  (has happened to me many times!).
 """
 function processes_to_mtkmodel(args...;
         type = ODESystem, name = nameof(type), independent = t, kw...,
@@ -78,9 +82,10 @@ processes_to_mtkeqs(procs, default_dict(v); kw...)
 # The main implementation has the defaults to be a map from variable to process
 # because this simplifies a bit the code
 function processes_to_mtkeqs(_processes::Vector, default::Dict{Num, Any};
-        type = ODESystem, name = nameof(type), independent = t, warn_default::Bool = true,
+        warn_default::Bool = true, check_rhs::Bool = true,
     )
     processes = expand_multi_processes(_processes)
+    check_rhs && check_rhs_validity(processes)
     # Setup: obtain lhs-variables so we can track new variables that are not
     # in this vector. The vector has to be of type `Num`
     lhs_vars = Num[lhs_variable(p) for p in processes]
@@ -154,6 +159,18 @@ function expand_multi_processes(procs::Vector)
     return expanded
 end
 
+function check_rhs_validity(processes)
+    for p in processes
+        if rhs(p) isa Equation
+            lvar = lhs_variable(p)
+            throw(ArgumentError("Process assigned to variable $(lvar) is ill defined. "*
+                "The RHS is an `<: Equation` type but it shouldn't be."
+            ))
+        end
+    end
+    return
+end
+
 function default_dict(processes)
     default = Dict{Num, Any}()
     for proc in processes

test/Project.toml

@@ -1,6 +1,6 @@
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/runtests.jl

@@ -1,7 +1,18 @@
 using ProcessBasedModelling
 using Test
-using OrdinaryDiffEq
+using OrdinaryDiffEqTsit5
 
+# This module is used in one of the tests
+module TestDefault
+    using ProcessBasedModelling
+    @variables x(t) = 0.5 y(t) = 0.2
+    register_default_process!.([
+        Differential(t)(x) ~ 0.2y - x,
+        y ~ x^2,
+    ], Ref(TestDefault))
+end
+
+@testset "ProcessBasedModelling" begin
 @testset "construction + evolution" begin
     # The model, as defined below, is bistable due to ice albedo feedback
     # so two initial conditions should go to two attractors
@@ -61,7 +72,7 @@ using OrdinaryDiffEq
     ufs = []
     for u0 in u0s
         p = ODEProblem(sys, u0, (0.0, 1000.0*365*24*60*60.0))
-        sol = solve(p, Tsit5())
+        sol = solve(p, Tsit5(); abstol = 1e-9, reltol = 1e-9)
         push!(ufs, sol.u[end])
     end
 
@@ -215,13 +226,16 @@ end
     @test sort(ModelingToolkit.getname.(unknowns(sys2))) == [:w, :x, :y, :z]
 end
 
-module TestDefault
-    using ProcessBasedModelling
-    @variables x(t) = 0.5 y(t) = 0.2
-    register_default_process!.([
-        Differential(t)(x) ~ 0.2y - x,
-        y ~ x^2,
-    ], Ref(TestDefault))
+@testset "equation in RHS" begin
+    @variables z(t) = 0.0
+    @variables x(t) = 0.0
+    @variables y(t) = 0.0
+    procs = [
+        ExpRelaxation(z, x^2, 1.0), # introduces x and y variables
+        y ~ z-x,                    # is an equation, not a process!
+        z ~ (z ~ x^2),
+    ]
+    @test_throws ["an `<: Equation` type"] processes_to_mtkeqs(procs)
 end
 
 @testset "registering default" begin
@@ -233,3 +247,6 @@ end
     @test has_symbolic_var(mtk, z)
     @test has_symbolic_var(mtk, TestDefault.x)
 end
+
+
+end # @testset