diff --git a/Dockerfile b/Dockerfile index 157129f7c223252d23f973daa26544175d660beb..a713bdfa5b73018cf248ea372f346b1846d0d423 100644 --- a/Dockerfile +++ b/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/* diff --git a/nifty5/library/bench_gridder.py b/nifty5/library/bench_gridder.py index 88fd8b3df68b4dff4b6008dc76b86533ddc5eb5d..ff8cb183863b84e78b045b8b59822fc90ce08065 100644 --- a/nifty5/library/bench_gridder.py +++ b/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) diff --git a/nifty5/library/gridder.py b/nifty5/library/gridder.py index 80d168395adf0270f1af1ded65eed175ccfe1f3b..8baaa1cc7a26b461669f6c0087398aa6d1a58a4e 100644 --- a/nifty5/library/gridder.py +++ b/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) diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py index b4834a768a363d806f312a3e5bb5446406f14267..29da0e5b8869d2631947e1c623b19e29d236e473 100644 --- a/test/test_operators/test_nft.py +++ b/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)