Implement AbstractNFFTs

[tknopp]

Hi,

over at NFFT.jl I have been working on an AbstractNFFTs.jl interface package that allows to use different NFFT implementations. The design is similar to that of AbstractFFTs.jl, which also allows for different implementations of the ordinary FFT.

Some more background can be found here:

JuliaMath/NFFT.jl#59

The question is if

• FINUFFT.jl would be interested in implementing AbstractNFFTs.jl
• If there are any missing things in the interface package that would make it difficult to wrap FINUFFT.jl.

Currently we have one CPU (NFFT.jl) and one GPU (CuNFFT.jl) implementation. Further I have a wrapper for NFFT3.jl (https://github.com/JuliaMath/NFFT.jl/blob/master/NFFT3/NFFT3.jl) that will hopefully be incorporated in NFFT3.jl soon.

The wrapper would not need to be your primary interface. It's just an additional interface being compatible with the common interface.

Best,

Tobias

[tknopp]

(Side comment: I just read through the papers of FINUFFT and this is awesome work that you have done! This makes me wanting the wrapper for myself so that we can use it over in MRIReco.jl )

And one more direct comment to the interface: While AbstractNFFTs.jl does not require a plan to have specific parameters (like the oversampling factor and the kernel size), it is actually nice of different implementation share also the same parameters since it allows for exchanging the plan. Does FINUFFT.jl provides a way to set the oversampling factor and the kernel size or can you just pass the desired accuracy? In fact, motivated by your paper I would like to get an interface into AbstractNFFTs.jl where one just passes the desired accuracy and the oversampling factor and the kernel size are chosen automatically.

[ahbarnett]

Dear Tobias,
Thanks for reaching out. Yes, we do wipe the floor with other known CPU implementations :)
(But since you claim to be faster than NFFT3 which we beat, it would be good to benchmark FINUFFT vs NFFT.jl)
We would love to have FINUFFT.jl wrapped in this way so that more Julia users can access it.
We also have cuFINUFFT which does not yet have a jl wrapper, but that would be nice. (On the Py side, people are starting to wrap our CPU and GPU codes into Tensorflow, etc.)

Our devs on the Julia side is @ludvigak @lu1and10 and myself, and I don't know how much time people have right now, but it seems like a short self-contained project.

Comments:

0. A good/clean application-agnostic interface design is crucial. (Ie, one that is not just good for MRI community.) The NFFT is used in many disjoint areas of science, just like the FFT. It would be great if you reach out eg to those users cited in the FINUFFT docs, in case they want to chime in on your abstract interface. It will be more work, but much better in the end.

1. As you say, I encourage you to make your abstract interface be as method-independent as possible. This means that there should be no required references to kernel parameters (ie, width, kernel type, etc), since those are "under the hood". The basic user should simply set only the transform sizes and desired accuracy. Even the upsampling should be an option - FINUFFT does let you set upsampling as an option (to two values: 2 and 5/4). You will notice that our default is to auto-choose upsampling (eg type 3 prefer 5/4 but most type 1,2 prefer upsampling of 2), but having it as an option lets the user experiment.

2. I think it's crucial to remove all MRI-specific notation from AbstractNFFTs, and stick to the mathematical definitions only.
   I like the fact you have proper math definitions in the docs for NFFT.jl. But currently you refer in
   https://github.com/JuliaMath/NFFT.jl/blob/master/AbstractNFFTs/src/AbstractNFFTs.jl
   to "apodization" and "convolve" which few outside MRI would understand. One of our goals with FINUFFT (and NFFT3 does this too) was to make it application-agnostic, as a pure math library should be.

3. Multiple transforms with the same NU pts. I see your "directional NFFT" definition in your docs.
   This is similar to our "vectorized" (multiple data) except we always transform along the innermost dims.
   We don't allow arbitrary "middle" dims to be transformed. Is there a reason you need this?
   (I don't know of any NFFT implementations that do this, eg does FFTW let you do this in the FFT case?) Your interface could be realised by us by sandwiching FINUFFT between two multi-dim transposes, but you'd probably have to help us.

4. Type 3. This comes up in applications, and shaped our interface. AbstractNFFT should be designed to allow it in future.

5. I'm glad you have proper mathematical definitions! But there are lots of math typos (bold, etc) - I will file an Issue there.

6. in AbstractNFFT you may want to think about something like our plan vs setpts interface. This allows NU points
   to move without re-planning the FFT. This was influenced by NFFT3's design... Ah, I see you have nodes! which does this, good.

7. I don't see the plan p object defined in the AbstractNFFT. This seems crucial: does it contain the tolerance, and NU pts arrays?

Best wishes, Alex

[tknopp]

Thanks so much for you positive feedback! In the meantime I was greedy and gave it a try:
https://github.com/JuliaMath/NFFT.jl/blob/master/Wrappers/FINUFFT.jl
And it works! Accuracy tests pass. I think at this point its actually quite nice to have the wrapper evolve within NFFT.jl and if we are satisfied we can integrate it into NFFT3.jl and FINUFFT.jl

Thanks for reaching out. Yes, we do wipe the floor with other known CPU implementations :)
(But since you claim to be faster than NFFT3 which we beat, it would be good to benchmark FINUFFT vs NFFT.jl)

we will see ... :-) I have the benchmarks already running but at this point they are not yet serious since my wrapper is certainly not ideal.

We also have cuFINUFFT which does not yet have a jl wrapper, but that would be nice. (On the Py side, people are starting to wrap our CPU and GPU codes into Tensorflow, etc.)

Actually I am very(!) interested in cuFINUFFT. So please, please make this available as a Julia package. We have a very preliminary CuNFFT but that is more a toy thing. Getting cuFINUFFT would thus be very very much appreciated.

Our devs on the Julia side is @ludvigak @lu1and10 and myself, and I don't know how much time people have right now, but it seems like a short self-contained project.

I am not completely in a hurry but would like to streamline the interface of AbstractNFFT.jl now. Still waiting on feedback of NFFT3 folks.

A good/clean application-agnostic interface design is crucial. (Ie, one that is not just good for MRI community.) The NFFT is used in many disjoint areas of science, just like the FFT. It would be great if you reach out eg to those users cited in the FINUFFT docs, in case they want to chime in on your abstract interface. It will be more work, but much better in the end.

Totally agree. I know the origins of NFFT and various applications. I implemented the NNFFT (NuFFT type 3) in NFFT3 when I studied computer science. But actually NFFT.jl should already be completely independent of applications.

As you say, I encourage you to make your abstract interface be as method-independent as possible. This means that there should be no required references to kernel parameters (ie, width, kernel type, etc), since those are "under the hood". The basic user should simply set only the transform sizes and desired accuracy. Even the upsampling should be an option - FINUFFT does let you set upsampling as an option (to two values: 2 and 5/4). You will notice that our default is to auto-choose upsampling (eg type 3 prefer 5/4 but most type 1,2 prefer upsampling of 2), but having it as an option lets the user experiment.

We will have "standard" parameters but it will always be possible to tweak a plan with implementation specific parts. Regarding the parameters my proposal right now is:

• default is the parameter reltol (that keyword is a standard Julia name for that) and if that is chosen, the other parts are derived automatically (no work on your side, since you do this already)
• if the user instead specifies kernel size and sigma that parameter is chosen. For FINUFFT this means I need to reverse engineer the window size, which should be doable with formula 3.2 of your paper.

I think it's crucial to remove all MRI-specific notation from AbstractNFFTs, and stick to the mathematical definitions only.
I like the fact you have proper math definitions in the docs for NFFT.jl. But currently you refer in
https://github.com/JuliaMath/NFFT.jl/blob/master/AbstractNFFTs/src/AbstractNFFTs.jl
to "apodization" and "convolve" which few outside MRI would understand. One of our goals with FINUFFT (and NFFT3 does this too) was to make it application-agnostic, as a pure math library should be.

I agree, no MRI speech! But still, these two names make also sense when looking at the NFFT purely from the signal processing / mathematical perspective. For the term "apodization" I am not 100% sure but the term "convolution" definitely fits. How would you call them otherwise? When teaching NFFT I always give the intuitive explanation via the convolution theorem before I make the formal derivation, which is more technical.

Multiple transforms with the same NU pts. I see your "directional NFFT" definition in your docs.
This is similar to our "vectorized" (multiple data) except we always transform along the innermost dims.
We don't allow arbitrary "middle" dims to be transformed. Is there a reason you need this?
(I don't know of any NFFT implementations that do this, eg does FFTW let you do this in the FFT case?) Your interface could be realised by us by sandwiching FINUFFT between two multi-dim transposes, but you'd probably have to help us.

:-) why we did this: because we can! But don't worry about this since it is certainly pretty exotic. We will simply implement this part only partially with FINUFFT. Actually NFFT3 does not provide this at all and thus is more restricted. But as my goal for the abstract interface is to be very general I came up with this solution. By the way: FINUFFT seems also not be capable of doing 4D transforms, which nobody is doing anyway. In that case we need also to raise an error but that is fine. So we will always have feature that one backend will not have.

Type 3. This comes up in applications, and shaped our interface. AbstractNFFT should be designed to allow it in future.

This is probably the most important point I need to think about. Right now my idea would be to have a dedicated plan for the type 3 NUFFT (/NNFFT). So it would be AbstractNFFTPlan{T,D,R} and AbstractNNFFTPlan{T,D,R}. The later will get two nodes in the planner. the remaining parts will be the same. We definitely want to have this as a static parameter and I don't see any benefit merging these two into a single type.

I'm glad you have proper mathematical definitions! But there are lots of math typos (bold, etc) - I will file an Issue there.

I am glad to get feedback and corrections! Just note that this is work in progress. So I am aware that the documentation is full of typos. This is the downside of developing things in the open. :-)
But regarding the mathematics I don't want to make anything special. I am certainly biased towards NFFT3 notation but I hope that this is not blocking us here to find common terminology.

in AbstractNFFT you may want to think about something like our plan vs setpts interface. This allows NU points
to move without re-planning the FFT. This was influenced by NFFT3's design...

Isn't that nodes!: https://github.com/JuliaMath/NFFT.jl/blob/master/AbstractNFFTs/src/interface.jl#L55 ? It allows to exchange the nodes of a plan.

Thank you very much already for all these very good points!

[ahbarnett]

Tobias and crew have wrapped FINUFFT in NFFT.jl here:
https://github.com/JuliaMath/NFFT.jl/blob/master/Wrappers/FINUFFT.jl
Which is great. I think type 3 is not wrapped, and directional transforms not done.

I'm not sure there's anything to do right now, until AbstractNFFTs.jl becomes a separate repo, includes mathematical definitions/conventions of transforms (like AbstractFFTs.jl does, although for them that task was very easy), and has various interface changes
(remove references to sigma and m which are specific to certain methods/kernels).

I like the idea of eventual intergration, so I'll leave this open.