# Issue with PDMP and Forwardiff - DifferentialEquation

Hi,

I repost here an issue posted in DifferentialEquations because I am not sure this issue concerns them entirely.

In the following example, Forwardiff which is called from Rosenbrock23() fails. I spent quite some time to debug it. Please note that Forwardiff works with PDMP.jl. An example is the here.

Hence, I am a bit stuck here. Maybe some of you have an idea.

``````using PDMP, LinearAlgebra, Random, DifferentialEquations

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] = -1
end
nothing
end

rate_tcp(x) = 1/x

function R_fd!(rate, xc, xd, t, parms, sum_rate::Bool)
#replacing xc[1] -> xd[1] in the next line removes the error
rate[1] = rate_tcp(xc[1]) * xd[1]
rate[2] = 0.0
if sum_rate==false
return 0., 0.
else
return sum(rate), 0.
end
end

xc0 = [ 1.0 ]
xd0 = [1, 1]

nu_fd = [[1 0];[0 -1]]
parms = [0.0]

# works:
res =  @time PDMP.pdmp!(xc0, xd0, F_fd!, R_fd!, nu_fd, parms, 0.0, 1.0, n_jumps = 3,   ode = Tsit5())
# fail because of autodiff
res =  @time PDMP.pdmp!(xc0, xd0, F_fd!, R_fd!, nu_fd, parms, 0.0, 1.0, n_jumps = 3,   ode = Rosenbrock23())
``````

The error is the following:

``````ERROR: LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PDMP, Symbol("##31#33")){PDMP.PDMPProblem{Float64,Int64,Array{Float64,1},Array{Int64,1},Array{Int64,2},Array{Float64,1},typeof(F_fd!),typeof(R_fd!),typeof(PDMP.Delta_dummy)}},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Nothing},Float64},Float64,2})
Closest candidates are:
Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:185
Float64(::T<:Number) where T<:Number at boot.jl:725
Float64(::Int8) at float.jl:60
...
``````

@ChrisRackauckas As I told you on gitter, I cannot understand why your trick allocates in the case of PDMP.jl.

To explain the issue, I have a MWE.

``````using Revise, PiecewiseDeterministicMarkovProcesses, LinearAlgebra, Random, DifferentialEquations, Sundials
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] = -1.0
end
nothing
end

rate_tcp(x) = 1/x

function R_fd!(rate, xc, xd, t, parms, sum_rate::Bool)
# @show typeof(rate)
rate[1] = 0.
rate[1] = 1.0 + rate_tcp(xd[1])
if sum_rate == false
return 0., 0.
else
return sum(rate), 0.
end
end

Dummy! = PDMP.Delta_dummy

xc0 = [ 1.0 ]
xd0 = [ 1  ]

nu_fd = [[1 0];[0 -1]]
parms = [0.0]
``````

We create a problem with

``````problem = PDMP.chv_diffeq!(xc0, xd0, F_fd!, R_fd!,Dummy!, nu_fd, parms, 0.0, 10000.0 ;n_jumps = 3000, ode = Tsit5(), save_positions=(false,false), return_pb = true)
``````

In order to show the issue, we have to change the callable struct. I did not put it on PDMP.jl in order to avoid changing too much the master branch.

``````# callable struct for the CHV method
function (prob::PDMP.PDMPProblem)(xdot, x, data, t)
# @show typeof(xdot) typeof(x)
if prob.chv # this is to simulate with the CHV method
tau = x[end]
# rate = similar(x,length(prob.rate)) #This is to use autodiff but it slows things down
tmp = PDMP.get_rate(prob.rateCache, eltype(x))
_tmp = reinterpret(eltype(x), tmp)
sr = prob.pdmpFunc.R(_tmp,x,prob.xd,tau,prob.parms,true)[1]
prob.pdmpFunc.F(xdot,x,prob.xd,tau,prob.parms)
xdot[end] = 1.0
@inbounds for i in eachindex(xdot)
xdot[i] = xdot[i] / sr
end
else # this is to simulate with rejection method
prob.pdmpFunc.F(xdot, x, prob.xd, t, prob.parms)
end
nothing
end
``````

We can then test the call:

``````using BenchmarkTools
dx = similar(xc0)
@btime problem(\$dx,\$xc0,\$parms,0.)
``````

The result is

``````44.763 ns (1 allocation: 32 bytes)
``````

whereas in the docs it gives 0 allocations. Does anybody have an idea?

`xd = [1]` which is an integer, and then you have `ẋ[1] = 1 + xd[1]` which adds an integer, while `ẋ[1] = -1.0` adds a float. If I had to guess that would be it, but other than that I don’t have a good guess.

Well, if I dont use your trick by invoking the code on master, ie which uses:

``````function (prob::PDMPProblem)(xdot, x, data, t)
# @show typeof(xdot) typeof(x)
if prob.chv # this is to simulate with the CHV method
tau = x[end]
# rate = similar(x,length(prob.rate)) #This is to use autodiff but it slows things down
# tmp = get_rate(prob.rateCache, eltype(x))
# _tmp = reinterpret(eltype(x), tmp)
sr = prob.pdmpFunc.R(prob.rateCache.rate,x,prob.xd,tau,prob.parms,true)[1]
prob.pdmpFunc.F(xdot,x,prob.xd,tau,prob.parms)
xdot[end] = 1.0
@inbounds for i in eachindex(xdot)
xdot[i] = xdot[i] / sr
end
else # this is to simulate with rejection method
prob.pdmpFunc.F(xdot, x, prob.xd, t, prob.parms)
end
nothing
end
``````

then

``````julia> @btime problem(\$dx,\$xc0,\$parms,0.)
26.434 ns (0 allocations: 0 bytes)
``````

I have tracked allocations `julia --track-allocation=user` using

``````problem(dx,xc0,parms,0.)

Profile.clear_malloc_data()
problem(dx,xc0,parms,0.)
problem(dx,xc0,parms,0.)

``````

and I get the following:

``````        - # callable struct for the CHV method
- function (prob::PDMPProblem)(xdot, x, data, t)
- 	# @show typeof(xdot) typeof(x)
- 	# if prob.chv # this is to simulate with the CHV method
13456 		tau = x[end]
0 		tmp = get_rate(prob.rateCache, eltype(x))
64 		_tmp = reinterpret(eltype(x), tmp)
0 		sr = prob.pdmpFunc.R(_tmp,x,prob.xd,tau,prob.parms,true)[1]
0 		prob.pdmpFunc.F(xdot,x,prob.xd,tau,prob.parms)
0 		xdot[end] = 1.0
0 		@inbounds for i in eachindex(xdot)
0 			xdot[i] = xdot[i] / sr
- 		end
- 	# else # this is to simulate with rejection method
- 	# 	prob.pdmpFunc.F(xdot, x, prob.xd, t, prob.parms)
- 	# end
0 	nothing
- end
``````

I don’t really get the number (like why `tau=x[end]` seems to generate allocations ) but the line `_tmp = reinterpret(eltype(x), tmp)` is interesting as it seems to hint where the allocations observed above is coming from. Beyond, this, I have no clue on how to remove it.

Try putting that in a let block.

1 Like

I tried many things like

``````function (prob::PDMPProblem)(xdot, x, data, t)
tau = x[end]
let
tmp = get_rate(prob.rateCache, eltype(x))
_tmp = reinterpret(eltype(x), tmp)
sr = prob.pdmpFunc.R(_tmp,x,prob.xd,tau,prob.parms,true)[1]
prob.pdmpFunc.F(xdot,x,prob.xd,tau,prob.parms)
xdot[end] = 1.0
@inbounds for i in eachindex(xdot)
xdot[i] = xdot[i] / sr
end
end
end

``````

but I could not remove the allocation

Coming back to the issue, I may have a more urgent problem than allocations for now. I am trying to see if the reinterpret thing works well. I have an example of ODE where it fails. It seems related to `ForwardDiff.pickchunksize`.

The code I am using is borrowed from @ChrisRackauckas, you can find it here. In my use case, I have `length(u) = 81` and I get `ForwardDiff.pickchunksize(length(u)) = 12`.

Running the simulation yields:

``````ERROR: ArgumentError: cannot reinterpret an `ForwardDiff.Dual{nothing,Float64,12}` array to `ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##33#36")){PiecewiseDeterministicMarkovProcesses.PDMPProblem{Float64,Int64,Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,12},1}},SparseMatrixCSC{Int64,Int64},Synapse_params,getfield(Main, Symbol("##873#882")){Float64},getfield(Main, Symbol("##871#880")){Float64},typeof(Delta_synapse)}},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Nothing},Float64},Float64,9}` whose first dimension has size `81`.
The resulting array would have non-integral first dimension.
``````

Does anyone have an idea by any chance?

It seems the mistake is that `mod(81,12) >0`

It seems linked to Memory allocation tracking

Hi,

It seems I have fixed the allocation issue in PDMP using the cache trick from `DiffEqBase` and a profound rewrite of the library. For scalar vector fields, autodiff seems to work fine. But then I still see the error `The resulting array would have non-integral first dimension.`

A MWE would be:

``````using PiecewiseDeterministicMarkovProcesses, Random, DifferentialEquations
const PDMP = PiecewiseDeterministicMarkovProcesses

function F_tcp!(ẋ, xc, xd, parms, t)
if mod(xd[1],2)==0
ẋ[1] =  1.0
ẋ[2] = -1.0 * xc[2]
else
ẋ[1] = -1.0 * xc[1]
ẋ[2] =  1.0
end
nothing
end

R(x) = x

function R_tcp!(rate, xc, xd, parms, t, issum::Bool)
if issum == false
rate[1] = R(xc[1]) +  R(xc[2])
rate[2] = parms[1]
return 0.
else
return R(xc[1]) +  R(xc[2]) + parms[1]
end
end

xc0 = vec([0.05,0.075])
xd0 = vec([0, 1])

nu_tcp = [[1 0];[0 -1]]
parms = vec([0.1])
tf = 10000.0
nj = 1000

Random.seed!(1234)
# define a PDMP problem `a la` ODEProblem
problem = PDMP.PDMPProblem(F_tcp!, R_tcp!, nu_tcp, xc0, xd0, parms, (0.0, tf))
``````

and then, the solve part which calls autodiff

``````result4 = PDMP.solve(problem, CHV(Rodas5()))
ERROR: ArgumentError: cannot reinterpret an `ForwardDiff.Dual{nothing,Float64,2}` array to `ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){PDMPProblem{Float64,Int64,Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),typeof(R_tcp!),PiecewiseDeterministicMarkovProcesses.RateJump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},CHV{Rodas5{0,true,DefaultLinSolve,DataType}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3}` whose first dimension has size `2`.
The resulting array would have non-integral first dimension.
``````

I can isolate the issue a bit more by extracting the function to be differentiated:

``````algo = CHV(Tsit5())

# vectors to be sent to the function vf
xd1 = zeros(Float64, length(xc0)+1)
xdd1 = similar(xd1)

# for hoilding the jacobian
J = zeros(3,3)

# this is the vector field, use with an ODEProblem
vf = (dx,x) -> algo(dx, x, problem.caract, 0.)

# it works!
vf(xdd1, xd1)

ForwardDiff.jacobian!(J, vf, xdd1, xd1)
``````

which gives

``````ERROR: ArgumentError: cannot reinterpret an `ForwardDiff.Dual{nothing,Float64,2}` array to `ForwardDiff.Dual{ForwardDiff.Tag{getfield(Main, Symbol("##9#10")),Float64},Float64,3}` whose first dimension has size `2`.
The resulting array would have non-integral first dimension.
``````

I have no idea why this pops out

If anybody has an idea, I am all ears

It looks like you need a call to `reshape` in the way you’re doing the dual cache? I’d have to dig in to find out why your integration is in 3D but the cache is 2D

Yep! here

Hopefully, I dont mess up with your code. I did not call directly `DiffEqBase` to have fields more meaningful to my problem.

I noticed that you put the following in `DiffEqBase` file:

``````@require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" begin
``````

Would it solve my issue?

I think I understand the issue!!

In the above example, the vector `rate` and `xc` do not have the same dimension. Hence, when calling
`rate = get_rate(prob.ratecache, x)` in the trick of @ChrisRackauckas, it goes there:

``````get_rate(dc::DiffCache, u::AbstractArray{T}) where T<:ForwardDiff.Dual = reinterpret(T, dc.dual_rate)
``````

and if the size of `` and `x` do not match, it errors. If I change it into:

``````get_rate(dc::DiffCache, u::AbstractArray{T}) where T <: ForwardDiff.Dual = (dc.dual_rate)
``````

It seems to work. This leads me to wonder why @ChrisRackauckas had to put a `reinterpret` in the first place…

If your dual cache has the right tag that’ll work. If the tag is incorrect then it’ll fail.

it will fail silently or an error will pop out?

It will throw an error. The ForwardDiff tagging system is made to prevent perturbation confusion by throwing explicit errors if the tags don’t match. The reinterpret is essentially saying “please just change the tag for me, since this cache isn’t going to cause perturbation confusion”, essentially opting out of the tag check.

So I guess the allocation is coming from the reinterpreted array not being elided?

(In ForwardDiff2 we’ll have a more direct way of opting out of tag checks for this reason.)

There is no allocation issue anymore. The issue is that `xc` and `rate` does not have the same length. If I remove the tag `T`, it works sometimes but may fail as you notice. I don’t know how to put the correct tag in `get_rate` for it to work. Somehow, I need to change the vector length in `T` for `reinterpret` to work.

@ChrisRackauckas You are right, if I remove the tag in here it fails. If I leave it, it fails.

I modified a bit the previous example to make AD fails

``````using PiecewiseDeterministicMarkovProcesses, Random, DifferentialEquations
const PDMP = PiecewiseDeterministicMarkovProcesses

function F_tcp!(ẋ, xc, xd, parms, t)
if mod(xd[1],2)==0
ẋ[1] =  1.0
ẋ[2] = -1.0 * xc[2]
else
ẋ[1] = -1.0 * xc[1]
ẋ[2] =  1.0
end
nothing
end

R(x) = x

function R_tcp!(rate, xc, xd, parms, t, issum::Bool)
rate[1] = R(xc[1]) +  R(xc[2])
rate[2] = parms[1] * xd[1] * xc[1]
if issum == false
return 0.
else
return sum(rate)
end
end

xc0 = [0.05, 0.075]
xd0 = [0, 1]

nu_tcp = [1 0;0 -1]
parms = [0.1]
tf = 10000.0
nj = 1000
``````

Without the tag, it errors as:

``````julia> result4 = @time PDMP.solve(problem, CHV(Rodas5()))
ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3})
Closest candidates are:
Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:194
Float64(::T<:Number) where T<:Number at boot.jl:718
Float64(::Int8) at float.jl:60
...
Stacktrace:
[1] convert(::Type{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3}) at ./number.jl:7
[2] convert(::Type{ForwardDiff.Dual{nothing,Float64,2}}, ::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3}) at /Users/rveltz/.julia/packages/ForwardDiff/N0wMF/src/dual.jl:357
[3] setindex!(::Array{ForwardDiff.Dual{nothing,Float64,2},1}, ::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3}, ::Int64) at ./array.jl:766
[4] R_tcp!(::Array{ForwardDiff.Dual{nothing,Float64,2},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3},1}, ::Array{Int64,1}, ::Array{Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3}, ::Bool) at /Users/rveltz/work/prog_gd/julia/dev/PiecewiseDeterministicMarkovProcesses.jl/examples/tcp2d.jl:18
[5] (::VariableRate{typeof(R_tcp!)})(::Array{ForwardDiff.Dual{nothing,Float64,2},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3},1}, ::Array{Int64,1}, ::Array{Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3}, ::Bool) at /Users/rveltz/work/prog_gd/julia/dev/PiecewiseDeterministicMarkovProcesses.jl/src/rate.jl:13
[6] (::CHV{Rodas5{0,true,DefaultLinSolve,DataType}})(::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3},1}, ::PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}, ::Float64) at /Users/rveltz/work/prog_gd/julia/dev/PiecewiseDeterministicMarkovProcesses.jl/src/chvdiffeq.jl:16
[7] (::getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}})(::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3},1}, ::DiffEqBase.NullParameters, ::Float64) at /Users/rveltz/work/prog_gd/julia/dev/PiecewiseDeterministicMarkovProcesses.jl/src/chvdiffeq.jl:118
[8] (::ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing})(::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3},1}, ::Vararg{Any,N} where N) at /Users/rveltz/.julia/packages/DiffEqBase/V40QH/src/diffeqfunction.jl:230
[9] (::DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters})(::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3},1}) at /Users/rveltz/.julia/packages/DiffEqDiffTools/qCskj/src/function_wrappers.jl:15
``````

and with the (wrong) tag, it errors like:

``````julia> result4 = @time PDMP.solve(problem, CHV(Rodas5()))
ERROR: ArgumentError: cannot reinterpret an `ForwardDiff.Dual{nothing,Float64,2}` array to `ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3}` whose first dimension has size `2`.
The resulting array would have non-integral first dimension.

Stacktrace:
[1] (::getfield(Base, Symbol("#thrownonint#204")){ForwardDiff.Dual{nothing,Float64,2},ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3}})(::Type{ForwardDiff.Dual{nothing,Float64,2}}, ::Type{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3}}, ::Int64) at ./reinterpretarray.jl:24
[2] reinterpret at ./reinterpretarray.jl:39 [inlined]
[3] get_rate at /Users/rveltz/work/prog_gd/julia/dev/PiecewiseDeterministicMarkovProcesses.jl/src/utilsforwarddiff.jl:17 [inlined]
[4] (::CHV{Rodas5{0,true,DefaultLinSolve,DataType}})(::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(PiecewiseDeterministicMarkovProcesses, Symbol("##34#37")){CHV{Rodas5{0,true,DefaultLinSolve,DataType}},PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64},Float64,3},1}, ::PiecewiseDeterministicMarkovProcesses.PDMPCaracteristics{typeof(F_tcp!),VariableRate{typeof(R_tcp!)},PiecewiseDeterministicMarkovProcesses.Jump{Int64,Array{Int64,2},typeof(PiecewiseDeterministicMarkovProcesses.Delta_dummy)},Array{Float64,1},Array{Int64,1},PiecewiseDeterministicMarkovProcesses.DiffCache{Array{Float64,1},Array{ForwardDiff.Dual{nothing,Float64,2},1}},Array{Float64,1}}, ::Float64) at /Users/rveltz/work/prog_gd/julia/dev/PiecewiseDeterministicMarkovProcesses.jl/src/chvdiffeq.jl:15
[5] #34 at /Users/rveltz/work/prog_gd/julia/dev/PiecewiseDeterministicMarkovProcesses.jl/src/chvdiffeq.jl:118 [inlined]
[6] ODEFunction at /Users/rveltz/.julia/packages/DiffEqBase/V40QH/src/diffeqfunction.jl:230 [inlined]
[7] UJacobianWrapper at /Users/rveltz/.julia/packages/DiffEqDiffTools/qCskj/src/function_wrappers.jl:15 [inlined]
``````

I thought I had solved this issue but it seems I did not. If anyone knows how to make this work with `ForwardDiff`, it would be very helpful to me.

Thank you

maybe @kristoffer.carlsson has an idea or a hint… I am a bit stuck…