diff --git a/Project.toml b/Project.toml index 08af165..3e63242 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,13 @@ 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" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" [compat] -ModelingToolkit = "8.73" +ModelingToolkit = "9.0" Reexport = "1.2" julia = "1.9.0" 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, ) diff --git a/docs/src/index.md b/docs/src/index.md index 9f20146..3e5c6e9 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 = [ @@ -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 ``` diff --git a/src/API.jl b/src/API.jl index 18fc153..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,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 diff --git a/src/ProcessBasedModelling.jl b/src/ProcessBasedModelling.jl index b1610d6..0c233d5 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") @@ -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 diff --git a/src/make.jl b/src/make.jl index 6e16905..c568494 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: @@ -60,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, 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 @@ -84,10 +86,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/src/processes_basic.jl b/src/processes_basic.jl index 3c9379f..fb4b195 100644 --- a/src/processes_basic.jl +++ b/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 \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl index 95983f4..c0078ea 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,3 +1,18 @@ +""" + 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 +# necessary for the macro +_literalvalue(x) = x +_literalvalue(p::LiteralParameter) = p.p + """ has_variable(eq, var) @@ -47,12 +62,23 @@ 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 +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 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 @@ -74,22 +100,34 @@ 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 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. + +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) @@ -97,8 +135,16 @@ 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( + $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. + # 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 d383ec0..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), @@ -52,10 +51,11 @@ 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 + @test has_variable(equations(sys), T) u0s = [[300.0], [100.0]] ufs = [] @@ -67,6 +67,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 @@ -86,7 +102,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,34 +111,73 @@ 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 -@testset "extending default processes" begin - # API not yet finished on this one +@testset "utility functions" begin + # Test an untested clause: + @test default_value(0.5) == 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 + 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 + @variables x(t) = 0.5 + 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 + end -@testset "utility functions" begin +@testset "default processes" 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 + @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