merge some improvements from minimal branch

gridder_cxx.h

@@ -837,9 +837,9 @@ template<class T, class T2> class VisServ
   {
   private:
     const Baselines &baselines;
-    const const_mav<uint32_t,1> &idx;
+    const_mav<uint32_t,1> idx;
     T2 vis;
-    const const_mav<T,1> &wgt;
+    const_mav<T,1> wgt;
     size_t nvis;
     bool have_wgt;
 
@@ -877,9 +877,9 @@ template<class T, class T2> class MsServ
   {
   private:
     const Baselines &baselines;
-    const const_mav<uint32_t,1> &idx;
+    const_mav<uint32_t,1> idx;
     T2 ms;
-    const const_mav<T,2> &wgt;
+    const_mav<T,2> wgt;
     size_t nvis;
     bool have_wgt;
 
@@ -1025,7 +1025,7 @@ template<typename T, typename Serv> void grid2x_c
       }
     if (flip) r=conj(r);
     if (do_w_gridding) r*=hlp.Wfac();
-    srv.setVis(ipart, r);
+    srv.addVis(ipart, r);
     }
 }
   }
@@ -1329,7 +1329,7 @@ template<typename T, typename Serv> void x2dirty(
         }
       gconf.grid2dirty_c(cmav(grid), tdirty);
       gconf.apply_wscreen(cmav(tdirty), tdirty, wcur, true);
-  #pragma omp parallel for num_threads(nthreads)
+#pragma omp parallel for num_threads(nthreads)
       for (size_t i=0; i<nx_dirty; ++i)
         for (size_t j=0; j<ny_dirty; ++j)
           dirty(i,j) += tdirty(i,j).real();
@@ -1389,7 +1389,7 @@ template<typename T, typename Serv> void dirty2x(
     apply_wcorr(gconf, tdirty, ES_Kernel_with_correction(supp, nthreads), dw);
 
     vector<complex<T>> vis_loc_;
-    vector<uint32_t> idx_loc_, mapper;
+    vector<uint32_t> idx_loc_, newidx;
     tmpStorage<complex<T>,2> grid_({gconf.Nu(),gconf.Nv()});
     auto grid=grid_.getMav();
     tmpStorage<complex<T>,2> tdirty2_({nx_dirty,ny_dirty});
@@ -1401,31 +1401,44 @@ template<typename T, typename Serv> void dirty2x(
         cout << "Working on plane " << iw << " containing " << nvis_plane[iw]
              << " visibilities" << endl;
       double wcur = wmin+iw*dw;
-      size_t cnt=0;
-      vis_loc_.resize(nvis_plane[iw]);
-      auto vis_loc = mav<complex<T>,1>(vis_loc_.data(), {nvis_plane[iw]});
-      idx_loc_.resize(nvis_plane[iw]);
-      auto idx_loc = mav<uint32_t,1>(idx_loc_.data(), {nvis_plane[iw]});
-      mapper.resize(nvis_plane[iw]);
-      for (size_t ipart=0; ipart<nvis; ++ipart)
-        if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+supp))
-          mapper[cnt++] = ipart;
-  #pragma omp parallel for num_threads(nthreads)
-      for (size_t i=0; i<nvis_plane[iw]; ++i)
-          idx_loc[i] = srv.getIdx(mapper[i]);
-
-  #pragma omp parallel for num_threads(nthreads)
+
+#pragma omp parallel for num_threads(nthreads)
       for (size_t i=0; i<nx_dirty; ++i)
         for (size_t j=0; j<ny_dirty; ++j)
           tdirty2(i,j) = tdirty(i,j);
       gconf.apply_wscreen(cmav(tdirty2), tdirty2, wcur, false);
       gconf.dirty2grid_c(cmav(tdirty2), grid);
-      grid2vis_c(srv.getBaselines(), gconf, cmav(idx_loc), cmav(grid), vis_loc,
-        nullmav<T,1>(), wcur, dw);
 
-  #pragma omp parallel for num_threads(nthreads)
-      for (size_t i=0; i<nvis_plane[iw]; ++i)
-        srv.addVis(mapper[i], vis_loc[i]);
+      newidx.resize(nvis_plane[iw]);
+      size_t cnt=0;
+      if (srv.isMsServ())
+        {
+        for (size_t ipart=0; ipart<nvis; ++ipart)
+          if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+supp))
+            newidx[cnt++] = srv.getIdx(ipart);
+        auto newidx2 = const_mav<uint32_t, 1>(newidx.data(), {newidx.size()});
+        Serv newsrv(srv, newidx2);
+        grid2x_c(gconf, cmav(grid), newsrv, wcur, dw);
+        }
+      else
+        {
+        for (size_t ipart=0; ipart<nvis; ++ipart)
+          if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+supp))
+            newidx[cnt++] = ipart;
+        vis_loc_.resize(nvis_plane[iw]);
+        auto vis_loc = mav<complex<T>,1>(vis_loc_.data(), {nvis_plane[iw]});
+        vis_loc.fill(0);
+        idx_loc_.resize(nvis_plane[iw]);
+        auto idx_loc = mav<uint32_t,1>(idx_loc_.data(), {nvis_plane[iw]});
+#pragma omp parallel for num_threads(nthreads)
+        for (size_t i=0; i<nvis_plane[iw]; ++i)
+          idx_loc[i] = srv.getIdx(newidx[i]);
+        grid2vis_c(srv.getBaselines(), gconf, cmav(idx_loc), cmav(grid), vis_loc,
+          nullmav<T,1>(), wcur, dw);
+#pragma omp parallel for num_threads(nthreads)
+        for (size_t i=0; i<nvis_plane[iw]; ++i)
+          srv.addVis(newidx[i], vis_loc[i]);
+        }
       }
     }
   else

nifty_gridder.cc

@@ -554,6 +554,7 @@ template<typename T> pyarr<complex<T>> Pygrid2vis_c(
   auto wgt2 = make_const_mav<1>(wgt);
   auto res = makeArray<complex<T>>({nvis});
   auto vis = make_mav<1>(res);
+  vis.fill(0);
   {
   py::gil_scoped_release release;
   grid2vis_c<T>(baselines, gconf, idx2, grid2, vis, wgt2);
@@ -751,6 +752,7 @@ template<typename T> pyarr<complex<T>> dirty2vis2(const PyBaselines &baselines,
   auto dirty2 = make_const_mav<2>(dirty_);
   auto vis = makeArray<complex<T>>({idx2.shape(0)});
   auto vis2 = make_mav<1>(vis);
+  vis2.fill(0);
   {
   py::gil_scoped_release release;
   vis2.fill(0);