Skip to content

Commit

Permalink
Finish docs! (#2)
Browse files Browse the repository at this point in the history
* typos

* correct incorrect differential in docs

* linebreaks

* fix docs problems

* add convert_to_parameters macro

* fix missing exports

* correct names in index

* finish readme
  • Loading branch information
Datseris authored Feb 21, 2024
1 parent e5a0bc7 commit 7203e4a
Show file tree
Hide file tree
Showing 8 changed files with 86 additions and 51 deletions.
7 changes: 4 additions & 3 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,15 +1,16 @@
# ProcessBasedModelling.jl

[![docsdev](https://img.shields.io/badge/docs-dev-lightblue.svg)](https://juliadynamics.github.io/ProcessBasedModelling,jl/dev/)
[![docsdev](https://img.shields.io/badge/docs-dev-lightblue.svg)](https://juliadynamics.github.io/ProcessBasedModelling.jl/dev/)
[![docsstable](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliadynamics.github.io/ProcessBasedModelling.jl/stable/)
[![CI](https://github.com/JuliaDynamics/ProcessBasedModelling.jl/workflows/CI/badge.svg)](https://github.com/JuliaDynamics/ProcessBasedModelling.jl/actions?query=workflow%3ACI)
[![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.
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 particular physical concept or observable and there are few (or any) duplicate variables to make the definition of MTK "factories" worthwhile, while, there plenty of different physical representations, or _processes_ to represent the given physical concept.
In many fields this approach parallels modelling reasoning line of the researcher more closely than the "components" approach.
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 many scientific fields this approach parallels the modelling reasoning of the researcher more closely than the "components" approach.

Beyond this reasoning style, the biggest strength of ProcessBasedModelling.jl is the **informative errors and automation** it provides regarding incorrect/incomplete equations. When building the MTK model via ProcessBasedModelling.jl the user provides a vector of "processes": equations or custom types that have a well defined and single left-hand-side variable.
This allows ProcessBasedModelling.jl to:
Expand Down
2 changes: 0 additions & 2 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,5 @@ pages = [

build_docs_with_style(pages, ProcessBasedModelling;
authors = "George Datseris <[email protected]>",
# We need to remove the cross references because we don't list here
# the whole `DynamicalSystem` API...
warnonly = true,
)
43 changes: 25 additions & 18 deletions docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ To make the equations we want, we can use MTK directly, and call
```@example MAIN
eqs = [
Differential(t)(z) ~ x^2 - z
Differential(x) ~ 0.1y
Differential(t)(x) ~ 0.1y
y ~ z - x
]
Expand All @@ -50,30 +50,32 @@ model = ODESystem(eqs, t; name = :example)
equations(model)
```

All good. Now, if we missed the process for one variable (because of our own error/sloppyness/very-large-codebase), MTK will not throw an error at model construction,
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
model = ODESystem(eqs[1:2], t; name = :example)
model = structural_simplify(model)
equations(model)
```

only at the construction of the "problem" (here the `ODEProblem`)

```@example MAIN
try
prob = ODEProblem(model)
model = structural_simplify(model)
catch e
return e.msg
end
```

Interestingly, the error is wrong. ``x`` is defined and has an equation, at least on the basis of our scientific reasoning. However ``y`` that ``x`` introduced 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.
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".
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.

**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).

**PBM** resolves these problems and always gives accurate error messages. 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).
Here is what the user defines to make the same system of equations:
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
Expand All @@ -82,14 +84,17 @@ processes = [
```

which is then given to

```@example MAIN
model = processes_to_mtkmodel(processes; name = :example)
equations(model)
```

Notice that the resulting **MTK** model is not `structural_simplify`-ed, to allow composing it with other models. By default `t` is taken as the independent variable.

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:
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]])
Expand All @@ -109,7 +114,7 @@ parameters(model)
```

Lastly, [`processes_to_mtkmodel`](@ref) also allows the concept of "default" processes, that can be used for introduced "process-less" variables.
Default processes like `processes` given as a 2nd argument to [`process_to_mtkmodel`](@ref).
Default processes are like `processes` and given as a 2nd argument to [`process_to_mtkmodel`](@ref).
For example,

```@example MAIN
Expand All @@ -119,14 +124,14 @@ equations(model)

does not throw any warnings as it obtained a process for ``y`` from the given default processes.

### Special handling of timescales
## Special handling of timescales

In dynamical systems modelling the timescale associated with a process is a special parameter. That is why, if a timescale is given for either the [`TimeDerivative`](@ref) or [`ExpRelaxation`](@ref) processes, it is converted to a named `@parameter` by default:

```@example MAIN
processes = [
ExpRelaxation(z, x^2, 2.0), # third argument is the timescale
TimeDerivative(x, 0.1*y, 0.5),
ExpRelaxation(z, x^2, 2.0), # third argument is the timescale
TimeDerivative(x, 0.1*y, 0.5), # third argument is the timescale
y ~ z-x,
]
Expand Down Expand Up @@ -163,6 +168,7 @@ This API describes how you can implement your own `Process` subtype, if the [exi
Process
ProcessBasedModelling.lhs_variable
ProcessBasedModelling.rhs
ProcessBasedModelling.timescale
ProcessBasedModelling.NoTimeDerivative
ProcessBasedModelling.lhs
```
Expand All @@ -173,4 +179,5 @@ ProcessBasedModelling.lhs
default_value
has_variable
new_derived_named_parameter
@convert_to_parameters
```
6 changes: 3 additions & 3 deletions src/API.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ A process subtype `p::Process` extends the following unexported functions:
(left-hand-side variable). There is a default implementation
`lhs_variable(p) = p.variable` if the field exists.
- `rhs(p)` which is the right-hand-side expression, i.e., the "actual" process.
- (optional) `timescale`, which defaults to [`NoTimeDerivative`](@ref).
- (optional) `timescale(p)`, which defaults to [`NoTimeDerivative`](@ref).
- (optional) `lhs(p)` which returns the left-hand-side. Let `τ = timescale(p)`.
Then default `lhs(p)` behaviour depends on `τ` as follows:
- Just `lhs_variable(p)` if `τ == NoTimeDerivative()`.
Expand All @@ -20,7 +20,7 @@ abstract type Process end
"""
ProcessBasedModelling.NoTimeDerivative()
Singleton value that is the default output of the [`timescale`](@ref) function
Singleton value that is the default output of the `timescale` function
for variables that do not vary in time autonomously, i.e., they have no d/dt derivative
and hence the concept of a "timescale" does not apply to them.
"""
Expand Down Expand Up @@ -50,7 +50,7 @@ timescale(::Process) = NoTimeDerivative()
ProcessBasedModelling.lhs(p::Process)
Return the left-hand-side of the equation that `p` represents as an `Expression`.
If [`timescale`](@ref) is implemented for `p`, typically `lhs` does not need to be as well.
If `timescale` is implemented for `p`, typically `lhs` does not need to be as well.
See [`Process`](@ref) for more.
"""
function lhs(p::Process)
Expand Down
3 changes: 2 additions & 1 deletion src/ProcessBasedModelling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ export t
export Process, ParameterProcess, TimeDerivative, ExpRelaxation
export processes_to_mtkmodel
export new_derived_named_parameter
export hs_variable, default_value
export has_variable, default_value
export @convert_to_parameters

end
12 changes: 8 additions & 4 deletions src/processes_basic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
The simplest process which equates a given `variable` to a constant value
that is encapsulated in a parameter. If `value isa Real`, then
hence, a named parameter with the name of `variable` and `_0` appended is created.
a named parameter with the name of `variable` and `_0` appended is created.
Else, if `valua isa Num` then it is taken as the paremeter directly.
Example:
Expand Down Expand Up @@ -58,15 +58,16 @@ given `expression`, with timescale `τ`. It creates the equation:
τn*Differential(t)(variable) ~ expression - variable
```
Where `τn` is a new named `@parameter` with the value of `τ`
and name `τ_(\$(variable))`. If instead `τ` is `nothing`, then 1 is used in its place.
and name `τ_(\$(variable))`. If instead `τ` is `nothing`, then 1 is used in its place
(this is the default behavior).
If `iszero(τ)`, then the equation `variable ~ expression` is created instead.
The convenience function
```julia
ExpRelaxation(process, τ)
```
allows converting an existing process (or equation) into an exponential relaxation
by using the `rhs` as the `expression` in the equation above.
by using the `rhs(process)` as the `expression` in the equation above.
"""
struct ExpRelaxation <: Process
variable
Expand All @@ -77,4 +78,7 @@ ExpRelaxation(v, e) = ExpRelaxation(v, e, nothing)
ExpRelaxation(proc::Union{Process,Equation}, τ) = ExpRelaxation(lhs_variable(proc), rhs(proc), τ)

timescale(e::ExpRelaxation) = e.timescale
rhs(e::ExpRelaxation) = iszero(e.timescale) ? e.expression : e.expression - e.variable
function rhs(e::ExpRelaxation)
dt = isnothing(e.timescale) || iszero(e.timescale)
dt ? e.expression : e.expression - e.variable
end
50 changes: 30 additions & 20 deletions src/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -68,28 +68,38 @@ function new_derived_named_parameter(newstring::String, value::Real)
return first(dummy)
end

# Macro thanks to Jonas Isensee,
# https://discourse.julialang.org/t/metaprogramming-macro-calling-another-macro-making-named-variables/109621/6
"""
@named_parameters vars...
@convert_to_parameters vars...
Convert all Julia variables `vars` into `@parameters` with name the same as `vars`
and default value the same as the value of `vars`.
Convert all variables `vars` into `@parameters` with name the same as `vars`
and default value the same as the value of `vars`. Example:
```
julia> A, B = 0.5, 0.5
(0.5, 0.5)
julia> @convert_to_parameters A B
2-element Vector{Num}:
A
B
julia> typeof(A) # `A` is not a number anymore!
Num
julia> default_value(A)
0.5
"""
macro named_parameters(vars...)
return quote
out = []
for v in vars
res = (ModelingToolkit.toparam)((Symbolics.wrap)((Symbolics.SymbolicUtils.setmetadata)((Symbolics.setdefaultval)((Sym){Real}($(QuoteNode(v))), $(esc(v))), Symbolics.VariableSource, (:parameters, $(QuoteNode(v))))))
push!(out, res)
end
$(out...)
macro convert_to_parameters(vars...)
expr = Expr(:block)
for var in vars
binding = esc(var)
varname = QuoteNode(var)
push!(expr.args,
:($binding = (ModelingToolkit.toparam)((Symbolics.wrap)((SymbolicUtils.setmetadata)((Symbolics.setdefaultval)((Sym){Real}($varname), $binding), Symbolics.VariableSource, (:parameters, $varname)))))
)
end
push!(expr.args, Expr(:vect, esc.(vars)...))
return expr
end

# This may help:

# julia> @macroexpand @parameters A=0.1 B=0.1
# quote
# A = (ModelingToolkit.toparam)((Symbolics.wrap)((SymbolicUtils.setmetadata)((Symbolics.setdefaultval)((Sym){Real}(:A), 0.1), Symbolics.VariableSource, (:parameters, :A))))
# B = (ModelingToolkit.toparam)((Symbolics.wrap)((SymbolicUtils.setmetadata)((Symbolics.setdefaultval)((Sym){Real}(:B), 0.1), Symbolics.VariableSource, (:parameters, :B))))
# [A, B]
# end
14 changes: 14 additions & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -112,3 +112,17 @@ end
@testset "extending default processes" begin
# API not yet finished on this one
end

@testset "utility functions" begin
@variables x(t) = 0.5
p = new_derived_named_parameter(x, 0.2, "t")
@test ModelingToolkit.getname(p) == :x_t
@test default_value(p) == 0.2
p = new_derived_named_parameter(x, p, "lala")
@test ModelingToolkit.getname(p) == :x_t

A, B = 0.5, 0.5
@convert_to_parameters A B
@test A isa Num
@test default_value(A) == 0.5
end

0 comments on commit 7203e4a

Please sign in to comment.