diff --git a/CHANGELOG.md b/CHANGELOG.md index e965a3d..4a72c46 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,4 +1,11 @@ # Changelog ProcessBasedModelling.jl follows semver 2.0. -Changelog is kept with respect to v1 release. \ No newline at end of file +Changelog is kept with respect to v1 release. + +## 1.1 + +- New keyword `warn_default` in `processes_to_mtkmodel`. +- Now `processes_to_mtkmodel` allows for `XDESystems` as elements of the input vector. +- Detection of LHS-variable in `Equation` has been improved significantly. + Now, LHS can be any of: `x`, `Differential(t)(x)`, `Differential(t)(x)*p`, `p*Differential(t)(x)`. \ No newline at end of file diff --git a/Project.toml b/Project.toml index c96ac7e..bcafa07 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ProcessBasedModelling" uuid = "ca969041-2cf3-4b10-bc21-86f4417093eb" authors = ["Datseris "] -version = "1.0.5" +version = "1.1.0" [deps] ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" diff --git a/docs/src/index.md b/docs/src/index.md index 876884d..b3fd084 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -76,22 +76,25 @@ More variables than equations, here are the potential extra variable(s): 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 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. +Moreover, in our experience these error messages become increasingly less accurate or helpful when a model has many equations and/or variables. This makes it difficult to quickly find out where the "mistake" happened in the equations. **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). +For the majority of cases, **PBM** can infer the LHS variable a process "defines" automatically, just by passing in a vector of `Equation`s, like in **MTK**. +For cases where this is not possible a dedicated `Process` type is provided, whose subtypes act as wrappers around equations providing some additional conveniences. + 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 - y ~ z - x, # can be an equation because LHS is single variable + ExpRelaxation(z, x^2), # defines z, introduces x; `Process` subtype (optional) + Differential(t)(x) ~ 0.1*y, # defines x, introduces y; normal `Equation` + y ~ z - x, # defines y; normal `Equation` ] ``` @@ -117,17 +120,20 @@ there is no default process for x(t), and (t)x 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 warn, and make ``y`` equal to a named parameter, since ``y`` has a default value: +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. +So, running: ```@example MAIN model = processes_to_mtkmodel(processes[1:2]) equations(model) ``` +Makes the named parameter: + ```@example MAIN parameters(model) ``` -and the warning thrown was: +and throws the warning: ```julia ┌ Warning: Variable y(t) was introduced in process of variable x(t). │ However, a process for y(t) was not provided, diff --git a/src/API.jl b/src/API.jl index 52520dc..3a462af 100644 --- a/src/API.jl +++ b/src/API.jl @@ -51,8 +51,9 @@ timescale(::Process) = NoTimeDerivative() """ ProcessBasedModelling.lhs(p::Process) + ProcessBasedModelling.lhs(eq::Equation) -Return the left-hand-side of the equation that `p` represents as an `Expression`. +Return the right-hand-side of the equation governing the process. If `timescale` is implemented for `p`, typically `lhs` does not need to be as well. See [`Process`](@ref) for more. """ @@ -70,10 +71,9 @@ function lhs(p::Process) end """ - ProcessBasedModelling.rhs(p::Process) + ProcessBasedModelling.rhs(process) -Return the right-hand-side of the equation that `p` represents as an `Expression`. -See [`Process`](@ref) for more. +Return the right-hand-side of the equation governing the process. """ function rhs(p::Process) error("Right-hand side (`rhs`) is not defined for process $(nameof(typeof(p))).") @@ -81,13 +81,29 @@ end # Extensions for `Equation`: rhs(e::Equation) = e.rhs -lhs(e::Equation) = lhs_variable(e) -function lhs_variable(e::Equation) - x = e.lhs - # we first check whether `x` is a variable - if !is_variable(x) - throw(ArgumentError("In given equation $(e), the left-hand-side does "* - "not represent a single variable.")) +lhs(e::Equation) = e.lhs +lhs_variable(e::Equation) = lhs_variable(lhs(e)) +lhs_variable(x::Num) = Num(lhs_variable(x.val)) +function lhs_variable(x) # basically x is SymbolicUtils.BasicSymbolic{Real} + # check whether `x` is a single variable already + if is_variable(x) + return x + end + # check Differential(t)(x) + if hasproperty(x, :f) + if x.f isa Differential + return x.arguments[1] + end + end + # check Differential(t)(x)*parameter + if hasproperty(x, :arguments) + args = x.arguments + di = findfirst(a -> hasproperty(a, :f) && a.f isa Differential, args) + if !isnothing(di) + return args[di].arguments[1] + end end - return x + # error if all failed + throw(ArgumentError("We analyzed the LHS `$(x)` "* + "but could not extract a single variable it represents.")) end diff --git a/src/make.jl b/src/make.jl index ad5da09..db85b4b 100644 --- a/src/make.jl +++ b/src/make.jl @@ -6,15 +6,22 @@ The model/system is _not_ structurally simplified. `processes` is a vector whose elements can be: -- Any instance of a subtype of [`Process`](@ref). -- An `Equation` which is of the form `variable ~ expression` with `variable` a single - variable resulting from an `@variables` call. -- 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. +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. + Let `x` be a `@variable` and `p` be a `@parameter`. Then, the LHS can only be one of: + `x`, `Differential(t)(x)`, `Differential(t)(x)*p`, `p*Differential(t)(x)`. + Anything else will either error or fail unexpectedly. +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. +3. A ModelingToolkit.jl `XDESystem`, 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 already existing systems. `default` is a vector that can contain the first two possibilities only -as it contains default processes that may be assigned to variables introduced in `processes` -but they don't themselves have an assigned process. +as it contains default processes that may be assigned to individual variables introduced in +`processes` but they don't themselves have an assigned process. It is expected that downstream packages that use ProcessBasedModelling.jl to make a field-specific library implement a 1-argument version of `processes_to_mtkmodel`, @@ -26,10 +33,10 @@ or provide a wrapper function for it, and add a default value for `default`. - `name = nameof(type)`: the name of the model - `independent = t`: the independent variable (default: `@variables t`). `t` is also exported by ProcessBasedModelling.jl for convenience. -- `warn_defaut::Bool = true`: if `true`, throw a warning when a variable does not +- `warn_default::Bool = true`: if `true`, throw a warning when a variable does not have an assigned process but it has a default value so that it becomes a parameter instead. """ -function processes_to_mtkmodel(_processes, _default = []; +function processes_to_mtkmodel(_processes::Vector, _default = []; type = ODESystem, name = nameof(type), independent = t, warn_default::Bool = true, ) processes = expand_multi_processes(_processes) @@ -86,19 +93,25 @@ function processes_to_mtkmodel(_processes, _default = []; end end end + # return eqs sys = type(eqs, independent; name) return sys end function expand_multi_processes(procs::Vector) + etypes = Union{Vector, ODESystem, SDESystem, PDESystem} + !any(p -> p isa etypes, procs) && return procs # Expand vectors of processes or ODESystems - !any(p -> p isa Vector, procs) && return procs - expanded = deepcopy(procs) - idxs = findall(p -> p isa Vector, procs) + expanded = Any[procs...] + idxs = findall(p -> p isa etypes, procs) multiprocs = expanded[idxs] deleteat!(expanded, idxs) for mp in multiprocs - append!(expanded, mp) + if mp isa Vector + append!(expanded, mp) + else # then it is XDE system + append!(expanded, equations(mp)) + end end return expanded end @@ -144,8 +157,8 @@ function nonunique(x::AbstractArray{T}) where T duplicatedset = Set{T}() duplicatedvector = Vector{T}() for i in x - if(i in uniqueset) - if !(i in duplicatedset) + if i ∈ uniqueset + if !(i ∈ duplicatedset) push!(duplicatedset, i) push!(duplicatedvector, i) end diff --git a/test/runtests.jl b/test/runtests.jl index abf519e..406e523 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -191,4 +191,21 @@ end @testset "addition process error" begin @variables x(t) y(t) z(t) @test_throws ArgumentError AdditionProcess(x ~ 0.1z, y ~ x^2) +end + +@testset "ODESystem as process" begin + @variables z(t) = 0.0 + @variables x(t) = 0.0 + @variables y(t) = 0.0 + @variables w(t) = 0.0 + procs = [ + ExpRelaxation(z, x^2, 1.0), # introduces x and y variables + TimeDerivative(x, 0.1*y), # introduces y variable! + y ~ z-x, # is an equation, not a process! + ] + + 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] end \ No newline at end of file