Catalyst.jl
Catalyst.jl copied to clipboard
DSL does not work for parameters/species/variables with complicated interdependencies
Sometimes parameters/species/variables depend on each other. In the DSL, we create a single block each, in order:
@parameters begin
end
@variables begin
end
@species begin
end
Here, something declared, which depends on a subsequent thing, will throw an error. Right now, this is pretty much only a problem for when default values depend on other quantities. I don't think there is any relevant metadata which includes other parameters/variables/species, but if that gets added, that could be a problem.
I think the only "natural" case is a species/variable with a parametric initial condition, which is fine. But someone wanting to implement a model where something like this happen is possible, so would be nice to fix some time.
Here is the full list of cases which generate errors.
# Parameters depending on variables.
rn = @reaction_network begin
@parameters a = A
@variables A
end
# Parameters depending on species.
rn = @reaction_network begin
@parameters a = X
@species X
end
# Variables depending on species.
rn = @reaction_network begin
@variables V = X
@species X
end
# Parameters depending on subsequent parameters.
rn = @reaction_network begin
@parameters a = b b
end
# Variables depending on subsequent variables.
rn = @reaction_network begin
@variables A = B B
end
# Species depending on subsequent species.
rn = @reaction_network begin
@species X = Y Y
end
Was just about to make an issue related to this, it's a particular problem when using remake
to update the initial conditions of an ODEProblem(rs, ...; remove_conserved = true)
.
- This is because the Γ conservation law constants are defined as parameters, which have symbolic dependencies on the initial conditions.
- Perhaps an overload is necessary for the internal
SciMLBase.remake
function here:
function _updated_u0_p_symmap(prob, u0, ::Val{true}, p, ::Val{true})
isu0dep = any(symbolic_type(v) !== NotSymbolic() for (_, v) in u0)
ispdep = any(symbolic_type(v) !== NotSymbolic() for (_, v) in p)
if !isu0dep && !ispdep
return remake_buffer(prob, state_values(prob), u0),
remake_buffer(prob, parameter_values(prob), p)
end
if !isu0dep
u0 = remake_buffer(prob, state_values(prob), u0)
return _updated_u0_p_symmap(prob, u0, Val(false), p, Val(true))
end
if !ispdep
p = remake_buffer(prob, parameter_values(prob), p)
return _updated_u0_p_symmap(prob, u0, Val(true), p, Val(false))
end
varmap = merge(u0, p)
u0 = anydict(k => symbolic_type(v) === NotSymbolic() ? v : symbolic_evaluate(v, varmap)
for (k, v) in u0)
p = anydict(k => symbolic_type(v) === NotSymbolic() ? v : symbolic_evaluate(v, varmap)
for (k, v) in p)
return remake_buffer(prob, state_values(prob), u0),
remake_buffer(prob, parameter_values(prob), p)
end
- Here the main issue I've found are that the
ispdep
check doesn't find that Γ in the parameters is a symbolic type, because it is evaluated to numerical values in the ODEProblem. The only place where the symbolic expressions for Γ are exposed are indefaults
field of theODESystem
. - This obstructs the ostensible allowance for interdependencies between u0 and p, seen with the
varmap = merge(u0, p)
line. - Currently my workaround is calling
symbolic_evaluate
prior to callingremake
using the symbolic expressions from the defaults map:
function compute_Γ(tunable_var_value_map, constraint_syms, sys_defs::T) where T
#* Merge the default system definitions with the new symbolic map
merged_defs = merge(sys_defs, Dict(tunable_var_value_map))
NTuple{4, eltype(tunable_var_value_map)}(constraint => symbolic_evaluate(merged_defs[constraint], merged_defs) for constraint in constraint_syms)
end
I can make a separate issue for this in this repo (not sure if it pertains most here, or ModelingToolkit.jl
or SciMLBase.jl
) if this is truly a problem and I'm not missing anything.
Sounds good, yes, feel free to make a separate issue as well.
Right now there are a couple of bugs related to remake
circulating within MTK (but I expect these to be sorted out soon). Once these are fixed, I will start having a look at your issue.