diff --git a/CHANGELOG.md b/CHANGELOG.md index fb540bd..a3074e7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,10 @@ ProcessBasedModelling.jl follows semver 2.0. Changelog is kept with respect to v1 release. +## 1.9 + +Updated to ModelingToolkit.jl v11, which also comes with a dependency to AGPL license. + ## 1.8 - Updated to ModelingToolkit.jl v10. `type` keyword in `processes_to_mtkmodel` is now no longer used. diff --git a/Project.toml b/Project.toml index e29a5f9..78d7bf0 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,15 @@ name = "ProcessBasedModelling" uuid = "ca969041-2cf3-4b10-bc21-86f4417093eb" authors = ["Datseris "] -version = "1.8.0" +version = "1.9.0" [deps] ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" [compat] -ModelingToolkit = "10.0" +ModelingToolkit = "11" +Symbolics = "7" Reexport = "1.2" julia = "1.9.0" diff --git a/README.md b/README.md index a734813..3722b6a 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ [![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. +ProcessBasedModelling.jl is an extension to [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) (and the wider MTK ecosystem) 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 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 equation form. diff --git a/docs/Project.toml b/docs/Project.toml index 077d0f9..7429447 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -4,4 +4,3 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" -ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" diff --git a/docs/make.jl b/docs/make.jl index b3b8ccb..7f80f02 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -15,4 +15,5 @@ pages = [ build_docs_with_style(pages, ProcessBasedModelling; authors = "George Datseris ", + warnonly = [:doctest, :linkcheck], ) diff --git a/docs/src/index.md b/docs/src/index.md index c517685..05a0f1d 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -31,10 +31,10 @@ Let's say we want to build the system of equations y = z - x ``` -symbolically using ModelingToolkit.jl (**MTK**). We define +symbolically using ModelingToolkit.jl (**MTKBase**). We define ```@example MAIN -using ModelingToolkit +using ProcessBasedModelling # re-exports ModelingToolkit @variables t # independent variable _without_ units @variables z(t) = 0.0 @@ -69,7 +69,7 @@ model = ODESystem(eqs[1:2], t; name = :example) ```julia # here is the error -model = structural_simplify(model) +model = mtkcompile(model) ``` ``` ERROR: ExtraVariablesSystemException: The system is unbalanced. diff --git a/src/API.jl b/src/API.jl index d086fcc..9d4851d 100644 --- a/src/API.jl +++ b/src/API.jl @@ -86,7 +86,7 @@ end rhs(e::Equation) = e.rhs lhs(e::Equation) = e.lhs lhs_variable(e::Equation) = lhs_variable(lhs(e)) -lhs_variable(x::Num) = Num(lhs_variable(x.val)) +lhs_variable(x::Num) = Num(lhs_variable(Symbolics.unwrap(x))) function lhs_variable(x) # basically x is SymbolicUtils.BasicSymbolic{Real} # check whether `x` is a single variable already if is_variable(x) @@ -95,15 +95,15 @@ function lhs_variable(x) # basically x is SymbolicUtils.BasicSymbolic{Real} # check Differential(t)(x) if hasproperty(x, :f) if x.f isa Differential - return x.arguments[1] + return x.args[1] end end # check Differential(t)(x)*parameter - if hasproperty(x, :arguments) - args = x.arguments + if hasproperty(x, :args) + args = x.args di = findfirst(a -> hasproperty(a, :f) && a.f isa Differential, args) if !isnothing(di) - return args[di].arguments[1] + return args[di].args[1] end end # error if all failed diff --git a/src/ProcessBasedModelling.jl b/src/ProcessBasedModelling.jl index 7ee3f7d..f727786 100644 --- a/src/ProcessBasedModelling.jl +++ b/src/ProcessBasedModelling.jl @@ -11,6 +11,8 @@ using Reexport using ModelingToolkit: t_nounits as t, D_nounits as D @reexport using ModelingToolkit +getname = ModelingToolkit.SymbolicIndexingInterface.getname + include("API.jl") include("utils.jl") include("default.jl") diff --git a/src/make.jl b/src/make.jl index e508c64..9d71311 100644 --- a/src/make.jl +++ b/src/make.jl @@ -18,7 +18,7 @@ During construction, the following automations improve user experience: `processes` is a `Vector` whose elements can be: -1. Any instance of a subtype of [`Process`](@ref). `Process` is like a +1. Any instance of a subtype of [`Process`](@ref). `Process` is a wrapper around `Equation` that provides some conveniences, e.g., handling of timescales or not having limitations on the left-hand-side (LHS) form. 1. An `Equation`. The LHS format of the equation is limited. @@ -28,7 +28,7 @@ During construction, the following automations improve user experience: 2. A `Vector` of the above two, which is then expanded. This allows the convenience of functions representing a physical process that may require many equations to be defined (because e.g., they may introduce more variables). -3. A ModelingToolkit.jl `XDESystem`, in which case the `equations` of the system are expanded +3. A ModelingToolkit.jl `System`, in which case the `equations` of the system are expanded as if they were given as a vector of equations like above. This allows the convenience of straightforwardly coupling with already existing `XDESystem`s. @@ -59,7 +59,7 @@ These registered default processes are used when `default` is a `Module`. (has happened to me many times!). """ function processes_to_mtkmodel(args...; - type = System, name = :model, independent = t, kw..., + name = :model, independent = t, kw..., ) eqs = processes_to_mtkeqs(args...; kw...) sys = System(eqs, independent; name) @@ -122,7 +122,7 @@ function processes_to_mtkeqs(_processes::Vector, default::Dict{Num, Any}; However, a process for $(added_var) 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($(ModelingToolkit.getname(added_var)))`. + `ParameterProcess($(getname(added_var)))`. """) end parproc = ParameterProcess(added_var) diff --git a/src/processes_basic.jl b/src/processes_basic.jl index 18c0851..e23cdc0 100644 --- a/src/processes_basic.jl +++ b/src/processes_basic.jl @@ -111,7 +111,7 @@ struct AdditionProcess <: Process function AdditionProcess(process, added::Vector) for add in added if typeof(add) <: Union{Process, Equation} - v1, v2 = ModelingToolkit.getname(lhs_variable(process)), ModelingToolkit.getname(lhs_variable(add)) + v1, v2 = getname(lhs_variable(process)), getname(lhs_variable(add)) if v1 ≠ v2 throw(ArgumentError( "Added processes do not have the same lhs variable. Got: $(v1), $(v2)" diff --git a/src/utils.jl b/src/utils.jl index 571af50..0992b08 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,3 +1,5 @@ +import Symbolics + """ LiteralParameter(p) @@ -17,8 +19,9 @@ _literalvalue(p::LiteralParameter) = p.p has_symbolic_var(eqs, var) Return `true` if symbolic variable `var` exists in the equation(s) `eq`, `false` otherwise. -This works for either `@parameters` or `@variables`. -If `var` is a `Symbol` isntead of a `Num`, all variables are converted to their names +This works for either `@parameters` or `@variables` and currently only compares versus these +objects (i.e., it ignores `Differentials` or other complex structures). +If `var` is a `Symbol` instead of a `Num`, all variables are converted to their names and equality is checked on the basis of the name only. has_symbolic_var(model, var) @@ -27,7 +30,8 @@ When given a MTK model (such as `ODESystem`) search in _all_ the equations of th including observed variables. """ function has_symbolic_var(eq::Equation, var) - vars = get_variables(eq) + vars = get_variables(eq) # this includes differentials + vars = filter(x -> is_parameter(x) || is_variable(x), vars) return _has_thing(var, vars) end has_symbolic_var(eqs::Vector{Equation}, var) = any(eq -> has_symbolic_var(eq, var), eqs) @@ -37,8 +41,8 @@ function _has_thing(var::Num, vars) return any(isequal(var), vars) end function _has_thing(var::Symbol, vars) - vars = ModelingToolkit.getname.(vars) - var = ModelingToolkit.getname(var) + vars = getname.(vars) + var = getname(var) return any(isequal(var), vars) end @@ -55,34 +59,30 @@ all_equations(model) = vcat(equations(model), observed(model)) Return the default value of a symbolic variable `x` or `nothing` if it doesn't have any. Return `x` if `x` is not a symbolic variable. The difference with `ModelingToolkit.getdefault` is that this function will -not error on the absence of a default value. +not error in the absence of a default value. """ default_value(x) = x -default_value(x::Num) = default_value(x.val) -function default_value(x::ModelingToolkit.SymbolicUtils.Symbolic) - if haskey(x.metadata, ModelingToolkit.Symbolics.VariableDefaultValue) - return x.metadata[ModelingToolkit.Symbolics.VariableDefaultValue] - else - @warn("No default value assigned to variable/parameter $(x).") - return nothing - end +default_value(x::Num) = default_value(Symbolics.unwrap(x)) +function default_value(x::Symbolics.SymbolicT) + val = Symbolics.getmetadata(x, Symbolics.VariableDefaultValue, nothing) + isnothing(val) && @warn("No default value assigned to variable/parameter $(x).") + return val end + # trick to get default values for state variables: # Base.Fix1(getindex, ModelingToolkit.defaults(ssys)).(states(ssys)) # while `defaults` returns all default assignments. -is_variable(x::Num) = is_variable(x.val) -function is_variable(x) - if x isa ModelingToolkit.SymbolicUtils.Symbolic - if isnothing(x.metadata) - return false - end - if haskey(x.metadata, ModelingToolkit.Symbolics.VariableSource) - src = x.metadata[ModelingToolkit.Symbolics.VariableSource] - return first(src) == :variables - end - end - return false +is_variable(x::Num) = is_variable(Symbolics.unwrap(x)) +function is_variable(x::Symbolics.SymbolicT) + value = getmetadata(x, Symbolics.VariableSource, nothing) + return value isa Tuple{Symbol, Symbol} && value[1] == :variables +end + +is_parameter(x::Num) = is_parameter(Symbolics.unwrap(x)) +function is_parameter(x::Symbolics.SymbolicT) + value = getmetadata(x, Symbolics.VariableSource, nothing) + return value isa Tuple{Symbol, Symbol} && value[1] == :parameters end """ @@ -111,7 +111,7 @@ Now `p` will be a parameter with name `:τ_x` and default value `0.5`. new_derived_named_parameter(v, value::Num, extra::String; kw...) = value new_derived_named_parameter(v, value::LiteralParameter, extra::String; kw...) = value.p function new_derived_named_parameter(v, value::Real, extra; connector = "_", prefix = true) - n = string(ModelingToolkit.getname(v)) + n = string(Symbolics.getname(v)) newstring = if prefix extra*connector*n else @@ -134,7 +134,7 @@ 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, +This macro is useful to convert e.g., keyword arguments into named parameters, while also allowing the user to give custom parameter names, or to leave some keywords as numeric literals. @@ -174,7 +174,7 @@ macro convert_to_parameters(vars...) $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)))) + (ModelingToolkit.toparam_validate)((Symbolics.wrap)((SymbolicUtils.setmetadata)((Symbolics.setdefaultval)((SymbolicUtils.Sym){SymbolicUtils.SymReal}($varname; type = Real), $binding), Symbolics.VariableSource, (:parameters, $varname)))) )) ) ) diff --git a/test/runtests.jl b/test/runtests.jl index e5cd1ae..97f0bfc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,7 +19,7 @@ end # If that's the case, we are sure model construction was valid # First, make some default processes - @variables T(t) = 300.0 # temperature, in Kelvin + @variables T(t) = 300.0 # temperature, in Kelvin @variables α(t) # albedo of ice, unitless @variables ε(t) # effective emissivity, unitless solar_constant = 340.25 # W/m^2, already divided by 4 @@ -66,7 +66,7 @@ end @test length(unknowns(sys)) == 1 @test has_symbolic_var(equations(sys), T) - u0s = [[300.0], [100.0]] + u0s = [[T => 300.0], [T => 100.0]] ufs = [] for u0 in u0s p = ODEProblem(sys, u0, (0.0, 1000.0*365*24*60*60.0)) @@ -128,7 +128,7 @@ end @test length(unknowns(sys)) == 3 sys = processes_to_mtkmodel(procs[1:3]) @test length(unknowns(sys)) == 3 - @test length(unknowns(structural_simplify(sys))) == 2 + @test length(unknowns(mtkcompile(sys))) == 2 end end @@ -139,10 +139,10 @@ end @testset "derived" begin @variables x(t) = 0.5 p = new_derived_named_parameter(x, 0.2, "t") - @test ModelingToolkit.getname(p) == :t_x + @test ModelingToolkit.SymbolicIndexingInterface.getname(p) == :t_x @test default_value(p) == 0.2 p = new_derived_named_parameter(x, 0.2, "t"; prefix = false, connector = "") - @test ModelingToolkit.getname(p) == :xt + @test ModelingToolkit.SymbolicIndexingInterface.getname(p) == :xt end @testset "convert" begin @@ -151,7 +151,7 @@ end @convert_to_parameters A B C @test A isa Num @test default_value(A) == 0.5 - @test ModelingToolkit.getname(C) == :X + @test ModelingToolkit.SymbolicIndexingInterface.getname(C) == :X end @testset "literal in derived" begin @@ -219,10 +219,10 @@ end sys = processes_to_mtkmodel(procs) sys2 = processes_to_mtkmodel([sys, w ~ x*y]) @test length(equations(sys2)) == 4 - @test sort(ModelingToolkit.getname.(unknowns(sys2))) == [:w, :x, :y, :z] + @test sort(ModelingToolkit.SymbolicIndexingInterface.getname.(unknowns(sys2))) == [:w, :x, :y, :z] end -@testset "equation in RHS" begin +@testset "duplcate equation in RHS" begin @variables z(t) = 0.0 @variables x(t) = 0.0 @variables y(t) = 0.0 @@ -231,7 +231,7 @@ end y ~ z-x, # is an equation, not a process! z ~ (z ~ x^2), ] - @test_throws ["an `<: Equation` type"] processes_to_mtkeqs(procs) + @test_throws ["more than one"] processes_to_mtkeqs(procs) end @testset "not actual process" begin