# DiffEq and SVector

**URL:** <https://discourse.julialang.org/t/diffeq-and-svector/51541>\
**Category:** Modelling & Simulations\
**Created:** [December 9, 2020, 5:55pm UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541 "2020-12-09T17:55:53Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [December 9, 2020, 5:55pm UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/1 "2020-12-09T17:55:53Z")

</div>

Hello,

I’ve been rereading some of the docs and tutorials in SciML and in particular this [Optimizing DiffEq Code](https://tutorials.sciml.ai/html/introduction/03-optimizing_diffeq_code.html), and tried using StaticArrays as in the lorenz example, but on a problem of my own. In my case, the execution time went from about 70 ms using the modifying form with normal vectors, to about 300 ms using SVector. My problem has only 4 ODEs, and is needing about 22000 points to solve (mostly driven by a PeriodicCallback that fires 6000 times).

The solutions appear to be the same, but I do see a difference in the number of points saved, 22241 (modifying form) vs 22533 (SVector) which seems suspicious.

I’m not specifying the solver, but both report CompositeAlgorithm of Vern7, Rodas4 in the sol.alg field.

Due to proprietary stuff in my file I can’t just post it, so before I try to make an MWE I’ll try asking a couple of questions:

Have changes in Julia made the SVector optimization less valuable or possibly counterproductive?

Does the heavy use of PeriodicCallback do something that would cause this to backfire?

I’m using Julia 1.5.3 on Windows.

EDIT: The time point sets are generally different (other than those forced by the DiscreteCallback). I confirmed by `all(isapprox(sol(t),sol2(t);atol=1e-6) for t ∈ 0.0:0.001:0.6)` that the answers are equivalent to within the abstol I specified, but why would the change of mutating or not cause the simulation to take a (slightly) different path?

Thanks

---

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [December 9, 2020, 11:58pm UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/2 "2020-12-09T23:58:33Z")

</div>

> [@klaff](#):
>
> to be the same, but I do see a difference in the number of points saved, 22241 (modifying form) vs 22533 (SVector) which seems suspicious.

Does SVector work for visualization?

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [December 10, 2020, 12:14am UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/3 "2020-12-10T00:14:04Z")

</div>

> [@brett\_knoss](#):
>
> Does SVector work for visualization?

I’m sorry, I don’t understand your question.

I am plotting the solution in either case, but I don’t see why changing the form of the ODE (modifying vector vs returning static vector) would affect the types in the solution object.

---

<div class="post-metadata">

**Author:** ![brett\_knoss](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brett_knoss/32/13050_2.png) [@brett\_knoss](https://discourse.julialang.org/u/brett_knoss)\
**Post date:** [December 10, 2020, 12:15am UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/4 "2020-12-10T00:15:18Z")

</div>

nevermind, I was confused by the term SVector. Thinking of .svg

---

<div class="post-metadata">

**Author:** ![WschW](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wschw/32/6575_2.png) [@WschW](https://discourse.julialang.org/u/WschW)\
**Post date:** [December 10, 2020, 2:39am UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/5 "2020-12-10T02:39:20Z")

</div>

Theres a lot of information that I can’t glean from your post but that could have performance implications.

How large are the `SVector`s? Beyond a certain size they perform less efficiently then `Array`.

Also how type stable is the code? If you put SVectors into containers or use them as feilds in structs it can be easy to end up in a type unstable situation, for example:

```julia
# This is an array type which's eltype is not concrete, and so requires dispatch when you try to do things with elements from it
Array{SVector{3,Float64}}

# This should be used instead
Array{SArray{Tuple{3}, Float64, 1, 3}}

# A different way to do the above which may be easier.
const Float64x3 = typeof(SVector(0,0,0))
Array{Float64x3}

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [December 10, 2020, 4:11am UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/6 "2020-12-10T04:11:57Z")

</div>

As usual in Julia, the best way to deal with types is to _not_ mention them, when possible:

```julia
julia> using StaticArrays

julia> x = SVector(1, 2, 3);

julia> typeof(x)
SVector{3, Int64} (alias for SArray{Tuple{3}, Int64, 1, 3})

julia> x
3-element SVector{3, Int64} with indices SOneTo(3):
 1
 2
 3

julia> v = [x]; # make a Vector containing just `x`

julia> typeof(v)
Vector{SVector{3, Int64}} (alias for Array{SArray{Tuple{3}, Int64, 1, 3}, 1})

```

(The nice printing is from Julia 1.6.)

---

<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:** [December 10, 2020, 4:46am UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/7 "2020-12-10T04:46:05Z")

</div>

> [@klaff](#):
>
> but why would the change of mutating or not cause the simulation to take a (slightly) different path?

Floating point is non-associative and SIMD changes the results, so SVector giving more SIMD can definitely lead to a different floating point result.

I think to understand what’s going on here I’d need to get an MWE. This seems really in the weeds and not a high level issue.

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [December 10, 2020, 4:51am UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/8 "2020-12-10T04:51:38Z")

</div>

It’s a rather simple system - 4 scalar ODEs, so the SVector is 4 elements long.

I was worried about how the `@SVector` macro knew about types, so I switched from the macro form and declared them explicitly to be Float64s and it performed the same. I think the code should be type stable, `u0`, `p`, and `t` are all Float64s and the ODEs are pure math. The ODEs are fairly simple, the complex code is in my callback which has a bunch of rounding, clamping, stuff done in smaller integer sizes, but then the result of the callback is always a Float64 that gets stuffed into one of the parameters. That should be irrelevant to the SVectors, I’ll cut it out and do more comparisons tomorrow.

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [December 10, 2020, 4:53am UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/9 "2020-12-10T04:53:51Z")

</div>

Ah, that makes sense. I’ll trim it down tomorrow.

---

<div class="post-metadata">

**Author:** ![WschW](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wschw/32/6575_2.png) [@WschW](https://discourse.julialang.org/u/WschW)\
**Post date:** [December 10, 2020, 5:05am UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/10 "2020-12-10T05:05:50Z")

</div>

The `@SVector` macro should be perfectly fine to use.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [December 10, 2020, 7:31am UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/11 "2020-12-10T07:31:09Z")

</div>

> [@WschW](#):
>
> ```julia
> # This is an array type which's eltype is not concrete, and so requires dispatch when you try to do things with elements from it
> Array{SVector{3,Float64}}
> 
> ```

```julia
julia> using StaticArrays

julia> isconcretetype(SVector{3,Float64})
true

```

So the `eltype` of `Array{SVector{3,Float64}}` is concrete, although `Array{SVector{3,Float64}}` itself isn’t, because it’s missing the dimensionality of the array.

```julia
julia> isconcretetype(Array{SVector{3,Float64}})
false

julia> isconcretetype(Array{SVector{3,Float64},5})
true

```

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [December 10, 2020, 8:14pm UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/12 "2020-12-10T20:14:49Z")

</div>

MWE:

```julia
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ 88d20350-3af7-11eb-0983-791f4579d925
begin
	using DifferentialEquations
	using StaticArrays
end

# ╔═╡ f2ce46b0-3af7-11eb-2a6a-4716734a17d5
function sim(odef)
	Δt = 0.1e-3
    u0 = [0.0, 0.0, 0.0, 0.0]
	p = [0.178,19.0,10.0,0.25,0.25,20.0,140e3,1.37e3,1.00e3,0.1e-6,0.5,10*1000e-12]
    tspan = (0.0, 0.1)
    prob = ODEProblem(odef, u0, tspan, p)
    sol = solve(prob, reltol=1e-6, abstol=1e-6)
end

# ╔═╡ e6f9b200-3afe-11eb-3df1-d5a083f77901
function ode(u,p,t)
    i, Vsa, Vcfb, Vrc = u
	Lv,Rv,Vv,Rs,Rds,Ga,fc,Ri,Rfb,Cfb,Vdiv,τrc = p
    di = 1/Lv*(Vv-u[1]*(Rv+Rs+Rds))
    dVsa = 2π*fc*(i*Rs*Ga-Vsa)
	Vv = Vdiv*Vsa
	dVcfb = 1/(Cfb*Rfb)*(Vv*Rfb/Ri-Vcfb)
	dVrc = 1/τrc*(Vcfb+Vv-Vrc)
	@SVector [di, dVsa, dVcfb, dVrc]
end

# ╔═╡ 0a24afa0-3aff-11eb-2a88-89b666495af8
@time sol_sv=sim(ode);

# ╔═╡ ec697740-3af7-11eb-2c48-f9c34e498151
function ode!(du,u,p,t)
    i, Vsa, Vcfb, Vrc = u
	Lv,Rv,Vv,Rs,Rds,Ga,fc,Ri,Rfb,Cfb,Vdiv,τrc = p
    di = 1/Lv*(Vv-u[1]*(Rv+Rs+Rds))
    dVsa = 2π*fc*(i*Rs*Ga-Vsa)
	Vv = Vdiv*Vsa
	dVcfb = 1/(Cfb*Rfb)*(Vv*Rfb/Ri-Vcfb)
	dVrc = 1/τrc*(Vcfb+Vv-Vrc)
	du .= [di, dVsa, dVcfb, dVrc]
end

# ╔═╡ 30761000-3afe-11eb-028c-9732640ebe6f
@time sol=sim(ode!);

# ╔═╡ 36b22c50-3aff-11eb-2120-efdff8eada1c
all(isapprox(sol(t),sol_sv(t);atol=1e-6) for t ∈ 0.0:0.001:sol.t[end])

# ╔═╡ Cell order:
# ╠═36b22c50-3aff-11eb-2120-efdff8eada1c
# ╠═30761000-3afe-11eb-028c-9732640ebe6f
# ╠═0a24afa0-3aff-11eb-2a88-89b666495af8
# ╠═f2ce46b0-3af7-11eb-2a6a-4716734a17d5
# ╠═e6f9b200-3afe-11eb-3df1-d5a083f77901
# ╠═ec697740-3af7-11eb-2c48-f9c34e498151
# ╠═88d20350-3af7-11eb-0983-791f4579d925

```

I’ve been testing within Pluto, but just in case there’s some odd interaction, I included this file at the repl and got similar timings:

```julia
julia> using BenchmarkTools

julia> @btime sim(ode!);
  814.061 μs (5164 allocations: 362.88 KiB)

julia> @btime sim(ode);
  2.698 ms (28782 allocations: 1.20 MiB)

```

---

<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:** [December 10, 2020, 9:29pm UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/13 "2020-12-10T21:29:41Z")

</div>

`ode!` can be easily improved too by doing `du .= (di, dVsa, dVcfb, dVrc)` so as to not allocate the array.

---

<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:** [December 10, 2020, 9:30pm UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/14 "2020-12-10T21:30:18Z")

</div>

> [@klaff](#):
>
> ```julia
> u0 = [0.0, 0.0, 0.0, 0.0]
> p = [0.178,19.0,10.0,0.25,0.25,20.0,140e3,1.37e3,1.00e3,0.1e-6,0.5,10*1000e-12]
> tspan = (0.0, 0.1)
> prob = ODEProblem(odef, u0, tspan, p)
> 
> ```

Your initial condition is not a static vector here.

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [December 10, 2020, 9:58pm UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/15 "2020-12-10T21:58:42Z")

</div>

New MWE with improved `ode!` and using `@SVector` for `u0`.

```julia
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ 88d20350-3af7-11eb-0983-791f4579d925
begin
	using DifferentialEquations
	using StaticArrays
end

# ╔═╡ eee69550-3b30-11eb-3592-076654279c24
function sim_sv(odef)
	Δt = 0.1e-3
    u0 = @SVector[0.0, 0.0, 0.0, 0.0]
	p = [0.178,19.0,10.0,0.25,0.25,20.0,140e3,1.37e3,1.00e3,0.1e-6,0.5,10*1000e-12]
    tspan = (0.0, 0.1)
    prob = ODEProblem(odef, u0, tspan, p)
    sol = solve(prob, reltol=1e-6, abstol=1e-6)
end

# ╔═╡ f2ce46b0-3af7-11eb-2a6a-4716734a17d5
function sim(odef)
	Δt = 0.1e-3
    u0 = [0.0, 0.0, 0.0, 0.0]
	p = [0.178,19.0,10.0,0.25,0.25,20.0,140e3,1.37e3,1.00e3,0.1e-6,0.5,10*1000e-12]
    tspan = (0.0, 0.1)
    prob = ODEProblem(odef, u0, tspan, p)
    sol = solve(prob, reltol=1e-6, abstol=1e-6)
end

# ╔═╡ e6f9b200-3afe-11eb-3df1-d5a083f77901
function ode(u,p,t)
    i, Vsa, Vcfb, Vrc = u
	Lv,Rv,Vv,Rs,Rds,Ga,fc,Ri,Rfb,Cfb,Vdiv,τrc = p
    di = 1/Lv*(Vv-u[1]*(Rv+Rs+Rds))
    dVsa = 2π*fc*(i*Rs*Ga-Vsa)
	Vv = Vdiv*Vsa
	dVcfb = 1/(Cfb*Rfb)*(Vv*Rfb/Ri-Vcfb)
	dVrc = 1/τrc*(Vcfb+Vv-Vrc)
	@SVector [di, dVsa, dVcfb, dVrc]
end

# ╔═╡ 0a24afa0-3aff-11eb-2a88-89b666495af8
@time sol_sv=sim_sv(ode);

# ╔═╡ ec697740-3af7-11eb-2c48-f9c34e498151
function ode!(du,u,p,t)
    i, Vsa, Vcfb, Vrc = u
	Lv,Rv,Vv,Rs,Rds,Ga,fc,Ri,Rfb,Cfb,Vdiv,τrc = p
    di = 1/Lv*(Vv-u[1]*(Rv+Rs+Rds))
    dVsa = 2π*fc*(i*Rs*Ga-Vsa)
	Vv = Vdiv*Vsa
	dVcfb = 1/(Cfb*Rfb)*(Vv*Rfb/Ri-Vcfb)
	dVrc = 1/τrc*(Vcfb+Vv-Vrc)
	du .= (di, dVsa, dVcfb, dVrc)
end

# ╔═╡ 30761000-3afe-11eb-028c-9732640ebe6f
@time sol=sim(ode!);

# ╔═╡ 36b22c50-3aff-11eb-2120-efdff8eada1c
all(isapprox(sol(t),sol_sv(t);atol=1e-6) for t ∈ 0.0:0.001:sol.t[end])

# ╔═╡ Cell order:
# ╠═36b22c50-3aff-11eb-2120-efdff8eada1c
# ╠═30761000-3afe-11eb-028c-9732640ebe6f
# ╠═0a24afa0-3aff-11eb-2a88-89b666495af8
# ╠═eee69550-3b30-11eb-3592-076654279c24
# ╠═f2ce46b0-3af7-11eb-2a6a-4716734a17d5
# ╠═e6f9b200-3afe-11eb-3df1-d5a083f77901
# ╠═ec697740-3af7-11eb-2c48-f9c34e498151
# ╠═88d20350-3af7-11eb-0983-791f4579d925

```

and results at REPL (`ode!` allocations reduced by replacing array w/ tuple):

```julia
julia> include("svector_mwe.jl")
 13.307219 seconds (37.73 M allocations: 1.982 GiB, 5.72% gc time)
 23.287191 seconds (70.86 M allocations: 3.400 GiB, 5.05% gc time)
true

julia> using BenchmarkTools

julia> @btime sim(ode!);
  525.576 μs (3527 allocations: 183.83 KiB)

julia> @btime sim_sv(ode);
  2.454 ms (28168 allocations: 1.15 MiB)

```

---

<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:** [December 10, 2020, 10:37pm UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/16 "2020-12-10T22:37:35Z")

</div>

Because you’re using the default algorithm, autodiff is turned off and thus you’re hitting [https://github.com/JuliaDiff/FiniteDiff.jl/issues/120](https://github.com/JuliaDiff/FiniteDiff.jl/issues/120) . If you specify `Rodas5()` as the algorithm (since the equations are stiff and do the switch), then you get ~200us which is the fastest.

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [December 10, 2020, 11:11pm UTC](https://discourse.julialang.org/t/diffeq-and-svector/51541/17 "2020-12-10T23:11:57Z")

</div>

Thank you!

I’ve applied the learnings to my original problem. Just for the record, I started this thread with what I thought was a fine solution needing 70 ms, tried to make it faster (which gave me 300 ms), then with changing algorithm and removing an array allocation got to 30 ms, and then adding SVector like I was trying to do am now down to 14 ms!