Move to MTK v9 and prepare v1 tag (#4)

* bump version of MTK

* rename states to unknowns

* make default choices more fitting for dynsys

* increase test coverage

* increase coverage more

* more informative make new parameter

* correct exponential relaxation

* explcitily display errors and warnings

* state value

* name, not symbol

* massive improvement on parameter macro

* typo

* don't use (t) in the wranings

* add LiteralParameter

* fix literal value

* fix wrong test in derived

* fix omitted x

* restore warnonly setting

* fix timescale reference

* use canonical MTK t and derivative as requested

* fix typo in docs

* bump version to 1.0

* add addition process and test for all processes

Project.toml

@@ -1,13 +1,13 @@
 name = "ProcessBasedModelling"
 uuid = "ca969041-2cf3-4b10-bc21-86f4417093eb"
 authors = ["Datseris <datseris.george@gmail.com>"]
-version = "0.1.0"
+version = "1.0.0"
 
 [deps]
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [compat]
-ModelingToolkit = "8.73"
+ModelingToolkit = "9.0"
 Reexport = "1.2"
 julia = "1.9.0"

docs/make.jl

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

docs/src/index.md

@@ -5,12 +5,17 @@ ProcessBasedModelling
 !!! note "Basic familiarity with ModelingToolkit.jl"
     These docs assume that you have some basic familiarity with ModelingToolkit.jl. If you don't going through the introductory tutorial of [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) should be enough to get you started!
 
+!!! note "Default `t` is unitless"
+    Like ModelingToolkit.jl, ProcessBasedModelling.jl also exports `t` as the independent variable representing time.
+    However, instead of the default `t` of ModelingToolkit.jl, here `t` is unitless.
+    Do `t = ModelingToolkit.t` to obtain the unitful version of `t`.
 
 ## Usage
 
 In ProcessBasedModelling.jl, each variable is governed by a "process".
 Conceptually this is just an equation that _defines_ the given variable.
-To couple the variable with the process it is governed by, a user either defines simple equations of the form "variable = expression", or creates an instance of [`Process`](@ref) if the left-hand-side of the equation needs to be anything more complex. In either case, the variable and the expression are both _symbolic expressions_ created via ModellingToolkit.jl (more specifically, via Symbolics.jl).
+To couple the variable with the process it is governed by, a user either defines simple equations of the form "variable = expression", or creates an instance of [`Process`](@ref) if the left-hand-side of the equation needs to be anything more complex (or, simply if you want to utilize the conveniences of predefined processes).
+In either case, the variable and the expression are both _symbolic expressions_ created via ModellingToolkit.jl (more specifically, via Symbolics.jl).
 
 Once all the processes about the physical system are collected, they are given as a `Vector` to the [`processes_to_mtkmodel`](@ref) central function, similarly to how one gives a `Vector` of `Equation`s to e.g., `ModelingToolkit.ODESystem`. This function also defines what quantifies as a "process" in more specificity.
 
@@ -30,13 +35,17 @@ symbolically using ModelingToolkit.jl (**MTK**). We define
 using ModelingToolkit
 using OrdinaryDiffEq: Tsit5
 
-@variables t # independent variable
+@variables t # independent variable _without_ units
 @variables z(t) = 0.0
 @variables x(t) # no default value
 @variables y(t) = 0.0
 ```
 ProcessBasedModelling.jl (**PBM**) strongly recommends that all defined variables have a default value at definition point. Here we didn't do this for ``x`` to illustrate what how such an "omission" will be treated by **PBM**.
 
+!!! note "ModelingToolkit.jl is re-exported"
+    ProcessBasedModelling.jl re-exports the whole `ModelingToolkit` package,
+    so you don't need to be `using` both of them, just `using ProcessBasedModelling`.
+
 To make the equations we want, we can use MTK directly, and call
 ```@example MAIN
 eqs = [
@@ -52,20 +61,23 @@ equations(model)
 
 All good. Now, if we missed the process for one variable (because of our own error/sloppyness/very-large-codebase), MTK will throw an error when we try to _structurally simplify_ the model (a step necessary before solving the ODE problem):
 
-```@example MAIN
+```julia
 model = ODESystem(eqs[1:2], t; name = :example)
-try
-  model = structural_simplify(model)
-catch e
-  return e.msg
-end
+model = structural_simplify(model)
+```
+```
+ERROR: ExtraVariablesSystemException: The system is unbalanced.
+There are 3 highest order derivative variables and 2 equations.
+More variables than equations, here are the potential extra variable(s):
+ z(t)
+ x(t)
+ y(t)
 ```
 
-As you can see, the error message is unhelpful even with such a trivial system of equations,
-as all variables are reported as "potentially missing".
+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 errors 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 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.
 
 **PBM** resolves these problems and always gives accurate error messages when it comes to
 the construction of the system of equations.
@@ -95,15 +107,17 @@ Notice that the resulting **MTK** model is not `structural_simplify`-ed, to allo
 Now, in contrast to before, if we "forgot" a process, **PBM** will react accordingly.
 For example, if we forgot the 2nd process, then the construction will error informatively,
 telling us exactly which variable is missing, and because of which processes it is missing:
-```@example MAIN
-try
-  model = processes_to_mtkmodel(processes[[1, 3]])
-catch e
-  return e.msg
-end
+```julia
+model = processes_to_mtkmodel(processes[[1, 3]])
+```
+```
+ERROR: ArgumentError: Variable x was introduced in process of variable z(t).
+However, a process for x was not provided,
+there is no default process for x, and x doesn't have a default value.
+Please provide a process for variable x.
 ```
 
-If instead we "forgot" the ``y`` process, **PBM** will not error, but instead warn, and make ``y`` equal to a named parameter:
+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:
 ```@example MAIN
 model = processes_to_mtkmodel(processes[1:2])
 equations(model)
@@ -113,8 +127,18 @@ equations(model)
 parameters(model)
 ```
 
+and the warning thrown was:
+```julia
+┌ Warning: Variable y was introduced in process of variable x(t).
+│ However, a process for y 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(y)`.
+└ @ ProcessBasedModelling ...\ProcessBasedModelling\src\make.jl:65
+```
+
 Lastly, [`processes_to_mtkmodel`](@ref) also allows the concept of "default" processes, that can be used for introduced "process-less" variables.
-Default processes are like `processes` and given as a 2nd argument to [`process_to_mtkmodel`](@ref).
+Default processes are like `processes` and given as a 2nd argument to [`processes_to_mtkmodel`](@ref).
 For example,
 
 ```@example MAIN
@@ -143,7 +167,8 @@ equations(model)
 parameters(model)
 ```
 
-This special handling is also why each process explicitly declares a timescale via the [`timescale`](@ref) function that one can optionally extend.
+This special handling is also why each process can declare a timescale via the [`ProcessBasedModelling.timescale`](@ref) function that one can optionally extend
+(although in our experience the default behaviour covers almost all cases).
 
 
 ## Main API function
@@ -158,6 +183,7 @@ processes_to_mtkmodel
 ParameterProcess
 TimeDerivative
 ExpRelaxation
+AdditionProcess
 ```
 
 ## `Process` API
@@ -180,4 +206,5 @@ default_value
 has_variable
 new_derived_named_parameter
 @convert_to_parameters
+LiteralParameter
 ```

src/API.jl

@@ -1,5 +1,8 @@
 """
-A process subtype `p::Process` extends the following unexported functions:
+    Process
+
+A new process must subtype `Process` and can be used in [`processes_to_mtkmodel`](@ref).
+The type must extend the following functions from the module `ProcessBasedModelling`:
 
 - `lhs_variable(p)` which returns the variable the process describes
   (left-hand-side variable). There is a default implementation
@@ -57,12 +60,12 @@ function lhs(p::Process)
     τ = timescale(p)
     v = lhs_variable(p)
     if isnothing(τ) # time variability exists but timescale is nonexistent (unity)
-        return Differential(t)(v)
+        return D(v) # `D` is the MTK canonical variable for time derivative
     elseif τ isa NoTimeDerivative || iszero(τ) # no time variability
         return v
     else # τ is either Num or Real
         τvar = new_derived_named_parameter(v, τ, "τ", false)
-        return τvar*Differential(t)(v)
+        return τvar*D(v)
     end
 end
 

src/ProcessBasedModelling.jl

@@ -8,9 +8,9 @@ module ProcessBasedModelling
 end ProcessBasedModelling
 
 using Reexport
+using ModelingToolkit: t_nounits as t, D_nounits as D
 @reexport using ModelingToolkit
-
-@variables t # independent variable (time)
+export t
 
 include("API.jl")
 include("utils.jl")
@@ -22,10 +22,11 @@ include("processes_basic.jl")
 
 # TODO: Perhaps not don't export `t`?
 export t
-export Process, ParameterProcess, TimeDerivative, ExpRelaxation
+export Process, ParameterProcess, TimeDerivative, ExpRelaxation, AdditionProcess
 export processes_to_mtkmodel
 export new_derived_named_parameter
 export has_variable, default_value
-export @convert_to_parameters
+export @convert_to_parameters, LiteralParameter
+export lhs_variable, rhs, lhs
 
 end

src/make.jl

@@ -1,7 +1,8 @@
 """
     processes_to_mtkmodel(processes::Vector, default::Vector = []; kw...)
 
-Construct a ModelingToolkit.jl model using the provided `processes` and `default` processes.
+Construct a ModelingToolkit.jl model/system using the provided `processes` and `default` processes.
+The model/system is _not_ structurally simplified.
 
 `processes` is a vector whose elements can be:
 
@@ -60,34 +61,31 @@ function processes_to_mtkmodel(_processes, _default = [];
             append_incomplete_variables!(incomplete, introduced, lhs_vars, def_proc)
         else
             def_val = default_value(added_var) # utilize default value (if possible)
+            varstr = ModelingToolkit.getname(added_var)
             if !isnothing(def_val)
                 @warn("""
-                Variable $(added_var) was introduced in process of variable $(introduced[added_var]).
-                However, a process for $(added_var) was not provided,
+                Variable $(varstr) was introduced in process of variable $(introduced[added_var]).
+                However, a process for $(varstr) 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($(added_var))`.
+                `ParameterProcess($(varstr))`.
                 """)
                 parproc = ParameterProcess(added_var)
                 push!(eqs, lhs(parproc) ~ rhs(parproc))
                 push!(lhs_vars, added_var)
             else
                 throw(ArgumentError("""
-                Variable $(added_var) was introduced in process of variable $(introduced[added_var]).
-                However, a process for $(added_var) was not provided,
-                there is no default process for, and it doesn't have a default value.
-                Please provide a process for variable $(added_var).
+                Variable $(varstr) was introduced in process of variable $(introduced[added_var]).
+                However, a process for $(varstr) was not provided,
+                there is no default process for $(varstr), and $(varstr) doesn't have a default value.
+                Please provide a process for variable $(varstr).
                 """))
             end
         end
     end
     sys = type(eqs, independent; name)
     return sys
 end
-# version without given processes
-function processes_to_mtkmodel(; kwargs...)
-    return processes_to_mtkmodel(collect(values(default_processes())); kwargs...)
-end
 
 function expand_multi_processes(procs::Vector)
     # Expand vectors of processes or ODESystems

src/processes_basic.jl

@@ -79,6 +79,46 @@ ExpRelaxation(proc::Union{Process,Equation}, τ) = ExpRelaxation(lhs_variable(pr
 
 timescale(e::ExpRelaxation) = e.timescale
 function rhs(e::ExpRelaxation)
-    dt = isnothing(e.timescale) || iszero(e.timescale)
-    dt ? e.expression : e.expression - e.variable
+    τ = timescale(e)
+    hasdt = if τ isa NoTimeDerivative
+        false
+    elseif isnothing(τ)
+        true
+    else
+        !iszero(τ)
+    end
+    hasdt ? e.expression - e.variable : e.expression
+end
+
+"""
+    AdditionProcess(process, added)
+
+A convenience process for adding `added` to the `rhs` of the given `process`.
+`added` can be a `Process` or `Equation`, in which case it is checked that
+the `lhs_variable` matches. Otherwise, it can be an arbitrary expression.
+"""
+struct AdditionProcess <: Process
+    process
+    added
+    function AdditionProcess(process, added)
+        if typeof(added) <: Union{Process, Equation}
+            if ModelingToolkit.getname(lhs_variable(process)) ≠ ModelingToolkit.getname(lhs_variable(added))
+                @show lhs_variable(process), lhs_variable(added)
+                throw(ArgumentError("Added component does not have the same lhs variable."))
+            end
+        end
+        return new(process, added)
+    end
+end
+
+lhs_variable(a::AdditionProcess) = lhs_variable(a.process)
+timescale(a::AdditionProcess) = timescale(a.process)
+
+function rhs(a::AdditionProcess)
+    if typeof(a.added) <: Union{Process, Equation}
+        plus = rhs(a.added)
+    else
+        plus = a.added
+    end
+    return rhs(a.process) + plus
 end