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

tweak

parent c7d0c24a
...@@ -998,6 +998,8 @@ template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &base ...@@ -998,6 +998,8 @@ template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &base
// determine w values for every visibility, and min/max w; // determine w values for every visibility, and min/max w;
T wmin=T(1e38), wmax=T(-1e38); T wmin=T(1e38), wmax=T(-1e38);
vector<T> wval(nvis); vector<T> wval(nvis);
#pragma omp parallel for num_threads(nthreads) reduction(min:wmin) reduction(max:wmax)
for (size_t ipart=0; ipart<nvis; ++ipart) for (size_t ipart=0; ipart<nvis; ++ipart)
{ {
wval[ipart] = srv.getCoord(ipart).w; wval[ipart] = srv.getCoord(ipart).w;
...@@ -1019,20 +1021,28 @@ template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &base ...@@ -1019,20 +1021,28 @@ template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &base
nplanes += w_supp-2; nplanes += w_supp-2;
vector<size_t> nvis_plane(nplanes,0); vector<size_t> nvis_plane(nplanes,0);
vector<int> minplane(nvis); vector<int> minplane(nvis);
#pragma omp parallel num_threads(nthreads)
{
vector<size_t> nvploc(nplanes,0);
#pragma omp for
for (size_t ipart=0; ipart<nvis; ++ipart) for (size_t ipart=0; ipart<nvis; ++ipart)
{ {
int plane0 = int((wval[ipart]-wmin)/dw+0.5*w_supp)-int(w_supp-1); int plane0 = int((wval[ipart]-wmin)/dw+0.5*w_supp)-int(w_supp-1);
minplane[ipart]=plane0; minplane[ipart]=plane0;
for (int i=max<int>(0,plane0); i<min<int>(nplanes,plane0+w_supp); ++i) for (int i=max<int>(0,plane0); i<min<int>(nplanes,plane0+w_supp); ++i)
++nvis_plane[i]; ++nvploc[i];
} }
#pragma omp critical
for (size_t i=0; i<nplanes; ++i)
nvis_plane[i]+=nvploc[i];
}
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) = 0.; dirty(i,j) = 0.;
{ {
vector<complex<T>> vis_loc_; vector<complex<T>> vis_loc_;
vector<uint32_t> idx_loc_; vector<uint32_t> idx_loc_;
vector<uint32_t> mapper;
tmpStorage<complex<T>,2> grid_({nu,nv}); tmpStorage<complex<T>,2> grid_({nu,nv});
auto grid=grid_.getMav(); auto grid=grid_.getMav();
tmpStorage<complex<T>,2> tdirty_({nx_dirty,ny_dirty}); tmpStorage<complex<T>,2> tdirty_({nx_dirty,ny_dirty});
...@@ -1046,15 +1056,20 @@ template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &base ...@@ -1046,15 +1056,20 @@ template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &base
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]);
auto idx_loc = mav<uint32_t,1>(idx_loc_.data(), {nvis_plane[iw]}); auto idx_loc = mav<uint32_t,1>(idx_loc_.data(), {nvis_plane[iw]});
mapper.resize(nvis_plane[iw]);
for (size_t ipart=0; ipart<nvis; ++ipart) for (size_t ipart=0; ipart<nvis; ++ipart)
if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+w_supp)) if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+w_supp))
{ mapper[cnt++] = ipart;
double x=min(1.,2./(w_supp*dw)*abs(wcur-wval[ipart]));
vis_loc[cnt] = srv.getVis(ipart)*kernel(x);
idx_loc[cnt] = srv.getIdx(ipart);
++cnt;
}
myassert(cnt==nvis_plane[iw],"must not happen 2"); myassert(cnt==nvis_plane[iw],"must not happen 2");
#pragma omp parallel for num_threads(nthreads)
for (size_t i=0; i<nvis_plane[iw]; ++i)
{
auto ipart = mapper[i];
double x=min(1.,2./(w_supp*dw)*abs(wcur-wval[ipart]));
vis_loc[i] = srv.getVis(ipart)*kernel(x);
idx_loc[i] = srv.getIdx(ipart);
}
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}));
gconf.grid2dirty_c(grid, tdirty); gconf.grid2dirty_c(grid, tdirty);
...@@ -1091,6 +1106,8 @@ template<typename T, typename Serv> void dirty2x_wstack(const Baselines<T> &base ...@@ -1091,6 +1106,8 @@ template<typename T, typename Serv> void dirty2x_wstack(const Baselines<T> &base
// determine w values for every visibility, and min/max w; // determine w values for every visibility, and min/max w;
T wmin=T(1e38), wmax=T(-1e38); T wmin=T(1e38), wmax=T(-1e38);
vector<T> wval(nvis); vector<T> wval(nvis);
#pragma omp parallel for num_threads(nthreads) reduction(min:wmin) reduction(max:wmax)
for (size_t ipart=0; ipart<nvis; ++ipart) for (size_t ipart=0; ipart<nvis; ++ipart)
{ {
wval[ipart] = srv.getCoord(ipart).w; wval[ipart] = srv.getCoord(ipart).w;
...@@ -1112,13 +1129,21 @@ template<typename T, typename Serv> void dirty2x_wstack(const Baselines<T> &base ...@@ -1112,13 +1129,21 @@ template<typename T, typename Serv> void dirty2x_wstack(const Baselines<T> &base
nplanes += w_supp-2; nplanes += w_supp-2;
vector<size_t> nvis_plane(nplanes,0); vector<size_t> nvis_plane(nplanes,0);
vector<int> minplane(nvis); vector<int> minplane(nvis);
#pragma omp parallel num_threads(nthreads)
{
vector<size_t> nvploc(nplanes,0);
#pragma omp for
for (size_t ipart=0; ipart<nvis; ++ipart) for (size_t ipart=0; ipart<nvis; ++ipart)
{ {
int plane0 = int((wval[ipart]-wmin)/dw+0.5*w_supp)-int(w_supp-1); int plane0 = int((wval[ipart]-wmin)/dw+0.5*w_supp)-int(w_supp-1);
minplane[ipart]=plane0; minplane[ipart]=plane0;
for (int i=max<int>(0,plane0); i<min<int>(nplanes,plane0+w_supp); ++i) for (int i=max<int>(0,plane0); i<min<int>(nplanes,plane0+w_supp); ++i)
++nvis_plane[i]; ++nvploc[i];
} }
#pragma omp critical
for (size_t i=0; i<nplanes; ++i)
nvis_plane[i]+=nvploc[i];
}
for (size_t i=0; i<nvis; ++i) for (size_t i=0; i<nvis; ++i)
srv.setVis(i, 0.); srv.setVis(i, 0.);
...@@ -1131,37 +1156,43 @@ template<typename T, typename Serv> void dirty2x_wstack(const Baselines<T> &base ...@@ -1131,37 +1156,43 @@ template<typename T, typename Serv> void dirty2x_wstack(const Baselines<T> &base
// correct for w gridding // correct for w gridding
apply_wcorr(gconf, tdirty, kernel, dw); apply_wcorr(gconf, tdirty, kernel, dw);
vector<complex<T>> vis_loc_;
vector<uint32_t> idx_loc_;
vector<uint32_t> mapper;
tmpStorage<complex<T>,2> grid_({nu,nv}); tmpStorage<complex<T>,2> grid_({nu,nv});
auto grid=grid_.getMav(); auto grid=grid_.getMav();
tmpStorage<complex<T>,2> tdirty2_({nx_dirty,ny_dirty}); tmpStorage<complex<T>,2> tdirty2_({nx_dirty,ny_dirty});
auto tdirty2=tdirty2_.getMav(); auto tdirty2=tdirty2_.getMav();
vector<complex<T>> vis_loc_;
vector<uint32_t> idx_loc_;
for (size_t iw=0; iw<nplanes; ++iw) for (size_t iw=0; iw<nplanes; ++iw)
{ {
if (nvis_plane[iw]==0) continue; if (nvis_plane[iw]==0) continue;
double wcur = wmin+iw*dw; double wcur = wmin+iw*dw;
gconf.apply_wscreen(tdirty, tdirty2, wcur, true);
gconf.dirty2grid_c(tdirty2, grid);
size_t cnt=0; size_t cnt=0;
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]);
auto idx_loc = mav<uint32_t,1>(idx_loc_.data(), {nvis_plane[iw]}); auto idx_loc = mav<uint32_t,1>(idx_loc_.data(), {nvis_plane[iw]});
mapper.resize(nvis_plane[iw]);
for (size_t ipart=0; ipart<nvis; ++ipart) for (size_t ipart=0; ipart<nvis; ++ipart)
if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+w_supp)) if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+w_supp))
idx_loc[cnt++] = srv.getIdx(ipart); mapper[cnt++] = ipart;
myassert(cnt==nvis_plane[iw],"must not happen 2"); myassert(cnt==nvis_plane[iw],"must not happen 2");
#pragma omp parallel for num_threads(nthreads)
for (size_t i=0; i<nvis_plane[iw]; ++i)
idx_loc[i] = srv.getIdx(mapper[i]);
gconf.apply_wscreen(tdirty, tdirty2, wcur, true);
gconf.dirty2grid_c(tdirty2, grid);
grid2vis_c(baselines, gconf, const_mav<uint32_t,1>(idx_loc), const_mav<complex<T>,2>(grid), vis_loc, const_mav<T,1>(nullptr,{0})); grid2vis_c(baselines, gconf, const_mav<uint32_t,1>(idx_loc), const_mav<complex<T>,2>(grid), vis_loc, const_mav<T,1>(nullptr,{0}));
cnt=0;
for (size_t ipart=0; ipart<nvis; ++ipart) #pragma omp parallel for num_threads(nthreads)
if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+w_supp)) for (size_t i=0; i<nvis_plane[iw]; ++i)
{ {
double x=min(1.,2./(w_supp*dw)*abs(wcur-wval[ipart])); auto ipart = mapper[i];
srv.addVis(ipart, vis_loc[cnt]*kernel(x)); double x=min(1.,2./(w_supp*dw)*abs(wcur-wval[ipart]));
++cnt; srv.addVis(ipart, vis_loc[i]*kernel(x));
} }
} }
} }
template<typename T> void dirty2vis_wstack(const Baselines<T> &baselines, template<typename T> void dirty2vis_wstack(const Baselines<T> &baselines,
......
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