Merge branch 'test_mask_feature' into 'ducc0'

Test mask feature

See merge request !34

python/gridder_cxx.h

@@ -607,7 +607,7 @@ template<size_t supp, bool wgrid, typename T> class HelperX2g2
       auto iv0old = iv0;
       gconf.getpix(in.u, in.v, u, v, iu0, iv0);
       T x0 = (iu0-T(u))*2+(supp-1);
-      T y0 = (iv0-T(v))*2+(supp-1); 
+      T y0 = (iv0-T(v))*2+(supp-1);
       if constexpr(wgrid)
         krn.eval2s(x0, y0, T(xdw*(w0-in.w)), &buf.simd[0]);
       else
@@ -695,7 +695,7 @@ template<size_t supp, bool wgrid, typename T> class HelperG2x2
       auto iv0old = iv0;
       gconf.getpix(in.u, in.v, u, v, iu0, iv0);
       T x0 = (iu0-T(u))*2+(supp-1);
-      T y0 = (iv0-T(v))*2+(supp-1); 
+      T y0 = (iv0-T(v))*2+(supp-1);
       if constexpr(wgrid)
         krn.eval2s(x0, y0, T(xdw*(w0-in.w)), &buf.simd[0]);
       else
@@ -840,7 +840,7 @@ template<bool wgrid, typename T, typename Serv> void x2grid_c
   gconf.timers.push("gridding proper");
   checkShape(grid.shape(), {gconf.Nu(), gconf.Nv()});
 
-  if constexpr (is_same<T, float>::value) 
+  if constexpr (is_same<T, float>::value)
     switch(gconf.Supp())
       {
       case  4: x2grid_c_helper< 4, wgrid>(gconf, srv, grid, w0, dw); break;
@@ -923,7 +923,7 @@ template<bool wgrid, typename T, typename Serv> void grid2x_c
   gconf.timers.push("degridding proper");
   checkShape(grid.shape(), {gconf.Nu(), gconf.Nv()});
 
-  if constexpr (is_same<T, float>::value) 
+  if constexpr (is_same<T, float>::value)
     switch(gconf.Supp())
       {
       case  4: grid2x_c_helper< 4, wgrid>(gconf, grid, srv, w0, dw); break;
@@ -1193,7 +1193,7 @@ template<typename T> void report(const GridderConfig<T> &gconf, size_t nvis,
   if (x0*x0+y0*y0>1.)
     nmin = -sqrt(abs(1.-x0*x0-y0*y0))-1.;
   double dw = 0.5/gconf.Ofactor()/abs(nmin);
-  cout << "  w=[" << wmin << "; " << wmax << "], nmin=" << nmin << ", dw=" << dw
+  cout << "  w=[" << wmin << "; " << wmax << "], min(n-1)=" << nmin << ", dw=" << dw
        << ", wmax/dw=" << wmax/dw << endl;
   }
 
@@ -1442,7 +1442,7 @@ template<typename T> auto scanData(const Baselines &baselines, const mav<complex
     });
   timers.pop();
   return make_tuple(wmin, wmax, nvis, mask_out);
-  } 
+  }
 
 // Note to self: divide_by_n should always be true when doing Bayesian imaging,
 // but wsclean needs it to be false, so this must be kept as a parameter.

python/test/test_wgridder.py

@@ -39,7 +39,7 @@ def _l2error(a, b):
 
 
 def explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, xpixsize, ypixsize,
-                     apply_w):
+                     apply_w, mask):
     speedoflight = 299792458.
     x, y = np.meshgrid(*[-ss/2 + np.arange(ss) for ss in [nxdirty, nydirty]],
                        indexing='ij')
@@ -55,6 +55,8 @@ def explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, xpixsize, ypixsize,
         n = 1.
     for row in range(ms.shape[0]):
         for chan in range(ms.shape[1]):
+            if mask is not None and mask[row, chan] == 0:
+                continue
             phase = (freq[chan]/speedoflight *
                      (x*uvw[row, 0] + y*uvw[row, 1] - uvw[row, 2]*nm1))
             if wgt is None:
@@ -74,9 +76,10 @@ def explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, xpixsize, ypixsize,
 @pmp("singleprec", (True, False))
 @pmp("wstacking", (True, False))
 @pmp("use_wgt", (True, False))
+@pmp("use_mask", (False, True))
 @pmp("nthreads", (1, 2, 7))
 def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
-                              singleprec, wstacking, use_wgt, nthreads):
+                              singleprec, wstacking, use_wgt, nthreads, use_mask):
     if singleprec and epsilon < 5e-5:
         pytest.skip()
     rng = np.random.default_rng(42)
@@ -87,6 +90,7 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
     uvw = (rng.random((nrow, 3))-0.5)/(pixsizey*f0/speedoflight)
     ms = rng.random((nrow, nchan))-0.5 + 1j*(rng.random((nrow, nchan))-0.5)
     wgt = rng.uniform(0.9, 1.1, (nrow, nchan)) if use_wgt else None
+    mask = (rng.uniform(0, 1, (nrow, nchan)) > 0.5).astype(np.uint8) if use_mask else None
     dirty = rng.random((nxdirty, nydirty))-0.5
     nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
     if nu & 1:
@@ -101,9 +105,9 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
         if wgt is not None:
             wgt = wgt.astype("f4")
     dirty2 = ng.ms2dirty(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex,
-                         pixsizey, nu, nv, epsilon, wstacking, nthreads, 0).astype("f8")
+                         pixsizey, nu, nv, epsilon, wstacking, nthreads, 0, mask).astype("f8")
     ms2 = ng.dirty2ms(uvw, freq, dirty, wgt, pixsizex, pixsizey, nu, nv, epsilon,
-                      wstacking, nthreads+1, 0).astype("c16")
+                      wstacking, nthreads+1, 0, mask).astype("c16")
     ref = max(my_vdot(ms,ms).real, my_vdot(ms2,ms2).real,
               my_vdot(dirty, dirty).real, my_vdot(dirty2, dirty2).real)
     tol = 1e-5*ref if singleprec else 1e-11*ref
@@ -119,9 +123,10 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
 @pmp("singleprec", (False,))
 @pmp("wstacking", (False, True))
 @pmp("use_wgt", (True,))
+@pmp("use_mask", (True,))
 @pmp("nthreads", (1, 2, 7))
 @pmp("fov", (1., 20.))
-def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, fov, nthreads):
+def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, use_mask, fov, nthreads):
     if singleprec and epsilon < 5e-5:
         pytest.skip()
     rng = np.random.default_rng(42)
@@ -132,6 +137,7 @@ def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
     uvw = (rng.random((nrow, 3))-0.5)/(pixsizex*f0/speedoflight)
     ms = rng.random((nrow, nchan))-0.5 + 1j*(rng.random((nrow, nchan))-0.5)
     wgt = rng.uniform(0.9, 1.1, (nrow, 1)) if use_wgt else None
+    mask = (rng.uniform(0, 1, (nrow, nchan)) > 0.5).astype(np.uint8) if use_mask else None
     wgt = np.broadcast_to(wgt, (nrow, nchan)) if use_wgt else None
     nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
     if nu & 1:
@@ -145,6 +151,8 @@ def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
         if wgt is not None:
             wgt = wgt.astype("f4")
     dirty = ng.ms2dirty(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex,
-                        pixsizey, nu, nv, epsilon, wstacking, nthreads, 0).astype("f8")
-    ref = explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex, pixsizey, wstacking)
+                        pixsizey, nu, nv, epsilon, wstacking, nthreads,
+                        0, mask).astype("f8")
+    ref = explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex,
+                           pixsizey, wstacking, mask)
     assert_allclose(_l2error(dirty, ref), 0, atol=epsilon)