latest updates from gridder project

Dockerfile

@@ -14,6 +14,7 @@ RUN apt-get update && apt-get install -y \
   # more optional NIFTy dependencies
   && pip3 install pyfftw \
   && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \
+  && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git \
   && pip3 install jupyter \
   && rm -rf /var/lib/apt/lists/*
 

nifty5/library/bench_gridder.py

@@ -24,7 +24,7 @@ for ii in range(1, 8):
 
     visspace = ift.UnstructuredDomain(N)
 
-    img = np.random.randn(nu*nv) + 1j*np.random.randn(nu*nv)
+    img = np.random.randn(nu*nv)
     img = img.reshape((nu, nv))
     img = ift.from_global_data(uvspace, img)
 

nifty5/library/gridder.py

@@ -20,7 +20,7 @@ import numpy as np
 from ..domain_tuple import DomainTuple
 from ..domains.rg_space import RGSpace
 from ..domains.unstructured_domain import UnstructuredDomain
-from ..fft import fftn, ifftn
+from ..fft import hartley
 from ..operators.linear_operator import LinearOperator
 from ..sugar import from_global_data, makeDomain
 
@@ -50,7 +50,7 @@ class GridderMaker(object):
         self._rest = _RestOperator(domain, oversampled_domain, r2lamb)
 
     def getReordering(self, uv):
-        from testgridder import peanoindex
+        from nifty_gridder import peanoindex
         nu2, nv2 = self._rest._domain.shape
         return peanoindex(uv, nu2, nv2)
 
@@ -79,27 +79,28 @@ class _RestOperator(LinearOperator):
         rng = np.arange(nv)
         k = np.minimum(rng, nv-rng)
         c = np.pi*r2lamb/nv2**2
-        self._deconv_v = np.roll(np.exp(c*k**2)/r2lamb, -nv//2).reshape((1, -1))
+        self._deconv_v = np.roll(
+            np.exp(c*k**2)/r2lamb, -nv//2).reshape((1, -1))
         self._capability = self.TIMES | self.ADJOINT_TIMES
 
     def apply(self, x, mode):
         self._check_input(x, mode)
         nu, nv = self._target.shape
-        res = x.to_global_data_rw()
+        res = x.to_global_data()
         if mode == self.TIMES:
-            res = ifftn(res)*res.size
+            res = hartley(res)
             res = np.roll(res, (nu//2, nv//2), axis=(0, 1))
             res = res[:nu, :nv]
             res *= self._deconv_u
             res *= self._deconv_v
         else:
-            res *= self._deconv_u
+            res = res*self._deconv_u
             res *= self._deconv_v
             nu2, nv2 = self._domain.shape
-            res = np.pad(res, ((0, nu2-nu), (0, nv2-nv)), 'constant',
+            res = np.pad(res, ((0, nu2-nu), (0, nv2-nv)), mode='constant',
                          constant_values=0)
             res = np.roll(res, (-nu//2, -nv//2), axis=(0, 1))
-            res = fftn(res)
+            res = hartley(res)
         return from_global_data(self._tgt(mode), res)
 
 
@@ -113,14 +114,16 @@ class RadioGridder(LinearOperator):
         self._uv = uv  # FIXME: should we write-protect this?
 
     def apply(self, x, mode):
-        from testgridder import to_grid, from_grid
+        from nifty_gridder import to_grid, from_grid
         self._check_input(x, mode)
         nu2, nv2 = self._target.shape
         x = x.to_global_data()
         if mode == self.TIMES:
-            res = to_grid(self._uv, x, nu2, nv2, self._nspread,
-                          self._r2lamb)
+            res = to_grid(self._uv, x, nu2, nv2, self._nspread, self._r2lamb)
+            res += np.conj(np.roll(res[::-1, ::-1], (1, 1), axis=(0, 1)))
+            res = 0.5*(res.real+res.imag)
         else:
-            res = from_grid(self._uv, x, nu2, nv2, self._nspread,
-                            self._r2lamb)
+            mirr = np.roll(x[::-1, ::-1], (1, 1), axis=(0, 1))
+            x = 0.5*(x+mirr + 1j*(x-mirr))
+            res = from_grid(self._uv, x, nu2, nv2, self._nspread, self._r2lamb)
         return from_global_data(self._tgt(mode), res)

test/test_operators/test_nft.py

@@ -47,7 +47,7 @@ def test_gridding(nu, nv, N):
         *[-ss/2 + np.arange(ss) for ss in [nu, nv]], indexing='ij')
     dft = pynu*0.
     for i in range(N):
-        dft += vis.val[i]*np.exp(2j*np.pi*(x*uv[i, 0] + y*uv[i, 1]))
+        dft += (vis.val[i]*np.exp(2j*np.pi*(x*uv[i, 0] + y*uv[i, 1]))).real
     assert_allclose(dft, pynu)
 
 
@@ -68,8 +68,9 @@ def test_build(nu, nv, N, eps):
     RF = GM.getFull(uv)
 
     # Consistency checks
-    dt = np.complex
-    ift.extra.consistency_check(R0, domain_dtype=dt, target_dtype=dt)
-    ift.extra.consistency_check(R1, domain_dtype=dt, target_dtype=dt)
-    ift.extra.consistency_check(R, domain_dtype=dt, target_dtype=dt)
-    ift.extra.consistency_check(RF, domain_dtype=dt, target_dtype=dt)
+    flt = np.float64
+    cmplx = np.complex128
+    ift.extra.consistency_check(R0, cmplx, flt, only_r_linear=True)
+    ift.extra.consistency_check(R1, flt, flt)
+    ift.extra.consistency_check(R, cmplx, flt, only_r_linear=True)
+    ift.extra.consistency_check(RF, cmplx, flt, only_r_linear=True)