MWE from docs is not working because of undefined `plan_nfft`

[roflmaostc]

I'm using the recent master branch and the code from the docs returns:

(tmp) pkg> st
      Status `/tmp/Project.toml`
  [052768ef] CUDA v2.6.3
  [7a1cc6ca] FFTW v1.4.1
  [efe261a4] NFFT v0.6.1 `https://github.com/tknopp/NFFT.jl#master`
  [e88e6eb3] Zygote v0.6.10 `~/.julia/dev/Zygote`

julia> p = plan_nfft(x, N)                  # create plan. m and sigma are optional parameters
ERROR: UndefVarError: plan_nfft not defined
Stacktrace:
 [1] top-level scope
   @ REPL[32]:1

[roflmaostc]

Hm, not sure what was wrong, but it works after REPL restart

[tknopp]

did you using NFFT?

[tknopp]

if some Pkg operation has errored (or precompilation) then you would need to restart the REPL.

[roflmaostc]

Yeah, I did, probably Revise failed. But thanks for the quick reply. Full log was:

(tmp) pkg> add https://github.com/tknopp/NFFT.jl
     Cloning git-repo `https://github.com/tknopp/NFFT.jl`
    Updating git-repo `https://github.com/tknopp/NFFT.jl`
   Resolving package versions...
   Installed ColorTypes ─ v0.11.0
    Updating `/tmp/Project.toml`
  [efe261a4] ~ NFFT v0.5.2 ⇒ v0.6.1 `https://github.com/tknopp/NFFT.jl#master`
    Updating `/tmp/Manifest.toml`
  [3da002f7] + ColorTypes v0.11.0
  [5ae59095] + Colors v0.12.8
  [53c48c17] + FixedPointNumbers v0.8.4
  [a2bd30eb] + Graphics v1.1.0
  [efe261a4] ~ NFFT v0.5.2 ⇒ v0.6.1 `https://github.com/tknopp/NFFT.jl#master`
Precompiling project...
  4 dependencies successfully precompiled in 11 seconds (44 already precompiled)

(tmp) pkg> 

julia> using NFFT

julia> M, N = 1024, 512
(1024, 512)

julia> x = range(-0.4, stop=0.4, length=M)  # nodes at which the NFFT is evaluated
-0.4:0.0007820136852394917:0.4

julia> fHat = randn(ComplexF64,M)     
1024-element Vector{ComplexF64}:
 -0.40073711423571046 + 0.8080661714789623im
 [...]

julia> p = plan_nfft(x, N)                  # create plan. m and sigma are optional parameters
ERROR: UndefVarError: plan_nfft not defined
Stacktrace:
 [1] top-level scope
   @ REPL[24]:1

[tknopp]

could be. Revise would fail if we changed the struct. But anyway: Import is that it now works! Hope you find NFFT.jl useful. I have worked some time in the field of OCT (Thorlabs) and know that there one also requires an NFFT :-)

[roflmaostc]

I'm encountering the issue that I either a need a fully CUDA+Zygote compatible irregular cubic spline interpolation or a CUDA+Zygote compatible NFFT.

But probably NFFT doesn't work with Zygote efficiently, does it?

[tknopp]

Could you elaborate what you mean with compatible? Is this about package versions? Or are you talking about fancy Automatic Differentiation stuff? Unfortunately I have no experience with Zygote.

[roflmaostc]

My work is related to inverse modeling like deconvolution. The forward model usually contains some FFT operations and to invert the forward model we want to optimize a loss function with respect to the reconstruction (the sample). Zygote helps us figuring out the gradients so that we can put that in an optimizer.

Zygote has a adjoint rule for the FFT defined (it's an IFFT) so that we can efficiently use reverse automatic differentiation. For the NFFT we probably need to define an adjoint rule (maybe its something like a INFFT :D?) Usually Zygote cannot differentiate through naive Julia code, so one really has to put emphasis on that to get it working.

[roflmaostc]

You already provide an adjoint method, maybe that's already enough. Honestly, I first need to understand more in depth what the full consequences of a NFFT are.

[tknopp]

yes adjoint is no problem. But probably you are missing that we overload A*x and adjoint(A)*x. If you need that, you can use NFFTOp from https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/master/src/Operators/NFFTOp.jl

right now its in MRIReco.jl but we could move that outside.

[roflmaostc]

I'm pretty confused how to understand the output actually From the 2D example: The shapes of adjoint and output of the nfft are somehow different

But I probably need to read the NFFT documentation to know what's going on behind that.

As you might know, the shapes of fft and ifft are the same, and hence also the adjoints do have same shape and sizes.

julia> f = nfft_adjoint(p, fHat)
16×16 Matrix{ComplexF64}:
   -49.362+24.6656im    3.85684+16.3056im  -14.1196+10.9864im   14.9051-34.903im   …    15.6725+19.5867im    -1.9155-6.42269im    -8.12476-17.0358im   17.2532-3.94638im
  -20.7768-14.2628im   -8.47732-8.63997im    22.433-4.52048im  -22.7806-1.39246im      -52.7173+11.3855im   -9.46708+27.7343im     -5.9962+30.8146im  -18.1814-25.159im
 -0.390413-27.0455im   -22.2497-15.5196im  -8.70234+8.2614im   -10.3551-19.6801im       15.6014+32.534im    -14.8827+20.7465im     3.84553-13.4409im    4.6665-15.1238im
   18.4431-3.49755im    10.7214+25.2727im  -14.0833+23.6791im  -41.8363-4.5671im        19.3307-60.428im     24.3737+40.8844im      2.3045-32.4964im   28.8032-13.157im
   6.79058-37.8686im   -9.12399+4.57613im   24.6994-21.0399im   38.3865-6.82765im        -16.56+6.7572im    -34.7619+21.971im     -35.8184+15.3437im    36.996+20.0195im
   14.4025+28.0988im    15.1117+12.4556im  -6.82204+17.2285im   3.29595+3.27664im  …  -0.460226+8.82571im   -21.0588-0.2144im     -56.5957-29.4727im  -21.5305-16.6351im
   8.79032-19.3376im    -16.601+18.2259im  -3.09162+45.1923im   4.57518-31.635im        3.79025-29.2063im   -8.78163-0.886927im  -0.181176-14.1378im   3.82555+6.65829im
   54.9912-6.47767im   -17.2924-7.2874im    2.64169+24.497im    7.13744-16.3165im        6.5193+26.2084im   -23.4799+6.33826im    -29.6393-15.6063im   19.6957+22.1706im
  -17.7861-73.9137im   -8.51353-56.7743im   9.49438+4.38144im  -10.7461-39.2762im       17.8828-3.78593im   -36.5708-1.51967im     37.2458+19.6917im   11.7919+7.52873im
 -0.413541-21.4387im     58.409-32.1434im   1.21211+14.1222im   8.76805+4.92245im       21.6865+0.172431im   16.6209+22.8114im    -11.6854-18.5881im   24.0139+8.06563im
   36.3227+5.33871im    7.33636-27.6286im   -3.1544-20.1562im   12.4236+14.1936im  …    78.6123-20.3949im    1.48099-20.1871im     25.2775-30.6724im  -3.32101+19.8425im
  -1.15296+58.8136im    14.5882+16.8595im   27.1711-28.6956im   9.12951+11.0161im       4.76799+33.4249im   -7.87194-18.2194im      1.2279-37.7843im  -7.85924-27.4377im
   3.45147+14.9375im   -22.9858+12.7173im   29.2509-13.9024im   29.5163-11.0595im       24.8401+1.09892im   -11.4318-8.84855im    -13.6129-20.4596im    27.831+15.4057im
  -21.3067+14.5054im   -10.4823-8.62603im  -12.7915+14.8526im   7.84331+18.9571im       21.6462+4.10781im   -2.88284+9.38693im     15.5868+15.8525im   -8.0779+35.0505im
  -21.7436+7.64438im    31.9306-12.3024im   -20.039-25.0644im   -4.6185-16.98im        -18.4446-7.26463im    7.02414-16.3395im    -48.5688+35.5891im   16.5613-6.0935im
   27.4556+0.627738im   -1.6565-10.1982im   20.8955-11.961im    16.4241-10.7803im  …   -3.63112-16.0762im    11.3541-7.7232im      3.70002-6.54362im   5.49114+38.4053im

julia> g = nfft(p, f)
1024-element Vector{ComplexF64}:
    81.9971800350672 + 320.68810772246013im
   627.2095068668118 + 27.52098734455273im
  -53.99511315282497 - 97.19165017625885im
[...]

[tknopp]

yes the NFFT matrix is rectangular (MxN) and the N dimension is actually something like N_1 x N_2 x ... For the M dimension we cannot do this since it is irregular. If in your case N=M this is a special case and you can just reshape.

[roflmaostc]

Just to clarify again: NFFT can handle data sampled on a 2D irregular grid? So basically a vector of positions and a vector of data?

Thanks for your help so far! I'll definitely keep that in mind and to have CUDA support sounds already nice.

[tknopp]

Just to clarify again: NFFT can handle data sampled on a 2D irregular grid? So basically a vector of positions and a vector of data?

Exactly. It is ND where N can be 2. And as you say you need to provide both, the nodes and the actual values at the positions.

[kiranshila]

@roflmaostc Did you ever get this working with Zygote, or any other AD engine? I'm also working on inverse modeling and would really like to use Flux with this. I opened an issue about it here. #50

[roflmaostc]

No I haven't

Usually the adjoint definition of an fft is the ifft with some scaling. But I'm confused about the "adjoint" given here and whether it has really the same meaning we are interested in.

[kiranshila]

Yeah I see what you mean. If I'm taking the NFFT of an image, I get a vector of fourier components at each spatial frequency of interest. If the adjoint is the ifft, what does that mean? Especially if in this case those are two different sizes.

[tknopp]

I am not sure what you are struggling with. The NFFT is just a matrix vector multiplication. The adjoint is the conjugated transposed of the matrix. In case of a rectangular matrix, it is clear the the dimensions need to switch.

[JeffFessler]

i think many people wish there was an "inverse NFFT" but of course there is not one in general.

[roflmaostc]

You say "in general". What are the conditions to have one?

[JeffFessler]

well, first of all only a square matrix is invertible, so a necessary condition is that M=N. after that, the only case where i know how to invert analytically is when the samples lie on the cartesian grid and then you are back to the DFT special case. i am sure there are other cases with M=N that are invertible, but i'll bet that no fast algorithms exist for their inverses in general.

[tknopp]

Exactly. One can approximate the inverse NFFT using an adjoint NFFT but one needs to account for the sampling density, which can be nicely seen when considering the continuous Fourier transformation pair and applying a (rectangular) quadrature rule. When doing this it also becomes apparent that it is an approximation and no analytical solution.