Skip to content
GitLab
Commit 09b8682a authored by Martin Reinecke's avatar Martin Reinecke
Browse files

polishing

parent 8dbc0f5f
Hide whitespace changes
Inline Side-by-side
gridder_cxx.h
View file @ 09b8682a
......@@ -224,9 +224,6 @@ template<size_t ndim> void checkShape
myassert(shp1[i]==shp2[i], "shape mismatch");
}
/*! Returns the remainder of the division \a v1/v2.
The result is non-negative.
\a v1 can be positive or negative; \a v2 must be positive. */
template<typename T> inline T fmod1 (T v)
{ return v-floor(v); }
......@@ -295,19 +292,6 @@ void legendre_prep(int n, vector<double> &x, vector<double> &w, size_t nthreads)
// Start of real gridder functionality
//
size_t get_supp(double epsilon)
{
static const vector<double> maxmaperr { 1e8, 0.32, 0.021, 6.2e-4,
1.08e-5, 1.25e-7, 8.25e-10, 5.70e-12, 1.22e-13, 2.48e-15, 4.82e-17,
6.74e-19, 5.41e-21, 4.41e-23, 7.88e-25, 3.9e-26 };
double epssq = epsilon*epsilon;
for (size_t i=1; i<maxmaperr.size(); ++i)
if (epssq>maxmaperr[i]) return i;
myfail("requested epsilon too small - minimum is 2e-13");
}
struct RowChan
{
size_t row, chan;
......@@ -352,7 +336,8 @@ template<typename T> void hartley2complex
}
}
template<typename T> void hartley2_2D(const const_mav<T,2> &in, const mav<T,2> &out, size_t nthreads)
template<typename T> void hartley2_2D(const const_mav<T,2> &in,
const mav<T,2> &out, size_t nthreads)
{
myassert(in.shape()==out.shape(), "shape mismatch");
size_t nu=in.shape(0), nv=in.shape(1);
......@@ -379,17 +364,28 @@ template<typename T> void hartley2_2D(const const_mav<T,2> &in, const mav<T,2> &
}
class EC_Kernel
class ES_Kernel
{
protected:
double beta;
public:
EC_Kernel(size_t supp) : beta(2.3*supp) {}
ES_Kernel(size_t supp) : beta(2.3*supp) {}
double operator()(double v) const { return exp(beta*(sqrt(1.-v*v)-1.)); }
static size_t get_supp(double epsilon)
{
static const vector<double> maxmaperr { 1e8, 0.32, 0.021, 6.2e-4,
1.08e-5, 1.25e-7, 8.25e-10, 5.70e-12, 1.22e-13, 2.48e-15, 4.82e-17,
6.74e-19, 5.41e-21, 4.41e-23, 7.88e-25, 3.9e-26 };
double epssq = epsilon*epsilon;
for (size_t i=1; i<maxmaperr.size(); ++i)
if (epssq>maxmaperr[i]) return i;
myfail("requested epsilon too small - minimum is 2e-13");
}
};
class EC_Kernel_with_correction: public EC_Kernel
class ES_Kernel_with_correction: public ES_Kernel
{
protected:
static constexpr double pi = 3.141592653589793238462643383279502884197;
......@@ -398,8 +394,8 @@ class EC_Kernel_with_correction: public EC_Kernel
size_t supp;
public:
EC_Kernel_with_correction(size_t supp_, size_t nthreads)
: EC_Kernel(supp_), p(int(1.5*supp_+2)), supp(supp_)
ES_Kernel_with_correction(size_t supp_, size_t nthreads)
: ES_Kernel(supp_), p(int(1.5*supp_+2)), supp(supp_)
{
legendre_prep(2*p,x,wgt,nthreads);
psi=x;
......@@ -421,7 +417,7 @@ class EC_Kernel_with_correction: public EC_Kernel
vector<double> correction_factors(size_t n, size_t nval, size_t supp,
size_t nthreads)
{
EC_Kernel_with_correction kernel(supp, nthreads);
ES_Kernel_with_correction kernel(supp, nthreads);
vector<double> res(nval);
double xn = 1./n;
#pragma omp parallel for schedule(static) num_threads(nthreads)
......@@ -560,7 +556,7 @@ template<typename T> class GridderConfig
double pixsize_x, double pixsize_y, size_t nthreads_)
: nx_dirty(nxdirty), ny_dirty(nydirty), eps(epsilon),
psx(pixsize_x), psy(pixsize_y),
supp(get_supp(epsilon)), nsafe((supp+1)/2),
supp(ES_Kernel::get_supp(epsilon)), nsafe((supp+1)/2),
nu(max(2*nsafe,2*nx_dirty)), nv(max(2*nsafe,2*ny_dirty)),
beta(2.3*supp),
cfu(nx_dirty), cfv(ny_dirty), nthreads(nthreads_)
......@@ -788,7 +784,8 @@ template<typename T, typename T2=complex<T>> class Helper
T2 *p0w;
T kernel[64] ALIGNED(64);
Helper(const GridderConfig<T> &gconf_, const T2 *grid_r_, T2 *grid_w_, T w0_=-1, T dw_=-1)
Helper(const GridderConfig<T> &gconf_, const T2 *grid_r_, T2 *grid_w_,
T w0_=-1, T dw_=-1)
: gconf(gconf_), nu(gconf.Nu()), nv(gconf.Nv()), nsafe(gconf.Nsafe()),
supp(gconf.Supp()), beta(gconf.Beta()), grid_r(grid_r_),
grid_w(grid_w_), su(2*nsafe+(1<<logsquare)), sv(2*nsafe+(1<<logsquare)),
......@@ -979,8 +976,8 @@ template<typename T> void ms2grid_c
{ x2grid_c(gconf, makeMsServ(baselines, idx, ms, wgt), grid); }
template<typename T, typename Serv> void grid2x_c
(const GridderConfig<T> &gconf, const const_mav<complex<T>,2> &grid, const Serv &srv,
T w0=-1, T dw=-1)
(const GridderConfig<T> &gconf, const const_mav<complex<T>,2> &grid,
const Serv &srv, T w0=-1, T dw=-1)
{
checkShape(grid.shape(), {gconf.Nu(), gconf.Nv()});
myassert(grid.contiguous(), "grid is not contiguous");
......@@ -1165,7 +1162,7 @@ template<typename T> void get_correlations
template<typename T, typename T2> void apply_wcorr(const GridderConfig<T> &gconf,
const mav<T2,2> &dirty, const EC_Kernel_with_correction &kernel, double dw)
const mav<T2,2> &dirty, const ES_Kernel_with_correction &kernel, double dw)
{
auto nx_dirty=gconf.Nxdirty();
auto ny_dirty=gconf.Nydirty();
......@@ -1205,7 +1202,8 @@ template<typename T, typename T2> void apply_wcorr(const GridderConfig<T> &gconf
}
template<typename T, typename Serv> void x2dirty(
const GridderConfig<T> &gconf, const Serv &srv, const mav<T,2> &dirty, size_t verbosity=0)
const GridderConfig<T> &gconf, const Serv &srv, const mav<T,2> &dirty,
size_t verbosity=0)
{
if (verbosity>0)
cout << "Gridding without w-stacking: " << srv.Nvis()
......@@ -1242,7 +1240,8 @@ template<typename T, typename Serv> void dirty2x(
template<typename T, typename Serv> void wstack_common(
const GridderConfig<T> &gconf, const Serv &srv, T &wmin,
double &dw, size_t &nplanes, vector<size_t> &nvis_plane, vector<int> &minplane, size_t verbosity)
double &dw, size_t &nplanes, vector<size_t> &nvis_plane,
vector<int> &minplane, size_t verbosity)
{
size_t nvis = srv.Nvis();
size_t nthreads = gconf.Nthreads();
......@@ -1265,14 +1264,14 @@ template<typename T, typename Serv> void wstack_common(
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 = max<size_t>(2,((wmax-wmin)/dw+2));
dw=(wmax-wmin)/(nplanes-1);
auto w_supp = gconf.Supp();
wmin -= (0.5*w_supp-1)*dw;
wmax += (0.5*w_supp-1)*dw;
nplanes += w_supp-2;
if (verbosity>0) cout << "Kernel support: " << w_supp << endl;
nplanes = size_t((wmax-wmin)/dw+2);
dw = (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;
nvis_plane.resize(nplanes,0);
minplane.resize(nvis);
......@@ -1282,9 +1281,9 @@ template<typename T, typename Serv> void wstack_common(
#pragma omp for
for (size_t ipart=0; ipart<nvis; ++ipart)
{
int plane0 = int((abs(srv.getCoord(ipart).w)-wmin)/dw+0.5*w_supp)-int(w_supp-1);
int plane0 = int((abs(srv.getCoord(ipart).w)-wmin)/dw+0.5*supp)-int(supp-1);
minplane[ipart]=plane0;
for (int i=max<int>(0,plane0); i<min<int>(nplanes,plane0+w_supp); ++i)
for (int i=max<int>(0,plane0); i<min<int>(nplanes,plane0+supp); ++i)
++nvploc[i];
}
#pragma omp critical
......@@ -1294,7 +1293,8 @@ template<typename T, typename Serv> void wstack_common(
}
template<typename T, typename Serv> void x2dirty_wstack(
const GridderConfig<T> &gconf, const Serv &srv, const mav<T,2> &dirty, size_t verbosity=0)
const GridderConfig<T> &gconf, const Serv &srv, const mav<T,2> &dirty,
size_t verbosity=0)
{
auto nx_dirty=gconf.Nxdirty();
auto ny_dirty=gconf.Nydirty();
......@@ -1341,7 +1341,8 @@ template<typename T, typename Serv> void x2dirty_wstack(
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, T(dw));
vis2grid_c(srv.getBaselines(), gconf, cmav(idx_loc), cmav(vis_loc), grid,
nullmav<T,1>(), wcur, T(dw));
gconf.grid2dirty_c(cmav(grid), tdirty);
gconf.apply_wscreen(cmav(tdirty), tdirty, wcur, false);
#pragma omp parallel for num_threads(nthreads)
......@@ -1350,16 +1351,18 @@ template<typename T, typename Serv> void x2dirty_wstack(
dirty(i,j) += tdirty(i,j).real();
}
// correct for w gridding
apply_wcorr(gconf, dirty, EC_Kernel_with_correction(supp, nthreads), dw);
apply_wcorr(gconf, dirty, ES_Kernel_with_correction(supp, nthreads), dw);
}
template<typename T> void vis2dirty_wstack(const Baselines<T> &baselines,
const GridderConfig<T> &gconf, const const_mav<uint32_t,1> &idx,
const const_mav<complex<T>,1> &vis, mav<T,2> &dirty)
{
x2dirty_wstack(gconf, VisServ<T,const_mav<complex<T>,1>,const_mav<T,1>>(baselines, idx, vis, const_mav<T,1>(nullptr,{0})), dirty);
x2dirty_wstack(gconf, VisServ<T,const_mav<complex<T>,1>,
const_mav<T,1>>(baselines, idx, vis, const_mav<T,1>(nullptr,{0})), dirty);
}
template<typename T, typename Serv> void x2dirty_wstack2(
const GridderConfig<T> &gconf, const Serv &srv, const mav<T,2> &dirty, size_t verbosity=0)
const GridderConfig<T> &gconf, const Serv &srv, const mav<T,2> &dirty,
size_t verbosity=0)
{
auto nx_dirty=gconf.Nxdirty();
auto ny_dirty=gconf.Nydirty();
......@@ -1404,13 +1407,13 @@ template<typename T, typename Serv> void x2dirty_wstack2(
dirty(i,j) += tdirty(i,j).real();
}
// correct for w gridding
apply_wcorr(gconf, dirty, EC_Kernel_with_correction(supp, nthreads), dw);
apply_wcorr(gconf, dirty, ES_Kernel_with_correction(supp, nthreads), dw);
}
template<typename T> void vis2dirty_wstack2(const Baselines<T> &baselines,
const GridderConfig<T> &gconf, const const_mav<uint32_t,1> &idx,
const const_mav<complex<T>,1> &vis, mav<T,2> &dirty)
{
x2dirty_wstack2(gconf, VisServ<T,const_mav<complex<T>,1>,const_mav<T,1>>(baselines, idx, vis, const_mav<T,1>(nullptr,{0})), dirty);
x2dirty_wstack2(gconf, makeVisServ(baselines, idx, vis, nullmav<T,1>()), dirty);
}
template<typename T, typename Serv> void dirty2x_wstack(
......@@ -1436,7 +1439,7 @@ template<typename T, typename Serv> void dirty2x_wstack(
for (size_t j=0; j<ny_dirty; ++j)
tdirty(i,j) = dirty(i,j);
// correct for w gridding
apply_wcorr(gconf, tdirty, EC_Kernel_with_correction(supp, nthreads), dw);
apply_wcorr(gconf, tdirty, ES_Kernel_with_correction(supp, nthreads), dw);
vector<complex<T>> vis_loc_;
vector<uint32_t> idx_loc_, mapper;
......@@ -1470,7 +1473,8 @@ template<typename T, typename Serv> void dirty2x_wstack(
tdirty2(i,j) = tdirty(i,j);
gconf.apply_wscreen(cmav(tdirty2), tdirty2, wcur, true);
gconf.dirty2grid_c(cmav(tdirty2), grid);
grid2vis_c(srv.getBaselines(), gconf, cmav(idx_loc), cmav(grid), vis_loc, const_mav<T,1>(nullptr,{0}), wcur, T(dw));
grid2vis_c(srv.getBaselines(), gconf, cmav(idx_loc), cmav(grid), vis_loc,
const_mav<T,1>(nullptr,{0}), wcur, T(dw));
#pragma omp parallel for num_threads(nthreads)
for (size_t i=0; i<nvis_plane[iw]; ++i)
......@@ -1481,12 +1485,13 @@ template<typename T> void dirty2vis_wstack(const Baselines<T> &baselines,
const GridderConfig<T> &gconf, const const_mav<uint32_t,1> &idx,
const const_mav<T,2> &dirty, mav<complex<T>,1> &vis)
{
dirty2x_wstack(gconf, dirty, VisServ<T,mav<complex<T>,1>,const_mav<T,1>>(baselines, idx, vis, const_mav<T,1>(nullptr,{0})));
dirty2x_wstack(gconf, dirty, makeVisServ(baselines, idx, vis, nullmav<T,1>()));
}
template<typename T> size_t getIdxSize(const Baselines<T> &baselines,
const const_mav<bool,2> &flags, int chbegin, int chend, T wmin, T wmax, size_t nthreads)
const const_mav<bool,2> &flags, int chbegin, int chend, T wmin, T wmax,
size_t nthreads)
{
size_t nrow=baselines.Nrows(), nchan=baselines.Nchannels();
if (chbegin<0) chbegin=0;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment