# Issue with StaticArrays and PDMP
**URL:** https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392
**Category:** General Usage
**Tags:** diffeq
**Created:** [May 20, 2019, 7:54am UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392 "2019-05-20T07:54:28Z")
**Posts on this page:** 17
**Page:** 1
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 20, 2019, 7:54am UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/1 "2019-05-20T07:54:28Z")
Dear All,
I am trying to make [PiecewiseDeterministicMarkovProcesses.jl](https://github.com/rveltz/PiecewiseDeterministicMarkovProcesses.jl) (I am the creator) work with [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl). To this end, I try running the example [file](https://github.com/rveltz/PiecewiseDeterministicMarkovProcesses.jl/blob/master/examples/issueForwardDiff.jl) by appending
```julia
using StaticArrays
sxc0 = @MVector [1.0]
sxd0 = @MVector [1]
res = @time PDMP.chv_diffeq!(sxc0, sxd0, F_fd!, R_fd!,Dummy!, nu_fd, parms, 0.0, 10.0,false; n_jumps = 30, ode = Tsit5(),save_positions = (false, false))
```
However, I am unable to even create a problem as I get the error:
```julia
ERROR: DimensionMismatch("expected input array of length 1, got length 2")
Stacktrace:
[1] dimension_mismatch_fail(::Type, ::Array{Float64,1}) at /Users/rveltz/.julia/packages/StaticArrays/VyRz3/src/convert.jl:18
[2] convert at /Users/rveltz/.julia/packages/StaticArrays/VyRz3/src/convert.jl:23 [inlined]
[3] PiecewiseDeterministicMarkovProcesses.PDMPProblem{Float64,Int64,MArray{Tuple{1},Float64,1,1},MArray{Tuple{1},Int64,1,1},Array{Int64,2},Array{Float64,1},typeof(F_fd!),typeof(R_fd!),typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)}(::MArray{Tuple{1},Float64,1,1}, ::MArray{Tuple{1},Int64,1,1}, ::typeof(F_fd!), ::typeof(R_fd!), ::typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy), ::Array{Int64,2}, ::Array{Float64,1}, ::Float64, ::Float64, ::Bool, ::Bool, ::Bool) at /Users/rveltz/work/prog_gd/julia/dev/PiecewiseDeterministicMarkovProcesses.jl/src/utils.jl:44
[4] #chv_diffeq!#70 at /Users/rveltz/work/prog_gd/julia/dev/PiecewiseDeterministicMarkovProcesses.jl/src/chvdiffeq.jl:68 [inlined]
[5] (::getfield(PiecewiseDeterministicMarkovProcesses, Symbol("#kw##chv_diffeq!")))(::NamedTuple{(:n_jumps, :ode, :save_positions),Tuple{Int64,Tsit5,Tuple{Bool,Bool}}}, ::typeof(chv_diffeq!), ::MArray{Tuple{1},Float64,1,1}, ::MArray{Tuple{1},Int64,1,1}, ::typeof(F_fd!), ::typeof(R_fd!), ::typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy), ::Array{Int64,2}, ::Array{Float64,1}, ::Float64, ::Float64, ::Bool) at none:0
[6] top-level scope at none:0
```
Please note that the following works!! So the issue seems to be related to `sxc0`.
```julia
res = @time PDMP.chv_diffeq!(xc0, sxd0, F_fd!, R_fd!,Dummy!, nu_fd, parms, 0.0, 10.0,false; n_jumps = 30, ode = Tsit5(),save_positions = (false, false))
```
I spend some time trying to debug it but I am stuck. I would appreciate if one can give me a little push.
Thank you for your help,
Best regards
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 26, 2019, 6:43pm UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/2 "2019-05-26T18:43:08Z")
Using the `Debugger`, I had difficulties setting up a breakpoint. I don’t really understand what calls the `convert` method yet.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 26, 2019, 8:05pm UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/3 "2019-05-26T20:05:53Z")
Your issue is that `rate` is set to `zeros(Tc,size(nu,1))` but `xc::vectype_xc ` and `rate::vectype_xc` so it’ll attempt to convert to the size 1 MArray to match the first arg, while `nu` is length 2.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 26, 2019, 8:34pm UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/4 "2019-05-26T20:34:32Z")
Thank you a lot!
May I ask you how you discovered it?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 26, 2019, 8:42pm UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/5 "2019-05-26T20:42:49Z")
I had a hunch
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 26, 2019, 8:57pm UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/6 "2019-05-26T20:57:45Z")
it is a good one!
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 27, 2019, 10:54am UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/7 "2019-05-27T10:54:18Z")
I posted a fix for this but it seems far from efficient. Basically, to bypass the lack of methods such as `similar(xc, 10)` for static vectors, I let the user pass the required static vectors and operate inplace. This can be seen in the following example.
However, running everything with StaticArrays seems type unstable as it generates a lot of allocations. From using `@profiler` or the option `--track-allocation=user`, I have the impression that the issue lies with `broadcast` in `OrdinaryDiffEq/src/perform_step/low_order_rk_perform_step.jl`.
Can anyone reproduce my timings please (see below)?
```julia
using PiecewiseDeterministicMarkovProcesses, LinearAlgebra, Random, DifferentialEquations
const PDMP = PiecewiseDeterministicMarkovProcesses
function F_fd!(ẋ, xc, xd, t, parms)
# vector field used for the continuous variable
if mod(xd[1], 2) == 0
ẋ[1] = 1 + xd[1]
else
ẋ[1] = -xc[1]
end
nothing
end
rate_tcp(x) = 1/x
function R_fd!(rate, xc, xd, t, parms, sum_rate::Bool)
rate[1] = 1.0 + rate_tcp(xd[1]) * xc[1]
if sum_rate == false
return 0., 0.
else
return sum(rate), 0.
end
end
Dummy! = PDMP.Delta_dummy
using StaticArrays
sxc0 = @MVector [1.0]
sxd0 = @MVector [1]
ratevec = similar(sxc0, Size(2))
sxc0_e = similar(sxc0, Size(2))
```
Here is the comparison:
```julia
# The following runs with generic vectors
# it runs in 0.006755 seconds (458 allocations: 31.281 KiB)
ress = @time PDMP.chv_diffeq!(sxc0, sxd0, F_fd!, R_fd!,Dummy!, nu_fd, parms, 0.0, 10000.0,false; n_jumps = 3000, ode = Tsit5(),save_positions = (false, false), rate = ratevec)#, xc0_extended = sxc0_e)
# The following runs with static vectors
# it runs 0.505335 seconds (7.58 M allocations: 130.156 MiB, 7.50% gc time)
ress = @time PDMP.chv_diffeq!(sxc0, sxd0, F_fd!, R_fd!,Dummy!, nu_fd, parms, 0.0, 10000.0,false; n_jumps = 3000, ode = Tsit5(),save_positions = (false, false), rate = ratevec, xc0_extended = sxc0_e)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 27, 2019, 1:09pm UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/8 "2019-05-27T13:09:18Z")
Is there an SVector version?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 27, 2019, 1:23pm UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/9 "2019-05-27T13:23:55Z")
You mean an out of place version?
In this case, no.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 27, 2019, 3:16pm UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/10 "2019-05-27T15:16:30Z")
It would be helpful if you could either come up with a DiffEq MWE of this, or show that it’s just MArray broadcast that’s unstable. I think it’s possible to just show with the latter, but I haven’t found out how yet… but I’ve mentioned to @andyferris this behavior before. We just need an MWE to clarify what’s going on.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 27, 2019, 3:57pm UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/11 "2019-05-27T15:57:54Z")
Well, this is a diffEq (M)WE: it is just an example with iterator. I will try to simplify it if possible.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 31, 2019, 7:43am UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/12 "2019-05-31T07:43:32Z")
OK, I think it has to do with the iterator and `MVectors`. First I show the result of `@profiler`
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [May 31, 2019, 7:48am UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/13 "2019-05-31T07:48:44Z")
Then here is a MWE
```julia
function F_fd!(ẋ, xc, xd, t, parms)
# vector field used for the continuous variable
if mod(xd[1], 2) == 0
ẋ[1] = 1 + xd[1]
else
ẋ[1] = -xc[1]
end
nothing
end
```
These 2 have the same behaviour, allocate ~400 things:
```julia
using DifferentialEquations
prob = ODEProblem((du,u,p,t)->F_fd!(du,u,xd0,t,p),xc0,(0.0,10.0),parms)
sol = @time solve(prob, Tsit5())
probs = ODEProblem((du,u,p,t)->F_fd!(du,u,sxd0,t,p),sxc0,(0.0,10.0),parms)
sols = @time solve(prob, Tsit5())
```
Now using iterators
```julia
cb = DiscreteCallback((u,t,integrator)->false, integrator -> nothing, save_positions = (false, false))
integrator = init(prob, Tsit5(), callback = cb, save_everystep = false )
integrators = init(probs, Tsit5(), callback = cb, save_everystep = false )
function solveODE!(integ)
while integ.t<10.0
step!(integ)
# @show integ.t
end
end
@time solveODE!(integrator) # 0.000026 seconds (112 allocations: 1.844 KiB)
@time solveODE!(integrators) #0.000505 seconds (5.13 k allocations: 90.156 KiB)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 31, 2019, 10:01am UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/14 "2019-05-31T10:01:43Z")
I made an MWE out of it.
```julia
using DifferentialEquations, StaticArrays
function F_fd!(ẋ, xc, xd, t, parms)
# vector field used for the continuous variable
if mod(xd[1], 2) == 0
ẋ[1] = 1 + xd[1]
else
ẋ[1] = -xc[1]
end
nothing
end
sxc0 = @MVector [1.0]
xc0 = [1.0]
xd0 = [1.0]
sxd0 = @MVector [1.0]
parms = nothing
prob = ODEProblem((du,u,p,t)->F_fd!(du,u,xd0,t,p),xc0,(0.0,10.0),parms)
sol = @time solve(prob, Tsit5())
probs = ODEProblem((du,u,p,t)->F_fd!(du,u,sxd0,t,p),sxc0,(0.0,10.0),parms)
sols = @time solve(prob, Tsit5())
cb = DiscreteCallback((u,t,integrator)->false, integrator -> nothing, save_positions = (false, false))
integrator = init(prob, Tsit5(), callback = cb, save_everystep = false )
integrators = init(probs, Tsit5(), callback = cb, save_everystep = false )
function solveODE!(integ)
while integ.t<10.0
step!(integ)
# @show integ.t
end
end
@time solveODE!(integrator) # 0.000026 seconds (112 allocations: 1.844 KiB)
@time solveODE!(integrators) #0.000505 seconds (5.13 k allocations: 90.156 KiB)
```
I gotta say, I have no idea what this is. This seems like compiler nonsense since they are running almost the same code…
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [June 1, 2019, 8:44am UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/15 "2019-06-01T08:44:48Z")
Note that this shows no difference:
```julia
using DifferentialEquations
prob = ODEProblem((du,u,p,t)->F_fd!(du,u,xd0,t,p),xc0,(0.0,10.0),parms)
sol = @time solve(prob, Tsit5())
probs = ODEProblem((du,u,p,t)->F_fd!(du,u,sxd0,t,p),sxc0,(0.0,10.0),parms)
sols = @time solve(prob, Tsit5())
```
which is probably a hint to the issue. Also, the profiler suggests to look at the function `low_order_rk_perform...`
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [June 1, 2019, 11:02am UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/16 "2019-06-01T11:02:00Z")
They are hitting the same piece of code in `low_order_rk_perform` so it’s some “compiler optimizes one time but not the other” thing. I am not sure how to debug that.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [June 1, 2019, 12:34pm UTC](https://discourse.julialang.org/t/issue-with-staticarrays-and-pdmp/24392/17 "2019-06-01T12:34:14Z")
it seems to go away for this
```julia
integrator = init(prob, Euler(), callback = cb, save_everystep = false, dt = 0.01 )
integrators = init(probs, Euler(), callback = cb, save_everystep = false, dt = 0.01 )
```