PreallocationTools.jl icon indicating copy to clipboard operation
PreallocationTools.jl copied to clipboard

Complex + dualcache

Open moble opened this issue 2 years ago • 2 comments

Nowadays ForwardDiff seems to work nicely when using Complex (internally, at least), so I was trying to use dualcache with a complex array. Unfortunately, it seems like the more interesting methods of get_tmp base the returned eltype on the type of u (which is always some kind of Real) rather than on the type underlying the DiffCache. I wonder if it might be possible to generalize that.

It is possible to work around this by constructing the cache as zd = dualcache(reinterpret(reshape, T, z)) and then unpacking it inside the function as

    ztmp = get_tmp(zd, t)
    z = reinterpret(reshape, Complex{eltype(z)}, ztmp)

Obviously, this is ugly, as well as hard to remember to do. I've fallen into this trap a few times, getting weird error messages far from the source of the problem. Automatic handling of Complex would really make things nicer for me. :)

moble avatar Nov 01 '22 20:11 moble

Sorry, I am not usually dealing with models where complex numbers appear (unless I am doing something wrong :)); can you provide a MWE that illustrates the issue and what you are trying to accomplish?

thomvet avatar Nov 03 '22 13:11 thomvet

Sorry. Here's an MWE. The uncommented code is what I expected to be able to do; the commented code is what I have to do to get it to work.

using PreallocationTools
using ForwardDiff

z = zeros(ComplexF64, 20)
zd = dualcache(z)
#zd = dualcache(reinterpret(reshape, real(eltype(z)), z))

function sum_cis(θ)
    z = get_tmp(zd, θ)
    #ztmp = get_tmp(zd, θ)
    #z = reinterpret(reshape, Complex{eltype(ztmp)}, ztmp)
    for i ∈ eachindex(z)
        z[i] = cis(i*θ)
    end
    abs(sum(z))
end

ForwardDiff.derivative(sum_cis, 1.1)

moble avatar Nov 03 '22 14:11 moble