diff --git a/gridder_cxx.h b/gridder_cxx.h index 6c4500fc0e82a16ec83b5ab647504e8c2341e801..01a0c0da91db5954a15b95c1cef2e49ddda81a1e 100644 --- a/gridder_cxx.h +++ b/gridder_cxx.h @@ -296,11 +296,6 @@ void legendre_prep(int n, vector &x, vector &w, size_t nthreads) // Start of real gridder functionality // -struct RowChan - { - size_t row, chan; - }; - template void complex2hartley (const const_mav, 2> &grid, const mav &grid2, size_t nthreads) { @@ -431,6 +426,13 @@ vector correction_factors(size_t n, size_t nval, size_t supp, return res; } +using idx_t = uint32_t; + +struct RowChan + { + idx_t row, chan; + }; + struct UVW { double u, v, w; @@ -452,8 +454,8 @@ class Baselines protected: vector coord; vector f_over_c; - size_t nrows, nchan; - size_t shift, mask; + idx_t nrows, nchan; + idx_t shift, mask; public: template Baselines(const const_mav &coord_, @@ -461,12 +463,16 @@ class Baselines { constexpr double speedOfLight = 299792458.; myassert(coord_.shape(1)==3, "dimension mismatch"); + auto hugeval = size_t(~(idx_t(0))); + myassert(coord_.shape(0)>shift, index&mask}; } UVW effectiveCoord(const RowChan &rc) const { return coord[rc.row]*f_over_c[rc.chan]; } - UVW effectiveCoord(uint32_t index) const + UVW effectiveCoord(idx_t index) const { return effectiveCoord(getRowChan(index)); } size_t Nrows() const { return nrows; } size_t Nchannels() const { return nchan; } - uint32_t getIdx(size_t irow, size_t ichan) const + idx_t getIdx(idx_t irow, idx_t ichan) const { return ichan+(irow< void effectiveUVW(const mav &idx, + template void effectiveUVW(const mav &idx, mav &res) const { size_t nvis = idx.shape(0); @@ -505,7 +511,7 @@ class Baselines } template void ms2vis(const mav &ms, - const mav &idx, mav &vis, size_t nthreads) const + const mav &idx, mav &vis, size_t nthreads) const { checkShape(ms.shape(), {nrows, nchan}); size_t nvis = idx.shape(0); @@ -519,7 +525,7 @@ class Baselines } template void vis2ms(const mav &vis, - const mav &idx, mav &ms, size_t nthreads) const + const mav &idx, mav &ms, size_t nthreads) const { size_t nvis = vis.shape(0); checkShape(idx.shape(), {nvis}); @@ -872,10 +878,10 @@ template class SubServ { private: const Serv &srv; - const_mav subidx; + const_mav subidx; public: - SubServ(const Serv &orig, const const_mav &subidx_) + SubServ(const Serv &orig, const const_mav &subidx_) : srv(orig), subidx(subidx_){} size_t Nvis() const { return subidx.size(); } const Baselines &getBaselines() const { return srv.getBaselines(); } @@ -883,7 +889,7 @@ template class SubServ { return srv.getCoord(subidx[i]); } complex getVis(size_t i) const { return srv.getVis(subidx[i]); } - uint32_t getIdx(size_t i) const { return srv.getIdx(subidx[i]); } + idx_t getIdx(size_t i) const { return srv.getIdx(subidx[i]); } void setVis (size_t i, const complex &v) const { srv.setVis(subidx[i], v); } void addVis (size_t i, const complex &v) const @@ -894,7 +900,7 @@ template class VisServ { private: const Baselines &baselines; - const_mav idx; + const_mav idx; T2 vis; const_mav wgt; size_t nvis; @@ -902,14 +908,14 @@ template class VisServ public: VisServ(const Baselines &baselines_, - const const_mav &idx_, T2 vis_, const const_mav &wgt_) + const const_mav &idx_, T2 vis_, const const_mav &wgt_) : baselines(baselines_), idx(idx_), vis(vis_), wgt(wgt_), nvis(vis.shape(0)), have_wgt(wgt.size()!=0) { checkShape(idx.shape(), {nvis}); if (have_wgt) checkShape(wgt.shape(), {nvis}); } - SubServ getSubserv(const const_mav &subidx) const + SubServ getSubserv(const const_mav &subidx) const { return SubServ(*this, subidx); } size_t Nvis() const { return nvis; } const Baselines &getBaselines() const { return baselines; } @@ -917,7 +923,7 @@ template class VisServ { return baselines.effectiveCoord(idx[i]); } complex getVis(size_t i) const { return have_wgt ? vis(i)*wgt(i) : vis(i); } - uint32_t getIdx(size_t i) const { return idx[i]; } + idx_t getIdx(size_t i) const { return idx[i]; } void setVis (size_t i, const complex &v) const { vis(i) = have_wgt ? v*wgt(i) : v; } void addVis (size_t i, const complex &v) const @@ -925,14 +931,14 @@ template class VisServ }; template VisServ makeVisServ (const Baselines &baselines, - const const_mav &idx, const T2 &vis, const const_mav &wgt) + const const_mav &idx, const T2 &vis, const const_mav &wgt) { return VisServ(baselines, idx, vis, wgt); } template class MsServ { private: const Baselines &baselines; - const_mav idx; + const_mav idx; T2 ms; const_mav wgt; size_t nvis; @@ -940,14 +946,14 @@ template class MsServ public: MsServ(const Baselines &baselines_, - const const_mav &idx_, T2 ms_, const const_mav &wgt_) + const const_mav &idx_, T2 ms_, const const_mav &wgt_) : baselines(baselines_), idx(idx_), ms(ms_), wgt(wgt_), nvis(idx.shape(0)), have_wgt(wgt.size()!=0) { checkShape(ms.shape(), {baselines.Nrows(), baselines.Nchannels()}); if (have_wgt) checkShape(wgt.shape(), ms.shape()); } - SubServ getSubserv(const const_mav &subidx) const + SubServ getSubserv(const const_mav &subidx) const { return SubServ(*this, subidx); } size_t Nvis() const { return nvis; } const Baselines &getBaselines() const { return baselines; } @@ -959,7 +965,7 @@ template class MsServ return have_wgt ? ms(rc.row, rc.chan)*wgt(rc.row, rc.chan) : ms(rc.row, rc.chan); } - uint32_t getIdx(size_t i) const { return idx[i]; } + idx_t getIdx(size_t i) const { return idx[i]; } void setVis (size_t i, const complex &v) const { auto rc = baselines.getRowChan(idx(i)); @@ -973,7 +979,7 @@ template class MsServ }; template MsServ makeMsServ (const Baselines &baselines, - const const_mav &idx, const T2 &ms, const const_mav &wgt) + const const_mav &idx, const T2 &ms, const const_mav &wgt) { return MsServ(baselines, idx, ms, wgt); } template void x2grid_c @@ -1027,13 +1033,13 @@ template void x2grid_c template void vis2grid_c (const Baselines &baselines, const GridderConfig &gconf, - const const_mav &idx, const const_mav,1> &vis, + const const_mav &idx, const const_mav,1> &vis, mav,2> &grid, const const_mav &wgt, double w0=-1, double dw=-1) { x2grid_c(gconf, makeVisServ(baselines, idx, vis, wgt), grid, w0, dw); } template void ms2grid_c (const Baselines &baselines, const GridderConfig &gconf, - const const_mav &idx, const const_mav,2> &ms, + const const_mav &idx, const const_mav,2> &ms, mav,2> &grid, const const_mav &wgt) { x2grid_c(gconf, makeMsServ(baselines, idx, ms, wgt), grid); } @@ -1087,19 +1093,19 @@ template void grid2x_c } template void grid2vis_c (const Baselines &baselines, const GridderConfig &gconf, - const const_mav &idx, const const_mav,2> &grid, + const const_mav &idx, const const_mav,2> &grid, mav,1> &vis, const const_mav &wgt, double w0=-1, double dw=-1) { grid2x_c(gconf, grid, makeVisServ(baselines, idx, vis, wgt), w0, dw); } template void grid2ms_c (const Baselines &baselines, const GridderConfig &gconf, - const const_mav &idx, const const_mav,2> &grid, + const const_mav &idx, const const_mav,2> &grid, mav,2> &ms, const const_mav &wgt) { grid2x_c(gconf, grid, makeMsServ(baselines, idx, ms, wgt)); } template void apply_holo (const Baselines &baselines, const GridderConfig &gconf, - const const_mav &idx, const const_mav,2> &grid, + const const_mav &idx, const const_mav,2> &grid, mav,2> &ogrid, const const_mav &wgt) { checkShape(grid.shape(), {gconf.Nu(), gconf.Nv()}); @@ -1169,7 +1175,7 @@ template void apply_holo template void get_correlations (const Baselines &baselines, const GridderConfig &gconf, - const const_mav &idx, int du, int dv, + const const_mav &idx, int du, int dv, mav &ogrid, const const_mav &wgt) { size_t nvis = idx.shape(0); @@ -1314,7 +1320,7 @@ template void update_idx(vector &v, vector &vold, template void wstack_common( const GridderConfig &gconf, const Serv &srv, double &wmin, double &dw, size_t &nplanes, vector &nvis_plane, - vector> &minplane, size_t verbosity) + vector> &minplane, size_t verbosity) { size_t nvis = srv.Nvis(); size_t nthreads = gconf.Nthreads(); @@ -1369,7 +1375,7 @@ template void wstack_common( for (size_t j=0; j void wstack_common( for (size_t ipart=0; ipart void x2dirty( double dw; size_t nplanes; vector nvis_plane; - vector> minplane; + vector> minplane; if (verbosity>0) cout << "Gridding using improved w-stacking" << endl; wstack_common(gconf, srv, wmin, dw, nplanes, nvis_plane, minplane, verbosity); dirty.fill(0); - vector subidx, subidx_old; + vector subidx, subidx_old; tmpStorage,2> grid_({gconf.Nu(),gconf.Nv()}); auto grid=grid_.getMav(); tmpStorage,2> tdirty_(dirty.shape()); @@ -1418,8 +1424,8 @@ template void x2dirty( << " visibilities" << endl; double wcur = wmin+iw*dw; - update_idx(subidx, subidx_old, minplane[iw], iw>=supp ? minplane[iw-supp] : vector()); - auto subidx2 = const_mav(subidx.data(), {subidx.size()}); + update_idx(subidx, subidx_old, minplane[iw], iw>=supp ? minplane[iw-supp] : vector()); + auto subidx2 = const_mav(subidx.data(), {subidx.size()}); auto subsrv = srv.getSubserv(subidx2); grid.fill(0); x2grid_c(gconf, subsrv, grid, wcur, dw); @@ -1451,7 +1457,7 @@ template void x2dirty( } } template void vis2dirty(const Baselines &baselines, - const GridderConfig &gconf, const const_mav &idx, + const GridderConfig &gconf, const const_mav &idx, const const_mav,1> &vis, const const_mav &wgt, mav &dirty, bool do_wstacking, size_t verbosity) { @@ -1473,7 +1479,7 @@ template void dirty2x( double dw; size_t nplanes; vector nvis_plane; - vector> minplane; + vector> minplane; if (verbosity>0) cout << "Degridding using improved w-stacking" << endl; wstack_common(gconf, srv, wmin, dw, nplanes, nvis_plane, minplane, verbosity); @@ -1485,7 +1491,7 @@ template void dirty2x( // correct for w gridding apply_wcorr(gconf, tdirty, ES_Kernel(supp, nthreads), dw); - vector subidx, subidx_old; + vector subidx, subidx_old; tmpStorage,2> grid_({gconf.Nu(),gconf.Nv()}); auto grid=grid_.getMav(); tmpStorage,2> tdirty2_({nx_dirty,ny_dirty}); @@ -1505,8 +1511,8 @@ template void dirty2x( gconf.apply_wscreen(cmav(tdirty2), tdirty2, wcur, false); gconf.dirty2grid_c(cmav(tdirty2), grid); - update_idx(subidx, subidx_old, minplane[iw], iw>=supp ? minplane[iw-supp] : vector()); - auto subidx2 = const_mav(subidx.data(), {subidx.size()}); + update_idx(subidx, subidx_old, minplane[iw], iw>=supp ? minplane[iw-supp] : vector()); + auto subidx2 = const_mav(subidx.data(), {subidx.size()}); auto subsrv = srv.getSubserv(subidx2); grid2x_c(gconf, cmav(grid), subsrv, wcur, dw); } @@ -1528,7 +1534,7 @@ template void dirty2x( } } template void dirty2vis(const Baselines &baselines, - const GridderConfig &gconf, const const_mav &idx, + const GridderConfig &gconf, const const_mav &idx, const const_mav &dirty, const const_mav &wgt, mav,1> &vis, bool do_wstacking, size_t verbosity) { @@ -1550,11 +1556,11 @@ size_t getIdxSize(const Baselines &baselines, checkShape(flags.shape(), {nrow, nchan}); size_t res=0; #pragma omp parallel for num_threads(nthreads) reduction(+:res) - for (size_t irow=0; irow=wmin) && (w &flags, int chbegin, - int chend, double wmin, double wmax, mav &res) + int chend, double wmin, double wmax, mav &res) { size_t nrow=baselines.Nrows(), nchan=baselines.Nchannels(), @@ -1577,13 +1583,13 @@ void fillIdx(const Baselines &baselines, constexpr int side=1<> logsquare, nbv = (gconf.Nv()+1+side-1) >> logsquare; - vector acc(nbu*nbv+1, 0); - vector tmp(nrow*(chend-chbegin)); - for (size_t irow=0, idx=0; irow acc(nbu*nbv+1, 0); + vector tmp(nrow*(chend-chbegin)); + for (idx_t irow=0, idx=0; irow=wmin) && (uvw.w=wmin) && (w vector getWgtIndices(const Baselines &baselines, +template vector getWgtIndices(const Baselines &baselines, const GridderConfig &gconf, const const_mav &wgt) { size_t nrow=baselines.Nrows(), @@ -1621,14 +1627,14 @@ template vector getWgtIndices(const Baselines &baselines, constexpr int side=1<> logsquare, nbv = (gconf.Nv()+1+side-1) >> logsquare; - vector acc(nbu*nbv+1, 0); - vector tmp(nrow*nchan); + vector acc(nbu*nbv+1, 0); + vector tmp(nrow*nchan); - for (size_t irow=0, idx=0; irow vector getWgtIndices(const Baselines &baselines, for (size_t i=1; i res(acc.back()); + vector res(acc.back()); for (size_t irow=0, idx=0; irow void ms2dirty(const const_mav &uvw, Baselines baselines(uvw, freq); GridderConfig gconf(dirty.shape(0), dirty.shape(1), epsilon, pixsize_x, pixsize_y, nthreads); auto idx = getWgtIndices(baselines, gconf, wgt); - auto idx2 = const_mav(idx.data(),{idx.size()}); + auto idx2 = const_mav(idx.data(),{idx.size()}); x2dirty(gconf, makeMsServ(baselines,idx2,ms,wgt), dirty, do_wstacking, verbosity); } @@ -1670,7 +1676,7 @@ template void dirty2ms(const const_mav &uvw, Baselines baselines(uvw, freq); GridderConfig gconf(dirty.shape(0), dirty.shape(1), epsilon, pixsize_x, pixsize_y, nthreads); auto idx = getWgtIndices(baselines, gconf, wgt); - auto idx2 = const_mav(idx.data(),{idx.size()}); + auto idx2 = const_mav(idx.data(),{idx.size()}); ms.fill(0); dirty2x(gconf, dirty, makeMsServ(baselines,idx2,ms,wgt), do_wstacking, verbosity); }