Merge branch 'adjust_gridder' into 'NIFTy_7'

Switch to ducc0.wgridder

See merge request !549

Dockerfile

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

README.md

@@ -50,9 +50,7 @@ Installation
 
 Optional dependencies:
 - [DUCC0](https://gitlab.mpcdf.mpg.de/mtr/ducc) for faster FFTs, spherical
-    harmonic transforms, and non-uniform Fourier transforms
-- [nifty_gridder](https://gitlab.mpcdf.mpg.de/ift/nifty_gridder) (for radio
-    interferometry responses)
+    harmonic transforms, and radio interferometry gridding support
 - [mpi4py](https://mpi4py.scipy.org) (for MPI-parallel execution)
 - [matplotlib](https://matplotlib.org/) (for field plotting)
 
@@ -86,10 +84,6 @@ If this library is present, NIFTy will detect it automatically and prefer
 DUCC's FFT is compiled with optimizations for the host CPU and can provide
 significantly faster transforms.
 
-Support for the radio interferometry gridder is added via:
-
-    pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
-
 MPI support is added via:
 
     sudo apt-get install python3-mpi4py

demos/bench_gridder.py

-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-# Copyright(C) 2019-2020 Max-Planck-Society
-# Author: Martin Reinecke
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-
-from time import time
-
-import matplotlib.pyplot as plt
-import numpy as np
-
-import nifty7 as ift
-
-
-def main():
-    N0s, a0s, b0s, c0s = [], [], [], []
-
-    for ii in range(10, 26):
-        nu = 1024
-        nv = 1024
-        N = int(2**ii)
-        print('N = {}'.format(N))
-
-        rng = ift.random.current_rng()
-        uv = rng.uniform(-.5, .5, (N, 2))
-        vis = rng.normal(0., 1., N) + 1j*rng.normal(0., 1., N)
-
-        uvspace = ift.RGSpace((nu, nv))
-
-        visspace = ift.UnstructuredDomain(N)
-
-        img = rng.standard_normal((nu, nv))
-        img = ift.makeField(uvspace, img)
-
-        t0 = time()
-        GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv)
-        vis = ift.makeField(visspace, vis)
-        op = GM.getFull().adjoint
-        t1 = time()
-        op(img).val
-        t2 = time()
-        op.adjoint(vis).val
-        t3 = time()
-        print(t2-t1, t3-t2)
-        N0s.append(N)
-        a0s.append(t1 - t0)
-        b0s.append(t2 - t1)
-        c0s.append(t3 - t2)
-
-    print('Measure rest operator')
-    sc = ift.StatCalculator()
-    op = GM.getRest().adjoint
-    for _ in range(10):
-        t0 = time()
-        res = op(img)
-        sc.add(time() - t0)
-    print('FFT shape', res.shape)
-
-    plt.scatter(N0s, a0s, label='Gridder mr')
-    plt.legend()
-    # no idea why this is necessary, but if it is omitted, the range is wrong
-    plt.ylim(min(a0s), max(a0s))
-    plt.ylabel('time [s]')
-    plt.title('Initialization')
-    plt.loglog()
-    plt.savefig('bench0.png')
-    plt.close()
-
-    plt.scatter(N0s, b0s, color='k', marker='^', label='Gridder mr times')
-    plt.scatter(N0s, c0s, color='k', label='Gridder mr adjoint times')
-    plt.axhline(sc.mean, label='FFT')
-    plt.axhline(sc.mean + np.sqrt(sc.var))
-    plt.axhline(sc.mean - np.sqrt(sc.var))
-    plt.legend()
-    plt.ylabel('time [s]')
-    plt.title('Apply')
-    plt.loglog()
-    plt.savefig('bench1.png')
-    plt.close()
-
-
-if __name__ == '__main__':
-    main()

docs/source/installation.rst

@@ -23,10 +23,6 @@ If this library is present, NIFTy will detect it automatically and prefer
 DUCC's FFT is compiled with optimizations for the host CPU and can provide
 significantly faster transforms.
 
-Support for the radio interferometry gridder is added via::
-
-    pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
-
 MPI support is added via::
 
     sudo apt-get install python3-mpi4py

src/__init__.py

@@ -86,7 +86,7 @@ from .library.light_cone_operator import LightConeOperator
 from .library.wiener_filter_curvature import WienerFilterCurvature
 from .library.adjust_variances import (make_adjust_variances_hamiltonian,
                                        do_adjust_variances)
-from .library.gridder import GridderMaker
+from .library.gridder import Gridder
 from .library.correlated_fields import CorrelatedFieldMaker
 
 from . import extra

src/library/gridder.py

@@ -24,88 +24,42 @@ from ..operators.linear_operator import LinearOperator
 from ..sugar import makeDomain, makeField
 
 
-class GridderMaker(object):
-    def __init__(self, dirty_domain, uv, eps=2e-13):
-        import nifty_gridder
-        dirty_domain = makeDomain(dirty_domain)
-        if (len(dirty_domain) != 1 or not isinstance(dirty_domain[0], RGSpace)
-                or not len(dirty_domain.shape) == 2):
-            raise ValueError("need dirty_domain with exactly one 2D RGSpace")
+class Gridder(LinearOperator):
+    def __init__(self, target, uv, eps=2e-10, nthreads=1):
+        self._capability = self.TIMES | self.ADJOINT_TIMES
+        self._target = makeDomain(target)
+        if (len(self._target) != 1 or not isinstance(self._target[0], RGSpace)
+                or not len(self._target.shape) == 2):
+            raise ValueError("need target with exactly one 2D RGSpace")
         if uv.ndim != 2:
             raise ValueError("uv must be a 2D array")
         if uv.shape[1] != 2:
             raise ValueError("second dimension of uv must have length 2")
-        dstx, dsty = dirty_domain[0].distances
-        # wasteful hack to adjust to shape required by nifty_gridder
-        uvw = np.empty((uv.shape[0], 3), dtype=np.float64)
-        uvw[:, 0:2] = uv
-        uvw[:, 2] = 0.
-        # Scale uv such that 0<uv<=1 which is assumed by nifty_gridder
-        uvw[:, 0] = uvw[:, 0]*dstx
-        uvw[:, 1] = uvw[:, 1]*dsty
-        speedOfLight = 299792458.
-        bl = nifty_gridder.Baselines(uvw, np.array([speedOfLight]))
-        nxdirty, nydirty = dirty_domain.shape
-        gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.)
-        nu, nv = gconf.Nu(), gconf.Nv()
-        self._idx = nifty_gridder.getIndices(
-            bl, gconf, np.zeros((uv.shape[0], 1), dtype=np.bool))
-        self._bl = bl
-
-        du, dv = 1./(nu*dstx), 1./(nv*dsty)
-        grid_domain = RGSpace([nu, nv], distances=[du, dv], harmonic=True)
-
-        self._rest = _RestOperator(dirty_domain, grid_domain, gconf)
-        self._gridder = RadioGridder(grid_domain, bl, gconf, self._idx)
-
-    def getGridder(self):
-        return self._gridder
-
-    def getRest(self):
-        return self._rest
-
-    def getFull(self):
-        return self.getRest() @ self._gridder
-
-    def ms2vis(self, x):
-        return self._bl.ms2vis(x, self._idx)
-
-
-class _RestOperator(LinearOperator):
-    def __init__(self, dirty_domain, grid_domain, gconf):
-        self._domain = makeDomain(grid_domain)
-        self._target = makeDomain(dirty_domain)
-        self._gconf = gconf
-        self._capability = self.TIMES | self.ADJOINT_TIMES
-
-    def apply(self, x, mode):
-        self._check_input(x, mode)
-        res = x.val
-        if mode == self.TIMES:
-            res = self._gconf.grid2dirty(res)
-        else:
-            res = self._gconf.dirty2grid(res)
-        return makeField(self._tgt(mode), res)
-
-
-class RadioGridder(LinearOperator):
-    def __init__(self, grid_domain, bl, gconf, idx):
         self._domain = DomainTuple.make(
-            UnstructuredDomain((bl.Nrows())))
-        self._target = DomainTuple.make(grid_domain)
-        self._bl = bl
-        self._gconf = gconf
-        self._idx = idx
-        self._capability = self.TIMES | self.ADJOINT_TIMES
+            UnstructuredDomain((uv.shape[0])))
+        # wasteful hack to adjust to shape required by ducc0.wgridder
+        self._uvw = np.empty((uv.shape[0], 3), dtype=np.float64)
+        self._uvw[:, 0:2] = uv
+        self._uvw[:, 2] = 0.
+        self._eps = float(eps)
+        self._nthreads = int(nthreads)
+
 
     def apply(self, x, mode):
-        import nifty_gridder
         self._check_input(x, mode)
+        speedOfLight = 299792458.
+        freq = np.array([speedOfLight])
+        x = x.val
+        nxdirty, nydirty = self._target[0].shape
+        nu, nv = max(2*nxdirty, 16), max(2*nydirty, 16)
+        dstx, dsty = self._target[0].distances
+        from ducc0.wgridder import ms2dirty, dirty2ms
         if mode == self.TIMES:
-            x = self._bl.ms2vis(x.val.reshape((-1, 1)), self._idx)
-            res = nifty_gridder.vis2grid(self._bl, self._gconf, self._idx, x)
+            res = ms2dirty(self._uvw, freq, x.reshape((-1,1)), None, nxdirty,
+                           nydirty, dstx, dsty, nu, nv,
+                           self._eps, False, self._nthreads, 0)
         else:
-            res = nifty_gridder.grid2vis(self._bl, self._gconf, self._idx,
-                                         x.val)
-            res = self._bl.vis2ms(res, self._idx).reshape((-1,))
+            res = dirty2ms(self._uvw, freq, x, None, dstx, dsty, nu, nv,
+                           self._eps, False, self._nthreads, 0)
+            res = res.reshape((-1,))
         return makeField(self._tgt(mode), res)

test/test_operators/test_nft.py

@@ -46,10 +46,9 @@ def test_gridding(nu, nv, N, eps):
     dstx, dsty = dom.distances
     uv[:, 0] = uv[:, 0]/dstx
     uv[:, 1] = uv[:, 1]/dsty
-    GM = ift.GridderMaker(dom, uv=uv, eps=eps)
+    Op = ift.Gridder(dom, uv=uv, eps=eps)
     vis2 = ift.makeField(ift.UnstructuredDomain(vis.shape), vis)
 
-    Op = GM.getFull()
     pynu = Op(vis2).val
     # DFT
     x, y = np.meshgrid(
@@ -72,8 +71,7 @@ def test_cartesian():
     tmp = np.vstack([uu[None, :], vv[None, :]])
     uv = np.transpose(tmp, (2, 1, 0)).reshape(-1, 2)
 
-    GM = ift.GridderMaker(dom, uv=uv)
-    op = GM.getFull().adjoint
+    op = ift.Gridder(dom, uv=uv).adjoint
 
     fld = ift.from_random(dom, 'normal')
     arr = fld.val
@@ -95,16 +93,9 @@ def test_cartesian():
 def test_build(nu, nv, N, eps):
     dom = ift.RGSpace([nu, nv])
     uv = ift.random.current_rng().random((N, 2)) - 0.5
-    GM = ift.GridderMaker(dom, uv=uv, eps=eps)
-    R0 = GM.getGridder()
-    R1 = GM.getRest()
-    R = R1@R0
-    RF = GM.getFull()
+    RF = ift.Gridder(dom, uv=uv, eps=eps)
 
     # Consistency checks
     flt = np.float64
     cmplx = np.complex128
-    ift.extra.check_linear_operator(R0, cmplx, flt, only_r_linear=True)
-    ift.extra.check_linear_operator(R1, flt, flt)
-    ift.extra.check_linear_operator(R, cmplx, flt, only_r_linear=True)
     ift.extra.check_linear_operator(RF, cmplx, flt, only_r_linear=True)