I didn't ask for plan_bfft, but got MethodError: no method matching plan_bfft(::Vector{Number}, ::UnitRange{Int64})

Specific Domains Signal and Image Processing
1

While coding an iterative FFT function, I came upon a strange bug that depends on the arguments to the function in a complicated way. I’ve managed to reduce the problematic code to the following:

using FFTW
unitize(x::Complex) = x==0 ? 0 : x/abs(x)
function problematic(ys::Tuple,phis::Tuple,nit)
    n = length(ys)
    x = sum(ifft(ys[i]) for i=1:n)/n
    for j=1:nit
        thetas = ((unitize.(fft(x.*phi)) for phi in phis)...,)
        x = sum(ifft(ys[k].*thetas[k]) for k=1:n)/n
    end
    return x
end
x = [exp(-j^2) for j=-10:.1:10]
quad = exp.(im*[-10:.1:10;].^2)
y0 = abs.(fft(x))
y1 = abs.(fft(x.*quad))
x1 = problematic((y0,y1),(ones(length(x)),quad),4);              # This works
x1 = problematic((y0,y1),(ones(length(x)),quad),5);              # But this does not work

The various parameters and the auxiliary function unitize are, unfortunately, necessary for this demonstration, as slight tweaks to the function eliminates the error. The final parameter to problematic determines the number of iterations, and changing it slightly affects whether or not an error is thrown. With slightly different inputs y0, y1, quad, etc. the exact number of iterations at which the error shows up will change also. Under some conditions, I also find that the number of iterations required for an error differs between the REPL and a Jupyter notebook.

At any rate, sometimes this code throws a MethodError: no method matching plan_bfft(::Vector{Number}, ::UnitRange{Int64}). I find this strange, because I didn’t ask for plan_bfft at all. Is the compiler automatically planning an FFT when it sees many FFTs in a loop? Did I do something silly, or is this a bug?

I’m using Julia v1.7.2 and FFTW v1.4.6.

2

I’m not certain about why this only fails for nit = 5 and not 4, but it might be something about inference giving up beyond a certain tuple length (just a guess.)

However, I think the problem originates from this type instability:

Try for example:

unitize(x::Complex) = x==0 ? x/one(x) : x/abs(x)

and see if that helps.

1 Like
3

This does indeed seem to fix the problem (thanks!), as does this modification:
unitize(x::Complex) = x==0 ? zero(x) : x/abs(x)
which I believe also resolves the type instability–correct? I don’t understand yet why this is important. Are you able to give me some intuition as to what significance this change has? Is there a reason you used the expression x/one(x) instead of zero(x)?

4

It’s ok for Complex{<:AbstractFloat}:

julia> @code_warntype unitize(1.0 + 2.0im)
...
Body::ComplexF64

but not for Complex{<:Integer} (or Complex{<:Rational} etc.)

julia> @code_warntype unitize(1 + 2im)
...
Body::Union{Complex{Int64}, ComplexF64}

Because:

julia> unitize(1 + 2im)
0.4472135954999579 + 0.8944271909999159im

julia> unitize(0 + 0im)
0 + 0im

Just in case, I would also use iszero(x) instead of x == 0, though it’s not a problem for any cases I can think of off the top of my head.

1 Like