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

Merge remote-tracking branch 'origin/ducc0' into new_thread_pool

parents 8d3286da ce292ce2
Pipeline #81185 passed with stages
in 17 minutes and 32 seconds
......@@ -607,7 +607,7 @@ template<size_t supp, bool wgrid, typename T> class HelperX2g2
auto iv0old = iv0;
gconf.getpix(in.u, in.v, u, v, iu0, iv0);
T x0 = (iu0-T(u))*2+(supp-1);
T y0 = (iv0-T(v))*2+(supp-1);
T y0 = (iv0-T(v))*2+(supp-1);
if constexpr(wgrid)
krn.eval2s(x0, y0, T(xdw*(w0-in.w)), &buf.simd[0]);
else
......@@ -695,7 +695,7 @@ template<size_t supp, bool wgrid, typename T> class HelperG2x2
auto iv0old = iv0;
gconf.getpix(in.u, in.v, u, v, iu0, iv0);
T x0 = (iu0-T(u))*2+(supp-1);
T y0 = (iv0-T(v))*2+(supp-1);
T y0 = (iv0-T(v))*2+(supp-1);
if constexpr(wgrid)
krn.eval2s(x0, y0, T(xdw*(w0-in.w)), &buf.simd[0]);
else
......@@ -840,7 +840,7 @@ template<bool wgrid, typename T, typename Serv> void x2grid_c
gconf.timers.push("gridding proper");
checkShape(grid.shape(), {gconf.Nu(), gconf.Nv()});
if constexpr (is_same<T, float>::value)
if constexpr (is_same<T, float>::value)
switch(gconf.Supp())
{
case 4: x2grid_c_helper< 4, wgrid>(gconf, srv, grid, w0, dw); break;
......@@ -923,7 +923,7 @@ template<bool wgrid, typename T, typename Serv> void grid2x_c
gconf.timers.push("degridding proper");
checkShape(grid.shape(), {gconf.Nu(), gconf.Nv()});
if constexpr (is_same<T, float>::value)
if constexpr (is_same<T, float>::value)
switch(gconf.Supp())
{
case 4: grid2x_c_helper< 4, wgrid>(gconf, grid, srv, w0, dw); break;
......@@ -1193,7 +1193,7 @@ template<typename T> void report(const GridderConfig<T> &gconf, size_t nvis,
if (x0*x0+y0*y0>1.)
nmin = -sqrt(abs(1.-x0*x0-y0*y0))-1.;
double dw = 0.5/gconf.Ofactor()/abs(nmin);
cout << " w=[" << wmin << "; " << wmax << "], nmin=" << nmin << ", dw=" << dw
cout << " w=[" << wmin << "; " << wmax << "], min(n-1)=" << nmin << ", dw=" << dw
<< ", wmax/dw=" << wmax/dw << endl;
}
......@@ -1442,7 +1442,7 @@ template<typename T> auto scanData(const Baselines &baselines, const mav<complex
});
timers.pop();
return make_tuple(wmin, wmax, nvis, mask_out);
}
}
// Note to self: divide_by_n should always be true when doing Bayesian imaging,
// but wsclean needs it to be false, so this must be kept as a parameter.
......
......@@ -39,7 +39,7 @@ def _l2error(a, b):
def explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, xpixsize, ypixsize,
apply_w):
apply_w, mask):
speedoflight = 299792458.
x, y = np.meshgrid(*[-ss/2 + np.arange(ss) for ss in [nxdirty, nydirty]],
indexing='ij')
......@@ -55,6 +55,8 @@ def explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, xpixsize, ypixsize,
n = 1.
for row in range(ms.shape[0]):
for chan in range(ms.shape[1]):
if mask is not None and mask[row, chan] == 0:
continue
phase = (freq[chan]/speedoflight *
(x*uvw[row, 0] + y*uvw[row, 1] - uvw[row, 2]*nm1))
if wgt is None:
......@@ -74,9 +76,10 @@ def explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, xpixsize, ypixsize,
@pmp("singleprec", (True, False))
@pmp("wstacking", (True, False))
@pmp("use_wgt", (True, False))
@pmp("use_mask", (False, True))
@pmp("nthreads", (1, 2, 7))
def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
singleprec, wstacking, use_wgt, nthreads):
singleprec, wstacking, use_wgt, nthreads, use_mask):
if singleprec and epsilon < 5e-5:
pytest.skip()
rng = np.random.default_rng(42)
......@@ -87,6 +90,7 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
uvw = (rng.random((nrow, 3))-0.5)/(pixsizey*f0/speedoflight)
ms = rng.random((nrow, nchan))-0.5 + 1j*(rng.random((nrow, nchan))-0.5)
wgt = rng.uniform(0.9, 1.1, (nrow, nchan)) if use_wgt else None
mask = (rng.uniform(0, 1, (nrow, nchan)) > 0.5).astype(np.uint8) if use_mask else None
dirty = rng.random((nxdirty, nydirty))-0.5
nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
if nu & 1:
......@@ -101,9 +105,9 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
if wgt is not None:
wgt = wgt.astype("f4")
dirty2 = ng.ms2dirty(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex,
pixsizey, nu, nv, epsilon, wstacking, nthreads, 0).astype("f8")
pixsizey, nu, nv, epsilon, wstacking, nthreads, 0, mask).astype("f8")
ms2 = ng.dirty2ms(uvw, freq, dirty, wgt, pixsizex, pixsizey, nu, nv, epsilon,
wstacking, nthreads+1, 0).astype("c16")
wstacking, nthreads+1, 0, mask).astype("c16")
ref = max(my_vdot(ms,ms).real, my_vdot(ms2,ms2).real,
my_vdot(dirty, dirty).real, my_vdot(dirty2, dirty2).real)
tol = 1e-5*ref if singleprec else 1e-11*ref
......@@ -119,9 +123,10 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
@pmp("singleprec", (False,))
@pmp("wstacking", (False, True))
@pmp("use_wgt", (True,))
@pmp("use_mask", (True,))
@pmp("nthreads", (1, 2, 7))
@pmp("fov", (1., 20.))
def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, fov, nthreads):
def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, use_mask, fov, nthreads):
if singleprec and epsilon < 5e-5:
pytest.skip()
rng = np.random.default_rng(42)
......@@ -132,6 +137,7 @@ def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
uvw = (rng.random((nrow, 3))-0.5)/(pixsizex*f0/speedoflight)
ms = rng.random((nrow, nchan))-0.5 + 1j*(rng.random((nrow, nchan))-0.5)
wgt = rng.uniform(0.9, 1.1, (nrow, 1)) if use_wgt else None
mask = (rng.uniform(0, 1, (nrow, nchan)) > 0.5).astype(np.uint8) if use_mask else None
wgt = np.broadcast_to(wgt, (nrow, nchan)) if use_wgt else None
nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
if nu & 1:
......@@ -145,6 +151,8 @@ def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
if wgt is not None:
wgt = wgt.astype("f4")
dirty = ng.ms2dirty(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)
pixsizey, nu, nv, epsilon, wstacking, nthreads,
0, mask).astype("f8")
ref = explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex,
pixsizey, wstacking, mask)
assert_allclose(_l2error(dirty, ref), 0, atol=epsilon)
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