# Getting Catalyst and Oscar to work together

**URL:** <https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299>\
**Category:** Modelling & Simulations\
**Tags:** algebra, symbolics, catalyst\
**Created:** [October 14, 2024, 5:06pm UTC](https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299 "2024-10-14T17:06:21Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![adhalanay](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adhalanay/32/7371_2.png) [@adhalanay](https://discourse.julialang.org/u/adhalanay)\
**Post date:** [October 14, 2024, 5:06pm UTC](https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299/1 "2024-10-14T17:06:22Z")

</div>

This more of an extension of the previous topic [https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/18](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/18).

First some background: commutative algebra is of great use for studying mass-action chemical reaction networks, more specifically questions about. There is a large body of literature about it and a great place to start is Chapter 5 of [this book](https://bookstore.ams.org/cbms-134).

The most feature-rich CAS package seems to me to be Oscar.jl (I think it is on par with Singular, Macaulay2 etc.) I was wondering if there is a way to convert a `NonlinearSystem` coming from a reaction network to a set of polynomials to be fed to Oscar e.g. as generators of some ideal etc? I thing that`NonlinearFunction()` must play some role, but I was unable to find an example of how to use it practically.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [October 14, 2024, 5:15pm UTC](https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299/2 "2024-10-14T17:15:30Z")

</div>

This is a work-in-progress repo that is adding some further network analysis tooling for Catalyst:

> **[GitHub - SciML/CatalystNetworkAnalysis.jl: Network analysis algorithms for reaction networks...](https://github.com/SciML/CatalystNetworkAnalysis.jl)**
>
> Network analysis algorithms for reaction networks modeled using Catalyst.jl

It has some more general Catalyst to polynomial conversion / analysis methods (but I’m not that familiar with the code or if it would work for converting to the Oscar representation you need).

---

<div class="post-metadata">

**Author:** ![adhalanay](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adhalanay/32/7371_2.png) [@adhalanay](https://discourse.julialang.org/u/adhalanay)\
**Post date:** [October 14, 2024, 5:46pm UTC](https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299/3 "2024-10-14T17:46:58Z")

</div>

The package looks very close to what I’m looking for. It seams that is uses Oscar (at least it is imported).

---

<div class="post-metadata">

**Author:** ![vydu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vydu/32/210874_2.png) [@vydu](https://discourse.julialang.org/u/vydu)\
**Post date:** [October 14, 2024, 7:56pm UTC](https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299/4 "2024-10-14T19:56:51Z")

</div>

I’m working on building out some of this functionality in the package that Sam linked. The easiest way to get the species formation rate polynomials as Oscar polynomials, as far as I know, is to build a symbolic function from the species formation rate function. Catalyst has a function called `assemble_oderhs` that essentially gives an array of symbolic expressions that correspond to the right hand side of the chemical reaction network’s ODE. Then, using `Symbolics.build_function`, one can pass an array of these symbolic expressions and variables (which can be accessed using `species(rn)`), and then output a Julia function that will be able to take other types, like Oscar polynomial variables. And then downstream you can do things like build ideals and such.

I’m taking this approach in trying to write a concentration robustness check, though it’s very work in progress. Would be very interested in hearing what other kind of functionality related to this would be useful, and would certainly welcome contributions.

---

<div class="post-metadata">

**Author:** ![adhalanay](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adhalanay/32/7371_2.png) [@adhalanay](https://discourse.julialang.org/u/adhalanay)\
**Post date:** [October 16, 2024, 11:17am UTC](https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299/5 "2024-10-16T11:17:40Z")

</div>

It is a really interesting. The function `assemble_oderhs()` doesn’t seam to be documented but I looked in the source-code. Unfortunately I’m stuck at `Symbolics.build_function`. I couldn’t call the created function. Here is my code

> using Catalyst  
> using Oscar  
> using Symbolics  
> rn=@reaction\_network begin  
> k12, A+A → A+B  
> k21, A+B → A+A  
> k23, A+B → B+B  
> k32, B+B → A+B  
> k13, B+B → A+B  
> k31, A+B → B+B  
> end  
> in1=Catalyst.assemble\_oderhs(rn,species(rn))  
> f=Symbolics.build\_function(in1,species(rn))  
> f(species(rn))

which fails with

```julia
MethodError: objects of type Tuple{Expr, Expr} are not callable
The object of type `Tuple{Expr, Expr}` exists, but no method is defined for this combination of argument types when trying to treat it as a callable object.

Stacktrace:
 [1] top-level scope
   @ In[10]:2

```

Maybe there is something obvious, but I’m not get it.

---

<div class="post-metadata">

**Author:** ![vydu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vydu/32/210874_2.png) [@vydu](https://discourse.julialang.org/u/vydu)\
**Post date:** [October 16, 2024, 12:43pm UTC](https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299/6 "2024-10-16T12:43:37Z")

</div>

Ah, `build_function` will create two expressions that must be evaluated to then get callable functions, see [here](https://symbolics.juliasymbolics.org/v3.5/tutorials/symbolic_functions/#Building-Functions-1) (one of the functions just evaluates the input, the other evaluates and then updates the input array in-place). To get the callable function you’d need

```julia
f1_expr, f2_expr = Symbolics.build_function(in1,species(rn)...)
f = eval(f1_expr)
f(species(rn)...)

```

Note that this would return a Symbolic. To get the output in the form of an Oscar polynomial it would then need to construct a ring and its polynomial variables, and then pass the variables as the input to the function:

```julia
R, polyvars = polynomial_ring(QQ, map(s -> Symbolics.tosymbol(s; escape=false), species(rn)) ) # this second argument just gets symbols from each of the species A(t), B(t) -> :A, :B
f(polyvars...) # this would output an array of Oscar polynomials 

```

---

<div class="post-metadata">

**Author:** ![vydu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vydu/32/210874_2.png) [@vydu](https://discourse.julialang.org/u/vydu)\
**Post date:** [October 16, 2024, 1:10pm UTC](https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299/7 "2024-10-16T13:10:08Z")

</div>

Just realized the rate constant values would still be undefined in this case so the function would error. I think depending on your use case you can either substitute the rate constants directly into the symbolic expression if they are known using something like

```julia
pmap = Dict([:k12 => 1//2, :k13 => 2//1, :k23 => 1//1, :k21 => 1//1, :k32 => 1//1, :k31 => 1//1]) # give rational values for the polynomial ring
pmap = symmap_to_varmap(pmap) # changes the keys to Symbolics variables
in1 = [substitute(eq, pmap) for eq in in1] # substituted expression

```

and then pass this into `build_function`. Or you could add the parameters as variables to the polynomial rings directly (which might be more useful for symbolic solving and stuff).

```julia
R, polyvars = polynomial_ring(QQ, vcat(
    map(s -> Symbolics.tosymbol(s; escape=false), species(rn)),
    map(p -> Symbolics.tosymbol(p; escape=false), parameters(rn))
))

```

---

<div class="post-metadata">

**Author:** ![adhalanay](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adhalanay/32/7371_2.png) [@adhalanay](https://discourse.julialang.org/u/adhalanay)\
**Post date:** [November 18, 2024, 9:34pm UTC](https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299/8 "2024-11-18T21:34:42Z")

</div>

Sorry to reopen the thread so late. But when I try to run  
`f(polyvars...)` (because I want to build the steady-state ideal) I get the error:

```julia
MethodError: no method matching (::var"#31#32")(::QQMPolyRingElem, ::QQMPolyRingElem, ::QQMPolyRingElem, ::QQMPolyRingElem, ::QQMPolyRingElem, ::QQMPolyRingElem, ::QQMPolyRingElem, ::QQMPolyRingElem)
The function `#31` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  (::var"#31#32")(::Any, ::Any)
   @ Main ~/.julia/packages/SymbolicUtils/jf8aQ/src/code.jl:385

Stacktrace:
 [1] top-level scope
   @ In[26]:5

```

---

<div class="post-metadata">

**Author:** ![vydu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vydu/32/210874_2.png) [@vydu](https://discourse.julialang.org/u/vydu)\
**Post date:** [November 20, 2024, 5:43pm UTC](https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299/9 "2024-11-20T17:43:59Z")

</div>

Oh I guess the function is only expecting two arguments because it called `build_function` with just the `species(rn)`. My bad for that. If you want the rate constants as variables in the polynomial ring you’d have to call `build_function(in1, vcat(species(rn), parameters(rn))...)` so it expects them as arguments.

Basically the choices are:

1. build a symbolic function in just the species from the expression that has numerical values substituted for the parameters (in which case call `build_function` with just the `species(rn)`). Then construct the polynomial ring with just the species
2. build a symbolic function in both the species and parameters from the original purely symbolic expression (in which case call `build_function` like above with both `species(rn)`, `parameters(rn)`. Then construct the polynomial ring with both the species and parameters

---

<div class="post-metadata">

**Author:** ![adhalanay](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/adhalanay/32/7371_2.png) [@adhalanay](https://discourse.julialang.org/u/adhalanay)\
**Post date:** [November 20, 2024, 9:14pm UTC](https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299/10 "2024-11-20T21:14:05Z")

</div>

I am trying to replicate the commutative algebraic approach. This is my approach

```julia
R1,vars1=polynomial_ring(QQ,map(p -> Symbolics.tosymbol(p;escape=false),parameters(rn)))
   #K1=fraction_field(R1)
   #R2,vars2=polynomial_ring(K1,map(s -> Symbolics.tosymbol(s;escape=false),species(rn)))

```

Then `f` should define an ideal in `R2`.

---

<div class="post-metadata">

**Author:** ![vydu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vydu/32/210874_2.png) [@vydu](https://discourse.julialang.org/u/vydu)\
**Post date:** [November 20, 2024, 9:34pm UTC](https://discourse.julialang.org/t/getting-catalyst-and-oscar-to-work-together/121299/11 "2024-11-20T21:34:06Z")

</div>

I think then that

```julia
in1=Catalyst.assemble_oderhs(rn,species(rn))
f1, f2 = Symbolics.build_function(in1,vcat(species(rn), parameters(rn))...)
f = eval(f1)

R1,vars1=polynomial_ring(QQ,map(p -> Symbolics.tosymbol(p;escape=false),parameters(rn)))
K1=fraction_field(R1)
R2,vars2=polynomial_ring(K1,map(s -> Symbolics.tosymbol(s;escape=false),species(rn)))
f(vars2..., vars1...)

```

should work to build the ideal that you want.