Commit 6a2c68ec authored by Martin Reinecke's avatar Martin Reinecke
Browse files

cleanup

parent 3d2ae5e6
...@@ -89,7 +89,7 @@ template<typename T, size_t ndim> class mav ...@@ -89,7 +89,7 @@ template<typename T, size_t ndim> class mav
for (auto v: shp) res*=v; for (auto v: shp) res*=v;
return res; return res;
} }
size_t stride(size_t i) const { return str[i]; } ptrdiff_t stride(size_t i) const { return str[i]; }
T *data() T *data()
{ return d; } { return d; }
const T *data() const const T *data() const
...@@ -454,8 +454,8 @@ template<typename T> class GridderConfig ...@@ -454,8 +454,8 @@ template<typename T> class GridderConfig
checkShape(grid.shape(), {nu,nv}); checkShape(grid.shape(), {nu,nv});
checkShape(dirty.shape(), {nx_dirty,ny_dirty}); checkShape(dirty.shape(), {nx_dirty,ny_dirty});
vector<complex<T>> tmpdat(nu*nv); vector<complex<T>> tmpdat(nu*nv);
auto tmp=mav<complex<T>,2>(tmpdat.data(),{nu,nv},{nv,1}); auto tmp=mav<complex<T>,2>(tmpdat.data(),{nu,nv},{ptrdiff_t(nv),ptrdiff_t(1)});
constexpr auto sc = sizeof(complex<T>); constexpr auto sc = ptrdiff_t(sizeof(complex<T>));
pocketfft::c2c({nu,nv},{grid.stride(0)*sc,grid.stride(1)*sc}, pocketfft::c2c({nu,nv},{grid.stride(0)*sc,grid.stride(1)*sc},
{tmp.stride(0)*sc, tmp.stride(1)*sc}, {0,1}, pocketfft::BACKWARD, {tmp.stride(0)*sc, tmp.stride(1)*sc}, {0,1}, pocketfft::BACKWARD,
grid.data(), tmp.data(), T(1), nthreads); grid.data(), tmp.data(), T(1), nthreads);
...@@ -487,7 +487,7 @@ template<typename T> class GridderConfig ...@@ -487,7 +487,7 @@ template<typename T> class GridderConfig
if (j2>=nv) j2-=nv; if (j2>=nv) j2-=nv;
grid(i2,j2) = dirty(i,j)*cfu[i]*cfv[j]; grid(i2,j2) = dirty(i,j)*cfu[i]*cfv[j];
} }
constexpr auto sc = sizeof(complex<T>); constexpr auto sc = ptrdiff_t(sizeof(complex<T>));
pocketfft::stride_t strides{grid.stride(0)*sc,grid.stride(1)*sc}; pocketfft::stride_t strides{grid.stride(0)*sc,grid.stride(1)*sc};
pocketfft::c2c({nu,nv}, strides, strides, {0,1}, pocketfft::FORWARD, pocketfft::c2c({nu,nv}, strides, strides, {0,1}, pocketfft::FORWARD,
grid.data(), grid.data(), T(1), nthreads); grid.data(), grid.data(), T(1), nthreads);
...@@ -874,8 +874,8 @@ cout << "nplanes: " << nplanes << endl; ...@@ -874,8 +874,8 @@ cout << "nplanes: " << nplanes << endl;
++nvis_plane[i]; ++nvis_plane[i];
} }
for (ptrdiff_t i=0; i<nx_dirty; ++i) for (size_t i=0; i<nx_dirty; ++i)
for (ptrdiff_t j=0; j<ny_dirty; ++j) for (size_t j=0; j<ny_dirty; ++j)
dirty(i,j) = 0.; dirty(i,j) = 0.;
for (size_t iw=0; iw<nplanes; ++iw) for (size_t iw=0; iw<nplanes; ++iw)
{ {
...@@ -905,8 +905,8 @@ myassert(cnt==nvis_plane[iw],"must not happen 2"); ...@@ -905,8 +905,8 @@ myassert(cnt==nvis_plane[iw],"must not happen 2");
auto tdirty=mav<complex<T>,2>(tdirty_.data(),{nx_dirty, ny_dirty}); auto tdirty=mav<complex<T>,2>(tdirty_.data(),{nx_dirty, ny_dirty});
gconf.grid2dirty_c(grid, tdirty); gconf.grid2dirty_c(grid, tdirty);
gconf.apply_wscreen(tdirty, tdirty, wcur, false); gconf.apply_wscreen(tdirty, tdirty, wcur, false);
for (ptrdiff_t i=0; i<nx_dirty; ++i) for (size_t i=0; i<nx_dirty; ++i)
for (ptrdiff_t j=0; j<ny_dirty; ++j) for (size_t j=0; j<ny_dirty; ++j)
dirty(i,j) += tdirty(i,j); dirty(i,j) += tdirty(i,j);
} }
// correct for w gridding // correct for w gridding
......
...@@ -489,7 +489,6 @@ template<typename T> pyarr<complex<T>> Pyms2grid_c( ...@@ -489,7 +489,6 @@ template<typename T> pyarr<complex<T>> Pyms2grid_c(
auto nchan = baselines.Nchannels(); auto nchan = baselines.Nchannels();
auto ms2 = make_const_mav<2>(ms_); auto ms2 = make_const_mav<2>(ms_);
auto idx2 = make_const_mav<1>(idx_); auto idx2 = make_const_mav<1>(idx_);
size_t nvis = idx2.shape(0);
pyarr<T> wgt = providePotentialArray<T>(wgt_, {nrows, nchan}); pyarr<T> wgt = providePotentialArray<T>(wgt_, {nrows, nchan});
auto wgt2 = make_const_mav<2>(wgt); auto wgt2 = make_const_mav<2>(wgt);
auto res = provideCArray<complex<T>>(grid_in, {gconf.Nu(), gconf.Nv()}); auto res = provideCArray<complex<T>>(grid_in, {gconf.Nu(), gconf.Nv()});
...@@ -551,7 +550,6 @@ template<typename T> pyarr<complex<T>> Pygrid2ms_c( ...@@ -551,7 +550,6 @@ template<typename T> pyarr<complex<T>> Pygrid2ms_c(
auto nchan = baselines.Nchannels(); auto nchan = baselines.Nchannels();
auto grid2 = make_const_mav<2>(grid_); auto grid2 = make_const_mav<2>(grid_);
auto idx2 = make_const_mav<1>(idx_); auto idx2 = make_const_mav<1>(idx_);
size_t nvis = idx2.shape(0);
pyarr<T> wgt = providePotentialArray<T>(wgt_, {nrows, nchan}); pyarr<T> wgt = providePotentialArray<T>(wgt_, {nrows, nchan});
auto wgt2 = make_const_mav<2>(wgt); auto wgt2 = make_const_mav<2>(wgt);
auto res = provideCArray<complex<T>>(ms_in, {nrows, nchan}); auto res = provideCArray<complex<T>>(ms_in, {nrows, nchan});
...@@ -562,65 +560,6 @@ template<typename T> pyarr<complex<T>> Pygrid2ms_c( ...@@ -562,65 +560,6 @@ template<typename T> pyarr<complex<T>> Pygrid2ms_c(
} }
return res; return res;
} }
template<typename T> pyarr<complex<T>> grid2ms_c(const PyBaselines<T> &baselines,
const PyGridderConfig<T> &gconf, const pyarr<uint32_t> &idx_,
const pyarr<complex<T>> &grid_, py::object &ms_in, const py::object &wgt_)
{
auto nrows = baselines.Nrows();
auto nchan = baselines.Nchannels();
size_t nu=gconf.Nu(), nv=gconf.Nv();
checkArray(idx_, "idx", {0});
auto grid = grid_.data();
checkArray(grid_, "grid", {nu, nv});
size_t nvis = size_t(idx_.shape(0));
auto idx = idx_.template unchecked<1>();
auto wgt2 = providePotentialArray<T>(wgt_, {nrows, nchan});
bool have_wgt = wgt2.size()>0;
auto wgt = wgt2.template unchecked<2>();
auto res = provideArray<complex<T>>(ms_in, {nrows, nchan});
auto ms = res.template mutable_unchecked<2>();
{
py::gil_scoped_release release;
T beta = gconf.Beta();
size_t supp = gconf.Supp();
// Loop over sampling points
#pragma omp parallel num_threads(nthreads)
{
Helper<T> hlp(gconf, grid, nullptr);
T emb = exp(-2*beta);
int jump = hlp.lineJump();
const T * ku = hlp.kernel.data();
const T * kv = hlp.kernel.data()+supp;
#pragma omp for schedule(guided,100)
for (size_t ipart=0; ipart<nvis; ++ipart)
{
auto tidx = idx(ipart);
auto row = tidx/nchan;
auto chan = tidx-row*nchan;
UVW<T> coord = baselines.effectiveCoord(tidx);
hlp.prep(coord.u, coord.v);
complex<T> r = 0;
const auto * ptr = hlp.p0r;
for (size_t cu=0; cu<supp; ++cu)
{
complex<T> tmp(0);
for (size_t cv=0; cv<supp; ++cv)
tmp += ptr[cv] * kv[cv];
r += tmp*ku[cu];
ptr += jump;
}
r*=emb;
if (have_wgt)
r*=wgt(row, chan);
ms(row,chan) += r*emb;
}
}
}
return res;
}
template<typename T> pyarr<complex<T>> apply_holo( template<typename T> pyarr<complex<T>> apply_holo(
const PyBaselines<T> &baselines, const PyGridderConfig<T> &gconf, const PyBaselines<T> &baselines, const PyGridderConfig<T> &gconf,
......
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