Commit 06c2054d authored by Martin Reinecke's avatar Martin Reinecke
Browse files

fix w plane spacing and positioning; adjust tests

parent 000144b1
...@@ -264,7 +264,7 @@ template<typename T> class EC_Kernel_with_correction: public EC_Kernel<T> ...@@ -264,7 +264,7 @@ template<typename T> class EC_Kernel_with_correction: public EC_Kernel<T>
protected: protected:
static constexpr T pi = T(3.141592653589793238462643383279502884197L); static constexpr T pi = T(3.141592653589793238462643383279502884197L);
int p; int p;
vector<T> x, wgt, psi; vector<double> x, wgt, psi;
using EC_Kernel<T>::supp; using EC_Kernel<T>::supp;
public: public:
...@@ -396,7 +396,7 @@ template<typename T> class GridderConfig ...@@ -396,7 +396,7 @@ template<typename T> class GridderConfig
{ {
protected: protected:
size_t nx_dirty, ny_dirty; size_t nx_dirty, ny_dirty;
double eps, psx, psy; T eps, psx, psy;
size_t supp, nsafe, nu, nv; size_t supp, nsafe, nu, nv;
T beta; T beta;
vector<T> cfu, cfv; vector<T> cfu, cfv;
...@@ -872,8 +872,7 @@ template<typename T> void apply_holo ...@@ -872,8 +872,7 @@ template<typename T> void apply_holo
#pragma omp for schedule(guided,100) #pragma omp for schedule(guided,100)
for (size_t ipart=0; ipart<nvis; ++ipart) for (size_t ipart=0; ipart<nvis; ++ipart)
{ {
auto tidx = idx(ipart); UVW<T> coord = baselines.effectiveCoord(idx(ipart));
UVW<T> coord = baselines.effectiveCoord(tidx);
hlp.prep(coord.u, coord.v); hlp.prep(coord.u, coord.v);
complex<T> r = 0; complex<T> r = 0;
const auto * ptr = hlp.p0r; const auto * ptr = hlp.p0r;
...@@ -1028,12 +1027,14 @@ cout << "data w range: " << wmin << " to " << wmax << endl; ...@@ -1028,12 +1027,14 @@ cout << "data w range: " << wmin << " to " << wmax << endl;
y0 = -0.5*ny_dirty*psy; y0 = -0.5*ny_dirty*psy;
double nmin = sqrt(max(1.-x0*x0-y0*y0,0.))-1.; double nmin = sqrt(max(1.-x0*x0-y0*y0,0.))-1.;
double dw = 0.25/abs(nmin); double dw = 0.25/abs(nmin);
size_t nplanes = max<size_t>(2,((wmax-wmin)/dw+1));
dw=(wmax-wmin)/(nplanes-1);
auto w_supp = gconf.Supp(); auto w_supp = gconf.Supp();
EC_Kernel_with_correction<T> kernel(w_supp); EC_Kernel_with_correction<T> kernel(w_supp);
wmin -= 0.5*w_supp*dw; wmin -= (0.5*w_supp-1)*dw;
wmax += 0.5*w_supp*dw; wmax += (0.5*w_supp-1)*dw;
size_t nplanes = size_t((wmax-wmin)/dw)+2; nplanes += w_supp-2;
cout << "nplanes: " << nplanes << endl; cout << "nplanes: " << nplanes << endl;
vector<size_t> nvis_plane(nplanes,0); vector<size_t> nvis_plane(nplanes,0);
vector<int> minplane(nvis); vector<int> minplane(nvis);
...@@ -1062,7 +1063,6 @@ cout << "working on w plane #" << iw << " containing " << nvis_plane[iw] ...@@ -1062,7 +1063,6 @@ cout << "working on w plane #" << iw << " containing " << nvis_plane[iw]
if (nvis_plane[iw]==0) continue; if (nvis_plane[iw]==0) continue;
double wcur = wmin+iw*dw; double wcur = wmin+iw*dw;
size_t cnt=0; size_t cnt=0;
cout << "blip0" << endl;
vis_loc_.resize(nvis_plane[iw]); vis_loc_.resize(nvis_plane[iw]);
auto vis_loc = mav<complex<T>,1>(vis_loc_.data(), {nvis_plane[iw]}); auto vis_loc = mav<complex<T>,1>(vis_loc_.data(), {nvis_plane[iw]});
idx_loc_.resize(nvis_plane[iw]); idx_loc_.resize(nvis_plane[iw]);
...@@ -1076,14 +1076,10 @@ cout << "blip0" << endl; ...@@ -1076,14 +1076,10 @@ cout << "blip0" << endl;
++cnt; ++cnt;
} }
myassert(cnt==nvis_plane[iw],"must not happen 2"); myassert(cnt==nvis_plane[iw],"must not happen 2");
cout << "blip1" << endl;
grid_.fill(0.); grid_.fill(0.);
vis2grid_c(baselines, gconf, const_mav<uint32_t,1>(idx_loc), const_mav<complex<T>,1>(vis_loc), grid, const_mav<T,1>(nullptr,{0})); vis2grid_c(baselines, gconf, const_mav<uint32_t,1>(idx_loc), const_mav<complex<T>,1>(vis_loc), grid, const_mav<T,1>(nullptr,{0}));
cout << "blip2" << endl;
gconf.grid2dirty_c(grid, tdirty); gconf.grid2dirty_c(grid, tdirty);
cout << "blip3" << endl;
gconf.apply_wscreen(tdirty, tdirty, wcur, false); gconf.apply_wscreen(tdirty, tdirty, wcur, false);
cout << "blip4" << endl;
for (size_t i=0; i<nx_dirty; ++i) for (size_t i=0; i<nx_dirty; ++i)
for (size_t j=0; j<ny_dirty; ++j) for (size_t j=0; j<ny_dirty; ++j)
dirty(i,j) += tdirty(i,j).real(); dirty(i,j) += tdirty(i,j).real();
...@@ -1118,18 +1114,20 @@ template<typename T> void dirty2vis_wstack(const Baselines<T> &baselines, ...@@ -1118,18 +1114,20 @@ template<typename T> void dirty2vis_wstack(const Baselines<T> &baselines,
wmin = min(wmin,wval[ipart]); wmin = min(wmin,wval[ipart]);
wmax = max(wmax,wval[ipart]); wmax = max(wmax,wval[ipart]);
} }
cout << "data w range: " << wmin << " to " << wmax << endl; //cout << "data w range: " << wmin << " to " << wmax << endl;
double x0 = -0.5*nx_dirty*psx, double x0 = -0.5*nx_dirty*psx,
y0 = -0.5*ny_dirty*psy; y0 = -0.5*ny_dirty*psy;
double nmin = sqrt(max(1.-x0*x0-y0*y0,0.))-1.; double nmin = sqrt(max(1.-x0*x0-y0*y0,0.))-1.;
double dw = 0.25/abs(nmin); double dw = 0.25/abs(nmin);
size_t nplanes = max<size_t>(2,((wmax-wmin)/dw+1));
dw=(wmax-wmin)/(nplanes-1);
auto w_supp = gconf.Supp(); auto w_supp = gconf.Supp();
EC_Kernel_with_correction<T> kernel(w_supp); EC_Kernel_with_correction<T> kernel(w_supp);
wmin -= 0.5*w_supp*dw; wmin -= (0.5*w_supp-1)*dw;
wmax += 0.5*w_supp*dw; wmax += (0.5*w_supp-1)*dw;
size_t nplanes = size_t((wmax-wmin)/dw)+2; nplanes += w_supp-2;
cout << "nplanes: " << nplanes << endl; cout << "nplanes: " << nplanes << endl;
vector<size_t> nvis_plane(nplanes,0); vector<size_t> nvis_plane(nplanes,0);
vector<int> minplane(nvis); vector<int> minplane(nvis);
...@@ -1150,7 +1148,7 @@ cout << "nplanes: " << nplanes << endl; ...@@ -1150,7 +1148,7 @@ cout << "nplanes: " << nplanes << endl;
for (size_t j=0; j<ny_dirty; ++j) for (size_t j=0; j<ny_dirty; ++j)
tdirty(i,j) = dirty(i,j); tdirty(i,j) = dirty(i,j);
// correct for w gridding // correct for w gridding
cout << "applying correction for gridding in w direction" << endl; //cout << "applying correction for gridding in w direction" << endl;
apply_wcorr(gconf, tdirty, kernel, dw); apply_wcorr(gconf, tdirty, kernel, dw);
tmpStorage<complex<T>,2> grid_({nu,nv}); tmpStorage<complex<T>,2> grid_({nu,nv});
......
...@@ -72,11 +72,11 @@ def test_adjointness(nxdirty, nydirty, nrow, nchan, epsilon): ...@@ -72,11 +72,11 @@ def test_adjointness(nxdirty, nydirty, nrow, nchan, epsilon):
def test_adjointness_wgridding(nxdirty, nydirty, nrow, nchan, epsilon): def test_adjointness_wgridding(nxdirty, nydirty, nrow, nchan, epsilon):
np.random.seed(42) np.random.seed(42)
bl, conf, idx = _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow) bl, conf, idx = _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow)
vis = np.random.rand(*idx.shape)-0.5 vis = np.random.rand(*idx.shape)-0.5 + 1j*(np.random.rand(*idx.shape)-0.5)
dirty = np.random.rand(conf.Nxdirty(), conf.Nydirty())-0.5 dirty = np.random.rand(conf.Nxdirty(), conf.Nydirty())-0.5
dirty2 = ng.vis2dirty_wstack(bl, conf, idx, vis) dirty2 = ng.vis2dirty_wstack(bl, conf, idx, vis)
vis2 = ng.dirty2vis_wstack(bl, conf, idx, dirty) vis2 = ng.dirty2vis_wstack(bl, conf, idx, dirty)
assert_allclose(np.vdot(vis, vis2), np.vdot(dirty2, dirty), rtol=epsilon) assert_allclose(np.vdot(vis, vis2).real, np.vdot(dirty2, dirty), rtol=epsilon)
@pmp("nxdirty", (128,)) @pmp("nxdirty", (128,))
......
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