From b2e63f34be0c4f3338e5aa446eda400b83b8a28e Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 10:27:42 +0000 Subject: [PATCH 01/23] bump version of MTK --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 08af165..646c9e0 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,6 @@ 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" From 51c6b15d90d190d91cf14b8bfb45651613f4bb97 Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 10:36:08 +0000 Subject: [PATCH 02/23] rename states to unknowns --- test/runtests.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index d383ec0..9bc8992 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -52,10 +52,10 @@ using OrdinaryDiffEq sys = processes_to_mtkmodel(processes) @test sys isa ODESystem - @test length(states(sys)) == 3 + @test length(unknowns(sys)) == 3 sys = structural_simplify(sys) - @test length(states(sys)) == 1 + @test length(unknowns(sys)) == 1 u0s = [[300.0], [100.0]] ufs = [] @@ -86,7 +86,7 @@ end @testset "first two, still missing y, but it has default" begin model = @test_logs (:warn, r"\W*((?i)Variable(?-i))\W*") processes_to_mtkmodel(procs[1:2]) - @test length(states(model)) == 3 + @test length(unknowns(model)) == 3 end @testset "first with default the third; missing x" begin @@ -95,17 +95,17 @@ end @testset "first with default the second; y gets contant value a different warning" begin model = @test_logs (:warn, r"\W*((?i)parameter(?-i))\W*") processes_to_mtkmodel(procs[1:1], procs[2:2]) - @test length(states(model)) == 3 + @test length(unknowns(model)) == 3 end @testset "all three processes given" begin sys = processes_to_mtkmodel(procs[1:1], procs[2:3]) - @test length(states(sys)) == 3 + @test length(unknowns(sys)) == 3 sys = processes_to_mtkmodel(procs[1:2], procs[3:3]) - @test length(states(sys)) == 3 + @test length(unknowns(sys)) == 3 sys = processes_to_mtkmodel(procs[1:3]) - @test length(states(sys)) == 3 - @test length(states(structural_simplify(sys))) == 2 + @test length(unknowns(sys)) == 3 + @test length(unknowns(structural_simplify(sys))) == 2 end end From e7cc8535feae6982aa09d211254d589432f5151e Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 11:06:14 +0000 Subject: [PATCH 03/23] make default choices more fitting for dynsys --- docs/src/index.md | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 9f20146..3f76a27 100644 --- a/docs/src/index.md +++ b/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 = [ From 24794303b790f1b6561b8e0690260fbcea91fba0 Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 11:39:53 +0000 Subject: [PATCH 04/23] increase test coverage --- src/make.jl | 3 ++- test/runtests.jl | 3 +++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/make.jl b/src/make.jl index 6e16905..96138b2 100644 --- a/src/make.jl +++ b/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: diff --git a/test/runtests.jl b/test/runtests.jl index 9bc8992..be35d92 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -56,6 +56,7 @@ using OrdinaryDiffEq sys = structural_simplify(sys) @test length(unknowns(sys)) == 1 + @test has_variable(equations(sys), T) u0s = [[300.0], [100.0]] ufs = [] @@ -125,4 +126,6 @@ end @convert_to_parameters A B @test A isa Num @test default_value(A) == 0.5 + + @test default_value(0.5) == 0.5 end From ae7cde8502d9306c755520dd5143a88a6736d116 Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 11:42:02 +0000 Subject: [PATCH 05/23] increase coverage more --- src/make.jl | 4 ---- test/runtests.jl | 16 ++++++++++++++++ 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/src/make.jl b/src/make.jl index 96138b2..59bf2f9 100644 --- a/src/make.jl +++ b/src/make.jl @@ -85,10 +85,6 @@ function processes_to_mtkmodel(_processes, _default = []; 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 diff --git a/test/runtests.jl b/test/runtests.jl index be35d92..61f0e29 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -68,6 +68,22 @@ using OrdinaryDiffEq @test ufs[1] ≈ [319] atol = 1 @test ufs[2] ≈ [245] atol = 1 + + # vector of processes + processes = [ + [TanhProcess(α, T, 0.7, 0.289, 10.0, 274.5), + TanhProcess(ε, T, 0.5, 0.41, 2.0, 288.0),], + HeatBalance() + ] + + sys = processes_to_mtkmodel(processes) + @test sys isa ODESystem + @test length(unknowns(sys)) == 3 + + sys = structural_simplify(sys) + @test length(unknowns(sys)) == 1 + @test has_variable(equations(sys), T) + end @testset "add missing processes" begin From 896f5633c3ede6ceabbb3eeaa261b15b925c09e1 Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 13:00:46 +0000 Subject: [PATCH 06/23] more informative make new parameter --- src/utils.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/utils.jl b/src/utils.jl index 95983f4..c5f0781 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -51,6 +51,14 @@ If `value isa Num`, return `value`. Otherwise, create a new MTK `@parameter` whose name is created from `variable` by adding the `extra` string. If `suffix = true` the extra is added at the end after a `_`. Otherwise it is added at the start, then a `_` and then the variable name. + +For example, + +``` +@variables x(t) +p = new_derived_named_parameter(x, 0.5, "τ") +``` +Now `p` will be a parameter with symbol `:τ_x` and default value `0.5`. """ new_derived_named_parameter(v, value::Num, extra, suffix = true) = value function new_derived_named_parameter(v, value::Real, extra, suffix = true) From e3460d69634694ec1496339021c9ef5df52048a5 Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 13:09:16 +0000 Subject: [PATCH 07/23] correct exponential relaxation --- src/processes_basic.jl | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/processes_basic.jl b/src/processes_basic.jl index 3c9379f..933a193 100644 --- a/src/processes_basic.jl +++ b/src/processes_basic.jl @@ -79,6 +79,13 @@ 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 \ No newline at end of file From 262aa497421f270bb3fe7e04183e77137ef1530a Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 13:09:31 +0000 Subject: [PATCH 08/23] explcitily display errors and warnings --- docs/src/index.md | 45 ++++++++++++++++++++++++++++++--------------- src/make.jl | 2 +- 2 files changed, 31 insertions(+), 16 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 3f76a27..1cca259 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -61,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. @@ -104,12 +107,14 @@ 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(t) was introduced in process of variable z(t). +However, a process for x(t) was not provided, +there is no default process for x(t), and x(t) 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 instead warn, and make ``y`` equal to a named parameter: @@ -122,6 +127,16 @@ equations(model) parameters(model) ``` +and the warning thrown was: +```julia +┌ Warning: Variable y(t) was introduced in process of variable x(t). +│ However, a process for y(t) 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(t))`. +└ @ 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). For example, diff --git a/src/make.jl b/src/make.jl index 59bf2f9..2bc0923 100644 --- a/src/make.jl +++ b/src/make.jl @@ -76,7 +76,7 @@ function processes_to_mtkmodel(_processes, _default = []; 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. + there is no default process for $(added_var), and $(added_var) doesn't have a default value. Please provide a process for variable $(added_var). """)) end From e80632b95a0bf24e4769fa5d139b7d90aeaabc4f Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 13:10:29 +0000 Subject: [PATCH 09/23] state value --- docs/src/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/index.md b/docs/src/index.md index 1cca259..edf8e88 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -117,7 +117,7 @@ there is no default process for x(t), and x(t) 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 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) From a0b704e82a7ec752a15281d881a7195eb54d455a Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 13:12:02 +0000 Subject: [PATCH 10/23] name, not symbol --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index c5f0781..e064a69 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -58,7 +58,7 @@ For example, @variables x(t) p = new_derived_named_parameter(x, 0.5, "τ") ``` -Now `p` will be a parameter with symbol `:τ_x` and default value `0.5`. +Now `p` will be a parameter with name `:τ_x` and default value `0.5`. """ new_derived_named_parameter(v, value::Num, extra, suffix = true) = value function new_derived_named_parameter(v, value::Real, extra, suffix = true) From e91e7cd7648c948ad0f7ed286a8b39ea8c56bbd0 Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 14:15:05 +0000 Subject: [PATCH 11/23] massive improvement on parameter macro --- src/utils.jl | 26 ++++++++++++++++++++++---- test/runtests.jl | 6 ++++-- 2 files changed, 26 insertions(+), 6 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index e064a69..397d109 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -82,22 +82,33 @@ end @convert_to_parameters vars... Convert all variables `vars` into `@parameters` with name the same as `vars` -and default value the same as the value of `vars`. Example: +and default value the same as the value of `vars`. The macro leaves unaltered +inputs that are alread `Num`, assumming they are already parameters. +This macro is extremely useful to convert e.g., keyword arguments into named parameters, +while also allowing the user to give custom parameter names. + +Example: ``` julia> A, B = 0.5, 0.5 (0.5, 0.5) -julia> @convert_to_parameters A B -2-element Vector{Num}: +julia> C = first(@parameters X = 0.5) + +julia> @convert_to_parameters A B C +3-element Vector{Num}: A B + X julia> typeof(A) # `A` is not a number anymore! Num julia> default_value(A) 0.5 + +julia> C # the binding `C` still corresponds to parameter named `:X`! + X """ macro convert_to_parameters(vars...) expr = Expr(:block) @@ -105,8 +116,15 @@ macro convert_to_parameters(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))))) + :($binding = ifelse( + # don't do anyting if this is already a Num + $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)))) + ) ) + ) end push!(expr.args, Expr(:vect, esc.(vars)...)) return expr diff --git a/test/runtests.jl b/test/runtests.jl index 61f0e29..70b996f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -139,9 +139,11 @@ end @test ModelingToolkit.getname(p) == :x_t A, B = 0.5, 0.5 - @convert_to_parameters A B + C = first(@parameters X = 0.5) + @convert_to_parameters A B C @test A isa Num @test default_value(A) == 0.5 - + @test ModelingToolkit.getname(C) == :X + # Test an untested clause: @test default_value(0.5) == 0.5 end From 7d1fe63bd0ce8c87d7201cdf98f12a362f497abf Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 14:18:31 +0000 Subject: [PATCH 12/23] typo --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index 397d109..76676eb 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -83,7 +83,7 @@ end 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 alread `Num`, assumming they are already parameters. +inputs that are of type `Num`, assumming they are already parameters. This macro is extremely useful to convert e.g., keyword arguments into named parameters, while also allowing the user to give custom parameter names. From 6a88a88a5f2065267ca5f8c5367a67dd0a6d711e Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 15:37:57 +0000 Subject: [PATCH 13/23] don't use (t) in the wranings --- docs/src/index.md | 14 +++++++------- src/ProcessBasedModelling.jl | 1 + src/make.jl | 15 ++++++++------- 3 files changed, 16 insertions(+), 14 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index edf8e88..82abb6a 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -111,10 +111,10 @@ telling us exactly which variable is missing, and because of which processes it model = processes_to_mtkmodel(processes[[1, 3]]) ``` ``` -ERROR: ArgumentError: Variable x(t) was introduced in process of variable z(t). -However, a process for x(t) was not provided, -there is no default process for x(t), and x(t) doesn't have a default value. -Please provide a process for variable x(t). +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 warn, and make ``y`` equal to a named parameter, since ``y`` has a default value: @@ -129,11 +129,11 @@ parameters(model) and the warning thrown was: ```julia -┌ Warning: Variable y(t) was introduced in process of variable x(t). -│ However, a process for y(t) was not provided, +┌ 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(t))`. +│ `ParameterProcess(y)`. └ @ ProcessBasedModelling ...\ProcessBasedModelling\src\make.jl:65 ``` diff --git a/src/ProcessBasedModelling.jl b/src/ProcessBasedModelling.jl index b1610d6..009b153 100644 --- a/src/ProcessBasedModelling.jl +++ b/src/ProcessBasedModelling.jl @@ -27,5 +27,6 @@ export processes_to_mtkmodel export new_derived_named_parameter export has_variable, default_value export @convert_to_parameters +export lhs_variable, rhs, lhs end diff --git a/src/make.jl b/src/make.jl index 2bc0923..c568494 100644 --- a/src/make.jl +++ b/src/make.jl @@ -61,23 +61,24 @@ 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 $(added_var), and $(added_var) 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 From 31e61d997bf88bb66360cf389bb784466b9ba47d Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 16:52:26 +0000 Subject: [PATCH 14/23] add LiteralParameter --- docs/src/index.md | 1 + src/ProcessBasedModelling.jl | 2 +- src/utils.jl | 23 ++++++++++++++++++++--- test/runtests.jl | 8 ++++++++ 4 files changed, 30 insertions(+), 4 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 82abb6a..44fe849 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -204,4 +204,5 @@ default_value has_variable new_derived_named_parameter @convert_to_parameters +LiteralParameter ``` diff --git a/src/ProcessBasedModelling.jl b/src/ProcessBasedModelling.jl index 009b153..2d54661 100644 --- a/src/ProcessBasedModelling.jl +++ b/src/ProcessBasedModelling.jl @@ -26,7 +26,7 @@ export Process, ParameterProcess, TimeDerivative, ExpRelaxation 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 diff --git a/src/utils.jl b/src/utils.jl index 76676eb..7fb1357 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,3 +1,15 @@ +""" + LiteralParameter(p) + +A wrapper around a value `p` to indicate to +[`new_derived_named_parameter`](@ref) or [`@convert_to_parameters`](@ref) +to _not_ convert the given parameter `p` into a named `@parameters` instance, +but rather keep it as a numeric literal in the generated equations. +""" +struct LiteralParameter{P} + p::P +end + """ has_variable(eq, var) @@ -47,7 +59,9 @@ end """ new_derived_named_parameter(variable, value, extra::String, suffix = true) -If `value isa Num`, return `value`. Otherwise, create a new MTK `@parameter` +If `value isa Num` return `value`. +If `value isa `[`LiteralParameter`](@ref), replace it with its literal value. +Otherwise, create a new MTK `@parameter` whose name is created from `variable` by adding the `extra` string. If `suffix = true` the extra is added at the end after a `_`. Otherwise it is added at the start, then a `_` and then the variable name. @@ -60,7 +74,8 @@ p = new_derived_named_parameter(x, 0.5, "τ") ``` Now `p` will be a parameter with name `:τ_x` and default value `0.5`. """ -new_derived_named_parameter(v, value::Num, extra, suffix = true) = value +new_derived_named_parameter(v, value::Num, args...) = value +new_derived_named_parameter(v, value::LiteralParameter, args...) = value.p function new_derived_named_parameter(v, value::Real, extra, suffix = true) n = string(ModelingToolkit.getname(v)) newstring = if suffix @@ -84,6 +99,7 @@ end 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, while also allowing the user to give custom parameter names. @@ -117,12 +133,13 @@ macro convert_to_parameters(vars...) varname = QuoteNode(var) push!(expr.args, :($binding = ifelse( + $binding isa LiteralParameter, $(binding).p, ifelse( # don't do anyting if this is already a Num $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)))) - ) + )) ) ) end diff --git a/test/runtests.jl b/test/runtests.jl index 70b996f..52305af 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -146,4 +146,12 @@ end @test ModelingToolkit.getname(C) == :X # Test an untested clause: @test default_value(0.5) == 0.5 + + p = LiteralParameter(0.5) + p = new_derived_named_parameter(x, p, "t") + @test p == 0.5 + + p = LiteralParameter(0.5) + @convert_to_parameters p + @test p == 0.5 end From 02c9fa681a486dc6a4c092f8ee4bfab9ed0886ce Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 17:02:40 +0000 Subject: [PATCH 15/23] fix literal value --- src/utils.jl | 7 +++++-- test/runtests.jl | 49 +++++++++++++++++++++++++++++------------------- 2 files changed, 35 insertions(+), 21 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 7fb1357..c0078ea 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -9,6 +9,9 @@ but rather keep it as a numeric literal in the generated equations. struct LiteralParameter{P} p::P end +# necessary for the macro +_literalvalue(x) = x +_literalvalue(p::LiteralParameter) = p.p """ has_variable(eq, var) @@ -63,7 +66,7 @@ If `value isa Num` return `value`. If `value isa `[`LiteralParameter`](@ref), replace it with its literal value. Otherwise, create a new MTK `@parameter` whose name is created from `variable` by adding the `extra` string. -If `suffix = true` the extra is added at the end after a `_`. Otherwise +If `suffix == true` the extra is added at the end after a `_`. Otherwise it is added at the start, then a `_` and then the variable name. For example, @@ -133,7 +136,7 @@ macro convert_to_parameters(vars...) varname = QuoteNode(var) push!(expr.args, :($binding = ifelse( - $binding isa LiteralParameter, $(binding).p, ifelse( + $binding isa LiteralParameter, _literalvalue($(binding)), ifelse( # don't do anyting if this is already a Num $binding isa Num, $binding, # Else, convert to modeling toolkit param. diff --git a/test/runtests.jl b/test/runtests.jl index 52305af..1c09681 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -131,27 +131,38 @@ end 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 - C = first(@parameters X = 0.5) - @convert_to_parameters A B C - @test A isa Num - @test default_value(A) == 0.5 - @test ModelingToolkit.getname(C) == :X # Test an untested clause: @test default_value(0.5) == 0.5 - p = LiteralParameter(0.5) - p = new_derived_named_parameter(x, p, "t") - @test p == 0.5 + @testset "derived" 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 + end + + @testset "convert" begin + A, B = 0.5, 0.5 + C = first(@parameters X = 0.5) + @convert_to_parameters A B C + @test A isa Num + @test default_value(A) == 0.5 + @test ModelingToolkit.getname(C) == :X + + end + + @testset "literal in derived" begin + p = LiteralParameter(0.5) + p = new_derived_named_parameter(x, p, "t") + @test p == 0.5 + end + + @testset "literal in covert" begin + p = LiteralParameter(0.5) + @convert_to_parameters p + @test p == 0.5 + end - p = LiteralParameter(0.5) - @convert_to_parameters p - @test p == 0.5 end From 4c4e96a4854b9b9f7444e10d1fe08d2474185d51 Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 18:06:36 +0000 Subject: [PATCH 16/23] fix wrong test in derived --- test/runtests.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 1c09681..292e377 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -139,8 +139,6 @@ end 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 end @testset "convert" begin @@ -150,7 +148,6 @@ end @test A isa Num @test default_value(A) == 0.5 @test ModelingToolkit.getname(C) == :X - end @testset "literal in derived" begin From c2036fb7d71dd123372f45bb80698f5b336c2e0d Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 18:22:32 +0000 Subject: [PATCH 17/23] fix omitted x --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 292e377..731ec50 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -151,6 +151,7 @@ end end @testset "literal in derived" begin + @variables x(t) = 0.5 p = LiteralParameter(0.5) p = new_derived_named_parameter(x, p, "t") @test p == 0.5 From 8628f5ac5010aa8b2a12830cb2d975af5d13669a Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 18:22:49 +0000 Subject: [PATCH 18/23] restore warnonly setting --- docs/make.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index ebee781..b3b8ccb 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -15,5 +15,4 @@ pages = [ build_docs_with_style(pages, ProcessBasedModelling; authors = "George Datseris ", - warnonly = true, ) From e4597196ca22d70f3f6f7613d74b22d46cbe9946 Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 18:52:25 +0000 Subject: [PATCH 19/23] fix timescale reference --- docs/src/index.md | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 44fe849..8aadd0e 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -138,7 +138,7 @@ and the warning thrown was: ``` 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 @@ -167,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 @@ -186,7 +187,7 @@ ExpRelaxation ## `Process` API -This API describes how you can implement your own `Process` subtype, if the [existing predefined subtypes](@ref predefined_processes) don't fit your bill! +This API describes how you can implement your own `Process` subtype, if the [existing predefined subtypes](@ref predefineProcessBasedModelling.timescaled_processes) don't fit your bill! ```@docs Process From 1061b33fbfd7cd56ef5cd5ca7f3ac0dd5e7d1bef Mon Sep 17 00:00:00 2001 From: Datseris Date: Thu, 22 Feb 2024 18:57:15 +0000 Subject: [PATCH 20/23] use canonical MTK t and derivative as requested --- src/API.jl | 4 ++-- src/ProcessBasedModelling.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/API.jl b/src/API.jl index 18fc153..f56c2c1 100644 --- a/src/API.jl +++ b/src/API.jl @@ -57,12 +57,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) 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 diff --git a/src/ProcessBasedModelling.jl b/src/ProcessBasedModelling.jl index 2d54661..46c6cfe 100644 --- a/src/ProcessBasedModelling.jl +++ b/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") From de0211c1f2b4b97c3c99481eba85a6813a92bfce Mon Sep 17 00:00:00 2001 From: Datseris Date: Fri, 23 Feb 2024 11:55:46 +0000 Subject: [PATCH 21/23] fix typo in docs --- docs/src/index.md | 2 +- src/API.jl | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 8aadd0e..75ce767 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -187,7 +187,7 @@ ExpRelaxation ## `Process` API -This API describes how you can implement your own `Process` subtype, if the [existing predefined subtypes](@ref predefineProcessBasedModelling.timescaled_processes) don't fit your bill! +This API describes how you can implement your own `Process` subtype, if the [existing predefined subtypes](@ref predefined_processes) don't fit your bill! ```@docs Process diff --git a/src/API.jl b/src/API.jl index f56c2c1..2e93e01 100644 --- a/src/API.jl +++ b/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,7 +60,7 @@ function lhs(p::Process) τ = timescale(p) v = lhs_variable(p) if isnothing(τ) # time variability exists but timescale is nonexistent (unity) - return D(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 From 6854a269800f296a86217394af20584e746ffd5b Mon Sep 17 00:00:00 2001 From: Datseris Date: Fri, 23 Feb 2024 12:40:22 +0000 Subject: [PATCH 22/23] bump version to 1.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 646c9e0..3e63242 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ProcessBasedModelling" uuid = "ca969041-2cf3-4b10-bc21-86f4417093eb" authors = ["Datseris "] -version = "0.1.0" +version = "1.0.0" [deps] ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" From 779256c3a32fa59b5078bd08d75001404817a606 Mon Sep 17 00:00:00 2001 From: Datseris Date: Fri, 23 Feb 2024 13:15:13 +0000 Subject: [PATCH 23/23] add addition process and test for all processes --- docs/src/index.md | 1 + src/ProcessBasedModelling.jl | 2 +- src/processes_basic.jl | 33 +++++++++++++++++++++++++++++++++ test/runtests.jl | 27 ++++++++++++++++++++++----- 4 files changed, 57 insertions(+), 6 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 75ce767..3e5c6e9 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -183,6 +183,7 @@ processes_to_mtkmodel ParameterProcess TimeDerivative ExpRelaxation +AdditionProcess ``` ## `Process` API diff --git a/src/ProcessBasedModelling.jl b/src/ProcessBasedModelling.jl index 46c6cfe..0c233d5 100644 --- a/src/ProcessBasedModelling.jl +++ b/src/ProcessBasedModelling.jl @@ -22,7 +22,7 @@ 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 diff --git a/src/processes_basic.jl b/src/processes_basic.jl index 933a193..fb4b195 100644 --- a/src/processes_basic.jl +++ b/src/processes_basic.jl @@ -88,4 +88,37 @@ function rhs(e::ExpRelaxation) !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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 731ec50..0c66a1a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -43,7 +43,6 @@ using OrdinaryDiffEq return left + (right - left)*(1 + tanh(2(T - reference)/(scale)))*0.5 end - processes = [ TanhProcess(α, T, 0.7, 0.289, 10.0, 274.5), TanhProcess(ε, T, 0.5, 0.41, 2.0, 288.0), @@ -126,10 +125,6 @@ end end end -@testset "extending default processes" begin - # API not yet finished on this one -end - @testset "utility functions" begin # Test an untested clause: @test default_value(0.5) == 0.5 @@ -164,3 +159,25 @@ end end end + +@testset "default processes" begin + @variables x(t) = 0.5 + @variables y(t) = 0.5 + @variables z(t) = 0.5 + @variables w(t) = 0.5 + @variables q(t) = 0.5 + processes = [ + TimeDerivative(x, x^2, 1.2), + ParameterProcess(y), + ExpRelaxation(z, x^2), + AdditionProcess(ParameterProcess(w), x^2), + AdditionProcess(TimeDerivative(q, x^2, 1.2), ExpRelaxation(q, x^2)) + ] + mtk = processes_to_mtkmodel(processes) + eqs = equations(mtk) + @test has_variable(eqs, x) + @test has_variable(eqs, y) + @test has_variable(eqs, z) + @test has_variable(eqs, w) + @test has_variable(eqs, q) +end \ No newline at end of file