Type 2 3D transform gives different result from (non-cu) finufft for some inputs

[JBlaschke]

trafficstars

This is a bit of a strange error, so I will need to unpack a little here:

For some inputs, the type 2 3D NUFFT implemented in cufinufft is giving different results from the finufft implementation. I am attaching a bitbucket repo containing the finufft output (stored) and the code the runs cufinufft: https://bitbucket.org/jpblaschke/cufinufft-error/src/master/ The get_data.sh script will download the input and output data (approx 3GB).

This is what I see so far:

1. output from our finufft install is stored here: https://bitbucket.org/jpblaschke/cufinufft-error/downloads/output-cpu.npy (we have a fairly old finufft install -- I am blowing off the cobwebs/api changes so that i can confirm that the latest finufft gives the same answer)
2. I tested 2 versions of cufinufft ("old" and "new") which in turn give different results. Use git submodule update --init to clone the latest versions. Especially the output from cufinufft-new.py surprizes me, at it is all zeros.
3. For my other unit tests, the results from finufft and cufinufft agree to within the specified tolerance.

I will add more info as I find out more, but I would appreciate any insights all y'all can give me.

My current theory is that this could have something to do with the input data size/density. I have tested both larger and smaller inputs, for which finufft and cufinufft agree, so it's not as simple as "for large inputs this breaks"

Also tagging @monarin and @irischang020 who brought this bug to my attention.

[JBlaschke]

Btw -- this is (a zoomed in) view of the "good" ouput from finufft:

This is the slightly older version of cufinufft (I didn't zoom in -- it's mainly noise, no spikes):

And this is for the latest cufinufft:

[garrettwrong]

I can reproduce the zeros output running the cufinufft-new.py script using a local cufinufft build. I added a few CUDA error checks to src/cufinufft.cu, but nothing erroneous was reported there.

Often all zeros suggests we've crashed along the way and never touched the array. I confirmed changing the input result array to ones yields output of ones as well. The next step might be to reproduce without the python layer just to rule it out.

[MelodyShih]

Hi Johannas, sorry for taking this long to get back to you. I take a look at it today and found that there is a mistake when you call cufinufft -- the third argument of cufinufft in the new interface is the number of transforms, which is now set to -1:

plan = cufinufft(
        2, shape, -1, tol, dtype=dtype, gpu_method=1, gpu_device_id=0
)

Also, the shape of ugrid and M are not compatible -- ugrid.shape is (151,151,151) while M is 81.

With the following changes:

shape         = ugrid.shape
...
plan = cufinufft(
    2, shape, 1, tol, isign=-1, dtype=dtype, gpu_method=1, gpu_device_id=0
)

, I am able to get the 1e-7 as the sum of difference between the gpu and cpu results: nuvect_cpu - nuvect = (-1.1938391513136895e-07+1.290784103912311e-10j)

[JBlaschke]

Hi Melody,

thanks for catching that -- I noticed something similar when I tested this by calling the C-interface directly. This saves me from having to find the "bug" in the python interface (me).

Let me verify that this works in our code and then close the issue.

@monarin and @irischang020 ☝️

[irischang020]

Hi Melody, I was wondering if there is an example script for 3D type 1 and type 2 NUFFT? I only saw example2d1many.py and example2d2many.py in examples directory. To clarify, (1) what is the number of transforms n_trans in cufinufft.py (the default value seems to be 1)? and (2) why is isign defaults to -1 for type 1 and +1 for type 2? From the change you made, shouldn't type 2 correspond to isign=1, according to the docstring in cufinufft.py?

shape         = ugrid.shape
...
plan = cufinufft(
    2, shape, 1, tol, isign=-1, dtype=dtype, gpu_method=1, gpu_device_id=0
)

Thanks!

[MelodyShih]

Hi Iris,

We unfortunately don't have a 3D python example code yet but it can be done similarly by

plan = cufinufft(1, (N1, N2, N3), n_transf, eps=eps, dtype=dtype)
plan.set_pts(to_gpu(kx), to_gpu(ky), to_gpu(kz)) #kz stores the z coord. of the nonuniform points

One can apply the transform to multiple data (c for type 1, f for type 2) with the same location of nonuniform points at one call of plan.execute, ntransf is the number of data you'd like to apply the transforms to. For example, in example2d1many.py, ntransf=2 and you can see how c is specified accordingly.

Here is the mathematical definition of type 1 and type 2 nonuniform FFTs: https://finufft.readthedocs.io/en/latest/math.html
How to assign the isign parameter depends on what exponential sum you'd like to calculate. I changed to isign=-1 because I vaguely remember from @JBlaschke that this is what the application want.

[irischang020]

Hi Melody, Thanks for the clarification. I'll try what you suggested. I just checked the definition of type 1 and type 2 in finufft, and figured that the docstring about isign in cufinufft.py seems to be swapped:

:param isign: +1 or -1, controls sign of imaginary component in
        complex exponential. Default is -1 for type 1 and +1 for type 2.
        :param dtype: Datatype for this plan (np.float32 or np.float64). \
        Defaults np.float32.

which might need to be changed to be consistent with the code.

        if isign is None:
            if nufft_type == 2:
                isign = -1
            else:
                isign = +1

[MelodyShih]

Hi Iris, Thanks for pointing this out. I fixes the doc in #108. Besides that, I am wondering if the changes fix the issue you observed on your end? I would like to close this issue first if that is the case. :)
And let's discuss other questions elsewhere -- feel free to email me regarding the usage of cufinufft.