Issue with StaticArrays and PDMP

General Usage
1

Dear All,

I am trying to make PiecewiseDeterministicMarkovProcesses.jl (I am the creator) work with StaticArrays.jl. To this end, I try running the example file by appending

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:

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.

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

2

Using the Debugger, I had difficulties setting up a breakpoint. I don’t really understand what calls the convert method yet.

3

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.

4

Thank you a lot!

May I ask you how you discovered it?

5

I had a hunch

6

it is a good one!

7

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)?

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:

# 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)
8

Is there an SVector version?

9

You mean an out of place version?

In this case, no.

10

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.

1 Like
11

Well, this is a diffEq (M)WE: it is just an example with iterator. I will try to simplify it if possible.

12

OK, I think it has to do with the iterator and MVectors. First I show the result of @profiler

00
002702Ă—1922 837 KB

13

Then here is a MWE

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:

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

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)
14

I made an MWE out of it.

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…

15

Note that this shows no difference:

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...

16

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.

17

it seems to go away for this

integrator = init(prob, Euler(), callback = cb,  save_everystep = false, dt = 0.01 )
integrators = init(probs, Euler(), callback = cb,  save_everystep = false, dt = 0.01 )