ENH: Added 2D type 1 function interfaces

Un-tested. Tests for the 2D case are not implemented yet.

Author: Thomas Arildsen

nufft/__init__.py

@@ -8,5 +8,5 @@
     __NUFFT_SETUP__ = False
 
 if not __NUFFT_SETUP__:
-    __all__ = ["nufft1freqs", "nufft1", "nufft2", "nufft3"]
-    from .nufft import nufft1freqs, nufft1, nufft2, nufft3
+    __all__ = ["nufft1d1freqs", "nufft1d1", "nufft1d2", "nufft1d3", "nufft2d1"]
+    from .nufft import nufft1d1freqs, nufft1d1, nufft1d2, nufft1d3, nufft2d1

nufft/nufft.py

@@ -2,21 +2,22 @@
 
 from __future__ import division, print_function
 
-__all__ = ["nufft1freqs", "nufft1", "nufft2", "nufft3"]
+__all__ = ["nufft1d1freqs", "nufft1d1", "nufft1d2", "nufft1d3", "nufft2d1"]
 
 import numpy as np
 from ._nufft import (
     dirft1d1, nufft1d1f90,
     dirft1d2, nufft1d2f90,
     dirft1d3, nufft1d3f90,
+    dirft2d1, nufft2d1f90,
 )
 
 
-def nufft1freqs(ms, df=1.0):
+def nufft1d1freqs(ms, df=1.0):
     return df * (np.arange(-ms // 2, ms // 2) + ms % 2)
 
 
-def nufft1(x, y, ms, df=1.0, eps=1e-15, iflag=1, direct=False):
+def nufft1d1(x, y, ms, df=1.0, eps=1e-15, iflag=1, direct=False):
     # Make sure that the data are properly formatted.
     x = np.ascontiguousarray(x, dtype=np.float64)
     y = np.ascontiguousarray(y, dtype=np.complex128)
@@ -34,7 +35,7 @@ def nufft1(x, y, ms, df=1.0, eps=1e-15, iflag=1, direct=False):
     return p
 
 
-def nufft2(x, p, df=1.0, eps=1e-15, iflag=1, direct=False):
+def nufft1d2(x, p, df=1.0, eps=1e-15, iflag=1, direct=False):
     # Make sure that the data are properly formatted.
     x = np.ascontiguousarray(x, dtype=np.float64)
     p = np.ascontiguousarray(p, dtype=np.complex128)
@@ -50,7 +51,7 @@ def nufft2(x, p, df=1.0, eps=1e-15, iflag=1, direct=False):
     return y
 
 
-def nufft3(x, y, f, eps=1e-15, iflag=1, direct=False):
+def nufft1d3(x, y, f, eps=1e-15, iflag=1, direct=False):
     # Make sure that the data are properly formatted.
     x = np.ascontiguousarray(x, dtype=np.float64)
     y = np.ascontiguousarray(y, dtype=np.complex128)
@@ -69,3 +70,22 @@ def nufft3(x, y, f, eps=1e-15, iflag=1, direct=False):
         if flag:
             raise RuntimeError("nufft1d3 failed with code {0}".format(flag))
     return p / len(x)
+
+
+def nufft2d1(x, y, z, ms, mt, df=1.0, eps=1e-15, iflag=1, direct=False):
+    # Make sure that the data are properly formatted.
+    x = np.ascontiguousarray(x, dtype=np.float64)
+    y = np.ascontiguousarray(x, dtype=np.float64)
+    z = np.ascontiguousarray(y, dtype=np.complex128)
+    if len(x) != len(y) or len(y) != len(z):
+        raise ValueError("Dimension mismatch")
+
+    # Run the Fortran code.
+    if direct:
+        p = dirft2d1(x * df, y * df, z, iflag, ms, mt)
+    else:
+        p, flag = nufft2d1f90(x * df, y * df, z, iflag, eps, ms, mt)
+        # Check the output and return.
+        if flag:
+            raise RuntimeError("nufft2d1 failed with code {0}".format(flag))
+    return p

nufft/nufft1d.pyf

@@ -65,5 +65,29 @@ python module _nufft ! in
             complex*16 dimension(nk),intent(out),depend(nk) :: fk
         end subroutine dirft1d3
 
-    end interface
+        subroutine nufft2d1f90(nj,xj,yj,cj, iflag,eps, ms,mt,fk,ier) ! in :_nufft:src/nufft1d/nufft2df90.f
+            integer, optional,check(len(xj)>=nj),depend(xj) :: nj=len(xj)
+            real*8 dimension(nj) :: xj
+            real*8 dimension(nj) :: yj
+            complex*16 dimension(nj),depend(nj) :: cj
+            integer :: iflag
+            real*8 :: eps
+            integer :: ms
+            integer :: mt
+            complex*16 dimension(ms, mt),intent(out),depend(ms, mt) :: fk
+            integer,intent(out) :: ier
+        end subroutine nufft2d1f90
+
+        subroutine dirft2d1(nj,xj,yj,cj, iflag, ms,mt,fk) ! in :_nufft:src/nufft1d/dirft2d.f
+            integer, optional,check(len(xj)>=nj),depend(xj) :: nj=len(xj)
+            real*8 dimension(nj) :: xj
+            real*8 dimension(nj) :: yj
+            complex*16 dimension(nj),depend(nj) :: cj
+            integer :: iflag
+            integer :: ms
+            integer :: mt
+            complex*16 dimension(ms, mt),intent(out),depend(ms, mt) :: fk
+        end subroutine dirft2d1
+
+     end interface
 end python module _nufft

nufft/test.py

@@ -5,7 +5,7 @@
 __all__ = ["test_type_1", "test_type_2", "test_type_3", "test_1_and_3"]
 
 import numpy as np
-from .nufft import nufft1freqs, nufft1, nufft2, nufft3
+from .nufft import nufft1d1freqs, nufft1d1, nufft1d2, nufft1d3, nufft2d1
 
 
 def _get_data():
@@ -20,6 +20,21 @@ def _get_data():
     f = 48 * np.cos((np.arange(ms) + 1) * np.pi / ms)
     return x, y, f
 
+def _get_data_2d():
+    ms = 20
+    mt = 20
+    nj = 128
+    k1 = np.arange(-0.5 * nj, 0.5 * nj)
+    j = k1 + 0.5 * nj + 1
+    x = np.pi * np.cos(-np.pi * j / nj)
+    y = np.pi * np.sin(-np.pi * j / nj)
+    z = np.empty_like(x, dtype=np.complex128)
+    z.real = np.sin(np.pi * j / nj)
+    z.imag = np.cos(np.pi * j / nj)
+    ms_val, mt_val = np.meshgrid(np.arange(ms)/ms, np.arange(mt)/mt, sparse=True)
+    f = 48 * np.cos((ms_val + mt_val + 1) * np.pi)
+    return x, y, f
+
 
 def _get_data_roundtrip():
     ms = 512

setup.py

@@ -28,7 +28,8 @@
 
 # Set up the compiled extension.
 sources = list(map(os.path.join("src", "nufft1d", "{0}").format,
-                   ["dfftpack.f", "dirft1d.f", "next235.f", "nufft1df90.f"]))
+                   ["dfftpack.f", "dirft1d.f", "dirft2d.f",
+                    "next235.f", "nufft1df90.f", "nufft2df90.f"]))
 sources += [os.path.join("nufft", "nufft1d.pyf")]
 extensions = [Extension("nufft._nufft", sources=sources)]