Commit a7869f66 authored by Martin Reinecke's avatar Martin Reinecke

update pocketfft; use ms values for flagging

parent 38f63ce9
......@@ -1619,7 +1619,7 @@ void fillIdx(const Baselines &baselines,
size_t lo, hi;
calc_share(nthr, id, nrow, lo, hi);
vector<idx_t> &lacc(acc[id]);
lacc.resize(nbu*nbv+1, 0);
lacc.assign(nbu*nbv+1, 0);
for (idx_t irow=lo, idx=lo*(chend-chbegin); irow<hi; ++irow)
for (int ichan=chbegin; ichan<chend; ++ichan, ++idx)
......@@ -1666,13 +1666,16 @@ void fillIdx(const Baselines &baselines,
}
template<typename T> vector<idx_t> getWgtIndices(const Baselines &baselines,
const GridderConfig &gconf, const const_mav<T,2> &wgt)
const GridderConfig &gconf, const const_mav<T,2> &wgt,
const const_mav<complex<T>,2> &ms)
{
size_t nrow=baselines.Nrows(),
nchan=baselines.Nchannels(),
nsafe=gconf.Nsafe();
bool have_wgt=wgt.size()!=0;
if (have_wgt) checkShape(wgt.shape(),{nrow,nchan});
bool have_ms=ms.size()!=0;
if (have_ms) checkShape(ms.shape(), {nrow,nchan});
constexpr int side=1<<logsquare;
size_t nbu = (gconf.Nu()+1+side-1) >> logsquare,
nbv = (gconf.Nv()+1+side-1) >> logsquare;
......@@ -1690,10 +1693,11 @@ template<typename T> vector<idx_t> getWgtIndices(const Baselines &baselines,
size_t lo, hi;
calc_share(nthr, id, nrow, lo, hi);
vector<idx_t> &lacc(acc[id]);
lacc.resize(nbu*nbv+1, 0);
lacc.assign(nbu*nbv+1, 0);
for (idx_t irow=lo, idx=lo*nchan; irow<hi; ++irow)
for (idx_t ichan=0; ichan<nchan; ++ichan, ++idx)
if ((!have_wgt) || (wgt(irow,ichan)!=0))
if (((!have_ms ) || (norm(ms (irow,ichan))!=0)) &&
((!have_wgt) || (wgt(irow,ichan)!=0)))
{
auto uvw = baselines.effectiveCoord(RowChan{irow,idx_t(ichan)});
if (uvw.w<0) uvw.Flip();
......@@ -1739,7 +1743,7 @@ template<typename T> void ms2dirty(const const_mav<double,2> &uvw,
{
Baselines baselines(uvw, freq);
GridderConfig gconf(dirty.shape(0), dirty.shape(1), epsilon, pixsize_x, pixsize_y, nthreads);
auto idx = getWgtIndices(baselines, gconf, wgt);
auto idx = getWgtIndices(baselines, gconf, wgt, ms);
auto idx2 = const_mav<idx_t,1>(idx.data(),{idx.size()});
x2dirty(gconf, makeMsServ(baselines,idx2,ms,wgt), dirty, do_wstacking, verbosity);
}
......@@ -1751,7 +1755,8 @@ template<typename T> void dirty2ms(const const_mav<double,2> &uvw,
{
Baselines baselines(uvw, freq);
GridderConfig gconf(dirty.shape(0), dirty.shape(1), epsilon, pixsize_x, pixsize_y, nthreads);
auto idx = getWgtIndices(baselines, gconf, wgt);
const_mav<complex<T>,2> null_ms(nullptr, {0,0});
auto idx = getWgtIndices(baselines, gconf, wgt, null_ms);
auto idx2 = const_mav<idx_t,1>(idx.data(),{idx.size()});
ms.fill(0);
dirty2x(gconf, dirty, makeMsServ(baselines,idx2,ms,wgt), do_wstacking, verbosity);
......
......@@ -3387,6 +3387,36 @@ template<typename T> void r2r_separable_hartley(const shape_t &shape,
false);
}
template<typename T> void r2r_genuine_hartley(const shape_t &shape,
const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes,
const T *data_in, T *data_out, T fct, size_t nthreads=1)
{
if (util::prod(shape)==0) return;
if (axes.size()==1)
return r2r_separable_hartley(shape, stride_in, stride_out, axes, data_in,
data_out, fct, nthreads);
util::sanity_check(shape, stride_in, stride_out, data_in==data_out, axes);
shape_t tshp(shape);
tshp[axes.back()] = tshp[axes.back()]/2+1;
arr<complex<T>> tdata(util::prod(tshp));
stride_t tstride(shape.size());
tstride.back()=sizeof(complex<T>);
for (size_t i=tstride.size()-1; i>0; --i)
tstride[i-1]=tstride[i]*ptrdiff_t(tshp[i]);
r2c(shape, stride_in, tstride, axes, true, data_in, tdata.data(), fct, nthreads);
cndarr<cmplx<T>> atmp(tdata.data(), tshp, tstride);
ndarr<T> aout(data_out, shape, stride_out);
simple_iter iin(atmp);
rev_iter iout(aout, axes);
while(iin.remaining()>0)
{
auto v = atmp[iin.ofs()];
aout[iout.ofs()] = v.r+v.i;
aout[iout.rev_ofs()] = v.r-v.i;
iin.advance(); iout.advance();
}
}
} // namespace detail
using detail::FORWARD;
......@@ -3398,6 +3428,7 @@ using detail::c2r;
using detail::r2c;
using detail::r2r_fftpack;
using detail::r2r_separable_hartley;
using detail::r2r_genuine_hartley;
using detail::dct;
using detail::dst;
......
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