-
Notifications
You must be signed in to change notification settings - Fork 0
/
make.jl
193 lines (178 loc) · 8.51 KB
/
make.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
"""
processes_to_mtkmodel(processes::Vector [, default]; kw...)
Construct a ModelingToolkit.jl model/system using the provided `processes` and `default` processes.
The model/system is _not_ structurally simplified. During construction, the following automations
improve user experience:
- Variable(s) introduced in `processes` that does not itself have a process obtain
a default process from `default`.
- If no default exists, but the variable(s) itself has a default numerical value,
a [`ParameterProcess`](@ref) is created for said variable and a warning is thrown.
- Else, an informative error is thrown.
- An error is also thrown if any variable has two or more processes assigned to it.
`processes` is a `Vector` whose elements can be:
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)`,
however, the versions with `p` may fail unexpectedly. Anything else will error.
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
as if they were given as a vector of equations like above. This allows the convenience
of straightforwardly coupling with already existing `XDESystem`s.
## Default processes
`processes_to_mtkmodel` allows for specifying default processes by giving `default`.
These default processes are assigned to variables introduced in the main input `processes`,
but themselves do not have an assigned process in the main input.
`default` can be a `Vector` of individual processes (`Equation` or `Process`).
Alternatively, `default` can be a `Module`. The recommended way to build field-specific
modelling libraries based on ProcessBasedModelling.jl is to define modules/submodules
that offer a pool of pre-defined variables and processes.
Modules may register their own default processes via the function
[`register_default_process!`](@ref).
These registered processes are used when `default` is a `Module`.
## Keyword arguments
- `type = ODESystem`: the model type to make.
- `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_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.
"""
processes_to_mtkmodel(procs::Vector; kw...) =
processes_to_mtkmodel(procs, Dict{Num, Any}(); kw...)
processes_to_mtkmodel(procs::Vector, m::Module; kw...) =
processes_to_mtkmodel(procs, default_processes(m); kw...)
processes_to_mtkmodel(procs::Vector, v::Vector; kw...) =
processes_to_mtkmodel(procs, default_dict(v); kw...)
# The main implementation has the defaults to be a map from variable to process
# because this simplifies a bit the code
function processes_to_mtkmodel(_processes::Vector, default::Dict{Num, Any};
type = ODESystem, name = nameof(type), independent = t, warn_default::Bool = true,
)
processes = expand_multi_processes(_processes)
# Setup: obtain lhs-variables so we can track new variables that are not
# in this vector. The vector has to be of type `Num`
lhs_vars = Num[lhs_variable(p) for p in processes]
ensure_unique_vars(lhs_vars)
# First pass: generate equations from given processes
# and collect variables without equations
incomplete = Num[]
introduced = Dict{Num, Num}()
eqs = Equation[]
for proc in processes
append_incomplete_variables!(incomplete, introduced, lhs_vars, proc)
# add the generated equation in the pool of equations
push!(eqs, lhs(proc) ~ rhs(proc))
end
# Second pass: attempt to add default processes to incomplete variables
# throw an error if default does not exist
while !isempty(incomplete) # using a `while` allows us check variables added by the defaults
added_var = popfirst!(incomplete)
if haskey(default, added_var)
# Then obtain default process
def_proc = default[added_var]
# add the equation to the equation pool
push!(eqs, lhs(def_proc) ~ rhs(def_proc))
# Also ensure that potentially new processes introduced by the default process
# are taken care of and have a default value
push!(lhs_vars, lhs_variable(def_proc))
append_incomplete_variables!(incomplete, introduced, lhs_vars, def_proc)
else
def_val = default_value(added_var) # utilize default value (if possible)
if !isnothing(def_val)
if warn_default == true
@warn("""
Variable $(added_var) was introduced in process of variable $(introduced[added_var]).
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)))`.
""")
end
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).
"""))
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
expanded = Any[procs...]
idxs = findall(p -> p isa etypes, procs)
multiprocs = expanded[idxs]
deleteat!(expanded, idxs)
for mp in multiprocs
if mp isa Vector
append!(expanded, mp)
else # then it is XDE system
append!(expanded, equations(mp))
end
end
return expanded
end
function default_dict(processes)
default = Dict{Num, Any}()
for proc in processes
key = lhs_variable(proc)
if haskey(default, key)
throw(ArgumentError("More than 1 processes for variable $(key) in default processes."))
end
default[key] = proc
end
return default
end
# Add variables to `incomplete` from expression `r`, provided they are not in `lhs_vars`
# already. Also record which variables added them
function append_incomplete_variables!(incomplete, introduced, lhs_vars, process)
proc_vars = filter(is_variable, get_variables(rhs(process)))
newvars = setdiff(proc_vars, lhs_vars)
# Store newly introduced variables without a lhs
for nv in newvars
# Note we can't use `(nv ∈ incomplete)`, that creates a symbolic expression
any(isequal(nv), incomplete) && continue # skip if is already recorded
numnv = (nv isa Num) ? nv : Num(nv)
Symbol(nv) == :t && continue # Time-dependence is not a new variable!
push!(incomplete, numnv)
# Also record which variable introduced it
introduced[numnv] = lhs_variable(process)
end
return
end
function ensure_unique_vars(lhs_vars)
nonun = nonunique(lhs_vars)
isempty(nonun) || error("The following variables have more than one processes assigned to them: $(nonun)")
return
end
function nonunique(x::AbstractArray{T}) where T
uniqueset = Set{T}()
duplicatedset = Set{T}()
duplicatedvector = Vector{T}()
for i in x
if i ∈ uniqueset
if !(i ∈ duplicatedset)
push!(duplicatedset, i)
push!(duplicatedvector, i)
end
else
push!(uniqueset, i)
end
end
duplicatedvector
end