@unpack from Catalyst reaction system and unused variables

Specific Domains Chemistry
,
1

Hi all. I have happily been using Catalyst.jl (and of course related packages, such as ModelingToolkit.jl and DifferentialEquations.jl), but when remaking problems while using EnsembleProblem I encountered some odd behavior. In general, when I am making my reaction systems, many of the parameters end up being unused. That means, when I use the @variables a[1:N] to create N variables, these come back when I use the @unpack macro, even though the parameters are not in the reaction system.

To illustrate this, consider the following simple MWE of logistic growth

using Catalyst

N = 4
@variables t
@species (x(t))[1:N]
@parameters a[1:N]

rxs = Array{Any}(undef, N)

for n in 1:N
    if n != 4
        rxs[n] = Reaction(a[n]*(1-x[n]), [x[n]], [x[n]], [1], [2])
    else
        rxs[n] = Reaction(1, [x[n]], [x[n]], [1], [2])
    end
end

@named rs = ReactionSystem(rxs)

In this example I have not used a[4], and indeed when doing parameters(rs) this parameter does not appear in the parameters.

julia> parameters(rs)
3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 a[1]
 a[2]
 a[3]

Yet, when using the @unpack macro, it does

julia> @unpack a = rs; a
a[1:4]

This becomes strange when using remake(...) as I am basically specify the value of (in my case, many) parameters that are not even in the reaction system. While I can potentially filter a by using the reaction system directly, e.g. something along the lines of

julia> p = Symbolics.Num.(ModelingToolkit.parameters(prob.f.sys))
julia> [_a => avalues for _a in a if any(isequal.(_a, p))]

this imposes a heavy computation time burden as any(...) of course needs to loop over the array until it finds some. As I have large systems with many parameters these loops slow down my code by a lot.

So,

  1. is there a way of extracting/unpacking only the parameters that are present in the reaction system?, and
  2. why does @unpack also unpack these variables even though they are not contained in the reaction system?

Many thanks.

2

Do you see the same issue with a directly defined ODESystem? Catalyst doesn’t really do anything different from ModelingToolkit as far as I am aware (in fact, @unpack is from ModelingToolkit where it is defined for general AbstractSystems), so I think this is likely an issue with how the macro is defined in ModelingToolkit.

Also, keep in mind that ModelingToolkit just released a major, new breaking version (V9), and Catalyst does not support it yet. So this behavior may now be completely different on V9. I’d suggest seeing if ModelingToolkit V9 has fixed this issue for a normal ODESystem, and if not, open an issue about it there. If it is fixed then Catalyst V14 should automatically pick up whatever changes ModelingToolkit has made.

3

The @unpack behavior you are seeing may actually be the desired behavior in V9 of MTK now – you are declaring an array parameter so it may no longer treat the components as separate individual parameters internally (I think there have been some changes in that regard).

4

Some other things you could try

  • Filter and save the list of a’s from parameters(rs) to later use in remake at system construction (so you only need to do it once).
  • Use namespaced indexing, i.e. right after building your system complete it, rs = complete(rs), and then use rs.a[2] to access a given variable.
  • Finally, remake itself may have some issues currently. The change to using SymbolicIndexingInterface led to many code changes in SciMLBase, and remake in particular still hasn’t been fully fixed I think. (But I’m not sure what does/doesn’t work currently.)
5

Yeah I believe when making an ODESystem the problem persists, e.g.

julia> odesys = ModelingToolkit.convert(ODESystem, rs)
Model rs with 4 equations
States (4):
  (x(t))[1]
  (x(t))[2]
  (x(t))[3]
  (x(t))[4]
Parameters (3):
  a[1]
  a[2]
  a[3]

julia> @unpack a = odesys;

julia> a
a[1:4]
6

I can indeed understand why this might be desired behavior. After all, one might want to set different parameters (or parameter indices) to 0 depending on the run, for example. When excluding these from the parameter list this becomes more troublesome I would say. Although I feel like one should be able to declare an array parameter where only specific elements of the array are non-zero. It happens in my systems a lot, where the array (or matrix) parameters of the systems are in reality relatively sparse.

I will try the suggested things. The first one was already on my list, but I am still not sure what to pass to the prob_func when declaring the EnsembleProblem. I assumed that, in general, “global” variables are not desirable so that I had to extract the relevant parameters from the ODEProblem (or prob.f.sys, or similar).

7

Make a parametric function-like object to store whatever data you need, and then define it as your prob_func, see Methods · The Julia Language for example. This avoids a global variable to store the as.

8

More generally, if you feel there are interface issues in updating parameters within an array parameter I would encourage you to open an issue on ModelingToolkit for further discussion.

1 Like
9

Will do in the future.

10

Sounds good! Array parameters have been clunky for a while. I think MTKv9 is / has addressed this by trying to make how they work more uniform and simple, but I’m not sure if all the planned updates with respect to array parameters are completed now.