ALlow XDESystems, expand possibilities for LHS of Equation (#7)

* allow DESystems in processes_to_mtkmodel

* allow Differential and differential*p in LHS of equations as processes

* add changelog entry

* add warning comment

* fix type of vector expansion container

* show that differential can be used in the docs

* add new info to docs

* correct version parse

Author: Datseris

CHANGELOG.md

@@ -1,4 +1,11 @@
 # Changelog
 
 ProcessBasedModelling.jl follows semver 2.0.
-Changelog is kept with respect to v1 release.
+Changelog is kept with respect to v1 release.
+
+## 1.1
+
+- New keyword `warn_default` in `processes_to_mtkmodel`.
+- Now `processes_to_mtkmodel` allows for `XDESystems` as elements of the input vector.
+- Detection of LHS-variable in `Equation` has been improved significantly.
+  Now, LHS can be any of: `x`, `Differential(t)(x)`, `Differential(t)(x)*p`, `p*Differential(t)(x)`.

Project.toml

@@ -1,7 +1,7 @@
 name = "ProcessBasedModelling"
 uuid = "ca969041-2cf3-4b10-bc21-86f4417093eb"
 authors = ["Datseris <[email protected]>"]
-version = "1.0.5"
+version = "1.1.0"
 
 [deps]
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"

docs/src/index.md

@@ -76,22 +76,25 @@ More variables than equations, here are the potential extra variable(s):
 The error message is unhelpful as all variables are reported as "potentially missing".
 At least on the basis of our scientific reasoning however, both ``x, z`` have an equation.
 It is ``y`` that ``x`` introduced that does not have an equation.
-Moreover, in our experience these error messages become increasingly less useful when a model has many equations and/or variables, as many variables get cited as "missing" from the variable map even when only one should be.
+Moreover, in our experience these error messages become increasingly less accurate or helpful when a model has many equations and/or variables.
 This makes it difficult to quickly find out where the "mistake" happened in the equations.
 
 **PBM** resolves these problems and always gives accurate error messages when
 it comes to the construction of the system of equations.
 This is because on top of the variable map that MTK constructs automatically, **PBM** requires the user to implicitly provide a map of variables to processes that govern said variables. **PBM** creates the map automatically, the only thing the user has to do is to define the equations in terms of what [`processes_to_mtkmodel`](@ref) wants (which are either [`Process`](@ref)es or `Equation`s as above).
 
+For the majority of cases, **PBM** can infer the LHS variable a process "defines" automatically, just by passing in a vector of `Equation`s, like in **MTK**.
+For cases where this is not possible a dedicated `Process` type is provided, whose subtypes act as wrappers around equations providing some additional conveniences.
+
 Here is what the user defines to make the same system of equations via **PBM**:
 
 ```@example MAIN
 using ProcessBasedModelling
 
 processes = [
-    ExpRelaxation(z, x^2),      # introduces x variable
-    TimeDerivative(x, 0.1*y),   # introduces y variable
-    y ~ z - x,                  # can be an equation because LHS is single variable
+    ExpRelaxation(z, x^2),      # defines z, introduces x; `Process` subtype (optional)
+    Differential(t)(x) ~ 0.1*y, # defines x, introduces y; normal `Equation`
+    y ~ z - x,                  # defines y; normal `Equation`
 ]
 ```
 
@@ -117,17 +120,20 @@ there is no default process for x(t), and (t)x doesn't have a default value.
 Please provide a process for variable x(t).
 ```
 
-If instead we "forgot" the ``y`` process, **PBM** will not error, but warn, and make ``y`` equal to a named parameter, since ``y`` has a default value:
+If instead we "forgot" the ``y`` process, **PBM** will not error, but warn, and make ``y`` equal to a named parameter, since ``y`` has a default value.
+So, running:
 ```@example MAIN
 model = processes_to_mtkmodel(processes[1:2])
 equations(model)
 ```
 
+Makes the named parameter:
+
 ```@example MAIN
 parameters(model)
 ```
 
-and the warning thrown was:
+and throws the warning:
 ```julia
 ┌ Warning: Variable y(t) was introduced in process of variable x(t).
 │ However, a process for y(t) was not provided,

src/API.jl

@@ -51,8 +51,9 @@ timescale(::Process) = NoTimeDerivative()
 
 """
     ProcessBasedModelling.lhs(p::Process)
+    ProcessBasedModelling.lhs(eq::Equation)
 
-Return the left-hand-side of the equation that `p` represents as an `Expression`.
+Return the right-hand-side of the equation governing the process.
 If `timescale` is implemented for `p`, typically `lhs` does not need to be as well.
 See [`Process`](@ref) for more.
 """
@@ -70,24 +71,39 @@ function lhs(p::Process)
 end
 
 """
-    ProcessBasedModelling.rhs(p::Process)
+    ProcessBasedModelling.rhs(process)
 
-Return the right-hand-side of the equation that `p` represents as an `Expression`.
-See [`Process`](@ref) for more.
+Return the right-hand-side of the equation governing the process.
 """
 function rhs(p::Process)
     error("Right-hand side (`rhs`) is not defined for process $(nameof(typeof(p))).")
 end
 
 # Extensions for `Equation`:
 rhs(e::Equation) = e.rhs
-lhs(e::Equation) = lhs_variable(e)
-function lhs_variable(e::Equation)
-    x = e.lhs
-    # we first check whether `x` is a variable
-    if !is_variable(x)
-        throw(ArgumentError("In given equation $(e), the left-hand-side does "*
-        "not represent a single variable."))
+lhs(e::Equation) = e.lhs
+lhs_variable(e::Equation) = lhs_variable(lhs(e))
+lhs_variable(x::Num) = Num(lhs_variable(x.val))
+function lhs_variable(x) # basically x is SymbolicUtils.BasicSymbolic{Real}
+    # check whether `x` is a single variable already
+    if is_variable(x)
+        return x
+    end
+    # check Differential(t)(x)
+    if hasproperty(x, :f)
+        if x.f isa Differential
+            return x.arguments[1]
+        end
+    end
+    # check Differential(t)(x)*parameter
+    if hasproperty(x, :arguments)
+        args = x.arguments
+        di = findfirst(a -> hasproperty(a, :f) && a.f isa Differential, args)
+        if !isnothing(di)
+            return args[di].arguments[1]
+        end
     end
-    return x
+    # error if all failed
+    throw(ArgumentError("We analyzed the LHS `$(x)` "*
+    "but could not extract a single variable it represents."))
 end

src/make.jl

@@ -6,15 +6,22 @@ The model/system is _not_ structurally simplified.
 
 `processes` is a vector whose elements can be:
 
-- Any instance of a subtype of [`Process`](@ref).
-- An `Equation` which is of the form `variable ~ expression` with `variable` a single
-  variable resulting from an `@variables` call.
-- 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.
+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.
+   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)`.
+   Anything else will either error or fail unexpectedly.
+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.
+3. A ModelingToolkit.jl `XDESystem`, 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 already existing systems.
 
 `default` is a vector that can contain the first two possibilities only
-as it contains default processes that may be assigned to variables introduced in `processes`
-but they don't themselves have an assigned process.
+as it contains default processes that may be assigned to individual variables introduced in
+`processes` but they don't themselves have an assigned process.
 
 It is expected that downstream packages that use ProcessBasedModelling.jl to make a
 field-specific library implement a 1-argument version of `processes_to_mtkmodel`,
@@ -26,10 +33,10 @@ or provide a wrapper function for it, and add a default value for `default`.
 - `name = nameof(type)`: the name of the model
 - `independent = t`: the independent variable (default: `@variables t`).
   `t` is also exported by ProcessBasedModelling.jl for convenience.
-- `warn_defaut::Bool = true`: if `true`, throw a warning when a variable does not
+- `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.
 """
-function processes_to_mtkmodel(_processes, _default = [];
+function processes_to_mtkmodel(_processes::Vector, _default = [];
         type = ODESystem, name = nameof(type), independent = t, warn_default::Bool = true,
     )
     processes = expand_multi_processes(_processes)
@@ -86,19 +93,25 @@ function processes_to_mtkmodel(_processes, _default = [];
             end
         end
     end
+    # return eqs
     sys = type(eqs, independent; name)
     return sys
 end
 
 function expand_multi_processes(procs::Vector)
+    etypes = Union{Vector, ODESystem, SDESystem, PDESystem}
+    !any(p -> p isa etypes, procs) && return procs
     # Expand vectors of processes or ODESystems
-    !any(p -> p isa Vector, procs) && return procs
-    expanded = deepcopy(procs)
-    idxs = findall(p -> p isa Vector, procs)
+    expanded = Any[procs...]
+    idxs = findall(p -> p isa etypes, procs)
     multiprocs = expanded[idxs]
     deleteat!(expanded, idxs)
     for mp in multiprocs
-        append!(expanded, mp)
+        if mp isa Vector
+            append!(expanded, mp)
+        else # then it is XDE system
+            append!(expanded, equations(mp))
+        end
     end
     return expanded
 end
@@ -144,8 +157,8 @@ function nonunique(x::AbstractArray{T}) where T
     duplicatedset = Set{T}()
     duplicatedvector = Vector{T}()
     for i in x
-        if(i in uniqueset)
-            if !(i in duplicatedset)
+        if i ∈ uniqueset
+            if !(i ∈ duplicatedset)
                 push!(duplicatedset, i)
                 push!(duplicatedvector, i)
             end

test/runtests.jl

@@ -191,4 +191,21 @@ end
 @testset "addition process error" begin
     @variables x(t) y(t) z(t)
     @test_throws ArgumentError AdditionProcess(x ~ 0.1z, y ~ x^2)
+end
+
+@testset "ODESystem as process" begin
+    @variables z(t) = 0.0
+    @variables x(t) = 0.0
+    @variables y(t) = 0.0
+    @variables w(t) = 0.0
+    procs = [
+        ExpRelaxation(z, x^2, 1.0), # introduces x and y variables
+        TimeDerivative(x, 0.1*y),   # introduces y variable!
+        y ~ z-x,                    # is an equation, not a process!
+    ]
+
+    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]
 end