simplifications

gridder_cxx.h

@@ -823,6 +823,28 @@ template<typename T, typename T2=complex<T>> class Helper
       }
   };
 
+template<class T, class Serv> class SubServ
+  {
+  private:
+    const Serv &srv;
+    const_mav<uint32_t,1> subidx;
+
+  public:
+    SubServ(const Serv &orig, const const_mav<uint32_t,1> &subidx_)
+      : srv(orig), subidx(subidx_){}
+    size_t Nvis() const { return subidx.size(); }
+    const Baselines &getBaselines() const { return srv.getBaselines(); }
+    UVW getCoord(size_t i) const
+      { return srv.getCoord(subidx[i]); }
+    complex<T> getVis(size_t i) const
+      { return srv.getVis(subidx[i]); }
+    uint32_t getIdx(size_t i) const { return srv.getIdx(subidx[i]); }
+    void setVis (size_t i, const complex<T> &v) const
+      { srv.setVis(subidx[i], v); }
+    void addVis (size_t i, const complex<T> &v) const
+      { srv.addVis(subidx[i], v); }
+  };
+
 template<class T, class T2> class VisServ
   {
   private:
@@ -842,9 +864,8 @@ template<class T, class T2> class VisServ
       checkShape(idx.shape(), {nvis});
       if (have_wgt) checkShape(wgt.shape(), {nvis});
       }
-    VisServ(const VisServ &orig, const const_mav<uint32_t,1> &/*newidx*/)
-      : baselines(orig.baselines), idx(orig.idx), vis(orig.vis), wgt(orig.wgt)
-      { myfail("must not happen"); }
+    SubServ<T, VisServ> getSubserv(const const_mav<uint32_t,1> &subidx) const
+      { return SubServ<T, VisServ>(*this, subidx); }
     size_t Nvis() const { return nvis; }
     const Baselines &getBaselines() const { return baselines; }
     UVW getCoord(size_t i) const
@@ -856,7 +877,6 @@ template<class T, class T2> class VisServ
       { vis(i) = have_wgt ? v*wgt(i) : v; }
     void addVis (size_t i, const complex<T> &v) const
       { vis(i) += have_wgt ? v*wgt(i) : v; }
-    bool isMsServ() const { return false; }
   };
 template<class T, class T2> VisServ<T, T2> makeVisServ
   (const Baselines &baselines,
@@ -882,8 +902,8 @@ template<class T, class T2> class MsServ
       checkShape(ms.shape(), {baselines.Nrows(), baselines.Nchannels()});
       if (have_wgt) checkShape(wgt.shape(), ms.shape());
       }
-    MsServ(const MsServ &orig, const const_mav<uint32_t,1> &newidx)
-      : MsServ(orig.baselines, newidx, orig.ms, orig.wgt) {}
+    SubServ<T, MsServ> getSubserv(const const_mav<uint32_t,1> &subidx) const
+      { return SubServ<T, MsServ>(*this, subidx); }
     size_t Nvis() const { return nvis; }
     const Baselines &getBaselines() const { return baselines; }
     UVW getCoord(size_t i) const
@@ -905,7 +925,6 @@ template<class T, class T2> class MsServ
       auto rc = baselines.getRowChan(idx(i));
       ms(rc.row, rc.chan) += have_wgt ? v*wgt(rc.row, rc.chan) : v;
       }
-    bool isMsServ() const { return true; }
   };
 template<class T, class T2> MsServ<T, T2> makeMsServ
   (const Baselines &baselines,
@@ -1270,9 +1289,7 @@ template<typename T, typename Serv> void x2dirty(
     wstack_common(gconf, srv, wmin, dw, nplanes, nvis_plane, minplane, verbosity);
 
     dirty.fill(0);
-    vector<complex<T>> vis_loc_;
-    vector<uint32_t> idx_loc_;
-    vector<uint32_t> newidx;
+    vector<uint32_t> subidx;
     tmpStorage<complex<T>,2> grid_({gconf.Nu(),gconf.Nv()});
     auto grid=grid_.getMav();
     tmpStorage<complex<T>,2> tdirty_(dirty.shape());
@@ -1286,37 +1303,14 @@ template<typename T, typename Serv> void x2dirty(
       double wcur = wmin+iw*dw;
 
       size_t cnt=0;
-      newidx.resize(nvis_plane[iw]);
-      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);
-        grid.fill(0);
-        auto newidx2 = const_mav<uint32_t, 1>(newidx.data(), {newidx.size()});
-        Serv newsrv(srv, newidx2);
-        x2grid_c(gconf, newsrv, grid, 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]});
-        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)
-          {
-          auto ipart = newidx[i];
-          vis_loc[i] = srv.getVis(ipart);
-          idx_loc[i] = srv.getIdx(ipart);
-          }
-        grid.fill(0);
-        vis2grid_c(srv.getBaselines(), gconf, cmav(idx_loc), cmav(vis_loc), grid,
-          nullmav<T,1>(), wcur, dw);
-        }
+      subidx.resize(nvis_plane[iw]);
+      for (size_t ipart=0; ipart<nvis; ++ipart)
+        if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+supp))
+          subidx[cnt++] = ipart;
+      auto subidx2 = const_mav<uint32_t, 1>(subidx.data(), {subidx.size()});
+      auto subsrv = srv.getSubserv(subidx2);
+      grid.fill(0);
+      x2grid_c(gconf, subsrv, grid, wcur, dw);
       gconf.grid2dirty_c(cmav(grid), tdirty);
       gconf.apply_wscreen(cmav(tdirty), tdirty, wcur, true);
 #pragma omp parallel for num_threads(nthreads)
@@ -1379,9 +1373,7 @@ template<typename T, typename Serv> void dirty2x(
     // correct for w gridding
     apply_wcorr(gconf, tdirty, ES_Kernel_with_correction(supp, nthreads), dw);
 
-    vector<complex<T>> vis_loc_;
-    vector<uint32_t> idx_loc_;
-    vector<uint32_t> newidx;
+    vector<uint32_t> subidx;
     tmpStorage<complex<T>,2> grid_({gconf.Nu(),gconf.Nv()});
     auto grid=grid_.getMav();
     tmpStorage<complex<T>,2> tdirty2_({nx_dirty,ny_dirty});
@@ -1401,36 +1393,14 @@ template<typename T, typename Serv> void dirty2x(
       gconf.apply_wscreen(cmav(tdirty2), tdirty2, wcur, false);
       gconf.dirty2grid_c(cmav(tdirty2), grid);
 
-      newidx.resize(nvis_plane[iw]);
+      subidx.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]);
-        }
+      for (size_t ipart=0; ipart<nvis; ++ipart)
+        if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+supp))
+          subidx[cnt++] = ipart;
+      auto subidx2 = const_mav<uint32_t, 1>(subidx.data(), {subidx.size()});
+      auto subsrv = srv.getSubserv(subidx2);
+      grid2x_c(gconf, cmav(grid), subsrv, wcur, dw);
       }
     }
   else