Issue with PDMP and Forwardiff - DifferentialEquation

rveltz · November 23, 2018, 8:17pm

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.

thank you for your help,

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

rveltz · June 27, 2019, 8:11pm

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

ChrisRackauckas · June 29, 2019, 5:39am

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.

rveltz · June 29, 2019, 7:49am

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)

rveltz · July 1, 2019, 8:59am

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.

ChrisRackauckas · July 1, 2019, 11:50am

Try putting that in a let block.

rveltz · July 1, 2019, 4:43pm

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

rveltz · July 27, 2019, 7:59am

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

rveltz · July 30, 2019, 12:21pm

It seems linked to Memory allocation tracking

rveltz · September 27, 2019, 9:22am

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

Thank you for your help

ChrisRackauckas · September 27, 2019, 11:24am

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

rveltz · September 27, 2019, 11:27am

Yep! here

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

rveltz · October 16, 2019, 1:50pm

I noticed that you put the following in DiffEqBase file:

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

Would it solve my issue?

rveltz · October 17, 2019, 8:42am

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…

ChrisRackauckas · October 17, 2019, 10:44am

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

rveltz · October 17, 2019, 11:59am

it will fail silently or an error will pop out?

ChrisRackauckas · October 17, 2019, 12:05pm

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

rveltz · October 17, 2019, 3:14pm

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.

rveltz · October 17, 2019, 3:37pm

@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

rveltz · October 23, 2019, 7:58pm

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

Topic		Replies	Views
How to cache dual numbers correctly /dual_cache reinterpretation error General Usage question , diffeq	1	840	November 10, 2019
ForwardDiff caches usage General Usage forwarddiff	5	998	July 5, 2021
Problem with ForwardDiff and pre-allocating vector Numerics	6	1422	December 14, 2017
How to make this function compatible with ForwardDiff? General Usage forwarddiff , autodiff	15	1474	September 22, 2021
ForwardDiff problem: Cannot `convert` an object of type ForwardDiff.Dual{ General Usage question	2	1258	May 27, 2018

Issue with PDMP and Forwardiff - DifferentialEquation

Related topics