Skip to content

GitLab

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

dirty images really are real-valued :)

parent 8b2ec7ba
Hide whitespace changes
Inline Side-by-side
demo_wstack.py
View file @ a4356ac6
......@@ -13,7 +13,7 @@ def explicit_gridder(uvw, freq, ms, nxdirty, nydirty, xpixsize, ypixsize):
indexing='ij')
x *= xpixsize
y *= ypixsize
res = np.zeros((nxdirty, nydirty))+0j
res = np.zeros((nxdirty, nydirty))
eps = x**2+y**2
nm1 = -eps/(np.sqrt(1.-eps)+1.)
......@@ -22,7 +22,7 @@ def explicit_gridder(uvw, freq, ms, nxdirty, nydirty, xpixsize, ypixsize):
for chan in range(ms.shape[1]):
phase = (freq[chan]/speedoflight *
(x*uvw[row, 0] + y*uvw[row, 1] + uvw[row, 2]*nm1))
res += (ms[row, chan]*np.exp(2j*np.pi*phase))
res += (ms[row, chan]*np.exp(2j*np.pi*phase)).real
return res/n
......@@ -43,8 +43,7 @@ def test_against_wdft(nrow, nchan, nxdirty, nydirty, fov, epsilon, nthreads,
freq = f0 + np.arange(nchan)*(f0/nchan)
uvw = (np.random.rand(nrow, 3)-0.5)/(xpixsize*f0/speedoflight)
ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5)
tdirty = (np.random.rand(nxdirty, nydirty)-0.5
+ 1j*(np.random.rand(nxdirty, nydirty)-0.5))
tdirty = np.random.rand(nxdirty, nydirty)-0.5
single = epsilon > 5e-6
if single:
......@@ -74,7 +73,7 @@ def test_against_wdft(nrow, nchan, nxdirty, nydirty, fov, epsilon, nthreads,
xpixsize, ypixsize, epsilon, nthreads, verbosity=2),
tdirty)
adj2 = np.vdot(ms, degridder(uvw, freq, tdirty, None, xpixsize, ypixsize,
epsilon, nthreads, verbosity=2))
epsilon, nthreads, verbosity=2)).real
print("adjointness test:", np.abs(adj1-adj2)/np.maximum(np.abs(adj1),
np.abs(adj2)))
......
gridder_cxx.h
View file @ a4356ac6
......@@ -1048,7 +1048,7 @@ template<typename T, typename Serv> void wstack_common(
}
template<typename T, typename Serv> void x2dirty_wstack(
const GridderConfig<T> &gconf, const Serv &srv, const mav<complex<T>,2> &dirty, size_t verbosity=0)
const GridderConfig<T> &gconf, const Serv &srv, const mav<T,2> &dirty, size_t verbosity=0)
{
auto nx_dirty=gconf.Nxdirty();
auto ny_dirty=gconf.Nydirty();
......@@ -1110,7 +1110,7 @@ template<typename T, typename Serv> void x2dirty_wstack(
#pragma omp parallel for num_threads(nthreads)
for (size_t i=0; i<nx_dirty; ++i)
for (size_t j=0; j<ny_dirty; ++j)
dirty(i,j) += tdirty(i,j);
dirty(i,j) += tdirty(i,j).real();
}
}
// correct for w gridding
......@@ -1118,13 +1118,13 @@ template<typename T, typename Serv> void x2dirty_wstack(
}
template<typename T> void vis2dirty_wstack(const Baselines<T> &baselines,
const GridderConfig<T> &gconf, const const_mav<uint32_t,1> &idx,
const const_mav<complex<T>,1> &vis, mav<complex<T>,2> &dirty)
const const_mav<complex<T>,1> &vis, mav<T,2> &dirty)
{
x2dirty_wstack(gconf, VisServ<T,const_mav<complex<T>,1>,const_mav<T,1>>(baselines, idx, vis, const_mav<T,1>(nullptr,{0})), dirty);
}
template<typename T, typename Serv> void dirty2x_wstack(
const GridderConfig<T> &gconf, const const_mav<complex<T>,2> &dirty,
const GridderConfig<T> &gconf, const const_mav<T,2> &dirty,
const Serv &srv, size_t verbosity=0)
{
auto nx_dirty=gconf.Nxdirty();
......@@ -1146,7 +1146,7 @@ template<typename T, typename Serv> void dirty2x_wstack(
for (size_t i=0; i<nvis; ++i)
srv.setVis(i, 0.);
tmpStorage<complex<T>,2> tdirty_({nx_dirty,ny_dirty});
tmpStorage<T,2> tdirty_({nx_dirty,ny_dirty});
auto tdirty=tdirty_.getMav();
for (size_t i=0; i<nx_dirty; ++i)
for (size_t j=0; j<ny_dirty; ++j)
......@@ -1181,7 +1181,11 @@ template<typename T, typename Serv> void dirty2x_wstack(
for (size_t i=0; i<nvis_plane[iw]; ++i)
idx_loc[i] = srv.getIdx(mapper[i]);
gconf.apply_wscreen(tdirty, tdirty2, wcur, true);
#pragma omp parallel for num_threads(nthreads)
for (size_t i=0; i<nx_dirty; ++i)
for (size_t j=0; j<ny_dirty; ++j)
tdirty2(i,j) = tdirty(i,j);
gconf.apply_wscreen(tdirty2, tdirty2, wcur, true);
gconf.dirty2grid_c(tdirty2, grid);
grid2vis_c(srv.getBaselines(), gconf, const_mav<uint32_t,1>(idx_loc), const_mav<complex<T>,2>(grid), vis_loc, const_mav<T,1>(nullptr,{0}));
......@@ -1197,7 +1201,7 @@ template<typename T, typename Serv> void dirty2x_wstack(
}
template<typename T> void dirty2vis_wstack(const Baselines<T> &baselines,
const GridderConfig<T> &gconf, const const_mav<uint32_t,1> &idx,
const const_mav<complex<T>,2> &dirty, mav<complex<T>,1> &vis)
const const_mav<T,2> &dirty, mav<complex<T>,1> &vis)
{
dirty2x_wstack(gconf, dirty, VisServ<T,mav<complex<T>,1>,const_mav<T,1>>(baselines, idx, vis, const_mav<T,1>(nullptr,{0})));
}
......@@ -1316,7 +1320,7 @@ template<typename T> vector<uint32_t> getWgtIndices(const Baselines<T> &baseline
template<typename T> void full_gridding(const const_mav<T,2> &uvw,
const const_mav<T,1> &freq, const const_mav<complex<T>,2> &ms,
const const_mav<T,2> &wgt, T pixsize_x, T pixsize_y, double epsilon,
size_t nthreads, const mav<complex<T>,2> &dirty, size_t verbosity)
size_t nthreads, const mav<T,2> &dirty, size_t verbosity)
{
Baselines<T> baselines(uvw, freq);
GridderConfig<T> gconf(dirty.shape(0), dirty.shape(1), epsilon, pixsize_x, pixsize_y, nthreads);
......@@ -1327,7 +1331,7 @@ template<typename T> void full_gridding(const const_mav<T,2> &uvw,
}
template<typename T> void full_degridding(const const_mav<T,2> &uvw,
const const_mav<T,1> &freq, const const_mav<complex<T>,2> &dirty,
const const_mav<T,1> &freq, const const_mav<T,2> &dirty,
const const_mav<T,2> &wgt, T pixsize_x, T pixsize_y, double epsilon,
size_t nthreads, const mav<complex<T>,2> &ms, size_t verbosity)
{
......
nifty_gridder.cc
View file @ a4356ac6
......@@ -610,7 +610,7 @@ template<typename T> pyarr<uint32_t> PygetIndices(const PyBaselines<T> &baseline
return res;
}
template<typename T> pyarr<complex<T>> Pyvis2dirty_wstack(const PyBaselines<T> &baselines,
template<typename T> pyarr<T> Pyvis2dirty_wstack(const PyBaselines<T> &baselines,
const PyGridderConfig<T> &gconf, const pyarr<uint32_t> &idx_,
const pyarr<complex<T>> &vis_)
{
......@@ -618,7 +618,7 @@ template<typename T> pyarr<complex<T>> Pyvis2dirty_wstack(const PyBaselines<T> &
auto ny_dirty=gconf.Nydirty();
auto idx2=make_const_mav<1>(idx_);
auto vis2=make_const_mav<1>(vis_);
auto dirty = makeArray<complex<T>>({nx_dirty, ny_dirty});
auto dirty = makeArray<T>({nx_dirty, ny_dirty});
auto dirty2=make_mav<2>(dirty);
{
py::gil_scoped_release release;
......@@ -629,7 +629,7 @@ template<typename T> pyarr<complex<T>> Pyvis2dirty_wstack(const PyBaselines<T> &
template<typename T> pyarr<complex<T>> Pydirty2vis_wstack(const PyBaselines<T> &baselines,
const PyGridderConfig<T> &gconf, const pyarr<uint32_t> &idx_,
const pyarr<complex<T>> &dirty_)
const pyarr<T> &dirty_)
{
auto idx2=make_const_mav<1>(idx_);
auto nvis = idx2.shape(0);
......@@ -643,7 +643,7 @@ template<typename T> pyarr<complex<T>> Pydirty2vis_wstack(const PyBaselines<T> &
return vis;
}
template<typename T> pyarr<complex<T>> Pyfull_gridding(const pyarr<T> &uvw,
template<typename T> pyarr<T> Pyfull_gridding(const pyarr<T> &uvw,
const pyarr<T> &freq, const pyarr<complex<T>> &ms, const py::object &wgt_,
size_t npix_x, size_t npix_y, T pixsize_x, T pixsize_y, double epsilon,
size_t nthreads, size_t verbosity)
......@@ -653,14 +653,14 @@ template<typename T> pyarr<complex<T>> Pyfull_gridding(const pyarr<T> &uvw,
auto ms2 = make_const_mav<2>(ms);
auto wgt=providePotentialArray<T>(wgt_,{ms2.shape(0),ms2.shape(1)});
auto wgt2 = make_const_mav<2>(wgt);
auto dirty=makeArray<complex<T>>({npix_x,npix_y});
auto dirty=makeArray<T>({npix_x,npix_y});
auto dirty2 =make_mav<2>(dirty);
full_gridding(uvw2,freq2,ms2,wgt2,pixsize_x,pixsize_y,epsilon,nthreads,dirty2,verbosity);
return dirty;
}
template<typename T> pyarr<complex<T>> Pyfull_degridding(const pyarr<T> &uvw,
const pyarr<T> &freq, const pyarr<complex<T>> &dirty, const py::object &wgt_,
const pyarr<T> &freq, const pyarr<T> &dirty, const py::object &wgt_,
T pixsize_x, T pixsize_y, double epsilon, size_t nthreads, size_t verbosity)
{
auto uvw2 = make_const_mav<2>(uvw);
......
test.py
View file @ a4356ac6
......@@ -36,7 +36,7 @@ def _dft(uvw, vis, conf, apply_w):
phase = x*uvw[ii, 0] + y*uvw[ii, 1]
if apply_w:
phase += uvw[ii, 2]*nm1
dirty += (vv*np.exp(2j*np.pi*phase))
dirty += (vv*np.exp(2j*np.pi*phase)).real
if apply_w:
return dirty/n
return dirty
......@@ -76,7 +76,7 @@ def test_adjointness_wgridding(nxdirty, nydirty, nrow, nchan, epsilon):
dirty = np.random.rand(conf.Nxdirty(), conf.Nydirty())-0.5
dirty2 = ng.vis2dirty_wstack(bl, conf, idx, vis)
vis2 = ng.dirty2vis_wstack(bl, conf, idx, dirty)
assert_allclose(np.vdot(vis, vis2), np.vdot(dirty2, dirty),
assert_allclose(np.vdot(vis, vis2).real, np.vdot(dirty2, dirty),
rtol=epsilon)
......@@ -242,7 +242,7 @@ def test_against_dft(nxdirty, nydirty, epsilon, nchan, nrow):
ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5)
vis = bl.ms2vis(ms, idx)
uvw = bl.effectiveuvw(idx)
res0 = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, idx, vis))
res0 = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, idx, vis)).real
res1 = _dft(uvw, vis, conf, False)
assert_(_l2error(res0, res1) < epsilon)
......@@ -281,7 +281,7 @@ def test_wstack_against_wdft(nxdirty, nydirty, nchan, nrow, fov):
if len(iind) == 0:
continue
dd = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, iind, bl.ms2vis(ms, iind)))
res0 += conf.apply_wscreen(dd, 0.5*(ws[ii+1]+ws[ii]), adjoint=False)
res0 += conf.apply_wscreen(dd, 0.5*(ws[ii+1]+ws[ii]), adjoint=False).real
res1 = _dft(uvw, vis, conf, True)
assert_(_l2error(res0, res1) < 1e-4) # Very inaccurate
......@@ -313,7 +313,7 @@ def test_wplane_against_wdft(nxdirty, nydirty, nchan, nrow, fov):
w = np.min(uvw[:, 2])
iind = ng.getIndices(bl, conf, flags)
dd = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, iind, bl.ms2vis(ms, iind)))
res0 = conf.apply_wscreen(dd, w, adjoint=False)
res0 = conf.apply_wscreen(dd, w, adjoint=False).real
res1 = _dft(uvw, vis, conf, True)
assert_(_l2error(res0, res1) < epsilon)
......
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