From 318cec49e4d41064386ba637d16fb6a5ca3dce4b Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Wed, 9 Oct 2019 18:08:23 +0200 Subject: [PATCH] slight restructuring --- gridder_cxx.h | 296 +++++++++++++++++++++++++------------------------- 1 file changed, 149 insertions(+), 147 deletions(-) diff --git a/gridder_cxx.h b/gridder_cxx.h index db5a7bd..17937ec 100644 --- a/gridder_cxx.h +++ b/gridder_cxx.h @@ -907,6 +907,8 @@ template class VisServ bool have_wgt; public: + using Tsub = SubServ; + VisServ(const Baselines &baselines_, const const_mav &idx_, T2 vis_, const const_mav &wgt_) : baselines(baselines_), idx(idx_), vis(vis_), wgt(wgt_), @@ -915,8 +917,8 @@ template class VisServ checkShape(idx.shape(), {nvis}); if (have_wgt) checkShape(wgt.shape(), {nvis}); } - SubServ getSubserv(const const_mav &subidx) const - { return SubServ(*this, subidx); } + Tsub getSubserv(const const_mav &subidx) const + { return Tsub(*this, subidx); } size_t Nvis() const { return nvis; } const Baselines &getBaselines() const { return baselines; } UVW getCoord(size_t i) const @@ -945,6 +947,8 @@ template class MsServ bool have_wgt; public: + using Tsub = SubServ; + MsServ(const Baselines &baselines_, const const_mav &idx_, T2 ms_, const const_mav &wgt_) : baselines(baselines_), idx(idx_), ms(ms_), wgt(wgt_), @@ -953,8 +957,8 @@ template class MsServ checkShape(ms.shape(), {baselines.Nrows(), baselines.Nchannels()}); if (have_wgt) checkShape(wgt.shape(), ms.shape()); } - SubServ getSubserv(const const_mav &subidx) const - { return SubServ(*this, subidx); } + Tsub getSubserv(const const_mav &subidx) const + { return Tsub(*this, subidx); } size_t Nvis() const { return nvis; } const Baselines &getBaselines() const { return baselines; } UVW getCoord(size_t i) const @@ -1272,119 +1276,147 @@ template void apply_wcorr(const GridderConfig &gconf, } } -template void wminmax(const GridderConfig &gconf, - const Serv &srv, double &wmin, double &wmax) +template class WgridHelper { - size_t nvis = srv.Nvis(); - size_t nthreads = gconf.Nthreads(); + private: + const Serv &srv; + double wmin, dw; + size_t nplanes, supp; + vector> minplane; + size_t verbosity; + + int curplane; + vector subidx; + + static void wminmax(const GridderConfig &gconf, + const Serv &srv, double &wmin, double &wmax) + { + size_t nvis = srv.Nvis(); + size_t nthreads = gconf.Nthreads(); - wmin= 1e38; - wmax=-1e38; -// FIXME maybe this can be done more intelligently + wmin= 1e38; + wmax=-1e38; + // FIXME maybe this can be done more intelligently #pragma omp parallel for num_threads(nthreads) reduction(min:wmin) reduction(max:wmax) - for (size_t ipart=0; ipart void update_idx(vector &v, const vector &add, - const vector &del) - { - myassert(v.size()>=del.size(), "must not happen"); - vector res; - res.reserve((v.size()+add.size())-del.size()); - auto iin=v.begin(), ein=v.end(); - auto iadd=add.begin(), eadd=add.end(); - auto irem=del.begin(), erem=del.end(); - - while(iin!=ein) - { - if ((irem!=erem) && (*iin==*irem)) - { ++irem; ++iin; } // skip removed entry - else if ((iadd!=eadd) && (*iadd<*iin)) - res.push_back(*(iadd++)); // add new entry - else - res.push_back(*(iin++)); - } - myassert(irem==erem, "must not happen"); - while(iadd!=eadd) - res.push_back(*(iadd++)); - myassert(res.size()==(v.size()+add.size())-del.size(), "must not happen"); - v.swap(res); - } + template static void update_idx(vector &v, const vector &add, + const vector &del) + { + myassert(v.size()>=del.size(), "must not happen"); + vector res; + res.reserve((v.size()+add.size())-del.size()); + auto iin=v.begin(), ein=v.end(); + auto iadd=add.begin(), eadd=add.end(); + auto irem=del.begin(), erem=del.end(); + + while(iin!=ein) + { + if ((irem!=erem) && (*iin==*irem)) + { ++irem; ++iin; } // skip removed entry + else if ((iadd!=eadd) && (*iadd<*iin)) + res.push_back(*(iadd++)); // add new entry + else + res.push_back(*(iin++)); + } + myassert(irem==erem, "must not happen"); + while(iadd!=eadd) + res.push_back(*(iadd++)); + myassert(res.size()==(v.size()+add.size())-del.size(), "must not happen"); + v.swap(res); + } -template void wstack_common( - const GridderConfig &gconf, const Serv &srv, double &wmin, double &dw, - size_t &nplanes, vector> &minplane, size_t verbosity) - { - size_t nvis = srv.Nvis(); - size_t nthreads = gconf.Nthreads(); - double wmax; + public: + WgridHelper(const GridderConfig &gconf, const Serv &srv_, size_t verbosity_) + : srv(srv_), verbosity(verbosity_), curplane(-1) + { + size_t nvis = srv.Nvis(); + size_t nthreads = gconf.Nthreads(); + double wmax; - wminmax(gconf, srv, wmin, wmax); + wminmax(gconf, srv, wmin, wmax); #ifdef _OPENMP - if (verbosity>0) cout << "Using " << nthreads << " threads" << endl; + if (verbosity>0) cout << "Using " << nthreads << " threads" << endl; #else - if (verbosity>0) cout << "Using single-threaded mode" << endl; + if (verbosity>0) cout << "Using single-threaded mode" << endl; #endif - if (verbosity>0) cout << "W range: " << wmin << " to " << wmax << endl; - double x0 = -0.5*gconf.Nxdirty()*gconf.Pixsize_x(), - y0 = -0.5*gconf.Nydirty()*gconf.Pixsize_y(); - double nmin = sqrt(max(1.-x0*x0-y0*y0,0.))-1.; - dw = 0.25/abs(nmin); - nplanes = size_t((wmax-wmin)/dw+2); - dw = (1.+1e-13)*(wmax-wmin)/(nplanes-1); - - auto supp = gconf.Supp(); - wmin -= (0.5*supp-1)*dw; - wmax += (0.5*supp-1)*dw; - nplanes += supp-2; - if (verbosity>0) cout << "Kernel support: " << supp << endl; - if (verbosity>0) cout << "nplanes: " << nplanes << endl; - - minplane.resize(nplanes); - vector tcnt(my_max_threads()*nplanes,0); + if (verbosity>0) cout << "W range: " << wmin << " to " << wmax << endl; + double x0 = -0.5*gconf.Nxdirty()*gconf.Pixsize_x(), + y0 = -0.5*gconf.Nydirty()*gconf.Pixsize_y(); + double nmin = sqrt(max(1.-x0*x0-y0*y0,0.))-1.; + dw = 0.25/abs(nmin); + nplanes = size_t((wmax-wmin)/dw+2); + dw = (1.+1e-13)*(wmax-wmin)/(nplanes-1); + + supp = gconf.Supp(); + wmin -= (0.5*supp-1)*dw; + wmax += (0.5*supp-1)*dw; + nplanes += supp-2; + if (verbosity>0) cout << "Kernel support: " << supp << endl; + if (verbosity>0) cout << "nplanes: " << nplanes << endl; + + minplane.resize(nplanes); + vector tcnt(my_max_threads()*nplanes,0); #pragma omp parallel num_threads(nthreads) { - vector mytcnt(nplanes,0); - vector nvp(nplanes,0); - auto mythread = my_thread_num(); + vector mytcnt(nplanes,0); + auto mythread = my_thread_num(); #pragma omp for schedule(static) - for (size_t ipart=0; ipart(nplanes,plane0+supp); ++i) - ++nvp[i]; - } + for (size_t ipart=0; ipart myofs(nplanes, 0); - for (size_t j=0; j myofs(nplanes, 0); + for (size_t j=0; j(subidx.data(), {subidx.size()}); + return srv.getSubserv(subidx2); + } + double W() const { return wmin+curplane*dw; } + size_t Nvis() const { return subidx.size(); } + double DW() const { return dw; } + bool advance() + { + if (++curplane>=int(nplanes)) return false; + update_idx(subidx, minplane[curplane], curplane>=int(supp) ? minplane[curplane-supp] : vector()); + if (verbosity>1) + cout << "Working on plane " << curplane << " containing " << subidx.size() + << " visibilities" << endl; + return true; + } + }; template void x2dirty( const GridderConfig &gconf, const Serv &srv, const mav &dirty, @@ -1392,44 +1424,29 @@ template void x2dirty( { if (do_wstacking) { - auto nx_dirty=gconf.Nxdirty(); - auto ny_dirty=gconf.Nydirty(); - auto supp = gconf.Supp(); size_t nthreads = gconf.Nthreads(); - double wmin; - double dw; - size_t nplanes; - vector> minplane; if (verbosity>0) cout << "Gridding using improved w-stacking" << endl; - wstack_common(gconf, srv, wmin, dw, nplanes, minplane, verbosity); + WgridHelper hlp(gconf, srv, verbosity); + double dw = hlp.DW(); dirty.fill(0); - vector subidx; tmpStorage,2> grid_({gconf.Nu(),gconf.Nv()}); auto grid=grid_.getMav(); tmpStorage,2> tdirty_(dirty.shape()); auto tdirty=tdirty_.getMav(); - for (size_t iw=0; iw=supp ? minplane[iw-supp] : vector()); - if (subidx.size()==0) continue; - if (verbosity>1) - cout << "Working on plane " << iw << " containing " << subidx.size() - << " visibilities" << endl; - double wcur = wmin+iw*dw; - - auto subidx2 = const_mav(subidx.data(), {subidx.size()}); - auto subsrv = srv.getSubserv(subidx2); + if (hlp.Nvis()==0) continue; grid.fill(0); - x2grid_c(gconf, subsrv, grid, wcur, dw); + x2grid_c(gconf, hlp.getSubserv(), grid, hlp.W(), dw); gconf.grid2dirty_c(cmav(grid), tdirty); - gconf.apply_wscreen(cmav(tdirty), tdirty, wcur, true); + gconf.apply_wscreen(cmav(tdirty), tdirty, hlp.W(), true); #pragma omp parallel for num_threads(nthreads) - for (size_t i=0; i void dirty2x( { if (do_wstacking) { - auto nx_dirty=gconf.Nxdirty(); - auto ny_dirty=gconf.Nydirty(); - auto supp = gconf.Supp(); + size_t nx_dirty=gconf.Nxdirty(), ny_dirty=gconf.Nydirty(); size_t nthreads = gconf.Nthreads(); - double wmin; - double dw; - size_t nplanes; - vector> minplane; if (verbosity>0) cout << "Degridding using improved w-stacking" << endl; - wstack_common(gconf, srv, wmin, dw, nplanes, minplane, verbosity); - + WgridHelper hlp(gconf, srv, verbosity); + double dw = hlp.DW(); tmpStorage tdirty_({nx_dirty,ny_dirty}); auto tdirty=tdirty_.getMav(); for (size_t i=0; i subidx; + apply_wcorr(gconf, tdirty, ES_Kernel(gconf.Supp(), nthreads), dw); tmpStorage,2> grid_({gconf.Nu(),gconf.Nv()}); auto grid=grid_.getMav(); - tmpStorage,2> tdirty2_({nx_dirty,ny_dirty}); + tmpStorage,2> tdirty2_({nx_dirty, ny_dirty}); auto tdirty2=tdirty2_.getMav(); - for (size_t iw=0; iw=supp ? minplane[iw-supp] : vector()); - if (subidx.size()==0) continue; - if (verbosity>1) - cout << "Working on plane " << iw << " containing " << subidx.size() - << " visibilities" << endl; - double wcur = wmin+iw*dw; + if (hlp.Nvis()==0) continue; #pragma omp parallel for num_threads(nthreads) for (size_t i=0; i(subidx.data(), {subidx.size()}); - auto subsrv = srv.getSubserv(subidx2); - grid2x_c(gconf, cmav(grid), subsrv, wcur, dw); + grid2x_c(gconf, cmav(grid), hlp.getSubserv(), hlp.W(), dw); } } else -- GitLab