Plans to support symbolic coefficients in function netstoichmat()?

[sallyseal]

Are there plans to support symbolic coefficients in stoichiometry matrix construction? For example, I have a birth-death process where $M$ is produced in bursts that are geometrically distributed:

@register_symbolic Distributions.Geometric(b)
m = rand(Distributions.Geometric(1 / (1 + b)))

rs = @reaction_network rs begin
    K, 0 --> $m * M
    δ, M --> 0
end

Currently I'm unable to call netstoichmat(rs) without the error "ERROR: Stoichiometry matrices are not supported with symbolic stoichiometry" which means I'm having to hard code m. I'm using Catalyst v13.5.1. Thanks!

[isaacsas]

We can probably modify the function to allow one to provide a keyword to override this error and return the symbolic matrix, but note that no higher-level Catalyst functionality that relies on stoichiometry matrices will work with non-integer matrices.

Out of curiousity, what do you want to use the matrix for?

[sallyseal]

We can probably modify the function to allow one to provide a keyword to override this error and return the symbolic matrix, but note that no higher-level Catalyst functionality that relies on stoichiometry matrices will work with non-integer matrices.

Out of curiousity, what do you want to use the matrix for?

I'm wanting to use FiniteStateProjection.jl where the chemical reaction network contains bursty reactions with geometrically distributed burst sizes. I can use a workaround though, no problem.

[isaacsas]

We can add this in eventually, but right now we are hung up from making new releases due to waiting for ModelingToolkit V9 (and then needing to update Catalyst to v14 for it). Once that is all resolved we can start adding new functionality again.

Something like this can likely do what you want:

function netstoichmat2(::Type{Matrix{T}}, rn::ReactionSystem) where T
           smap = speciesmap(rn)
           nmat = Matrix{T}(undef, numspecies(rn), numreactions(rn))
           nmat .= 0
           for (k, rx) in pairs(reactions(rn))
               for (spec, coef) in rx.netstoich
                   Catalyst.isconstant(spec) && continue
                   nmat[smap[spec], k] = coef
               end
           end
           nmat
       end
ns = netstoichmat2(Matrix{Any}, rs)

giving

julia> ns
1×2 Matrix{Any}:
 rand(Geometric(1 / (1 + b)))  -1