# ModelingToolkit: using intermediate calculation results multiple times

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-using-intermediate-calculation-results-multiple-times/39738>\
**Category:** Modelling & Simulations\
**Created:** [May 19, 2020, 4:42am UTC](https://discourse.julialang.org/t/modelingtoolkit-using-intermediate-calculation-results-multiple-times/39738 "2020-05-19T04:42:19Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![caryan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caryan/32/14984_2.png) [@caryan](https://discourse.julialang.org/u/caryan)\
**Post date:** [May 19, 2020, 4:42am UTC](https://discourse.julialang.org/t/modelingtoolkit-using-intermediate-calculation-results-multiple-times/39738/1 "2020-05-19T04:42:19Z")

</div>

Let’s say I have an expensive scalar that shows up multiple times in my set of differential equations. Is there a way to tell ModelingToolkit to hoist that calculation out and only do it once? If I use an intermediate calculation it seems to just get copied everywhere it is used. E.g. using cosine as an dummy expensive calculation:

```julia
@parameters ω, t
@variables ψ₁(t), ψ₂(t)
@derivatives D'~t

expensive_f = cos(ω*t)

eqs = [D(ψ₁) ~ expensive_f*ψ₂,
       D(ψ₂) ~ expensive_f*ψ₁]

de = ODESystem(eqs)
generate_function(de)[2]

:((var"##MTIIPVar#335", var"##MTKArg#331", var"##MTKArg#332", var"##MTKArg#333")->begin
          @inbounds begin
                  let (ψ₁, ψ₂, ω, t) = (var"##MTKArg#331"[1], var"##MTKArg#331"[2], var"##MTKArg#332"[1], var"##MTKArg#333")
                      var"##MTIIPVar#335"[1] = cos(ω * t) * ψ₂
                      var"##MTIIPVar#335"[2] = cos(ω * t) * ψ₁
                  end
              end
          nothing
      end)

```

Looking through the @code\_llvm there are still two cos calls and but by @code\_native there is only one. Perhaps the compiler is always clever enough to pick this up…

---

<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:** [May 19, 2020, 5:03am UTC](https://discourse.julialang.org/t/modelingtoolkit-using-intermediate-calculation-results-multiple-times/39738/2 "2020-05-19T05:03:58Z")

</div>

> [@caryan](#):
>
> Looking through the @code\_llvm there are still two cos calls and but by @code\_native there is only one. Perhaps the compiler is always clever enough to pick this up…

That’s common subexpression elimination and indeed the compiler will know to do this once.

---

<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:** [May 19, 2020, 5:04am UTC](https://discourse.julialang.org/t/modelingtoolkit-using-intermediate-calculation-results-multiple-times/39738/3 "2020-05-19T05:04:30Z")

</div>

BTW, we will have a way to force CSE in the near future though (via “output variables”)

---

<div class="post-metadata">

**Author:** ![caryan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caryan/32/14984_2.png) [@caryan](https://discourse.julialang.org/u/caryan)\
**Post date:** [May 19, 2020, 7:01pm UTC](https://discourse.julialang.org/t/modelingtoolkit-using-intermediate-calculation-results-multiple-times/39738/4 "2020-05-19T19:01:22Z")

</div>

That sounds promising. Poking around it seems that even changing `expensive_f` to a complex exponential the compiler doesn’t do CSE and call `exp` twice.

---

<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:** [May 19, 2020, 7:15pm UTC](https://discourse.julialang.org/t/modelingtoolkit-using-intermediate-calculation-results-multiple-times/39738/5 "2020-05-19T19:15:51Z")

</div>

Could you give an MWE? We can use that to figure out how to get it fixed in the compiler.

---

<div class="post-metadata">

**Author:** ![caryan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/caryan/32/14984_2.png) [@caryan](https://discourse.julialang.org/u/caryan)\
**Post date:** [May 20, 2020, 9:19pm UTC](https://discourse.julialang.org/t/modelingtoolkit-using-intermediate-calculation-results-multiple-times/39738/6 "2020-05-20T21:19:54Z")

</div>

Of course. Here’s a [notebook](https://gist.github.com/caryan/f4c9afb2dfe4921ca72950bb90b5cfd0) with some of the expressions printed out but an example of the gap is summarized with:

```julia
# repeat with 100 terms to see a much bigger gap

N = 100
@parameters ω, t
@variables ψ[1:N](t)
@derivatives D'~t

drive_amplitude = cos(ω*t)

eqs = D.(ψ) .~ drive_amplitude.*ψ

de = ODESystem(eqs)

println("Without manual CSE")
u = collect(range(0,1; length=N)); du = similar(u);
f = eval(generate_function(de)[2])
@btime f($du, $u, 5e9, 2.0)

println("\n With manual CSE")
ex = generate_function(de)[2]
ex = postwalk(x -> x == :(cos(ω*t)) ? :(drive_amplitude) : x, ex)
pushfirst!(ex.args[2].args[1].args[3].args[1].args[2].args, :(drive_amplitude = cos(ω*t)))
f = eval(ex)
@btime f($du, $u, 2π*5, 2.0)

Without manual CSE
  3.385 μs (0 allocations: 0 bytes)

 With manual CSE
  80.211 ns (0 allocations: 0 bytes)

```

---

<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:** [May 20, 2020, 9:40pm UTC](https://discourse.julialang.org/t/modelingtoolkit-using-intermediate-calculation-results-multiple-times/39738/7 "2020-05-20T21:40:58Z")

</div>

I think it’s because it cannot prove `cos` is pure? @sdanisch might recall the issue.

---

<div class="post-metadata">

**Author:** ![sdanisch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdanisch/32/1406_2.png) [@sdanisch](https://discourse.julialang.org/u/sdanisch)\
**Post date:** [May 20, 2020, 10:42pm UTC](https://discourse.julialang.org/t/modelingtoolkit-using-intermediate-calculation-results-multiple-times/39738/8 "2020-05-20T22:42:03Z")

</div>

I guess it’s not about figuring out the purity (which seems to work), but actually forwarding that information to the LICM optimization:

> <https://github.com/JuliaLang/julia/issues/29285>
>
> Julia 1.0:
> 
> I just realized, that this is a gotcha one easily runs into, espec…ially when using the \`@.\` macro:
> \`\`\`Julia
> a = rand(1000, 1000)
> c = 1.0
> out = similar(a)
> julia\> @btime $(out) .= $a .- sin.($c);
> 5.695 ms (0 allocations: 0 bytes)
> julia\> @btime $(out) .= $a .- sin($c);
> 495.669 μs (0 allocations: 0 bytes)
> \`\`\`
> This likely happens because the compiler can't infer that sin is pure.
> I realize, with having access to the call tree in the new lazy broadcast, we could solve this for a predefined set of functions.
> 
> First trick could be to just overload \`broadcasted\` for known signatures:
> \`\`\`Julia
> Base.Broadcast.broadcasted(::typeof(sin), x::Number) = sin(x)
> @btime $out .= $a .+ sin.($c)
> \`\`\`
> this solves the problem for a chosen set of functions.
> We could also consider, if we introduce a purity trait to make this easier for multiple argument functions:
> \`\`\`Julia
> broadcasted(f, args...) = broadcasted(IsPure(f), f, args...)
> broadcasted(::Pure{true}, f, args...) = f(args...) # should probably not get applied to arrays
> broadcasted(::Pure{false}, f, args...) = Broadcasted(f, args...)
> \`\`\`
> I guess this has been discussed before, but I couldn't really find an issue about it...

I guess we have a PR for this?

> <https://github.com/JuliaLang/julia/pull/32368>
>
> Adds ReadNone and NoUnwind (and thunk) to call sites marked \`@pure\` in the sourc…e code. This allows LLVM to perform LICM and hoisting on these functions, which can be highly profitable.