JuliaMath / FFTW.jl

Julia bindings to the FFTW library for fast Fourier transforms
https://juliamath.github.io/FFTW.jl/stable
MIT License
269 stars 54 forks source link

both FFTW and Base export "FFTW" #47

Open EthanAnderes opened 6 years ago

EthanAnderes commented 6 years ago

Been trying to get some code running on v0.7 and running into difficulties with getting FFTW to work in my modules. Here is a simplified module that gives me errors.

module Fields

export r𝔽, Flat

using FFTW
import Base: *, \

abstract type Flat{Θpix,nside} end

struct r𝔽{P<:Flat,T<:Real,F}
    Δx::T
    Δk::T
    Ωk::T
    Ωx::T
    period::T
    nyq::T
    k::Vector{Matrix{T}}
    x::Vector{Matrix{T}}
    FFT::F
end

# real FFT generated function constructor
@generated function r𝔽(::Type{P},::Type{T}) where T<:Real where P<:Flat{Θpix, nside}  where {Θpix, nside}
    Δx     = deg2rad(Θpix/60)
    period = Δx*nside
    Δk     = 2π/period
    Ωk     = Δk^2
    Ωx   = Δx^2
    nyq    = 2π / (2Δx)
    k_side = ifftshift(-nside÷2:(nside-1)÷2) * Δk
    x_side = ifftshift(-nside÷2:(nside-1)÷2) * Δx
    k      = [reshape(k_side, 1, nside), reshape(k_side[1:nside÷2+1], nside÷2+1, 1)]
    x      = [reshape(x_side, 1, nside), reshape(x_side, nside, 1)]
    FFT    =  Ωx / (2π) * plan_rfft(rand(T,nside,nside); flags=FFTW.PATIENT, timelimit=4)
    r𝔽{P,T,typeof(FFT)}(Δx, Δk, Ωk, Ωx, period, nyq, k, x, FFT)
end

(*)(g::r𝔽{P,T}, x) where P<:Flat where T = g.FFT * x
(\)(g::r𝔽{P,T}, x) where P<:Flat where T = g.FFT \ x

end

I want the following code block to run ...

using Main.Fields

nside  = 1024
Θpix   = 2.0
Px     = Flat{Θpix,nside}
Tx     = Float32
g      =  r𝔽(Px,Tx)

fx = rand(nside, nside)
fk = g * fx
y = g \ fk

.... but I get the following error in the line:

julia> g      =  r𝔽(Px,Tx)
WARNING: both FFTW and Base export "FFTW"; uses of it in module Fields must be qualified
ERROR: UndefVarError: FFTW not defined
Stacktrace:
 [1] #s1#1(::Any, ::Any, ::Any, ::Any, ::Any, ::Type, ::Any) at ./REPL[1]:34
 [2] (::Core.GeneratedFunctionStub)(::Any, ::Vararg{Any,N} where N) at ./boot.jl:443
 [3] In toplevel scope

Any ideas what is going on?

Here is my version info:

julia> versioninfo()
Julia Version 0.7.0-DEV.2454
Commit 00dafe28c0* (2017-11-07 15:53 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin17.2.0)
  CPU: Intel(R) Core(TM) i7-5557U CPU @ 3.10GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, broadwell)
stevengj commented 6 years ago

Hmm, I can reproduce your error, but this works fine:

julia> module Foo
       using FFTW
       end
Main.Foo

julia> Foo.fft([3,4,5])
3-element Array{Complex{Float64},1}:
 12.0 + 0.0im               
 -1.5 + 0.8660254037844386im
 -1.5 - 0.8660254037844386im

so it may be a problem with JuliaLang/julia#22763. Can we boil it down to a simpler test case?

(Not sure why you need @generated here.)

stevengj commented 6 years ago

Okay, here is a simpler test case:

julia> module Foo
       using FFTW
       foo() = plan_fft([0.0+0.0im], flags=FFTW.PATIENT)
       end
Main.Foo

julia> Foo.foo()
WARNING: both FFTW and Base export "FFTW"; uses of it in module Foo must be qualified
ERROR: UndefVarError: FFTW not defined
Stacktrace:
 [1] foo() at ./REPL[1]:3

And even simpler:

julia> module Foo
       using FFTW
       const foo = FFTW.PATIENT
       end
WARNING: both FFTW and Base export "FFTW"; uses of it in module Foo must be qualified
ERROR: UndefVarError: FFTW not defined
EthanAnderes commented 6 years ago

Thanks for looking into it so quick @stevengj.

In response to your question "Not sure why you need at-generated here" ....

Basically, I'm working with multiple fields that can be defined on different grids (encoded by Flat{Θpix,nside}). The beauty of the generated version is that I can write methods which don't require me to keep track and pass planned fft's to my methods. E.g. most of my code has the form

function foo(f::Field{P,T}) where {P<:Flat, T<:Real}
  g = r𝔽(P,T)
  ... then do stuff with g and f, etc .
end

The first time foo gets called on a field with a new grid geometry P<:Flat, a new plan is generated. Subsequent calls to g = r𝔽(P,T) are instantious. It also gives me access to grid geometry stuff in the body of the methods that don't need to be passed or pre-computed.

Also, I believe (though untested as of yet) that this works well when parallelizing across multiple machines (all of which should be generating their own plan for the architecture, I think).

It's really amazing having this generated stuff!

stevengj commented 6 years ago

I see. Note that you can just use Array{T}(nside,nside) instead of rand(T,nside,nside), since the FFTW.PATIENT planner re-initializes the array anyway.

ararslan commented 6 years ago

I believe the original workaround for this was to use importall FFTW, but importall is now deprecated in 0.7.