Commit 03de01dd authored by Martin Reinecke's avatar Martin Reinecke
Browse files

cleanup

parent 1cf25cf5
......@@ -404,6 +404,7 @@ template<typename T> class GridderConfig
size_t nx_dirty, ny_dirty;
T ucorr, vcorr;
size_t w, nsafe, nu, nv;
T beta;
vector<T> cfu, cfv;
public:
......@@ -413,6 +414,7 @@ template<typename T> class GridderConfig
ucorr(1./urange), vcorr(1./vrange),
w(get_w(epsilon)), nsafe((w+1)/2),
nu(max(2*nsafe,2*nx_dirty)), nv(max(2*nsafe,2*ny_dirty)),
beta(2.3*w),
cfu(nx_dirty), cfv(ny_dirty)
{
myassert((nx_dirty&1)==0, "nx_dirty must be even");
......@@ -436,8 +438,7 @@ template<typename T> class GridderConfig
size_t Nv() const { return nv; }
size_t W() const { return w; }
size_t Nsafe() const { return nsafe; }
T Ucorr() const { return ucorr; }
T Vcorr() const { return vcorr; }
T Beta() const { return beta; }
pyarr_c<T> grid2dirty(const pyarr_c<T> &grid) const
{
myassert(grid.ndim()==2, "grid must be a 2D array");
......@@ -528,17 +529,22 @@ template<typename T> class GridderConfig
tmp.data(), tmp.mutable_data(), T(1), 0);
return tmp;
}
inline void getpix(T u_in, T v_in, T &u, T &v, int &iu0, int &iv0) const
{
u=fmodulo(u_in*ucorr, T(1))*nu,
iu0 = int(u-w*0.5 + 1 + nu) - nu;
if (iu0+w>nu+nsafe) iu0 = nu+nsafe-w;
v=fmodulo(v_in*vcorr, T(1))*nv;
iv0 = int(v-w*0.5 + 1 + nv) - nv;
if (iv0+w>nv+nsafe) iv0 = nv+nsafe-w;
}
};
template<typename T> class Helper
{
protected:
int nu, nv;
public:
int w;
T beta;
protected:
int nsafe, su;
const GridderConfig<T> &gconf;
int su;
public:
int sv;
......@@ -548,16 +554,14 @@ template<typename T> class Helper
NOINLINE void update(T u_in, T v_in)
{
auto u = fmodulo(u_in, T(1))*nu;
iu0 = int(u-w*0.5 + 1 + nu) - nu;
if (iu0+w>nu+nsafe) iu0 = nu+nsafe-w;
auto v = fmodulo(v_in, T(1))*nv;
iv0 = int(v-w*0.5 + 1 + nv) - nv;
if (iv0+w>nv+nsafe) iv0 = nv+nsafe-w;
T u, v;
size_t w = gconf.W();
gconf.getpix(u_in, v_in, u, v, iu0, iv0);
T xw=T(2)/w;
auto x0 = xw*(iu0-u);
auto y0 = xw*(iv0-v);
for (int i=0; i<w; ++i)
T beta = gconf.Beta();
for (size_t i=0; i<w; ++i)
{
auto x = x0+i*xw;
kernel[i ] = beta*sqrt(T(1)-x*x);
......@@ -569,35 +573,30 @@ template<typename T> class Helper
}
bool need_to_move() const
{ return (iu0<bu0) || (iv0<bv0) || (iu0+w>bu0+su) || (iv0+w>bv0+sv); }
{
int w = int(gconf.W());
return (iu0<bu0) || (iv0<bv0) || (iu0+w>bu0+su) || (iv0+w>bv0+sv);
}
void update_position()
{
bu0=((((iu0+nsafe)>>logsquare)<<logsquare))-nsafe;
bv0=((((iv0+nsafe)>>logsquare)<<logsquare))-nsafe;
bu0=((((iu0+gconf.Nsafe())>>logsquare)<<logsquare))-gconf.Nsafe();
bv0=((((iv0+gconf.Nsafe())>>logsquare)<<logsquare))-gconf.Nsafe();
}
protected:
Helper(int nu_, int nv_, int w_)
: nu(nu_), nv(nv_), w(w_), beta(2.3*w), nsafe((w+1)/2),
su(2*nsafe+(1<<logsquare)), sv(2*nsafe+(1<<logsquare)),
kernel(2*w),
Helper(const GridderConfig<T> &gconf_)
: gconf(gconf_),
su(2*gconf.Nsafe()+(1<<logsquare)), sv(2*gconf.Nsafe()+(1<<logsquare)),
kernel(2*gconf.W()),
bu0(-1000000), bv0(-1000000)
{
if (min(nu,nv)<2*nsafe) throw runtime_error("grid dimensions too small");
}
{}
};
template<typename T> class WriteHelper: public Helper<T>
{
protected:
using Helper<T>::nu;
using Helper<T>::nv;
public:
using Helper<T>::w;
using Helper<T>::beta;
protected:
using Helper<T>::nsafe;
using Helper<T>::gconf;
using Helper<T>::su;
public:
using Helper<T>::sv;
......@@ -616,7 +615,10 @@ template<typename T> class WriteHelper: public Helper<T>
void dump()
{
if (bu0<-nsafe) return; // nothing written into buffer yet
if (bu0<-int(gconf.Nsafe())) return; // nothing written into buffer yet
auto nu = int(gconf.Nu());
auto nv = int(gconf.Nv());
#pragma omp critical
{
int idxu = (bu0+nu)%nu;
......@@ -636,8 +638,8 @@ template<typename T> class WriteHelper: public Helper<T>
public:
complex<T> *p0;
WriteHelper(int nu_, int nv_, int w, complex<T> *grid_)
: Helper<T>(nu_, nv_, w), data(su*sv, T(0)), grid(grid_) {}
WriteHelper(const GridderConfig<T> &gconf_, complex<T> *grid_)
: Helper<T>(gconf_), data(su*sv, T(0)), grid(grid_) {}
~WriteHelper() { dump(); }
void prep_write(T u_in, T v_in)
......@@ -656,13 +658,7 @@ template<typename T> class WriteHelper: public Helper<T>
template<typename T> class ReadHelper: public Helper<T>
{
protected:
using Helper<T>::nu;
using Helper<T>::nv;
public:
using Helper<T>::w;
using Helper<T>::beta;
protected:
using Helper<T>::nsafe;
using Helper<T>::gconf;
using Helper<T>::su;
public:
using Helper<T>::sv;
......@@ -681,6 +677,8 @@ template<typename T> class ReadHelper: public Helper<T>
void load()
{
auto nu = int(gconf.Nu());
auto nv = int(gconf.Nv());
int idxu = (bu0+nu)%nu;
int idxv0 = (bv0+nv)%nv;
for (int iu=0; iu<su; ++iu)
......@@ -697,8 +695,8 @@ template<typename T> class ReadHelper: public Helper<T>
public:
const complex<T> *p0;
ReadHelper(int nu_, int nv_, int w_, const complex<T> *grid_)
: Helper<T>(nu_, nv_, w_), data(su*sv,T(0)), grid(grid_), p0(nullptr) {}
ReadHelper(const GridderConfig<T> &gconf_, const complex<T> *grid_)
: Helper<T>(gconf_), data(su*sv,T(0)), grid(grid_), p0(nullptr) {}
void prep_read(T u_in, T v_in)
{
......@@ -727,28 +725,28 @@ template<typename T> pyarr_c<T> vis2grid(const Baselines<T> &baselines,
auto res = makearray<complex<T>>({nu, nv});
auto grid = res.mutable_data();
for (size_t i=0; i<nu*nv; ++i) grid[i] = 0.;
T ucorr = gconf.Ucorr(), vcorr=gconf.Vcorr();
T beta = gconf.Beta();
size_t w = gconf.W();
#pragma omp parallel
{
WriteHelper<T> hlp(nu, nv, gconf.W(), grid);
T emb = exp(-2*hlp.beta);
WriteHelper<T> hlp(gconf, grid);
T emb = exp(-2*beta);
const T * RESTRICT ku = hlp.kernel.data();
const T * RESTRICT kv = hlp.kernel.data()+hlp.w;
const T * RESTRICT kv = hlp.kernel.data()+w;
// Loop over sampling points
#pragma omp for schedule(guided,100)
for (size_t ipart=0; ipart<nvis; ++ipart)
{
UVW<T> coord = baselines.effectiveCoord(idx[ipart]);
hlp.prep_write(coord.u*ucorr, coord.v*vcorr);
hlp.prep_write(coord.u, coord.v);
auto * RESTRICT ptr = hlp.p0;
int w = hlp.w;
auto v(vis[ipart]*emb);
for (int cu=0; cu<w; ++cu)
for (size_t cu=0; cu<w; ++cu)
{
complex<T> tmp(v*ku[cu]);
for (int cv=0; cv<w; ++cv)
for (size_t cv=0; cv<w; ++cv)
ptr[cv] += tmp*kv[cv];
ptr+=hlp.sv;
}
......@@ -772,28 +770,28 @@ template<typename T> pyarr_c<complex<T>> vis2grid_c(const Baselines<T> &baseline
auto res = makearray<complex<T>>({nu, nv});
auto grid = res.mutable_data();
for (size_t i=0; i<nu*nv; ++i) grid[i] = 0.;
T ucorr = gconf.Ucorr(), vcorr=gconf.Vcorr();
T beta = gconf.Beta();
size_t w = gconf.W();
#pragma omp parallel
{
WriteHelper<T> hlp(nu, nv, gconf.W(), grid);
T emb = exp(-2*hlp.beta);
WriteHelper<T> hlp(gconf, grid);
T emb = exp(-2*beta);
const T * RESTRICT ku = hlp.kernel.data();
const T * RESTRICT kv = hlp.kernel.data()+hlp.w;
const T * RESTRICT kv = hlp.kernel.data()+w;
// Loop over sampling points
#pragma omp for schedule(guided,100)
for (size_t ipart=0; ipart<nvis; ++ipart)
{
UVW<T> coord = baselines.effectiveCoord(idx[ipart]);
hlp.prep_write(coord.u*ucorr, coord.v*vcorr);
hlp.prep_write(coord.u, coord.v);
auto * RESTRICT ptr = hlp.p0;
int w = hlp.w;
auto v(vis[ipart]*emb);
for (int cu=0; cu<w; ++cu)
for (size_t cu=0; cu<w; ++cu)
{
complex<T> tmp(v*ku[cu]);
for (int cv=0; cv<w; ++cv)
for (size_t cv=0; cv<w; ++cv)
ptr[cv] += tmp*kv[cv];
ptr+=hlp.sv;
}
......@@ -818,28 +816,28 @@ template<typename T> pyarr_c<complex<T>> grid2vis(const Baselines<T> &baselines,
auto res = makearray<complex<T>>({nvis});
auto vis = res.mutable_data();
T ucorr = gconf.Ucorr(), vcorr=gconf.Vcorr();
T beta = gconf.Beta();
size_t w = gconf.W();
// Loop over sampling points
#pragma omp parallel
{
ReadHelper<T> hlp(nu, nv, gconf.W(), grid);
T emb = exp(-2*hlp.beta);
ReadHelper<T> hlp(gconf, grid);
T emb = exp(-2*beta);
const T * RESTRICT ku = hlp.kernel.data();
const T * RESTRICT kv = hlp.kernel.data()+hlp.w;
const T * RESTRICT kv = hlp.kernel.data()+w;
#pragma omp for schedule(guided,100)
for (size_t ipart=0; ipart<nvis; ++ipart)
{
UVW<T> coord = baselines.effectiveCoord(idx[ipart]);
hlp.prep_read(coord.u*ucorr, coord.v*vcorr);
hlp.prep_read(coord.u, coord.v);
complex<T> r = 0;
auto * RESTRICT ptr = hlp.p0;
int w = hlp.w;
for (int cu=0; cu<w; ++cu)
for (size_t cu=0; cu<w; ++cu)
{
complex<T> tmp(0);
for (int cv=0; cv<w; ++cv)
for (size_t cv=0; cv<w; ++cv)
tmp += ptr[cv] * kv[cv];
r += tmp*ku[cu];
ptr += hlp.sv;
......@@ -865,28 +863,28 @@ template<typename T> pyarr_c<complex<T>> grid2vis_c(const Baselines<T> &baseline
auto res = makearray<complex<T>>({nvis});
auto vis = res.mutable_data();
T ucorr = gconf.Ucorr(), vcorr=gconf.Vcorr();
T beta = gconf.Beta();
size_t w = gconf.W();
// Loop over sampling points
#pragma omp parallel
{
ReadHelper<T> hlp(nu, nv, gconf.W(), grid);
T emb = exp(-2*hlp.beta);
ReadHelper<T> hlp(gconf, grid);
T emb = exp(-2*beta);
const T * RESTRICT ku = hlp.kernel.data();
const T * RESTRICT kv = hlp.kernel.data()+hlp.w;
const T * RESTRICT kv = hlp.kernel.data()+w;
#pragma omp for schedule(guided,100)
for (size_t ipart=0; ipart<nvis; ++ipart)
{
UVW<T> coord = baselines.effectiveCoord(idx[ipart]);
hlp.prep_read(coord.u*ucorr, coord.v*vcorr);
hlp.prep_read(coord.u, coord.v);
complex<T> r = 0;
auto * RESTRICT ptr = hlp.p0;
int w = hlp.w;
for (int cu=0; cu<w; ++cu)
for (size_t cu=0; cu<w; ++cu)
{
complex<T> tmp(0);
for (int cv=0; cv<w; ++cv)
for (size_t cv=0; cv<w; ++cv)
tmp += ptr[cv] * kv[cv];
r += tmp*ku[cu];
ptr += hlp.sv;
......@@ -903,12 +901,7 @@ template<typename T> pyarr_c<uint32_t> getIndices(const Baselines<T> &baselines,
{
size_t nrow=baselines.Nrows(),
nchan=baselines.Nchannels(),
nu=gconf.Nu(),
nv=gconf.Nv(),
nsafe=gconf.Nsafe(),
w=gconf.W();
T ucorr=gconf.Ucorr(),
vcorr=gconf.Vcorr();
nsafe=gconf.Nsafe();
if (chbegin<0) chbegin=0;
if (chend<0) chend=nchan;
myassert(flags_.ndim()==2, "flags must be 2D");
......@@ -927,16 +920,11 @@ template<typename T> pyarr_c<uint32_t> getIndices(const Baselines<T> &baselines,
auto uvw = baselines.effectiveCoord(idx);
if ((uvw.w>=wmin) && (uvw.w<wmax))
{
auto u=fmodulo(uvw.u*ucorr, T(1))*nu,
v=fmodulo(uvw.v*vcorr, T(1))*nv;
int iu0 = int(u-w*0.5 + 1 + nu) - nu;
if (iu0+w>nu+nsafe) iu0 = nu+nsafe-w;
iu0+=nsafe;
int iv0 = int(v-w*0.5 + 1 + nv) - nv;
if (iv0+w>nv+nsafe) iv0 = nv+nsafe-w;
iv0+=nsafe;
iu0>>=logsquare;
iv0>>=logsquare;
T u, v;
int iu0, iv0;
gconf.getpix(uvw.u, uvw.v, u, v, iu0, iv0);
iu0 = (iu0+nsafe)>>logsquare;
iv0 = (iv0+nsafe)>>logsquare;
++bincnt[nbv*iu0 + iv0];
}
}
......@@ -954,16 +942,11 @@ template<typename T> pyarr_c<uint32_t> getIndices(const Baselines<T> &baselines,
auto uvw = baselines.effectiveCoord(idx);
if ((uvw.w>=wmin) && (uvw.w<wmax))
{
auto u=fmodulo(uvw.u*ucorr, T(1))*nu,
v=fmodulo(uvw.v*vcorr, T(1))*nv;
int iu0 = int(u-w*0.5 + 1 + nu) - nu;
if (iu0+w>nu+nsafe) iu0 = nu+nsafe-w;
iu0+=nsafe;
int iv0 = int(v-w*0.5 + 1 + nv) - nv;
if (iv0+w>nv+nsafe) iv0 = nv+nsafe-w;
iv0+=nsafe;
iu0>>=logsquare;
iv0>>=logsquare;
T u, v;
int iu0, iv0;
gconf.getpix(uvw.u, uvw.v, u, v, iu0, iv0);
iu0 = (iu0+nsafe)>>logsquare;
iv0 = (iv0+nsafe)>>logsquare;
iout[acc[nbv*iu0 + iv0]] = idx;
++acc[nbv*iu0 + iv0];
}
......
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