# Detecting Function Composition in Symbolics.jl

**URL:** <https://discourse.julialang.org/t/detecting-function-composition-in-symbolics-jl/115885>\
**Category:** Modelling & Simulations\
**Tags:** modelingtoolkit, symbolics, sympy\
**Created:** [June 19, 2024, 8:40pm UTC](https://discourse.julialang.org/t/detecting-function-composition-in-symbolics-jl/115885 "2024-06-19T20:40:56Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![nicrummel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicrummel/32/6304_2.png) [@nicrummel](https://discourse.julialang.org/u/nicrummel)\
**Post date:** [June 19, 2024, 8:40pm UTC](https://discourse.julialang.org/t/detecting-function-composition-in-symbolics-jl/115885/1 "2024-06-19T20:40:56Z")

</div>

# Motivation

I have a system of differential equations  
\dot{U} = F(U,w)  
where U\in \mathbb{R}^D are dependent variables and w \in\mathbb{R}^J are parameters. And I can quickly build (in-place) Julia functions from [Symbolics.jl](https://docs.sciml.ai/Symbolics/stable/) and [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) of the following functions of U and w:  
F(U,w), \; \nabla\_uF(U,w), \;\nabla\_wF(U,W), and \nabla\_w[\nabla\_uF](U,w)  
I then have data \{u^{(m)}\}\_{m=1}^M, where M\in \mathbb{N} is large. I don’t want assume anything about F a priori, but I need to evaluate all four of these functions at every u^{(m)} for different w during an optimization routine. This is quickly becoming expensive…

## Goal

I would like to further speed up evaluating these functions by auto detecting  
\nabla\_uF[u\_m,w] = h(g(u), w)

### Concrete Example (Hindmarsh Rose)

This example currently is causing too much of slow down, but larger systems are very slow. Hopefully it can help ground the discussion:  
\dot{U} = F(U,w) = \left[\begin{array}{l} {w}\_1 U\_2+w\_2 U\_1^3+w\_3 U\_1^2+w\_4 U\_3 \\ w\_5+w\_6 U\_1^2+w\_7 U\_2 \\ w\_8 U\_1+w\_9+w\_{10} U\_3 \end{array}\right]  
\nabla\_U F(U,w) = \left[\begin{array}{ccc} 3 U\_1^2 w\_2 + 2U\_1w\_3 & {w}\_1 & w\_4 \\ 2U\_1 w\_6 & w\_7 & 0 \\ w\_8 & 0 & w\_{10} \end{array}\right]  
I want to programmatically recognize that  
g(U) = \left[\begin{array}{c} 3U\_1^2\\ 2U\_1 \end{array}\right]  
Then we can precompute \{y^{(m)}\}\_{m=1}^M := \{g(u^{(m)})\}\_{m=1}^M and then quickly compute  
\begin{align\*}\nabla\_UF(u^{(m)},w) &= h(g(u^{(m)},w) = h(y^{(m)},w) \\ &= \left[\begin{array}{ccc} y^{(m)}\_1 w\_2 + y\_2^{(m)}w\_3 & {w}\_1 & w\_4 \\ y\_2^{(m)} w\_6 & w\_7 & 0 \\ w\_8 & 0 & w\_{10} \end{array}\right] \end{align\*}

## Ideas on how to start

- There should be a way to isolate each expressions that contain just U and scalars. We would need to leverage when operators are/aren’t associative, commutative, etc.. In the example this would result in [3U\_1^2, 2U\_1,2U\_1].
- Then we define a vector of the unique expressions. In the example: [3U\_1^2, 2U\_1,] .
- Substitute the y\_i in the original equation for each expression to define h(y,w).

I am having a hard time getting started in the Symbolics API with what I can see from the documentation. Perhaps there is way to do this in sympy through [SymPy.jl](https://jverzani.github.io/SymPyCore.jl/dev/)? Any help or starting points would be a great help!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 20, 2024, 10:44am UTC](https://discourse.julialang.org/t/detecting-function-composition-in-symbolics-jl/115885/2 "2024-06-20T10:44:47Z")

</div>

This is just common subexpression elimination?

---

<div class="post-metadata">

**Author:** ![nicrummel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicrummel/32/6304_2.png) [@nicrummel](https://discourse.julialang.org/u/nicrummel)\
**Post date:** [June 22, 2024, 12:33am UTC](https://discourse.julialang.org/t/detecting-function-composition-in-symbolics-jl/115885/3 "2024-06-22T00:33:25Z")

</div>

@ChrisRackauckas I think so! This is readily available in [pysym](https://docs.sympy.org/latest/modules/rewriting.html) (link goes to relevant documentation). It looks like others have used [SymbolicUtils.jl](https://docs.sciml.ai/SymbolicUtils/stable/) and the (unmaintained) [SymbolicCodegen.jl](https://github.com/cscherrer/SymbolicCodegen.jl/tree/main) to do CSE. It can be found in this [show and tell](https://github.com/JuliaMath/MeasureTheory.jl/discussions/43).

My code already has symbolic expression from [Symbolics.jl](https://docs.sciml.ai/Symbolics/stable/), Is there a way to use the high level module to:

1. Find the common subexpressions
2. Build Julia functions of each subexpression
3. Substitute the expressions back into the original
4. Build a function out of the final expression

Or should I bring in some [SymPy.jl](https://jverzani.github.io/SymPyCore.jl/dev/)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 22, 2024, 3:08pm UTC](https://discourse.julialang.org/t/detecting-function-composition-in-symbolics-jl/115885/4 "2024-06-22T15:08:45Z")

</div>

If you use `build_function`, then you just add `cse=true`. We don’t do it by default though because there’s an unknown bug according to @shashi IIRC, but I don’t know what it is or how to trigger it so if you can help figure out what that is we can just make it standard (since it’s a standard codegen optimization, even Julia will run CSE on the generated code).

---

<div class="post-metadata">

**Author:** ![nicrummel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicrummel/32/6304_2.png) [@nicrummel](https://discourse.julialang.org/u/nicrummel)\
**Post date:** [July 3, 2024, 4:49pm UTC](https://discourse.julialang.org/t/detecting-function-composition-in-symbolics-jl/115885/5 "2024-07-03T16:49:53Z")

</div>

Hey @ChrisRackauckas,

Just following up on this CSE discussion. (I know that most of this explanation is not necessary for you, but I am putting it here so if someone else is trying to follow the discussion they can). Here is what I have done:

- I looked into turning on the key word argument `cse`, for `build_function`. The outputs are two functions, `(f,f!)` , where the latter is an in place function. The allocating function `f` appears to have the proper common subexpression elimination, while `f!` does not. This may be the “bug”, that was mentioned. I probably can go make the necessary changed in Symbolics.jl to have this function properly, but have not done so yet.

- For my use case, even if `f!` worked properly this is not the full solution. Turning `cse=true` would produce (using the notation above):

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 5, 2024, 6:06pm UTC](https://discourse.julialang.org/t/detecting-function-composition-in-symbolics-jl/115885/6 "2024-07-05T18:06:05Z")

</div>

> [@nicrummel](#):
>
> I looked into turning on the key word argument `cse`, for `build_function`. The outputs are two functions, `(f,f!)` , where the latter is an in place function. The allocating function `f` appears to have the proper common subexpression elimination, while `f!` does not. This may be the “bug”, that was mentioned. I probably can go make the necessary changed in [Symbolics.jl](https://juliahub.com/ui/Packages/General/Symbolics) to have this function properly, but have not done so yet.

Yeah that would be a bug. Open an issue?

> [@nicrummel](#):
>
> This not the desired result for my use case. I need something slightly different. Because I am trying to run optimization methods to find the “best” w^_w∗w^_ I want to evaluate fff at fixed uuu’s but many different www. I need the subexpression that only contain uuu and scalars, g(u)g(u)g(u), so I can produce {y^{(m)}}_{m=1}^M := {g(u^{(m)}}_{m=1}^M{y(m)}Mm=1:={g(u(m)}Mm=1{y^{(m)}}_{m=1}^M := {g(u^{(m)}}_{m=1}^M. With these precomputed y^{(m)}y(m)y^{(m)}, I can call h(y,w)h(y,w)h(y,w) which should have fewer allocations and runtime.
> 
> To accomplish this goal I have tried to use some prebuilt tools in [Symbolics.jl](https://juliahub.com/ui/Packages/General/Symbolics) and SymPy (through [SymPy.jl](https://juliahub.com/ui/Packages/General/SymPy)). SymPy is a pain for my application because I am taking [Symbolics.jl](https://juliahub.com/ui/Packages/General/Symbolics) types turning them into SymPy types then back. So there is a lot of plumbing which adds overhead and decreases readability. Furthermore, neither [Symbolilics.jl](https://juliahub.com/ui/Packages/General/Symbolilics) or SymPy gives me the result I am looking for. These tools are not good at isolating one variable from the other. I `dev`’d [Symbolics.jl](https://juliahub.com/ui/Packages/General/Symbolics) module to extend the `coeff` function to handle a broader class of expressions, but it honestly is patch work. I am looking for something more than just subexpressions that are coefficients.
> 
> Do you have any suggestions on how to proceed? Thanks again for your time.

That probably needs an issue. I don’t know if we cover that case right now.