Catalyst.jl

Catalyst.jl copied to clipboard

ported from my MTK PR

I'd like to know the best way to integrate this into the lowering code that exists. Should I add a kwarg in the convert(<:ODESystem, rn) code that defaults to the current ode rate law?

@anandijain I'll try to review carefully tomorrow. A few high level comments for now:

1. I'd suggest something like truncated_cme as the name since this isn't actually providing a representation for the (possibly infinite) full CME in general (and there are other approximations one might want to use).
2. Have you looked at https://github.com/kaandocal/FiniteStateProjection.jl? Is there anything from there that would be helpful to use / can be reused to avoid redundant code? (I don't think we should add a dependency on it, but perhaps @kaandocal would have some ideas on components that make sense for handling CMEs in Catalyst.)
3. Could you add some tests where you actually write out the truncated CME ODEs for one or more simple systems, hitting the major reaction types (zero-order, first-order, and second-order), and check they match what is generated from your code (or check that the corresponding ODE rhs functions have the same value for a given u). We want to make sure the generated equations are correct.

Thanks!

This seems similar to what I've been trying to do in FiniteStateProjection.jl, namely converting a reaction network into a (truncated) version of the CME. I did not opt to use ODESystems for this, mainly because I sometimes deal with very state spaces (e.g. 10k states) and thought this might be inefficient, but I never did any benchmarks. The package can convert RNs into Julia functions computing the CME directly using loops, or into matrices that represent the right-hand side of the CME dP/dt = A * P. I think it would be cool to incorporate this as well, ie. add support for conversion to ODESystems along the lines of this PR :)

This seems similar to what I've been trying to do in FiniteStateProjection.jl, namely converting a reaction network into a (truncated) version of the CME. I did not opt to use ODESystems for this, mainly because I sometimes deal with very state spaces (e.g. 10k states) and thought this might be inefficient, but I never did any benchmarks. The package can convert RNs into Julia functions computing the CME directly using loops, or into matrices that represent the right-hand side of the CME dP/dt = A * P. I think it would be cool to incorporate this as well, ie. add support for conversion to ODESystems along the lines of this PR :)

My guess is that with 10K states using loops or matrix-vector form will work better as we may run up against compiler limitations for the size of the MT generated rhs function. This is more what I had in mind; representing the CME as a (possibly infinite) differential-difference equation in many discrete variables. i.e. it is symbolically represented by one equation for an abstract differential-difference equation in the species population variables, and then we can custom generate appropriate functions from that depending on the desired approximation. It would be cool if there was a way to get the symbolic code generated into loops. I guess we need a version of build_function that is designed for difference-type equations/operators. @ChrisRackauckas is there something along those lines used by the PDE conversion tooling?