Commit ef5aea1d authored by Martin Reinecke's avatar Martin Reinecke

introduce flexible padding

parent a7869f66
......@@ -381,14 +381,17 @@ class ES_Kernel
size_t supp;
public:
ES_Kernel(size_t supp_, size_t nthreads)
: beta(get_beta(supp_)*supp_), p(int(1.5*supp_+2)), supp(supp_)
ES_Kernel(size_t supp_, double ofactor, size_t nthreads)
: beta(get_beta(supp_,ofactor)*supp_), p(int(1.5*supp_+2)), supp(supp_)
{
legendre_prep(2*p,x,wgt,nthreads);
psi=x;
for (auto &v:psi)
v=operator()(v);
}
ES_Kernel(size_t supp_, size_t nthreads)
: ES_Kernel(supp_, 2., nthreads){}
double operator()(double v) const { return exp(beta*(sqrt(1.-v*v)-1.)); }
/* Compute correction factors for the ES gridding kernel
This implementation follows eqs. (3.8) to (3.10) of Barnett et al. 2018 */
......@@ -399,33 +402,59 @@ class ES_Kernel
tmp += wgt[i]*psi[i]*cos(pi*supp*v*x[i]);
return 1./(supp*tmp);
}
static double get_beta(size_t supp)
static double get_beta(size_t supp, double ofactor=2)
{
static const vector<double> opt_beta {-1, 0.14, 1.70, 2.08, 2.205, 2.26,
2.29, 2.307, 2.316, 2.3265, 2.3324, 2.282, 2.294, 2.304, 2.3138, 2.317};
myassert(supp<opt_beta.size(), "bad support size");
return opt_beta[supp];
myassert((supp>=2) && (supp<=15), "unsupported support size");
if (ofactor>=2)
{
static const vector<double> opt_beta {-1, 0.14, 1.70, 2.08, 2.205, 2.26,
2.29, 2.307, 2.316, 2.3265, 2.3324, 2.282, 2.294, 2.304, 2.3138, 2.317};
myassert(supp<opt_beta.size(), "bad support size");
return opt_beta[supp];
}
if (ofactor>=1.2)
{
// empirical, but pretty accurate approximation
static const array<double,16> betacorr{0,0,-0.51,-0.21,-0.1,-0.05,-0.025,-0.0125,0,0,0,0,0,0,0,0};
auto x0 = 1./(2*ofactor);
auto bcstrength=1.+(x0-0.25)*2.5;
return 2.32+bcstrength*betacorr[supp]+(0.25-x0)*3.1;
}
myfail("oversampling factor is too small");
}
static size_t get_supp(double epsilon)
static size_t get_supp(double epsilon, double ofactor=2)
{
static const vector<double> maxmaperr { 1e8, 0.19, 2.98e-3, 5.98e-5,
1.11e-6, 2.01e-8, 3.55e-10, 5.31e-12, 8.81e-14, 1.34e-15, 2.17e-17,
2.12e-19, 2.88e-21, 3.92e-23, 8.21e-25, 7.13e-27 };
double epssq = epsilon*epsilon;
if (ofactor>=2)
{
static const vector<double> maxmaperr { 1e8, 0.19, 2.98e-3, 5.98e-5,
1.11e-6, 2.01e-8, 3.55e-10, 5.31e-12, 8.81e-14, 1.34e-15, 2.17e-17,
2.12e-19, 2.88e-21, 3.92e-23, 8.21e-25, 7.13e-27 };
for (size_t i=1; i<maxmaperr.size(); ++i)
if (epssq>maxmaperr[i]) return i;
myfail("requested epsilon too small - minimum is 1e-13");
for (size_t i=1; i<maxmaperr.size(); ++i)
if (epssq>maxmaperr[i]) return i;
myfail("requested epsilon too small - minimum is 1e-13");
}
if (ofactor>=1.2)
{
for (size_t w=2; w<16; ++w)
{
auto estimate = 12*exp(-2.*w*ofactor); // empirical, not very good approximation
if (epssq>estimate) return w;
}
myfail("requested epsilon too small");
}
myfail("oversampling factor is too small");
}
};
/* 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 supp,
vector<double> correction_factors(size_t n, double ofactor, size_t nval, size_t supp,
size_t nthreads)
{
ES_Kernel kernel(supp, nthreads);
ES_Kernel kernel(supp, ofactor, nthreads);
vector<double> res(nval);
double xn = 1./n;
#pragma omp parallel for schedule(static) num_threads(nthreads)
......@@ -550,9 +579,9 @@ class Baselines
class GridderConfig
{
protected:
size_t nx_dirty, ny_dirty;
double eps, psx, psy;
size_t supp, nsafe, nu, nv;
size_t nx_dirty, ny_dirty, nu, nv;
double ofactor, eps, psx, psy;
size_t supp, nsafe;
double beta;
vector<double> cfu, cfv;
size_t nthreads;
......@@ -572,34 +601,44 @@ class GridderConfig
}
public:
GridderConfig(size_t nxdirty, size_t nydirty, double epsilon,
double pixsize_x, double pixsize_y, size_t nthreads_)
: nx_dirty(nxdirty), ny_dirty(nydirty), eps(epsilon),
GridderConfig(size_t nxdirty, size_t nydirty, size_t nu_, size_t nv_,
double epsilon, double pixsize_x, double pixsize_y, size_t nthreads_)
: nx_dirty(nxdirty), ny_dirty(nydirty), nu(nu_), nv(nv_),
ofactor(min(double(nu)/nxdirty, double(nv)/nydirty)),
eps(epsilon),
psx(pixsize_x), psy(pixsize_y),
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(ES_Kernel::get_beta(supp)*supp),
supp(ES_Kernel::get_supp(epsilon, ofactor)), nsafe((supp+1)/2),
beta(ES_Kernel::get_beta(supp, ofactor)*supp),
cfu(nx_dirty), cfv(ny_dirty), nthreads(nthreads_),
ushift(supp*(-0.5)+1+nu), vshift(supp*(-0.5)+1+nv),
maxiu0((nu+nsafe)-supp), maxiv0((nv+nsafe)-supp)
{
myassert(nu>=2*nsafe, "nu too small");
myassert(nv>=2*nsafe, "nv too small");
myassert((nx_dirty&1)==0, "nx_dirty must be even");
myassert((ny_dirty&1)==0, "ny_dirty must be even");
myassert((nu&1)==0, "nu must be even");
myassert((nv&1)==0, "nv must be even");
myassert(epsilon>0, "epsilon must be positive");
myassert(pixsize_x>0, "pixsize_x must be positive");
myassert(pixsize_y>0, "pixsize_y must be positive");
myassert(ofactor>=1.2, "oversampling factor smaller than 1.2");
auto tmp = correction_factors(nu, nx_dirty/2+1, supp, nthreads);
auto tmp = correction_factors(nu, ofactor, nx_dirty/2+1, supp, nthreads);
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, supp, nthreads);
tmp = correction_factors(nv, ofactor, ny_dirty/2+1, supp, nthreads);
cfv[ny_dirty/2]=tmp[0];
cfv[0]=tmp[ny_dirty/2];
for (size_t i=1; i<ny_dirty/2; ++i)
cfv[ny_dirty/2-i] = cfv[ny_dirty/2+i] = tmp[i];
}
GridderConfig(size_t nxdirty, size_t nydirty,
double epsilon, double pixsize_x, double pixsize_y, size_t nthreads_)
: GridderConfig(nxdirty, nydirty, 2*nxdirty, 2*nydirty, epsilon,
pixsize_x, pixsize_y, nthreads_) {}
size_t Nxdirty() const { return nx_dirty; }
size_t Nydirty() const { return ny_dirty; }
double Epsilon() const { return eps; }
......@@ -611,6 +650,7 @@ class GridderConfig
size_t Nsafe() const { return nsafe; }
double Beta() const { return beta; }
size_t Nthreads() const { return nthreads; }
double Ofactor() const{ return ofactor; }
template<typename T> void apply_taper(const mav<const T,2> &img, mav<T,2> &img2, bool divide) const
{
......@@ -1466,7 +1506,7 @@ template<typename T, typename Serv> void x2dirty(
dirty(i,j) += tdirty(i,j).real();
}
// correct for w gridding
apply_wcorr(gconf, dirty, ES_Kernel(gconf.Supp(), nthreads), dw);
apply_wcorr(gconf, dirty, ES_Kernel(gconf.Supp(), gconf.Ofactor(), nthreads), dw);
}
else
{
......@@ -1511,7 +1551,7 @@ template<typename T, typename Serv> void dirty2x(
for (size_t j=0; j<ny_dirty; ++j)
tdirty(i,j) = dirty(i,j);
// correct for w gridding
apply_wcorr(gconf, tdirty, ES_Kernel(gconf.Supp(), nthreads), dw);
apply_wcorr(gconf, tdirty, ES_Kernel(gconf.Supp(), gconf.Ofactor(), nthreads), dw);
tmpStorage<complex<T>,2> grid_({gconf.Nu(),gconf.Nv()});
auto grid=grid_.getMav();
tmpStorage<complex<T>,2> tdirty2_({nx_dirty, ny_dirty});
......@@ -1736,31 +1776,49 @@ template<typename T> vector<idx_t> getWgtIndices(const Baselines &baselines,
return res;
}
template<typename T> void ms2dirty(const const_mav<double,2> &uvw,
template<typename T> void ms2dirty_general(const const_mav<double,2> &uvw,
const const_mav<double,1> &freq, const const_mav<complex<T>,2> &ms,
const const_mav<T,2> &wgt, double pixsize_x, double pixsize_y, double epsilon,
const const_mav<T,2> &wgt, double pixsize_x, double pixsize_y, size_t nu, size_t nv, double epsilon,
bool do_wstacking, size_t nthreads, const mav<T,2> &dirty, size_t verbosity)
{
Baselines baselines(uvw, freq);
GridderConfig gconf(dirty.shape(0), dirty.shape(1), epsilon, pixsize_x, pixsize_y, nthreads);
GridderConfig gconf(dirty.shape(0), dirty.shape(1), nu, nv, epsilon, pixsize_x, pixsize_y, nthreads);
auto idx = getWgtIndices(baselines, gconf, wgt, ms);
auto idx2 = const_mav<idx_t,1>(idx.data(),{idx.size()});
x2dirty(gconf, makeMsServ(baselines,idx2,ms,wgt), dirty, do_wstacking, verbosity);
}
template<typename T> void ms2dirty(const const_mav<double,2> &uvw,
const const_mav<double,1> &freq, const const_mav<complex<T>,2> &ms,
const const_mav<T,2> &wgt, double pixsize_x, double pixsize_y, double epsilon,
bool do_wstacking, size_t nthreads, const mav<T,2> &dirty, size_t verbosity)
{
ms2dirty_general(uvw, freq, ms, wgt, pixsize_x, pixsize_y,
2*dirty.shape(0), 2*dirty.shape(1), epsilon, do_wstacking, nthreads,
dirty, verbosity);
}
template<typename T> void dirty2ms(const const_mav<double,2> &uvw,
template<typename T> void dirty2ms_general(const const_mav<double,2> &uvw,
const const_mav<double,1> &freq, const const_mav<T,2> &dirty,
const const_mav<T,2> &wgt, double pixsize_x, double pixsize_y, double epsilon,
const const_mav<T,2> &wgt, double pixsize_x, double pixsize_y, size_t nu, size_t nv,double epsilon,
bool do_wstacking, size_t nthreads, const mav<complex<T>,2> &ms, size_t verbosity)
{
Baselines baselines(uvw, freq);
GridderConfig gconf(dirty.shape(0), dirty.shape(1), epsilon, pixsize_x, pixsize_y, nthreads);
GridderConfig gconf(dirty.shape(0), dirty.shape(1), nu, nv, epsilon, pixsize_x, pixsize_y, nthreads);
const_mav<complex<T>,2> null_ms(nullptr, {0,0});
auto idx = getWgtIndices(baselines, gconf, wgt, null_ms);
auto idx2 = const_mav<idx_t,1>(idx.data(),{idx.size()});
ms.fill(0);
dirty2x(gconf, dirty, makeMsServ(baselines,idx2,ms,wgt), do_wstacking, verbosity);
}
template<typename T> void dirty2ms(const const_mav<double,2> &uvw,
const const_mav<double,1> &freq, const const_mav<T,2> &dirty,
const const_mav<T,2> &wgt, double pixsize_x, double pixsize_y, double epsilon,
bool do_wstacking, size_t nthreads, const mav<complex<T>,2> &ms, size_t verbosity)
{
dirty2ms_general(uvw, freq, dirty, wgt, pixsize_x, pixsize_y,
2*dirty.shape(0), 2*dirty.shape(1), epsilon, do_wstacking, nthreads, ms,
verbosity);
}
} // namespace detail
......
......@@ -883,9 +883,9 @@ py::array Pydirty2vis(const PyBaselines &baselines,
myfail("type matching failed: 'dirty' has neither type 'f4' nor 'f8'");
}
template<typename T> py::array ms2dirty2(const py::array &uvw_,
template<typename T> py::array ms2dirty_general2(const py::array &uvw_,
const py::array &freq_, const py::array &ms_, const py::object &wgt_,
size_t npix_x, size_t npix_y, double pixsize_x, double pixsize_y, double epsilon,
size_t npix_x, size_t npix_y, double pixsize_x, double pixsize_y, size_t nu, size_t nv, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
{
auto uvw = getPyarr<double>(uvw_, "uvw");
......@@ -900,12 +900,24 @@ template<typename T> py::array ms2dirty2(const py::array &uvw_,
auto dirty2 = make_mav<2>(dirty);
{
py::gil_scoped_release release;
ms2dirty(uvw2,freq2,ms2,wgt2,pixsize_x,pixsize_y,epsilon,do_wstacking,
ms2dirty_general(uvw2,freq2,ms2,wgt2,pixsize_x,pixsize_y, nu, nv, epsilon,do_wstacking,
nthreads,dirty2,verbosity);
}
return dirty;
}
py::array Pyms2dirty_general(const py::array &uvw,
const py::array &freq, const py::array &ms, const py::object &wgt,
size_t npix_x, size_t npix_y, double pixsize_x, double pixsize_y, size_t nu, size_t nv, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
{
if (isPytype<complex<float>>(ms))
return ms2dirty_general2<float>(uvw, freq, ms, wgt, npix_x, npix_y,
pixsize_x, pixsize_y, nu, nv, epsilon, do_wstacking, nthreads, verbosity);
if (isPytype<complex<double>>(ms))
return ms2dirty_general2<double>(uvw, freq, ms, wgt, npix_x, npix_y,
pixsize_x, pixsize_y, nu, nv, epsilon, do_wstacking, nthreads, verbosity);
myfail("type matching failed: 'ms' has neither type 'c8' nor 'c16'");
}
constexpr auto ms2dirty_DS = R"""(
Converts an MS object to dirty image.
......@@ -948,18 +960,14 @@ py::array Pyms2dirty(const py::array &uvw,
size_t npix_x, size_t npix_y, double pixsize_x, double pixsize_y, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
{
if (isPytype<complex<float>>(ms))
return ms2dirty2<float>(uvw, freq, ms, wgt, npix_x, npix_y,
pixsize_x, pixsize_y, epsilon, do_wstacking, nthreads, verbosity);
if (isPytype<complex<double>>(ms))
return ms2dirty2<double>(uvw, freq, ms, wgt, npix_x, npix_y,
pixsize_x, pixsize_y, epsilon, do_wstacking, nthreads, verbosity);
myfail("type matching failed: 'ms' has neither type 'c8' nor 'c16'");
return Pyms2dirty_general(uvw, freq, ms, wgt, npix_x, npix_y,
pixsize_x, pixsize_y, 2*npix_x, 2*npix_y, epsilon, do_wstacking, nthreads,
verbosity);
}
template<typename T> py::array dirty2ms2(const py::array &uvw_,
template<typename T> py::array dirty2ms_general2(const py::array &uvw_,
const py::array &freq_, const py::array &dirty_, const py::object &wgt_,
double pixsize_x, double pixsize_y, double epsilon,
double pixsize_x, double pixsize_y, size_t nu, size_t nv, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
{
auto uvw = getPyarr<double>(uvw_, "uvw");
......@@ -974,12 +982,24 @@ template<typename T> py::array dirty2ms2(const py::array &uvw_,
auto ms2 = make_mav<2>(ms);
{
py::gil_scoped_release release;
dirty2ms(uvw2,freq2,dirty2,wgt2,pixsize_x,pixsize_y,epsilon,do_wstacking,
dirty2ms_general(uvw2,freq2,dirty2,wgt2,pixsize_x,pixsize_y,nu, nv,epsilon,do_wstacking,
nthreads,ms2,verbosity);
}
return ms;
}
py::array Pydirty2ms_general(const py::array &uvw,
const py::array &freq, const py::array &dirty, const py::object &wgt,
double pixsize_x, double pixsize_y, size_t nu, size_t nv, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
{
if (isPytype<float>(dirty))
return dirty2ms_general2<float>(uvw, freq, dirty, wgt,
pixsize_x, pixsize_y, nu, nv, epsilon, do_wstacking, nthreads, verbosity);
if (isPytype<double>(dirty))
return dirty2ms_general2<double>(uvw, freq, dirty, wgt,
pixsize_x, pixsize_y, nu, nv, epsilon, do_wstacking, nthreads, verbosity);
myfail("type matching failed: 'dirty' has neither type 'f4' nor 'f8'");
}
constexpr auto dirty2ms_DS = R"""(
Converts a dirty image to an MS object.
......@@ -1020,13 +1040,9 @@ py::array Pydirty2ms(const py::array &uvw,
double pixsize_x, double pixsize_y, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
{
if (isPytype<float>(dirty))
return dirty2ms2<float>(uvw, freq, dirty, wgt,
pixsize_x, pixsize_y, epsilon, do_wstacking, nthreads, verbosity);
if (isPytype<double>(dirty))
return dirty2ms2<double>(uvw, freq, dirty, wgt,
pixsize_x, pixsize_y, epsilon, do_wstacking, nthreads, verbosity);
myfail("type matching failed: 'dirty' has neither type 'f4' nor 'f8'");
return Pydirty2ms_general(uvw, freq, dirty, wgt, pixsize_x, pixsize_y,
2*dirty.shape(0), 2*dirty.shape(1), epsilon, do_wstacking, nthreads,
verbosity);
}
} // unnamed namespace
......@@ -1116,7 +1132,13 @@ PYBIND11_MODULE(nifty_gridder, m)
m.def("ms2dirty", &Pyms2dirty, ms2dirty_DS, "uvw"_a, "freq"_a, "ms"_a,
"wgt"_a=None, "npix_x"_a, "npix_y"_a, "pixsize_x"_a, "pixsize_y"_a,
"epsilon"_a, "do_wstacking"_a=false, "nthreads"_a=1, "verbosity"_a=0);
m.def("ms2dirty_general", &Pyms2dirty_general, "uvw"_a, "freq"_a, "ms"_a,
"wgt"_a=None, "npix_x"_a, "npix_y"_a, "pixsize_x"_a, "pixsize_y"_a, "nu"_a, "nv"_a,
"epsilon"_a, "do_wstacking"_a=false, "nthreads"_a=1, "verbosity"_a=0);
m.def("dirty2ms", &Pydirty2ms, dirty2ms_DS, "uvw"_a, "freq"_a, "dirty"_a,
"wgt"_a=None, "pixsize_x"_a, "pixsize_y"_a, "epsilon"_a,
"do_wstacking"_a=false, "nthreads"_a=1, "verbosity"_a=0);
m.def("dirty2ms_general", &Pydirty2ms_general, "uvw"_a, "freq"_a, "dirty"_a,
"wgt"_a=None, "pixsize_x"_a, "pixsize_y"_a, "nu"_a, "nv"_a, "epsilon"_a,
"do_wstacking"_a=false, "nthreads"_a=1, "verbosity"_a=0);
}
......@@ -89,6 +89,46 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, nrow, nchan, epsilon, singleprec
assert_allclose(np.vdot(ms, ms2).real, np.vdot(dirty2, dirty), rtol=tol)
@pmp("nxdirty", (30, 128))
@pmp("nydirty", (128, 250))
@pmp("ofactor", (1.2, 1.5, 1.7, 2.0))
@pmp("nrow", (2, 27))
@pmp("nchan", (1, 5))
@pmp("epsilon", (1e-1, 1e-3, 1e-5))
@pmp("singleprec", (True, False))
@pmp("wstacking", (True, False))
@pmp("use_wgt", (True, False))
@pmp("nthreads", (1, 2))
def test_adjointness_ms2dirty_general(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, nthreads):
if singleprec and epsilon < 5e-5:
return
np.random.seed(42)
pixsizex = np.pi/180/60/nxdirty*0.2398
pixsizey = np.pi/180/60/nxdirty
speedoflight, f0 = 299792458., 1e9
freq = f0 + np.arange(nchan)*(f0/nchan)
uvw = (np.random.rand(nrow, 3)-0.5)/(pixsizey*f0/speedoflight)
ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5)
wgt = np.random.rand(nrow, nchan) if use_wgt else None
dirty = np.random.rand(nxdirty, nydirty)-0.5
nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
if nu&1:
nu+=1
if nv&1:
nv+=1
if singleprec:
ms = ms.astype("c8")
dirty = dirty.astype("f4")
if wgt is not None:
wgt = wgt.astype("f4")
dirty2 = ng.ms2dirty_general(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex,
pixsizey, nu, nv, epsilon, wstacking, nthreads, 0).astype("f8")
ms2 = ng.dirty2ms_general(uvw, freq, dirty, wgt, pixsizex, pixsizey, nu, nv, epsilon,
wstacking, nthreads, 0).astype("c16")
tol = 5e-4 if singleprec else 5e-11
assert_allclose(np.vdot(ms, ms2).real, np.vdot(dirty2, dirty), rtol=tol)
@pmp('nxdirty', [16, 64])
@pmp('nydirty', [64])
@pmp("nrow", (2, 27))
......@@ -120,6 +160,43 @@ def test_ms2dirty_against_wdft(nxdirty, nydirty, nrow, nchan, epsilon, singlepre
assert_allclose(_l2error(dirty, ref), 0, atol=epsilon)
@pmp('nxdirty', [16, 64])
@pmp('nydirty', [64])
@pmp('ofactor', [1.2, 1.4, 1.7, 2])
@pmp("nrow", (2, 27))
@pmp("nchan", (1, 5))
@pmp("epsilon", (1e-2, 1e-3, 1e-4, 1e-7))
@pmp("singleprec", (True,))
@pmp("wstacking", (True,))
@pmp("use_wgt", (False, True))
@pmp("nthreads", (1, 2))
@pmp("fov", (1., 20.))
def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, fov, nthreads):
if singleprec and epsilon < 5e-5:
return
np.random.seed(40)
pixsizex = fov*np.pi/180/nxdirty
pixsizey = fov*np.pi/180/nydirty*1.1
speedoflight, f0 = 299792458., 1e9
freq = f0 + np.arange(nchan)*(f0/nchan)
uvw = (np.random.rand(nrow, 3)-0.5)/(pixsizex*f0/speedoflight)
ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5)
wgt = np.random.rand(nrow, nchan) if use_wgt else None
nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
if nu&1:
nu+=1
if nv&1:
nv+=1
if singleprec:
ms = ms.astype("c8")
if wgt is not None:
wgt = wgt.astype("f4")
dirty = ng.ms2dirty_general(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex,
pixsizey, nu, nv, epsilon, wstacking, nthreads, 0).astype("f8")
ref = explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex, pixsizey, wstacking)
assert_allclose(_l2error(dirty, ref), 0, atol=epsilon)
@pmp("nxdirty", (128, 300))
@pmp("nydirty", (128, 250))
@pmp("nrow", (1, 100))
......
Markdown is supported
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