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

rename kernel width from w to supp

parent 736bf3b3
......@@ -157,7 +157,7 @@ template<typename T>
template<typename T> pyarr_c<T> makeArray(const vector<size_t> &shape)
{ return pyarr_c<T>(shape); }
size_t get_w(double epsilon)
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,
......@@ -318,12 +318,12 @@ template<typename T> void hartley2_2D(const pyarr_c<T> &in, pyarr_c<T> &out)
/* Compute correction factors for the ES gridding kernel
This implementation follows eqs. (3.8) to (3.10) of Barnett et al. 2018 */
vector<double> correction_factors (size_t n, size_t nval, size_t w)
vector<double> correction_factors (size_t n, size_t nval, size_t supp)
{
constexpr double pi = 3.141592653589793238462643383279502884197;
auto beta = 2.3*w;
auto p = int(1.5*w+2);
double alpha = pi*w/n;
auto beta = 2.3*supp;
auto p = int(1.5*supp+2);
double alpha = pi*supp/n;
vector<double> x, wgt;
legendre_prep(2*p,x,wgt);
auto psi = x;
......@@ -336,7 +336,7 @@ vector<double> correction_factors (size_t n, size_t nval, size_t w)
double tmp=0;
for (int i=0; i<p; ++i)
tmp += wgt[i]*psi[i]*cos(alpha*k*x[i]);
res[k] = 1./(w*tmp);
res[k] = 1./(supp*tmp);
}
return res;
}
......@@ -587,7 +587,7 @@ template<typename T> class GridderConfig
private:
size_t nx_dirty, ny_dirty;
double eps, psx, psy;
size_t w, nsafe, nu, nv;
size_t supp, nsafe, nu, nv;
T beta;
vector<T> cfu, cfv;
......@@ -608,9 +608,9 @@ template<typename T> class GridderConfig
double pixsize_x, double pixsize_y)
: nx_dirty(nxdirty), ny_dirty(nydirty), eps(epsilon),
psx(pixsize_x), psy(pixsize_y),
w(get_w(epsilon)), nsafe((w+1)/2),
supp(get_supp(epsilon)), nsafe((supp+1)/2),
nu(max(2*nsafe,2*nx_dirty)), nv(max(2*nsafe,2*ny_dirty)),
beta(2.3*w),
beta(2.3*supp),
cfu(nx_dirty), cfv(ny_dirty)
{
{
......@@ -621,12 +621,12 @@ template<typename T> class GridderConfig
myassert(pixsize_x>0, "pixsize_x must be positive");
myassert(pixsize_y>0, "pixsize_y must be positive");
auto tmp = correction_factors(nu, nx_dirty/2+1, w);
auto tmp = correction_factors(nu, nx_dirty/2+1, supp);
cfu[nx_dirty/2]=tmp[0];
cfu[0]=tmp[nx_dirty/2];
for (size_t i=1; i<nx_dirty/2; ++i)
cfu[nx_dirty/2-i] = cfu[nx_dirty/2+i] = tmp[i];
tmp = correction_factors(nv, ny_dirty/2+1, w);
tmp = correction_factors(nv, ny_dirty/2+1, supp);
cfv[ny_dirty/2]=tmp[0];
cfv[0]=tmp[ny_dirty/2];
for (size_t i=1; i<ny_dirty/2; ++i)
......@@ -640,7 +640,7 @@ template<typename T> class GridderConfig
double Pixsize_y() const { return psy; }
size_t Nu() const { return nu; }
size_t Nv() const { return nv; }
size_t W() const { return w; }
size_t Supp() const { return supp; }
size_t Nsafe() const { return nsafe; }
T Beta() const { return beta; }
......@@ -761,11 +761,11 @@ template<typename T> class GridderConfig
inline void getpix(T u_in, T v_in, T &u, T &v, int &iu0, int &iv0) const
{
u=fmodulo(u_in*psx, T(1))*nu,
iu0 = int(u-w*0.5 + 1 + nu) - nu;
if (iu0+w>nu+nsafe) iu0 = nu+nsafe-w;
iu0 = int(u-supp*0.5 + 1 + nu) - nu;
if (iu0+supp>nu+nsafe) iu0 = nu+nsafe-supp;
v=fmodulo(v_in*psy, T(1))*nv;
iv0 = int(v-w*0.5 + 1 + nv) - nv;
if (iv0+w>nv+nsafe) iv0 = nv+nsafe-w;
iv0 = int(v-supp*0.5 + 1 + nv) - nv;
if (iv0+supp>nv+nsafe) iv0 = nv+nsafe-supp;
}
pyarr_c<complex<T>> apply_wscreen(const pyarr_c<complex<T>> &dirty_, double w, bool adjoint) const
{
......@@ -810,7 +810,7 @@ template<typename T, typename T2=complex<T>> class Helper
{
private:
const GridderConfig<T> &gconf;
int nu, nv, nsafe, w;
int nu, nv, nsafe, supp;
T beta;
const T2 *grid_r;
T2 *grid_w;
......@@ -864,12 +864,12 @@ template<typename T, typename T2=complex<T>> class Helper
Helper(const GridderConfig<T> &gconf_, const T2 *grid_r_, T2 *grid_w_)
: gconf(gconf_), nu(gconf.Nu()), nv(gconf.Nv()), nsafe(gconf.Nsafe()),
w(gconf.W()), beta(gconf.Beta()), grid_r(grid_r_), grid_w(grid_w_),
su(2*nsafe+(1<<logsquare)), sv(2*nsafe+(1<<logsquare)),
supp(gconf.Supp()), beta(gconf.Beta()), grid_r(grid_r_),
grid_w(grid_w_), su(2*nsafe+(1<<logsquare)), sv(2*nsafe+(1<<logsquare)),
bu0(-1000000), bv0(-1000000),
rbuf(su*sv*(grid_r!=nullptr),T(0)),
wbuf(su*sv*(grid_w!=nullptr),T(0)),
kernel(2*w)
kernel(2*supp)
{}
~Helper() { if (grid_w) dump(); }
......@@ -879,20 +879,20 @@ template<typename T, typename T2=complex<T>> class Helper
{
T u, v;
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 xsupp=T(2)/supp;
auto x0 = xsupp*(iu0-u);
auto y0 = xsupp*(iv0-v);
for (int i=0; i<supp; ++i)
{
auto x = x0+i*xw;
auto x = x0+i*xsupp;
kernel[i ] = beta*sqrt(T(1)-x*x);
auto y = y0+i*xw;
kernel[i+w] = beta*sqrt(T(1)-y*y);
auto y = y0+i*xsupp;
kernel[i+supp] = 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))
if ((iu0<bu0) || (iv0<bv0) || (iu0+supp>bu0+su) || (iv0+supp>bv0+sv))
{
if (grid_w) { dump(); fill(wbuf.begin(), wbuf.end(), T(0)); }
bu0=((((iu0+nsafe)>>logsquare)<<logsquare))-nsafe;
......@@ -948,7 +948,7 @@ template<typename T> pyarr_c<complex<T>> vis2grid_c(
{
py::gil_scoped_release release;
T beta = gconf.Beta();
size_t w = gconf.W();
size_t supp = gconf.Supp();
#pragma omp parallel num_threads(nthreads)
{
......@@ -956,7 +956,7 @@ template<typename T> pyarr_c<complex<T>> vis2grid_c(
T emb = exp(-2*beta);
int jump = hlp.lineJump();
const T * ku = hlp.kernel.data();
const T * kv = hlp.kernel.data()+w;
const T * kv = hlp.kernel.data()+supp;
// Loop over sampling points
#pragma omp for schedule(guided,100)
......@@ -968,10 +968,10 @@ template<typename T> pyarr_c<complex<T>> vis2grid_c(
auto v(vis(ipart)*emb);
if (have_wgt)
v*=wgt(ipart);
for (size_t cu=0; cu<w; ++cu)
for (size_t cu=0; cu<supp; ++cu)
{
complex<T> tmp(v*ku[cu]);
for (size_t cv=0; cv<w; ++cv)
for (size_t cv=0; cv<supp; ++cv)
ptr[cv] += tmp*kv[cv];
ptr+=jump;
}
......@@ -1056,7 +1056,7 @@ template<typename T> pyarr_c<complex<T>> ms2grid_c(
{
py::gil_scoped_release release;
T beta = gconf.Beta();
size_t w = gconf.W();
size_t supp = gconf.Supp();
#pragma omp parallel num_threads(nthreads)
{
......@@ -1064,7 +1064,7 @@ template<typename T> pyarr_c<complex<T>> ms2grid_c(
T emb = exp(-2*beta);
int jump = hlp.lineJump();
const T * ku = hlp.kernel.data();
const T * kv = hlp.kernel.data()+w;
const T * kv = hlp.kernel.data()+supp;
// Loop over sampling points
#pragma omp for schedule(guided,100)
......@@ -1079,10 +1079,10 @@ template<typename T> pyarr_c<complex<T>> ms2grid_c(
auto v(ms(row,chan)*emb);
if (have_wgt)
v*=wgt(row, chan);
for (size_t cu=0; cu<w; ++cu)
for (size_t cu=0; cu<supp; ++cu)
{
complex<T> tmp(v*ku[cu]);
for (size_t cv=0; cv<w; ++cv)
for (size_t cv=0; cv<supp; ++cv)
ptr[cv] += tmp*kv[cv];
ptr+=jump;
}
......@@ -1117,7 +1117,7 @@ template<typename T> pyarr_c<complex<T>> grid2vis_c(
{
py::gil_scoped_release release;
T beta = gconf.Beta();
size_t w = gconf.W();
size_t supp = gconf.Supp();
// Loop over sampling points
#pragma omp parallel num_threads(nthreads)
......@@ -1126,7 +1126,7 @@ template<typename T> pyarr_c<complex<T>> grid2vis_c(
T emb = exp(-2*beta);
int jump = hlp.lineJump();
const T * ku = hlp.kernel.data();
const T * kv = hlp.kernel.data()+w;
const T * kv = hlp.kernel.data()+supp;
#pragma omp for schedule(guided,100)
for (size_t ipart=0; ipart<nvis; ++ipart)
......@@ -1135,10 +1135,10 @@ template<typename T> pyarr_c<complex<T>> grid2vis_c(
hlp.prep(coord.u, coord.v);
complex<T> r = 0;
const auto * ptr = hlp.p0r;
for (size_t cu=0; cu<w; ++cu)
for (size_t cu=0; cu<supp; ++cu)
{
complex<T> tmp(0);
for (size_t cv=0; cv<w; ++cv)
for (size_t cv=0; cv<supp; ++cv)
tmp += ptr[cv] * kv[cv];
r += tmp*ku[cu];
ptr += jump;
......@@ -1199,7 +1199,7 @@ template<typename T> pyarr_c<complex<T>> grid2ms_c(const Baselines<T> &baselines
{
py::gil_scoped_release release;
T beta = gconf.Beta();
size_t w = gconf.W();
size_t supp = gconf.Supp();
// Loop over sampling points
#pragma omp parallel num_threads(nthreads)
......@@ -1208,7 +1208,7 @@ template<typename T> pyarr_c<complex<T>> grid2ms_c(const Baselines<T> &baselines
T emb = exp(-2*beta);
int jump = hlp.lineJump();
const T * ku = hlp.kernel.data();
const T * kv = hlp.kernel.data()+w;
const T * kv = hlp.kernel.data()+supp;
#pragma omp for schedule(guided,100)
for (size_t ipart=0; ipart<nvis; ++ipart)
......@@ -1220,10 +1220,10 @@ template<typename T> pyarr_c<complex<T>> grid2ms_c(const Baselines<T> &baselines
hlp.prep(coord.u, coord.v);
complex<T> r = 0;
const auto * ptr = hlp.p0r;
for (size_t cu=0; cu<w; ++cu)
for (size_t cu=0; cu<supp; ++cu)
{
complex<T> tmp(0);
for (size_t cv=0; cv<w; ++cv)
for (size_t cv=0; cv<supp; ++cv)
tmp += ptr[cv] * kv[cv];
r += tmp*ku[cu];
ptr += jump;
......@@ -1265,7 +1265,7 @@ template<typename T> pyarr_c<complex<T>> apply_holo(
for (size_t i=0; i<nu*nv; ++i) ogrid[i] = T(0);
T beta = gconf.Beta();
size_t w = gconf.W();
size_t supp = gconf.Supp();
// Loop over sampling points
#pragma omp parallel num_threads(nthreads)
......@@ -1274,7 +1274,7 @@ template<typename T> pyarr_c<complex<T>> apply_holo(
T emb = exp(-2*beta);
int jump = hlp.lineJump();
const T * ku = hlp.kernel.data();
const T * kv = hlp.kernel.data()+w;
const T * kv = hlp.kernel.data()+supp;
#pragma omp for schedule(guided,100)
for (size_t ipart=0; ipart<nvis; ++ipart)
......@@ -1284,10 +1284,10 @@ template<typename T> pyarr_c<complex<T>> apply_holo(
hlp.prep(coord.u, coord.v);
complex<T> r = 0;
const auto * ptr = hlp.p0r;
for (size_t cu=0; cu<w; ++cu)
for (size_t cu=0; cu<supp; ++cu)
{
complex<T> tmp(0);
for (size_t cv=0; cv<w; ++cv)
for (size_t cv=0; cv<supp; ++cv)
tmp += ptr[cv] * kv[cv];
r += tmp*ku[cu];
ptr += jump;
......@@ -1299,10 +1299,10 @@ template<typename T> pyarr_c<complex<T>> apply_holo(
r*=twgt*twgt;
}
auto * wptr = hlp.p0w;
for (size_t cu=0; cu<w; ++cu)
for (size_t cu=0; cu<supp; ++cu)
{
complex<T> tmp(r*ku[cu]);
for (size_t cv=0; cv<w; ++cv)
for (size_t cv=0; cv<supp; ++cv)
wptr[cv] += tmp*kv[cv];
wptr += jump;
}
......@@ -1317,9 +1317,9 @@ template<typename T> pyarr_c<T> get_correlations(
const pyarr<uint32_t> &idx_, int du, int dv, const py::object &wgt_)
{
size_t nu=gconf.Nu(), nv=gconf.Nv();
size_t w = gconf.W();
myassert(size_t(abs(du))<w, "|du| must be smaller than W");
myassert(size_t(abs(dv))<w, "|dv| must be smaller than W");
size_t supp = gconf.Supp();
myassert(size_t(abs(du))<supp, "|du| must be smaller than Supp");
myassert(size_t(abs(dv))<supp, "|dv| must be smaller than Supp");
checkArray(idx_, "idx", {0});
size_t nvis = size_t(idx_.shape(0));
auto idx = idx_.template unchecked<1>();
......@@ -1336,13 +1336,13 @@ template<typename T> pyarr_c<T> get_correlations(
size_t u0, u1, v0, v1;
if (du>=0)
{ u0=0; u1=w-du; }
{ u0=0; u1=supp-du; }
else
{ u0=-du; u1=w; }
{ u0=-du; u1=supp; }
if (dv>=0)
{ v0=0; v1=w-dv; }
{ v0=0; v1=supp-dv; }
else
{ v0=-dv; v1=w; }
{ v0=-dv; v1=supp; }
// Loop over sampling points
#pragma omp parallel num_threads(nthreads)
......@@ -1351,7 +1351,7 @@ template<typename T> pyarr_c<T> get_correlations(
T emb = exp(-2*beta);
int jump = hlp.lineJump();
const T * ku = hlp.kernel.data();
const T * kv = hlp.kernel.data()+w;
const T * kv = hlp.kernel.data()+supp;
#pragma omp for schedule(guided,100)
for (size_t ipart=0; ipart<nvis; ++ipart)
......@@ -1488,7 +1488,7 @@ PYBIND11_MODULE(nifty_gridder, m)
.def("Pixsize_y", &GridderConfig<double>::Pixsize_y)
.def("Nu", &GridderConfig<double>::Nu)
.def("Nv", &GridderConfig<double>::Nv)
.def("W", &GridderConfig<double>::W)
.def("Supp", &GridderConfig<double>::Supp)
.def("apply_taper", &GridderConfig<double>::apply_taper, apply_taper_DS,
"img"_a, "divide"_a=false)
.def("grid2dirty", &GridderConfig<double>::grid2dirty,
......
......@@ -220,7 +220,7 @@ def test_applyholo(nxdirty, nydirty, nrow, nchan, epsilon, weight):
def test_correlations(nxdirty, nydirty, nrow, nchan, epsilon, du, dv, weight):
np.random.seed(42)
bl, conf, idx = _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow)
if conf.W() <= abs(du) or conf.W() <= dv:
if conf.Supp() <= abs(du) or conf.Supp() <= dv:
return
wgt = np.random.rand(*idx.shape) if weight else None
y0 = ng.get_correlations(bl, conf, idx, du, dv, wgt)
......
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