[bug report] incorrect 2d type 2 nufft result

[tianrluo]

(Description edited to make assessment easier)

It is likely an error of rounding or mod operations related to coordinates. Play with the x_ofst, y_ofst below to see the error shows/disappears.

The following lines can reproduce it:

import numpy as np
import finufft, cufinufft, pycuda
import pycuda.autoinit

from pycuda import gpuarray

sx, sy = 4, 4
x_ofst, y_ofst = 0.00, 0.00  # change this to reproduce the error

x, y = np.mgrid[-sx//2:sx//2:1., -sy//2:sy//2:1.]
x, y = x+x_ofst, y+y_ofst
xr, yr = x.reshape(-1)/sx * 2 * np.pi, y.reshape(-1)/sy * 2 * np.pi

p = np.zeros((sx, sy), dtype=np.complex128)
p[0][0] = 1

xr_cu, yr_cu, p_cu = gpuarray.to_gpu(xr), gpuarray.to_gpu(yr), gpuarray.to_gpu(p)

plan = finufft.Plan(2, (sx, sy), 1, eps=1e-12)
plan.setpts(xr.reshape(-1), yr.reshape(-1))

plan_cu = cufinufft.cufinufft(2, (sx, sy), eps=1e-12, dtype=np.float64)
plan_cu.set_pts(xr_cu, yr_cu)

kp = plan.execute(p)

kp_cu = gpuarray.GPUArray((sx*sy,), dtype=np.complex128)
plan_cu.execute(kp_cu, p_cu)
kp_cu_np = kp_cu.get()

kp, kp_cu_np = kp.reshape((sx,sy)), kp_cu_np.reshape((sx,sy))

print(kp.real)
print(kp_cu_np.real)

and the outputs are:

# kp.real, correct
array([[ 1. -1.  1. -1.]
 [-1.  1. -1.  1.]
 [ 1. -1.  1. -1.]
 [-1.  1. -1.  1.]])

# kp_cu_np, incorrect
array([[ 1.  0.  1. nan]
 [nan  0.  0.  0.]
 [ 1.  0.  0. nan]
 [ 0. nan nan nan]])

[MelodyShih]

Hi tianrluo, thank you for reporting this. I am able to reproduce the error. I will investigate it and update with you what went wrong.

[tianrluo]

Thanks! Looking forward to a fix!

[MelodyShih]

Hi, I believe the pull request #113 fixes the issue. Would you mind double check that if it fixes on your side as well? Thanks.

[tianrluo]

I just compiled your branch for testing, and the error is solved. I guess you are also going to update the cufinufft package on pip soonish...? Thanks!

[tianrluo]

Hi, @garrettwrong @janden,

Hope things are going well!

I was just wondering if there's a plan to update the cufinufft on pypi regarding this bug fix, or if there are more updates to be waited for...? (Sorry if this is a spam, mentioning you two since you are listed package maintainers on pypi... )

Thanks!