Commit 11a10e5f authored by Martin Reinecke's avatar Martin Reinecke
Browse files

better names

parent f3d8625e
......@@ -21,7 +21,7 @@ def explicit_gridder(uvw, freq, ms, nxdirty, nydirty, xpixsize, ypixsize):
for row in range(ms.shape[0]):
for chan in range(ms.shape[1]):
phase = (freq[chan]/speedoflight *
(x*uvw[row, 0] + y*uvw[row, 1] + uvw[row, 2]*nm1))
(x*uvw[row, 0] + y*uvw[row, 1] - uvw[row, 2]*nm1))
res += (ms[row, chan]*np.exp(2j*np.pi*phase)).real
return res/n
......@@ -48,32 +48,29 @@ def test_against_wdft(nrow, nchan, nxdirty, nydirty, fov, epsilon, nthreads,
single = epsilon > 5e-6
if single:
print("\nCalling single-precision functions")
uvw = uvw.astype("f4")
freq = freq.astype("f4")
ms = ms.astype("c8")
tdirty = tdirty.astype("c8")
gridder, degridder = ng.full_gridding_f, ng.full_degridding_f
else:
print("\nCalling double-precision functions")
gridder, degridder = ng.full_gridding, ng.full_degridding
if test_against_explicit:
print("\nTesting against explicit transform "
"(potentially VERY slow!)...")
truth = explicit_gridder(uvw, freq, ms, nxdirty, nydirty,
xpixsize, ypixsize)
res = gridder(uvw, freq, ms, None, nxdirty, nydirty, xpixsize,
ypixsize, epsilon, nthreads)
res = ng.ms2dirty(uvw, freq, ms, None, nxdirty, nydirty, xpixsize,
ypixsize, epsilon, nthreads)
print("L2 error between explicit transform and gridder:",
_l2error(truth, res))
# test adjointness
print("\nTesting adjointness of the gridding/degridding operation")
adj1 = np.vdot(gridder(uvw, freq, ms, None, nxdirty, nydirty,
xpixsize, ypixsize, epsilon, nthreads, verbosity=2),
adj1 = np.vdot(ng.ms2dirty(uvw, freq, ms, None, nxdirty, nydirty,
xpixsize, ypixsize, epsilon, nthreads,
verbosity=2),
tdirty)
adj2 = np.vdot(ms, degridder(uvw, freq, tdirty, None, xpixsize, ypixsize,
epsilon, nthreads, verbosity=2)).real
adj2 = np.vdot(ms, ng.dirty2ms(uvw, freq, tdirty, None, xpixsize, ypixsize,
epsilon, nthreads, verbosity=2)).real
print("adjointness test:", np.abs(adj1-adj2)/np.maximum(np.abs(adj1),
np.abs(adj2)))
......
......@@ -26,11 +26,11 @@ from casacore.tables import table
# - Only one spectral window
# - Flag both LL and RR if one is flagged
ms = 'supernovashell.55.7+3.4.spw0.ms'
t = table(ms, readonly=True)
name = 'supernovashell.55.7+3.4.spw0.ms'
t = table(name, readonly=True)
uvw = t.getcol("UVW") # [:, uvw]
vis = np.array(t.getcol("DATA"), dtype=np.complex128) # [:, ch, corr]
wgt = t.getcol("WEIGHT")
ms = np.array(t.getcol("DATA"), dtype=np.complex128) # [:, ch, corr]
wgt = t.getcol("WEIGHT").astype("f8")
# Flag if one correlation is flagged
flags = np.any(np.array(t.getcol('FLAG'), np.bool), axis=2) # [:, ch]
if len(set(t.getcol('FIELD_ID'))) != 1:
......@@ -39,24 +39,24 @@ if len(set(t.getcol('DATA_DESC_ID'))) != 1:
raise RuntimeError
t.close()
print('# Rows: {}'.format(vis.shape[0]))
print('# Channels: {}'.format(vis.shape[1]))
print('# Correlations: {}'.format(vis.shape[2]))
print('# Rows: {}'.format(ms.shape[0]))
print('# Channels: {}'.format(ms.shape[1]))
print('# Correlations: {}'.format(ms.shape[2]))
print("{} % flagged".format(np.sum(flags)/flags.size*100))
t = table(join(ms, 'SPECTRAL_WINDOW'), readonly=True)
t = table(join(name, 'SPECTRAL_WINDOW'), readonly=True)
freq = t.getcol('CHAN_FREQ')[0]
t.close()
# Select either RR+LL or XX+YY
t = table(join(ms, 'POLARIZATION'), readonly=True)
t = table(join(name, 'POLARIZATION'), readonly=True)
pol = list(t.getcol('CORR_TYPE')[0])
t.close()
if set(pol) <= set([5, 6, 7, 8]):
ind = [pol.index(5), pol.index(8)]
else:
ind = [pol.index(9), pol.index(12)]
vis = np.sum(vis[:, :, ind], axis=2)
ms = np.sum(ms[:, :, ind], axis=2)
wgt = 1/np.sum(1/wgt, axis=1)
# WEIGHT -> WEIGHT_SPECTRUM
......@@ -72,15 +72,14 @@ nthreads = 4
epsilon = 6e-6
t0 = time()
print('Start gridding...')
if (epsilon > 5e-6):
dirty = ng.full_gridding_f(
uvw.astype("f4"), freq.astype("f4"), vis.astype("c8"),
wgt.astype("f4"), npixdirty, npixdirty, pixsize, pixsize, epsilon,
nthreads, verbosity=2)
else:
dirty = ng.full_gridding(
uvw, freq, vis, wgt, npixdirty, npixdirty, pixsize, pixsize, epsilon,
nthreads, verbosity=2)
if epsilon > 5e-6:
ms = ms.astype("c8")
wgt = wgt.astype("f4")
dirty = ng.ms2dirty(
uvw, freq, ms, wgt, npixdirty, npixdirty, pixsize, pixsize, epsilon,
do_wstacking=True, nthreads=nthreads, verbosity=2)
print('Done')
t = time() - t0
print("{} s".format(t))
......
......@@ -1565,7 +1565,7 @@ template<typename T> vector<uint32_t> getWgtIndices(const Baselines &baselines,
return res;
}
template<typename T> void gridding(const const_mav<double,2> &uvw,
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)
......@@ -1577,7 +1577,7 @@ template<typename T> void gridding(const const_mav<double,2> &uvw,
x2dirty(gconf, makeMsServ(baselines,idx2,ms,wgt), dirty, do_wstacking, verbosity);
}
template<typename T> void degridding(const const_mav<double,2> &uvw,
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)
......@@ -1603,8 +1603,8 @@ using detail::vis2grid_c;
using detail::grid2vis_c;
using detail::vis2dirty;
using detail::dirty2vis;
using detail::gridding;
using detail::degridding;
using detail::ms2dirty;
using detail::dirty2ms;
using detail::streamDump__;
using detail::CodeLocation;
} // namespace gridder
......
......@@ -769,7 +769,7 @@ py::array Pydirty2vis(const PyBaselines &baselines,
myfail("type matching failed: 'dirty' has neither type 'f4' nor 'f8'");
}
template<typename T> py::array gridding2(const py::array &uvw_,
template<typename T> py::array ms2dirty2(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,
bool do_wstacking, size_t nthreads, size_t verbosity)
......@@ -786,27 +786,27 @@ template<typename T> py::array gridding2(const py::array &uvw_,
auto dirty2 = make_mav<2>(dirty);
{
py::gil_scoped_release release;
gridding(uvw2,freq2,ms2,wgt2,pixsize_x,pixsize_y,epsilon,do_wstacking,
ms2dirty(uvw2,freq2,ms2,wgt2,pixsize_x,pixsize_y,epsilon,do_wstacking,
nthreads,dirty2,verbosity);
}
return dirty;
}
py::array Pygridding(const py::array &uvw,
py::array Pyms2dirty(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,
bool do_wstacking, size_t nthreads, size_t verbosity)
{
if (isPytype<complex<float>>(ms))
return gridding2<float>(uvw, freq, ms, wgt, npix_x, npix_y,
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 gridding2<double>(uvw, freq, ms, wgt, npix_x, npix_y,
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'");
}
template<typename T> py::array degridding2(const py::array &uvw_,
template<typename T> py::array dirty2ms2(const py::array &uvw_,
const py::array &freq_, const py::array &dirty_, const py::object &wgt_,
double pixsize_x, double pixsize_y, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
......@@ -823,22 +823,22 @@ template<typename T> py::array degridding2(const py::array &uvw_,
auto ms2 = make_mav<2>(ms);
{
py::gil_scoped_release release;
degridding(uvw2,freq2,dirty2,wgt2,pixsize_x,pixsize_y,epsilon,do_wstacking,
dirty2ms(uvw2,freq2,dirty2,wgt2,pixsize_x,pixsize_y,epsilon,do_wstacking,
nthreads,ms2,verbosity);
}
return ms;
}
py::array Pydegridding(const py::array &uvw,
py::array Pydirty2ms(const py::array &uvw,
const py::array &freq, const py::array &dirty, const py::object &wgt,
double pixsize_x, double pixsize_y, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
{
if (isPytype<float>(dirty))
return degridding2<float>(uvw, freq, dirty, wgt,
return dirty2ms2<float>(uvw, freq, dirty, wgt,
pixsize_x, pixsize_y, epsilon, do_wstacking, nthreads, verbosity);
if (isPytype<double>(dirty))
return degridding2<double>(uvw, freq, dirty, wgt,
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'");
}
......@@ -927,10 +927,10 @@ PYBIND11_MODULE(nifty_gridder, m)
m.def("dirty2vis",&Pydirty2vis, "baselines"_a, "gconf"_a,
"idx"_a, "dirty"_a, "do_wstacking"_a=false);
m.def("gridding",&Pygridding,"uvw"_a,"freq"_a,"ms"_a,
m.def("ms2dirty",&Pyms2dirty,"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("degridding",&Pydegridding,"uvw"_a,"freq"_a,"dirty"_a,
m.def("dirty2ms",&Pydirty2ms,"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);
}
......@@ -103,10 +103,10 @@ def test_adjointness_wgridding_highlevel(nxdirty, nydirty, nrow, nchan, epsilon,
if singleprec:
vis = vis.astype("c8")
dirty = dirty.astype("f4")
dirty2 = ng.gridding(uvw, freq, vis, None, nxdirty, nydirty, pixsizex,
dirty2 = ng.ms2dirty(uvw, freq, vis, None, nxdirty, nydirty, pixsizex,
pixsizey, epsilon, wstacking, 1, 0).astype("f8")
vis2 = ng.degridding(uvw, freq, dirty, None, pixsizex, pixsizey, epsilon,
wstacking, 1, 0).astype("c16")
vis2 = ng.dirty2ms(uvw, freq, dirty, None, pixsizex, pixsizey, epsilon,
wstacking, 1, 0).astype("c16")
tol = 5e-5 if singleprec else 5e-13
assert_allclose(np.vdot(vis, vis2).real, np.vdot(dirty2, dirty),
rtol=tol)
......
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