Hi, have you found a solution to this problem? I’m stuck with the exact same issue.
No.
To me, the easiest way to see the issue without the PDMP burden is to look at my post above where vf is exposed
Thanks, actually after investigation my issue was slightly different and was related to the computation of the Hessian in the optim.jl package.
It produces a Dual{Tag, Dual{Tag, T}} type. I solved it using the modified code below:
struct DiffCache2{T<:AbstractArray, S1<:AbstractArray, S2<:AbstractArray}
du::T
dual_du::S1
dual_dual_du::S2
end
function DiffCache2(u::AbstractArray{T}, siz, ::Type{Val{chunk_size}}) where {T, chunk_size}
DiffCache2(u, zeros(Dual{nothing,T,chunk_size}, siz...), zeros(Dual{nothing,Dual{nothing,T,chunk_size},chunk_size}, siz...))
end
dualcache2(u::AbstractArray, N=Val{pickchunksize(length(u))}) = DiffCache2(u, size(u), N)
get_tmp2(dc::DiffCache2, u::AbstractArray{T1}) where T1<:Dual{S,T2} where S<:Tag where T2<:Dual = reinterpret(T1, dc.dual_dual_du)
get_tmp2(dc::DiffCache2, u::AbstractArray{T}) where T<:Dual = reinterpret(T, dc.dual_du)
get_tmp2(dc::DiffCache2, u::AbstractArray) = dc.du
What was your issue and how did you debug it?
I was trying to use dualcache for an optimization problem involving large arrays (trying to replicate this), and I kept getting reinterpret errors when calling get_tmp.
It turns out that because I was using optim.jl with a Newton algorithm, it needs to calculate the Hessian (second derivative) of the objective function, and instead of an Array of Dual numbers, I got an Array of Dual{Dual}, which couldn’t be reinterpreted in the get_tmp function because dual_du was not large enough.
The code mentioned above solves the issue by having another cache array dual_dual_du which is of Dual{Dual{…}} type.
Sorry to bring up this old thread, but when I ran into the The resulting array would have non-integral first dimension. error thrown by reinterpret in the context above, this thread popped up and it was disappointing that there was no solution to it (or at least not reported). So, I thought it makes sense to post the solution here and spread the word.
TL;DR version:
After some tinkering together with @ChrisRackauckas , stuff like this is now fixed in PreallocationTools.jl. The library now allows flexible changing of chunk sizes (which is the underlying source of the error, see below) for cache arrays that contain dual numbers.
By solving that problem, it also became possible to define dual number caches for multiple nested AD levels (essentially a more flexible version of what @touste showed above). While that’s now possible, it’s (often) not advisable in larger scale problems (for example calculation of “large” Hessians with forward over forward AD). So, users should still do appropriate analysis and performance testing on the specific problem at hand.
Some background to the error messages and what the way to figuring this out was:
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.
gets thrown because we ask to change the type-interpretation of a block of memory, but the size (in this context: number of bytes) of the original type, the size of the final type, the number of elements in the array to be re-interpreted, and the output dimensions of the array are not consistent.
In the context of automatic differentiation with ForwardDiff.jl, which uses chunked dual numbers, this occurs when defining a cache with one chunk size and then using it with another.
To make a specific example:
#define Dual numbers with different chunk sizes
origintype = ForwardDiff.Dual{nothing,Float64,2} #nothing is the tag here
finaltype = ForwardDiff.Dual{typeof(something),Float64,3} #now with a tag and different chunk size
finaltype2 = ForwardDiff.Dual{typeof(something),Float64,1}
#check sizes of these types
sizeof(origintype) #24; chunk size two means one value and two partials, hence three Float64 (each is 8 bytes)
sizeof(finaltype) #32
sizeof(finaltype2) #16
#define a matrix of dual numbers to be used as cache.
cache_dual = zeros(origintype, 2,2)
#use re-interpret to change type (tag and chunk size)
R1 = reinterpret(finaltype, A) #non-integer dimension error
R2 = reinterpret(finaltype2, A) #results in a 3 x 2 array (!) of Duals with chunksize 1 and the correct tag.
So, in the case of R1 we attempt to reinterpret a block of memory of 2 x 2 x 24 bytes = 96 bytes into a new array with an eltype size of 32 bytes. We ask for a Y x 2 array and since 96/32 = 3, we obviously get the non-integer dimension error.
In R2, we reinterpret 96 bytes into elements of 16 bytes and get 6 elements, rearranged into a 3 x 2 array. Since we really wanted to have a 2 x 2 array and the rest of our code probably relied on that fact (matrix multiplication or the like), that outcome is equally bad and will likely throw a dimension mismatch error somewhere else.
To resolve this problem, PreallocationTools.jl is now using a mix of views into the cache array and array reinterpretation (for example, four elements in a float64 array are reinterpreted into one dual number with three partials). A little hacky, but it works and it was quite straightforward to extend this to nested dual numbers. The issue is not anymore to type them correctly, but just to provide a large enough cache array that gets reinterpreted into an array of the right eltype and the desired size.