Commit d3a37d82 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

compactification

parent f14b7c1d
......@@ -497,78 +497,20 @@ template<typename T> class GridderConfig
template<typename T> class Helper
{
protected:
private:
const GridderConfig<T> &gconf;
complex<T> *grid;
bool write;
int su;
public:
int sv;
vector<T> kernel;
private:
int iu0, iv0; // start index of the current visibility
int bu0, bv0; // start index of the current buffer
NOINLINE void update(T u_in, T v_in)
{
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);
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);
auto y = y0+i*xw;
kernel[i+w] = beta*sqrt(T(1)-y*y);
}
for (auto &k : kernel)
k = exp(k);
}
bool need_to_move() const
{
int w = int(gconf.W());
return (iu0<bu0) || (iv0<bv0) || (iu0+w>bu0+su) || (iv0+w>bv0+sv);
}
void update_position()
{
bu0=((((iu0+gconf.Nsafe())>>logsquare)<<logsquare))-gconf.Nsafe();
bv0=((((iv0+gconf.Nsafe())>>logsquare)<<logsquare))-gconf.Nsafe();
}
protected:
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)
{}
};
template<typename T> class WriteHelper: public Helper<T>
{
protected:
using Helper<T>::gconf;
using Helper<T>::su;
public:
using Helper<T>::sv;
using Helper<T>::kernel;
using Helper<T>::iu0;
using Helper<T>::iv0;
using Helper<T>::bu0;
using Helper<T>::bv0;
using Helper<T>::need_to_move;
using Helper<T>::update_position;
using Helper<T>::update;
private:
vector<complex<T>> data;
complex<T> *grid;
void dump()
void dump() const
{
if (bu0<-int(gconf.Nsafe())) return; // nothing written into buffer yet
auto nu = int(gconf.Nu());
......@@ -591,45 +533,6 @@ template<typename T> class WriteHelper: public Helper<T>
}
}
public:
complex<T> *p0;
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)
{
update(u_in, v_in);
if (need_to_move())
{
dump();
update_position();
fill(data.begin(), data.end(), T(0));
}
p0 = data.data() + sv*(iu0-bu0) + iv0-bv0;
}
};
template<typename T> class ReadHelper: public Helper<T>
{
protected:
using Helper<T>::gconf;
using Helper<T>::su;
public:
using Helper<T>::sv;
using Helper<T>::kernel;
using Helper<T>::iu0;
using Helper<T>::iv0;
using Helper<T>::bu0;
using Helper<T>::bv0;
using Helper<T>::need_to_move;
using Helper<T>::update_position;
using Helper<T>::update;
private:
vector<complex<T>> data;
const complex<T> *grid;
void load()
{
auto nu = int(gconf.Nu());
......@@ -649,17 +552,44 @@ template<typename T> class ReadHelper: public Helper<T>
}
public:
const complex<T> *p0;
ReadHelper(const GridderConfig<T> &gconf_, const complex<T> *grid_)
: Helper<T>(gconf_), data(su*sv,T(0)), grid(grid_), p0(nullptr) {}
complex<T> *p0;
vector<T> kernel;
void prep_read(T u_in, T v_in)
Helper(const GridderConfig<T> &gconf_, const complex<T> *grid_, bool write_)
: gconf(gconf_), grid(const_cast<complex<T> *>(grid_)), write(write_),
su(2*gconf.Nsafe()+(1<<logsquare)), sv(2*gconf.Nsafe()+(1<<logsquare)),
bu0(-1000000), bv0(-1000000), data(su*sv,T(0)), kernel(2*gconf.W())
{}
~Helper() { if (write) dump(); }
void prep(T u_in, T v_in)
{
update(u_in, v_in);
if (need_to_move())
T u, v;
int 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);
T beta = gconf.Beta();
for (int i=0; i<w; ++i)
{
auto x = x0+i*xw;
kernel[i ] = beta*sqrt(T(1)-x*x);
auto y = y0+i*xw;
kernel[i+w] = beta*sqrt(T(1)-y*y);
}
for (auto &k : kernel)
k = exp(k);
if ((iu0<bu0) || (iv0<bv0) || (iu0+w>bu0+su) || (iv0+w>bv0+sv))
{
update_position();
load();
if (write) dump();
bu0=((((iu0+gconf.Nsafe())>>logsquare)<<logsquare))-gconf.Nsafe();
bv0=((((iv0+gconf.Nsafe())>>logsquare)<<logsquare))-gconf.Nsafe();
if (!write)
load();
else
fill(data.begin(), data.end(), T(0));
}
p0 = data.data() + sv*(iu0-bu0) + iv0-bv0;
}
......@@ -685,7 +615,7 @@ template<typename T> pyarr_c<complex<T>> vis2grid_c(const Baselines<T> &baseline
#pragma omp parallel
{
WriteHelper<T> hlp(gconf, grid);
Helper<T> hlp(gconf, grid, true);
T emb = exp(-2*beta);
const T * RESTRICT ku = hlp.kernel.data();
const T * RESTRICT kv = hlp.kernel.data()+w;
......@@ -695,7 +625,7 @@ template<typename T> pyarr_c<complex<T>> vis2grid_c(const Baselines<T> &baseline
for (size_t ipart=0; ipart<nvis; ++ipart)
{
UVW<T> coord = baselines.effectiveCoord(idx[ipart]);
hlp.prep_write(coord.u, coord.v);
hlp.prep(coord.u, coord.v);
auto * RESTRICT ptr = hlp.p0;
auto v(vis[ipart]*emb);
for (size_t cu=0; cu<w; ++cu)
......@@ -736,7 +666,7 @@ template<typename T> pyarr_c<complex<T>> grid2vis_c(const Baselines<T> &baseline
// Loop over sampling points
#pragma omp parallel
{
ReadHelper<T> hlp(gconf, grid);
Helper<T> hlp(gconf, grid, false);
T emb = exp(-2*beta);
const T * RESTRICT ku = hlp.kernel.data();
const T * RESTRICT kv = hlp.kernel.data()+w;
......@@ -745,7 +675,7 @@ template<typename T> pyarr_c<complex<T>> grid2vis_c(const Baselines<T> &baseline
for (size_t ipart=0; ipart<nvis; ++ipart)
{
UVW<T> coord = baselines.effectiveCoord(idx[ipart]);
hlp.prep_read(coord.u, coord.v);
hlp.prep(coord.u, coord.v);
complex<T> r = 0;
auto * RESTRICT ptr = hlp.p0;
for (size_t cu=0; cu<w; ++cu)
......
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