HI. This question is a continuation of a previous topic here, on constructing anonymous functions dynamically. In my case, these functions are intended to form the dynamical rule used for initialising CoupledODEs from a set of stock, flow and variable specifications entered by the user.
Here is the sample code, as far as I have been able to strip it down for this question. Please note that the SD (system dynamics) model rabbits (created in the method demo() at the bottom) could contain arbitrarily many stocks, variables and flows. Also variables could contain arithmetic expressions, as well as (in this case) the numerical value 0.3:
module SD
using DynamicalSystems
#-----------------------------------------------------------------------------------------
# Module types:
#-----------------------------------------------------------------------------------------
# VExpr represents anything that can be eval()ed to a Number: Exprs, identifiers and Numbers.
VExpr = Union{Expr,Symbol,Number}
mutable struct SDModel
stocks::Dict{Symbol,Float64}
flows::Dict{Symbol,Vector{VExpr}}
variables::Dict{Symbol,VExpr}
ode_model::Union{Nothing,CoupledODEs}
function SDModel( title="An Exponential SD model")
new( Dict{Symbol,Float64}(),
Dict{Symbol,Vector{VExpr}}(),
Dict{Symbol,VExpr}(),
nothing
)
end
end
#-----------------------------------------------------------------------------------------
# Module methods:
#-----------------------------------------------------------------------------------------
"""
Build the CoupledODEs required for the SDModel to run.
"""
function build!( sdm)
p_keys = collect(keys(sdm.variables))
p_vals = deepcopy(collect(values(sdm.variables)))
u_keys = collect(keys(sdm.stocks))
u_vals = collect(values(sdm.stocks))
du_vals = deepcopy(collect(values(sdm.flows)))
p_builder = Vector{Function}(undef,length(p_keys))
du_builder = Vector{Function}(undef,length(u_keys))
# Build dynamical rules for parameters:
for (j,expr) in enumerate(p_vals)
for (i,k) in enumerate(p_keys)
substitute!(expr, k, :(p[$i]))
end
for (i,k) in enumerate(u_keys)
substitute!(expr, k, :(u[$i]))
end
println("!!! Dump of function to update p:")
dump(Expr(:(->),:(p,u),expr))
p_builder[j] = eval(Expr(:(->),:(p,u),expr))
end
# Build dynamical rules for derivatives:
for (j,flow_set) in enumerate(du_vals)
for (i,k) in enumerate(p_keys)
for flow_term in flow_set
substitute!(flow_term,k,:(p[$i]))
end
end
for (i,k) in enumerate(u_keys)
for flow_term in flow_set
substitute!(flow_term,k,:(u[$i]))
end
end
println("\n!!! Dump of function to update du:")
dump(Expr(:(->),:(p,u),Expr(:call,:(+),flow_set...)))
du_builder[j] = eval(Expr(:(->),:(p,u),Expr(:call,:(+),flow_set...)))
end
# Build general dynamical rule for CoupledODEs :
dynamical_rule! = function (du,u,p,t)
for i in eachindex(p)
p[i] = Float64(p_builder[i](p,u))
end
for i in 1:length(u_keys)
du[i] = Float64(eval(du_builder[i](p,u)))
end
nothing
end
# Alternative dynamical_rule! that is capable of running:
#= dynamical_rule! = function (du,u,p,t)
p[1] = 0.3
du[1] = p[1] * u[1]
nothing
end=#
sdm.ode_model = CoupledODEs(dynamical_rule!, u_vals, p_vals)
end
#-----------------------------------------------------------------------------------------
"""
Substitute IN PLACE each occurrence of old_argument by the new_arg in the given expression.
"""
function substitute!( expr::Expr, old_arg, new_arg)
for (i,arg) in enumerate(expr.args)
if arg==old_arg
expr.args[i] = new_arg
elseif arg isa Expr
substitute!(arg, old_arg, new_arg)
end
end
nothing
end
function substitute!( sym::Symbol, old_arg, new_arg)
if old_arg == sym
sym = new_arg
end
nothing
end
substitute!( expr::Number, old_arg, new_arg) = nothing
#-----------------------------------------------------------------------------------------
"""
Use-case for building a simple exponential model.
"""
function demo()
rabbits = SDModel()
rabbits.stocks[:Rabbits] = 2.0
rabbits.variables[:mu] = Meta.quot(0.3)
rabbits.flows[:Rabbits] = [Meta.parse("mu * Rabbits")]
println()
println( "!!! Dump of pre-build rabbits model:")
display(rabbits)
println()
build!( rabbits)
println()
trajectory(rabbits.ode_model, 3; Δt=0.25)
end
end
If this code is run with the alternative dynamical_rule! commented out in build!(), it runs correctly. If, however, I run it with the currently used dynamical_rule! constructed within the build!(), this rule method appears to do its job correctly in setting the values of p and du, yet when I pass it to the CoupledODEs constructor, this throws errors connected to the fact that I use the union type VExpr = Union{Expr, Number, Symbol} to manage expression syntax.
My question is this: I feel I could solve these issues if only I could get a handle on what precisely is the nature of the internal problem for CoupledODEs(). However, I’m finding it very difficult to obtain that information. Any ideas?
Yup, fair enough. I’d already scrapped a ton of code, but you’re right - what’s left still isn’t minimal. Just note for the moment that dynamical_rule! is the assembled result of a load of substitution and manipulation of Exprs, and I’ll make sure I get out a minimal example by the end of the day. Thanks very much! 