Unpacking of variables that do not exist

General Usage
, ,
1

I am trying to optimise Catalyst code by omitting many reactions as their rates will, eventually, be put to zero. Briefly, I am considering systems of which many possible reactions will have a rate of 0. To avoid creating these reactions and thereby including (potentially many) reactions/terms that do not change the system dynamics, I check if their numerical value is non-zero and only create the Reaction if it is.

However, when inserting the parameter values when creating the ODEProblem from the ReactionSystem, I now have the issue that the variables are “hard-coded”. This is because, in general, one does not know beforehand (as in, when writing the function) what rates will be 0. As a result, the main issue is that I would like to @unpack the symbols, but that the symbols now do not exist as they will have a rate of 0 and are thus left out.

To illustrate this, consider a super simple system with (logarithmic) growth and death, i.e.

\dot{x} = \eta x(1-x) - \mu x

where \eta and \mu the growth and death rate. Now let us say that the death rate is 0, and thus the reaction x \rightarrow_\mu 0 is left out. However, when I now do

julia> @unpack η, μ = reaction_system
ERROR: ArgumentError: System reaction_system: variable ÎĽ does not exist

I obviously get the error that the variable μ does not exist; as expected. However, is it possible to “conditionally” @unpack, only if the variable exists? As in, can I create a macro/function to do something like:

julia> @unpack_ifexists η, μ = reaction_system; η
η

I know that the @unpack macro can be customized, but I could not find a customization that only unpacks when the variable actually exists.

If I could, then I can write a very general function that, depending on if variables exist in the ReactionSystem or not, assigns them a (non-zero) value when defining the ODEProblem.

I hope my problem and question are clear.

2

Assuming reaction_system is a named tuple, you can merge it with another named tuple to provide defaults for some parameters. For example:

julia> params1 = (a=3, b=4)
(a = 3, b = 4)

julia> (; a, b) = merge((a=0,), params1) # default a=0 is ignored
(a = 3, b = 4)

julia> params2 = (b=4,) # missing a
(b = 4,)

julia> (; a, b) = merge((a=0,), params2) # default a=0 is used
(a = 0, b = 4)

(Here, I’m using the built-in property-destructuring syntax rather than @unpack from UnPack.jl.)

3

The catalyst reaction system has parameters function that returns a vector with the parameters. Maybe that’s enough to directly check if your paramter is there with a simple conditional if my_par in parameters(reaction_system) ... .

4

When property-destructuring syntax (; a, b) = params was implemented in Julia (julia#39285), @simeonschaub commented that “We could think about whether we want to allow specifying default values for properties as well.”

Following up on this, it might be nice to have a syntax like

(; a=default_a, b) = params

as a (non-breaking) extension to this syntax. The merge trick I suggested above to implement defaults only works for named tuples, whereas destructuring works for more arbitrary types.

I filed an issue to track this: default values for destructuring syntax: (; x=default) = y · Issue #51940 · JuliaLang/julia · GitHub

1 Like
5

The problem with this is that it does not seem straightforward to know what my_par is. Consider for example a generalized Lotka-Volterra system without death, for which

julia> parameters(rs)
22-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 η[1]
 η[2]
 A[1, 1]
 A[1, 2]
 A[2, 1]
 A[2, 2]
julia> @unpack η = rs;
julia> η
η[1:2]
julia> typeof(η)
Symbolics.Arr{Num, 1}

It does not seem easy to check whether the variable my_var is in the reaction system as it is not defined before the @unpack command:

julia> A in parameters(rs)
ERROR: UndefVarError: `A` not defined
6

Unfortunately the reaction_system is not a named tuple, but a ReactionSystem from Catalyst. So I cannot simply merge it easily with another named tuple to assign default values and circumvent the issue.

7

I see. I the worst case, since probably you don’t need this in a perfomance-critical part of the code, use a try:

julia> r = @reaction_network Test begin
           k1, A --> B
           k2, B --> C
       end
Model Test
States (3):
  A(t)
  B(t)
  C(t)
Parameters (2):
  k1
  k2

julia> function my_pars(r)
           k1, k2, k3 = try
               @unpack k1, k2, k3 = r
               (k1, k2, k3)
           catch
               @unpack k1, k2 = r
               (k1, k2, 0.0)
           end
           return k1, k2, k3
       end
my_pars (generic function with 1 method)

julia> k1, k2, k3 = my_pars(r)
(k1, k2, 0.0)

ps: Seems that you can convert the output from parameters to an array of regular symbols:

julia> p = Symbol.(parameters(r))
2-element Vector{Symbol}:
 :k1
 :k2

julia> :k1 in p
true

Than you may be able to use a conditional syntax.

1 Like
8

Does :my_var in propertynames(rs) work?

(If my_var, = @unpack rs or (; my_var) = rs works, then rs.my_var should work, in which case :my_var should be a symbol in the propertynames list.)

1 Like
9

Yeah I though about this, but in my (general) case I have many parameters, e.g. \eta, \mu, A, etc. So a try catch for all different combinations of parameters that might be 0 is not really feasible…

But I was also struggling with the scope of variables when using try catch, and your suggestion made me realize I “just” need to my_var = try @unpack my_var = rs catch nothing end, to let my_var be the variable if it exists or not.

10

Yes this works indeed, thanks a lot. This will omit a bunch of try catch clauses that I was thinking of implementing.

11

Better yet, use hasproperty(rs, :my_var) (see hasproperty).

3 Likes
12

Yes the standard property functions are defined for MTK’s system types, so hasproperty and propertynames should work. Feel free to open an issue if they do not.

We should probably remove the scalarization there.