diff --git a/gridder_cxx.h b/gridder_cxx.h
index 6ddec362398c2660436af5782efdf328da2a50b6..8c23d529ad6de5b428d665b716cc0d52d4e18487 100644
--- a/gridder_cxx.h
+++ b/gridder_cxx.h
@@ -31,35 +31,32 @@
#include "pocketfft_hdronly.h"
#include <array>
-template<typename T, size_t ndim, bool c_contiguous> class mav
+namespace gridder {
+
+using namespace std;
+
+template<typename T, size_t ndim> class mav
{
static_assert((ndim>0) && (ndim<3), "only supports 1D and 2D arrays");
private:
T *d;
- std::array<size_t, ndim> shp;
- std::array<ptrdiff_t, ndim> str;
+ array<size_t, ndim> shp;
+ array<ptrdiff_t, ndim> str;
public:
- mav(T *d_, const std::array<size_t,ndim> &shp_,
- const std::array<ptrdiff_t,ndim> &str_)
- : d(d_), shp(shp_), str(str_)
- { static_assert(!c_contiguous, "this type does not accept strides"); }
- mav(T *d_, const std::array<size_t,ndim> &shp_)
+ mav(T *d_, const array<size_t,ndim> &shp_,
+ const array<ptrdiff_t,ndim> &str_)
+ : d(d_), shp(shp_), str(str_) {}
+ mav(T *d_, const array<size_t,ndim> &shp_)
: d(d_), shp(shp_)
{
- static_assert(c_contiguous, "this type requires strides");
str[ndim-1]=1;
for (size_t d=2; d<=ndim; ++d)
str[ndim-d] = str[ndim-d+1]*shp[ndim-d+1];
}
- operator mav<const T, ndim, true>() const
- {
- static_assert(c_contiguous);
- return mav<const T, ndim, true>(d,shp);
- }
- operator mav<const T, ndim, false>() const
- { return mav<const T, ndim, false>(d,shp, str); }
+ operator mav<const T, ndim>() const
+ { return mav<const T, ndim>(d,shp,str); }
T &operator[](size_t i)
{ return operator()(i); }
const T &operator[](size_t i) const
@@ -67,25 +64,25 @@ template<typename T, size_t ndim, bool c_contiguous> class mav
T &operator()(size_t i)
{
static_assert(ndim==1, "ndim must be 1");
- return c_contiguous ? d[i] : d[str[0]*i];
+ return d[str[0]*i];
}
const T &operator()(size_t i) const
{
static_assert(ndim==1, "ndim must be 1");
- return c_contiguous ? d[i] : d[str[0]*i];
+ return d[str[0]*i];
}
T &operator()(size_t i, size_t j)
{
static_assert(ndim==2, "ndim must be 2");
- return c_contiguous ? d[str[0]*i + j] : d[str[0]*i + str[1]*j];
+ return d[str[0]*i + str[1]*j];
}
const T &operator()(size_t i, size_t j) const
{
static_assert(ndim==2, "ndim must be 2");
- return c_contiguous ? d[str[0]*i + j] : d[str[0]*i + str[1]*j];
+ return d[str[0]*i + str[1]*j];
}
size_t shape(size_t i) const { return shp[i]; }
- const std::array<size_t,ndim> &shape() const { return shp; }
+ const array<size_t,ndim> &shape() const { return shp; }
size_t size() const
{
size_t res=1;
@@ -94,20 +91,11 @@ template<typename T, size_t ndim, bool c_contiguous> class mav
}
size_t stride(size_t i) const { return str[i]; }
T *data()
- {
-// static_assert(c_contiguous, "type is not C contiguous");
- return d;
- }
+ { return d; }
const T *data() const
- {
-// static_assert(c_contiguous, "type is not C contiguous");
- return d;
- }
+ { return d; }
};
-template<typename T, size_t ndim> using cmav = mav<T,ndim,true>;
-template<typename T, size_t ndim> using smav = mav<T,ndim,false>;
-
//
// basic utilities
//
@@ -115,11 +103,11 @@ template<typename T, size_t ndim> using smav = mav<T,ndim,false>;
void myassert(bool cond, const char *msg)
{
if (cond) return;
- throw std::runtime_error(msg);
+ throw runtime_error(msg);
}
template<size_t ndim> void checkShape
- (const std::array<size_t, ndim> &shp1, const std::array<size_t, ndim> &shp2)
+ (const array<size_t, ndim> &shp1, const array<size_t, ndim> &shp2)
{
for (size_t i=0; i<ndim; ++i)
myassert(shp1[i]==shp2[i], "shape mismatch");
@@ -138,27 +126,12 @@ template<typename T> inline T fmodulo (T v1, T v2)
static size_t nthreads = 1;
-constexpr auto set_nthreads_DS = R"""(
-Specifies the number of threads to be used by the module
-
-Parameters
-==========
-nthreads: int
- the number of threads to be used. Must be >=1.
-)""";
void set_nthreads(size_t nthreads_)
{
myassert(nthreads_>=1, "nthreads must be >= 1");
nthreads = nthreads_;
}
-constexpr auto get_nthreads_DS = R"""(
-Returns the number of threads used by the module
-
-Returns
-=======
-int : the number of threads used by the module
-)""";
size_t get_nthreads()
{ return nthreads; }
@@ -169,7 +142,7 @@ size_t get_nthreads()
static inline double one_minus_x2 (double x)
{ return (fabs(x)>0.1) ? (1.+x)*(1.-x) : 1.-x*x; }
-void legendre_prep(int n, std::vector<double> &x, std::vector<double> &w)
+void legendre_prep(int n, vector<double> &x, vector<double> &w)
{
constexpr double pi = 3.141592653589793238462643383279502884197;
constexpr double eps = 3e-14;
@@ -229,7 +202,7 @@ void legendre_prep(int n, std::vector<double> &x, std::vector<double> &w)
size_t get_supp(double epsilon)
{
- static const std::vector<double> maxmaperr { 1e8, 0.32, 0.021, 6.2e-4,
+ static const vector<double> maxmaperr { 1e8, 0.32, 0.021, 6.2e-4,
1.08e-5, 1.25e-7, 8.25e-10, 5.70e-12, 1.22e-13, 2.48e-15, 4.82e-17,
6.74e-19, 5.41e-21, 4.41e-23, 7.88e-25, 3.9e-26 };
@@ -237,70 +210,7 @@ size_t get_supp(double epsilon)
for (size_t i=1; i<maxmaperr.size(); ++i)
if (epssq>maxmaperr[i]) return i;
- throw std::runtime_error("requested epsilon too small - minimum is 2e-13");
- }
-
-template<typename T> void complex2hartley
- (const smav<const std::complex<T>, 2> &grid, smav<T,2> &grid2)
- {
- myassert(grid.shape()==grid2.shape(), "shape mismatch");
- size_t nu=grid.shape(0), nv=grid.shape(1);
-
-#pragma omp parallel for num_threads(nthreads)
- for (size_t u=0; u<nu; ++u)
- {
- size_t xu = (u==0) ? 0 : nu-u;
- for (size_t v=0; v<nv; ++v)
- {
- size_t xv = (v==0) ? 0 : nv-v;
- grid2(u,v) += T(0.5)*(grid( u, v).real()+grid( u, v).imag()+
- grid(xu,xv).real()-grid(xu,xv).imag());
- }
- }
- }
-
-template<typename T> void hartley2complex
- (const smav<const T,2> &grid, smav<std::complex<T>,2> &grid2)
- {
- myassert(grid.shape()==grid2.shape(), "shape mismatch");
- size_t nu=grid.shape(0), nv=grid.shape(1);
-
-#pragma omp parallel for num_threads(nthreads)
- for (size_t u=0; u<nu; ++u)
- {
- size_t xu = (u==0) ? 0 : nu-u;
- for (size_t v=0; v<nv; ++v)
- {
- size_t xv = (v==0) ? 0 : nv-v;
- T v1 = T(0.5)*grid( u, v);
- T v2 = T(0.5)*grid(xu,xv);
- grid2(u,v) = std::complex<T>(v1+v2, v1-v2);
- }
- }
- }
-
-template<typename T> void hartley2_2D(const smav<const T,2> &in, smav<T,2> &out)
- {
- myassert(in.shape()==out.shape(), "shape mismatch");
- size_t nu=in.shape(0), nv=in.shape(1), sz=sizeof(T);
- pocketfft::stride_t str{ptrdiff_t(sz*nv), ptrdiff_t(sz)};
- auto d_i = in.data();
- auto ptmp = out.data();
- pocketfft::r2r_separable_hartley({nu, nv}, str, str, {0,1}, d_i, ptmp, T(1),
- nthreads);
-#pragma omp parallel for num_threads(nthreads)
- for(size_t i=1; i<(nu+1)/2; ++i)
- for(size_t j=1; j<(nv+1)/2; ++j)
- {
- T a = ptmp[i*nv+j];
- T b = ptmp[(nu-i)*nv+j];
- T c = ptmp[i*nv+nv-j];
- T d = ptmp[(nu-i)*nv+nv-j];
- ptmp[i*nv+j] = T(0.5)*(a+b+c-d);
- ptmp[(nu-i)*nv+j] = T(0.5)*(a+b+d-c);
- ptmp[i*nv+nv-j] = T(0.5)*(a+c+d-b);
- ptmp[(nu-i)*nv+nv-j] = T(0.5)*(b+c+d-a);
- }
+ throw runtime_error("requested epsilon too small - minimum is 2e-13");
}
template<typename T> class EC_Kernel
@@ -322,7 +232,7 @@ template<typename T> class EC_Kernel_with_correction: public EC_Kernel<T>
protected:
static constexpr T pi = T(3.141592653589793238462643383279502884197L);
int p;
- std::vector<T> x, wgt, psi;
+ vector<T> x, wgt, psi;
using EC_Kernel<T>::supp;
public:
@@ -347,10 +257,10 @@ template<typename T> class EC_Kernel_with_correction: public EC_Kernel<T>
};
/* Compute correction factors for the ES gridding kernel
This implementation follows eqs. (3.8) to (3.10) of Barnett et al. 2018 */
-std::vector<double> correction_factors (size_t n, size_t nval, size_t supp)
+vector<double> correction_factors (size_t n, size_t nval, size_t supp)
{
EC_Kernel_with_correction<double> kernel(supp);
- std::vector<double> res(nval);
+ vector<double> res(nval);
double xn = 1./n;
#pragma omp parallel for schedule(static) num_threads(nthreads)
for (size_t k=0; k<nval; ++k)
@@ -370,12 +280,12 @@ template<typename T> struct UVW
template<typename T> class Baselines
{
protected:
- std::vector<UVW<T>> coord;
- std::vector<T> f_over_c;
+ vector<UVW<T>> coord;
+ vector<T> f_over_c;
size_t nrows, nchan;
public:
- Baselines(const smav<const T,2> &coord_, const smav<const T,1> &freq)
+ Baselines(const mav<const T,2> &coord_, const mav<const T,1> &freq)
{
constexpr double speedOfLight = 299792458.;
myassert(coord_.shape(1)==3, "dimension mismatch");
@@ -401,7 +311,7 @@ template<typename T> class Baselines
size_t Nrows() const { return nrows; }
size_t Nchannels() const { return nchan; }
- void effectiveUVW(const smav<const uint32_t,1> &idx, smav<T,2> &res) const
+ void effectiveUVW(const mav<const uint32_t,1> &idx, mav<T,2> &res) const
{
size_t nvis = idx.shape(0);
myassert(res.shape(0)==nvis, "shape mismatch");
@@ -415,8 +325,8 @@ template<typename T> class Baselines
}
}
- template<typename T2> void ms2vis(const smav<const T2,2> &ms,
- const smav<const uint32_t,1> &idx, smav<T2,1> &vis) const
+ template<typename T2> void ms2vis(const mav<const T2,2> &ms,
+ const mav<const uint32_t,1> &idx, mav<T2,1> &vis) const
{
myassert(ms.shape(0)==nrows, "shape mismatch");
myassert(ms.shape(1)==nchan, "shape mismatch");
@@ -432,8 +342,8 @@ template<typename T> class Baselines
}
}
- template<typename T2> void vis2ms(const smav<const T2,1> &vis,
- const smav<const uint32_t,1> &idx, smav<T2,2> &ms) const
+ template<typename T2> void vis2ms(const mav<const T2,1> &vis,
+ const mav<const uint32_t,1> &idx, mav<T2,2> &ms) const
{
size_t nvis = vis.shape(0);
myassert(idx.shape(0)==nvis, "shape mismatch");
@@ -457,11 +367,10 @@ template<typename T> class GridderConfig
double eps, psx, psy;
size_t supp, nsafe, nu, nv;
T beta;
- std::vector<T> cfu, cfv;
+ vector<T> cfu, cfv;
- std::complex<T> wscreen(double x, double y, double w, bool adjoint) const
+ complex<T> wscreen(double x, double y, double w, bool adjoint) const
{
- using namespace std;
constexpr double pi = 3.141592653589793238462643383279502884197;
#if 1
double eps = sqrt(x+y);
@@ -482,7 +391,7 @@ template<typename T> class GridderConfig
: nx_dirty(nxdirty), ny_dirty(nydirty), eps(epsilon),
psx(pixsize_x), psy(pixsize_y),
supp(get_supp(epsilon)), nsafe((supp+1)/2),
- nu(std::max(2*nsafe,2*nx_dirty)), nv(std::max(2*nsafe,2*ny_dirty)),
+ nu(max(2*nsafe,2*nx_dirty)), nv(max(2*nsafe,2*ny_dirty)),
beta(2.3*supp),
cfu(nx_dirty), cfv(ny_dirty)
{
@@ -514,78 +423,30 @@ template<typename T> class GridderConfig
size_t Nsafe() const { return nsafe; }
T Beta() const { return beta; }
- void grid2dirty(const smav<const T,2> &grid, smav<T,2> &grid2) const
- {
- checkShape(grid.shape(), {nu,nv});
- std::vector<T> tmpdat(nu*nv);
- auto tmp=smav<T,2>(tmpdat.data(),{nu,nv},{nv,1});
- hartley2_2D<T>(grid, tmp);
- checkShape(grid2.shape(), {nx_dirty, ny_dirty});
- for (size_t i=0; i<nx_dirty; ++i)
- for (size_t j=0; j<ny_dirty; ++j)
- {
- size_t i2 = nu-nx_dirty/2+i;
- if (i2>=nu) i2-=nu;
- size_t j2 = nv-ny_dirty/2+j;
- if (j2>=nv) j2-=nv;
- grid2(i,j) = tmp(i2,j2)*cfu[i]*cfv[j];
- }
- }
-#if 0
- pyarr_c<T> apply_taper(const pyarr_c<T> &img, bool divide) const
- {
- checkArray(img, "img", {nx_dirty, ny_dirty});
- auto pin = img.data();
- auto res = makeArray<T>({nx_dirty, ny_dirty});
- auto pout = res.mutable_data();
+ void apply_taper(const mav<const T,2> &img, mav<T,2> &img2, bool divide) const
{
- py::gil_scoped_release release;
+ checkShape(img.shape(), {nx_dirty, ny_dirty});
+ checkShape(img2.shape(), {nx_dirty, ny_dirty});
if (divide)
for (size_t i=0; i<nx_dirty; ++i)
for (size_t j=0; j<ny_dirty; ++j)
- pout[ny_dirty*i + j] = pin[ny_dirty*i + j]/(cfu[i]*cfv[j]);
+ img2(i,j) = img(i,j)/(cfu[i]*cfv[j]);
else
for (size_t i=0; i<nx_dirty; ++i)
for (size_t j=0; j<ny_dirty; ++j)
- pout[ny_dirty*i + j] = pin[ny_dirty*i + j]*cfu[i]*cfv[j];
- }
- return res;
- }
- pyarr_c<complex<T>> grid2dirty_c(const pyarr_c<complex<T>> &grid) const
- {
- checkArray(grid, "grid", {nu, nv});
- auto tmp = makeArray<complex<T>>({nu, nv});
- auto ptmp = tmp.mutable_data();
- pocketfft::c2c({nu,nv},{grid.strides(0),grid.strides(1)},
- {tmp.strides(0), tmp.strides(1)}, {0,1}, pocketfft::BACKWARD,
- grid.data(), tmp.mutable_data(), T(1), nthreads);
- auto res = makeArray<complex<T>>({nx_dirty, ny_dirty});
- auto pout = res.mutable_data();
- {
- py::gil_scoped_release release;
- for (size_t i=0; i<nx_dirty; ++i)
- for (size_t j=0; j<ny_dirty; ++j)
- {
- size_t i2 = nu-nx_dirty/2+i;
- if (i2>=nu) i2-=nu;
- size_t j2 = nv-ny_dirty/2+j;
- if (j2>=nv) j2-=nv;
- pout[ny_dirty*i + j] = ptmp[nv*i2+j2]*cfu[i]*cfv[j];
- }
- }
- return res;
+ img2(i,j) = img(i,j)*cfu[i]*cfv[j];
}
- pyarr_c<T> dirty2grid(const pyarr_c<T> &dirty) const
- {
- checkArray(dirty, "dirty", {nx_dirty, ny_dirty});
- auto pdirty = dirty.data();
- auto tmp = makeArray<T>({nu, nv});
- auto ptmp = tmp.mutable_data();
+ void grid2dirty_c(const mav<const complex<T>,2> &grid, mav<complex<T>,2> &dirty) const
{
- py::gil_scoped_release release;
- for (size_t i=0; i<nu*nv; ++i)
- ptmp[i] = 0.;
+ checkShape(grid.shape(), {nu,nv});
+ checkShape(dirty.shape(), {nx_dirty,ny_dirty});
+ vector<complex<T>> tmpdat(nu*nv);
+ auto tmp=mav<complex<T>,2>(tmpdat.data(),{nu,nv},{nv,1});
+ constexpr auto sc = sizeof(complex<T>);
+ pocketfft::c2c({nu,nv},{grid.stride(0)*sc,grid.stride(1)*sc},
+ {tmp.stride(0)*sc, tmp.stride(1)*sc}, {0,1}, pocketfft::BACKWARD,
+ grid.data(), tmp.data(), T(1), nthreads);
for (size_t i=0; i<nx_dirty; ++i)
for (size_t j=0; j<ny_dirty; ++j)
{
@@ -593,12 +454,10 @@ template<typename T> class GridderConfig
if (i2>=nu) i2-=nu;
size_t j2 = nv-ny_dirty/2+j;
if (j2>=nv) j2-=nv;
- ptmp[nv*i2+j2] = pdirty[ny_dirty*i + j]*cfu[i]*cfv[j];
+ dirty(i,j) = tmp(i2,j2)*cfu[i]*cfv[j];
}
}
- hartley2_2D<T>(tmp, tmp);
- return tmp;
- }
+#if 0
pyarr_c<complex<T>> dirty2grid_c(const pyarr_c<complex<T>> &dirty) const
{
checkArray(dirty, "dirty", {nx_dirty, ny_dirty});
@@ -677,7 +536,7 @@ template<typename T> class GridderConfig
constexpr int logsquare=4;
-template<typename T, typename T2=std::complex<T>> class Helper
+template<typename T, typename T2=complex<T>> class Helper
{
private:
const GridderConfig<T> &gconf;
@@ -689,7 +548,7 @@ template<typename T, typename T2=std::complex<T>> class Helper
int iu0, iv0; // start index of the current visibility
int bu0, bv0; // start index of the current buffer
- std::vector<T2> rbuf, wbuf;
+ vector<T2> rbuf, wbuf;
void dump() const
{
@@ -731,7 +590,7 @@ template<typename T, typename T2=std::complex<T>> class Helper
public:
const T2 *p0r;
T2 *p0w;
- std::vector<T> kernel;
+ vector<T> kernel;
Helper(const GridderConfig<T> &gconf_, const T2 *grid_r_, T2 *grid_w_)
: gconf(gconf_), nu(gconf.Nu()), nv(gconf.Nv()), nsafe(gconf.Nsafe()),
@@ -775,4 +634,6 @@ template<typename T, typename T2=std::complex<T>> class Helper
}
};
+} // namespace gridder
+
#endif
diff --git a/nifty_gridder.cc b/nifty_gridder.cc
index fccb64f2b3cf19c5aa7082e3072b3e3f7c6af274..5afa07207aaf45cbcccb630dc82d09889f0205d3 100644
--- a/nifty_gridder.cc
+++ b/nifty_gridder.cc
@@ -25,6 +25,7 @@
#include "gridder_cxx.h"
using namespace std;
+using namespace gridder;
namespace py = pybind11;
@@ -39,15 +40,27 @@ auto None = py::none();
//
// Start of real gridder functionality
//
+constexpr auto set_nthreads_DS = R"""(
+Specifies the number of threads to be used by the module
+
+Parameters
+==========
+nthreads: int
+ the number of threads to be used. Must be >=1.
+)""";
+constexpr auto get_nthreads_DS = R"""(
+Returns the number of threads used by the module
+
+Returns
+=======
+int : the number of threads used by the module
+)""";
template<typename T>
using pyarr = py::array_t<T>;
-// The "_c" suffix here stands for "C memory order, contiguous"
-template<typename T>
- using pyarr_c = py::array_t<T, py::array::c_style | py::array::forcecast>;
-template<typename T> pyarr_c<T> makeArray(const vector<size_t> &shape)
- { return pyarr_c<T>(shape); }
+template<typename T> pyarr<T> makeArray(const vector<size_t> &shape)
+ { return pyarr<T>(shape); }
void checkArray(const py::array &arr, const char *aname,
const vector<size_t> &shape)
@@ -94,7 +107,7 @@ template<typename T> pyarr<T> providePotentialArray(const py::object &in,
return tmp_;
}
-template<typename T> pyarr_c<T> provideCArray(py::object &in,
+template<typename T> pyarr<T> provideCArray(py::object &in,
const vector<size_t> &shape)
{
if (in.is_none())
@@ -106,19 +119,12 @@ template<typename T> pyarr_c<T> provideCArray(py::object &in,
tmp[i] = T(0);
return tmp_;
}
- auto tmp_ = in.cast<pyarr_c<T>>();
+ auto tmp_ = in.cast<pyarr<T>>();
checkArray(tmp_, "temporary", shape);
return tmp_;
}
-template<size_t ndim, typename T> cmav<T,ndim> make_cmav(pyarr_c<T> &in)
- {
- myassert(ndim==in.ndim(), "dimension mismatch");
- array<size_t,ndim> dims;
- for (size_t i=0; i<ndim; ++i) dims[i]=in.shape(i);
- return cmav<T, ndim>(in.mutable_data(),dims);
- }
-template<size_t ndim, typename T> smav<T,ndim> make_smav(pyarr<T> &in)
+template<size_t ndim, typename T> mav<T,ndim> make_mav(pyarr<T> &in)
{
myassert(ndim==in.ndim(), "dimension mismatch");
array<size_t,ndim> dims;
@@ -129,9 +135,9 @@ template<size_t ndim, typename T> smav<T,ndim> make_smav(pyarr<T> &in)
str[i]=in.strides(i)/sizeof(T);
myassert(str[i]*ptrdiff_t(sizeof(T))==in.strides(i), "weird strides");
}
- return smav<T, ndim>(in.mutable_data(),dims,str);
+ return mav<T, ndim>(in.mutable_data(),dims,str);
}
-template<size_t ndim, typename T> smav<T,ndim> make_smav(pyarr_c<T> &in)
+template<size_t ndim, typename T> mav<const T,ndim> make_const_mav(const pyarr<T> &in)
{
myassert(ndim==in.ndim(), "dimension mismatch");
array<size_t,ndim> dims;
@@ -142,77 +148,7 @@ template<size_t ndim, typename T> smav<T,ndim> make_smav(pyarr_c<T> &in)
str[i]=in.strides(i)/sizeof(T);
myassert(str[i]*ptrdiff_t(sizeof(T))==in.strides(i), "weird strides");
}
- return smav<T, ndim>(in.mutable_data(),dims,str);
- }
-template<size_t ndim, typename T> cmav<const T,ndim> make_const_cmav(const pyarr_c<T> &in)
- {
- myassert(ndim==in.ndim(), "dimension mismatch");
- array<size_t,ndim> dims;
- for (size_t i=0; i<ndim; ++i) dims[i]=in.shape(i);
- return cmav<const T, ndim>(in.data(),dims);
- }
-template<size_t ndim, typename T> smav<const T,ndim> make_const_smav(const pyarr<T> &in)
- {
- myassert(ndim==in.ndim(), "dimension mismatch");
- array<size_t,ndim> dims;
- array<ptrdiff_t,ndim> str;
- for (size_t i=0; i<ndim; ++i)
- {
- dims[i]=in.shape(i);
- str[i]=in.strides(i)/sizeof(T);
- myassert(str[i]*ptrdiff_t(sizeof(T))==in.strides(i), "weird strides");
- }
- return smav<const T, ndim>(in.data(),dims,str);
- }
-template<size_t ndim, typename T> smav<const T,ndim> make_const_smav(const pyarr_c<T> &in)
- {
- myassert(ndim==in.ndim(), "dimension mismatch");
- array<size_t,ndim> dims;
- array<ptrdiff_t,ndim> str;
- for (size_t i=0; i<ndim; ++i)
- {
- dims[i]=in.shape(i);
- str[i]=in.strides(i)/sizeof(T);
- myassert(str[i]*ptrdiff_t(sizeof(T))==in.strides(i), "weird strides");
- }
- return smav<const T, ndim>(in.data(),dims,str);
- }
-template<typename T> pyarr_c<T> complex2hartley
- (const pyarr_c<complex<T>> &grid_, py::object &grid_in)
- {
- auto grid = make_const_smav<2>(grid_);
- size_t nu=grid.shape(0), nv=grid.shape(1);
-
- auto res = provideCArray<T>(grid_in, {nu, nv});
- auto grid2 = make_smav<2>(res);
- {
- py::gil_scoped_release release;
- complex2hartley(grid, grid2);
- }
- return res;
- }
-
-template<typename T> pyarr_c<complex<T>> hartley2complex
- (const pyarr_c<T> &grid_)
- {
- auto grid = make_const_smav<2>(grid_);
- size_t nu=grid.shape(0), nv=grid.shape(1);
-
- auto res = makeArray<complex<T>>({nu, nv});
- auto grid2 = make_smav<2>(res);
- {
- py::gil_scoped_release release;
- hartley2complex(grid, grid2);
- }
- return res;
- }
-
-template<typename T> void hartley2_2D(const pyarr_c<T> &in, pyarr_c<T> &out)
- {
- auto grid = make_const_smav<2>(in);
- auto grid2 = make_smav<2>(out);
- py::gil_scoped_release release;
- hartley2_2D(grid, grid2);
+ return mav<const T, ndim>(in.data(),dims,str);
}
constexpr auto PyBaselines_DS = R"""(
@@ -236,7 +172,7 @@ template<typename T> class PyBaselines: public Baselines<T>
public:
using Baselines<T>::Baselines;
PyBaselines(const pyarr<T> &coord, const pyarr<T> &freq)
- : Baselines<T>(make_const_smav<2>(coord), make_const_smav<1>(freq))
+ : Baselines<T>(make_const_mav<2>(coord), make_const_mav<1>(freq))
{}
using Baselines<T>::effectiveCoord;
@@ -260,12 +196,12 @@ template<typename T> class PyBaselines: public Baselines<T>
)""";
// using Baselines<T>::effectiveUVW;
- pyarr_c<T> effectiveuvw(const pyarr<uint32_t> &idx_) const
+ pyarr<T> effectiveuvw(const pyarr<uint32_t> &idx_) const
{
size_t nvis = size_t(idx_.shape(0));
- auto idx=make_const_smav<1>(idx_);
+ auto idx=make_const_mav<1>(idx_);
auto res_=makeArray<T>({nvis, 3});
- auto res=make_smav<2>(res_);
+ auto res=make_mav<2>(res_);
{
py::gil_scoped_release release;
Baselines<T>::effectiveUVW(idx,res);
@@ -273,14 +209,14 @@ template<typename T> class PyBaselines: public Baselines<T>
return res_;
}
- template<typename T2> pyarr_c<T2> ms2vis(const pyarr<T2> &ms_,
- const pyarr_c<uint32_t> &idx_) const
+ template<typename T2> pyarr<T2> ms2vis(const pyarr<T2> &ms_,
+ const pyarr<uint32_t> &idx_) const
{
- auto idx=make_const_smav<1>(idx_);
+ auto idx=make_const_mav<1>(idx_);
size_t nvis = size_t(idx.shape(0));
- auto ms = make_const_smav<2>(ms_);
+ auto ms = make_const_mav<2>(ms_);
auto res=makeArray<T2>({nvis});
- auto vis = make_smav<1>(res);
+ auto vis = make_mav<1>(res);
{
py::gil_scoped_release release;
Baselines<T>::ms2vis(ms, idx, vis);
@@ -305,13 +241,13 @@ template<typename T> class PyBaselines: public Baselines<T>
np.array((nrows, nchannels), dtype=np.complex)
the measurement set's visibility data (0 where not covered by idx)
)""";
- template<typename T2> pyarr_c<T2> vis2ms(const pyarr<T2> &vis_,
+ template<typename T2> pyarr<T2> vis2ms(const pyarr<T2> &vis_,
const pyarr<uint32_t> &idx_, py::object &ms_in) const
{
- auto vis=make_const_smav<1>(vis_);
- auto idx=make_const_smav<1>(idx_);
+ auto vis=make_const_mav<1>(vis_);
+ auto idx=make_const_mav<1>(idx_);
auto res = provideArray<T2>(ms_in, {nrows, nchan});
- auto ms = make_smav<2>(res);
+ auto ms = make_mav<2>(res);
{
py::gil_scoped_release release;
Baselines<T>::vis2ms(vis, idx, ms);
@@ -320,34 +256,6 @@ template<typename T> class PyBaselines: public Baselines<T>
}
};
-constexpr auto grid2dirty_DS = R"""(
-Converts from UV grid to dirty image (FFT, cropping, correction)
-
-Parameters
-==========
-grid: np.array((nu, nv), dtype=np.float64)
- gridded UV data
-
-Returns
-=======
-nd.array((nxdirty, nydirty), dtype=np.float64)
- the dirty image
-)""";
-
-constexpr auto dirty2grid_DS = R"""(
-Converts from a dirty image to a UV grid (correction, padding, FFT)
-
-Parameters
-==========
-dirty: nd.array((nxdirty, nydirty), dtype=np.float64)
- the dirty image
-
-Returns
-=======
-np.array((nu, nv), dtype=np.float64)
- gridded UV data
-)""";
-
constexpr auto apply_taper_DS = R"""(
Applies the taper (or its inverse) to an image
@@ -432,86 +340,30 @@ template<typename T> class PyGridderConfig: public GridderConfig<T>
using GridderConfig<T>::Nsafe;
using GridderConfig<T>::Beta;
- pyarr_c<T> grid2dirty(const pyarr_c<T> &grid) const
+ pyarr<T> apply_taper(const pyarr<T> &img, bool divide) const
{
- checkArray(grid, "grid", {nu, nv});
auto res = makeArray<T>({nx_dirty, ny_dirty});
- auto grid2=make_const_smav<2>(grid);
- auto res2=make_smav<2>(res);
+ auto img2 = make_const_mav<2>(img);
+ auto res2 = make_mav<2>(res);
{
py::gil_scoped_release release;
- GridderConfig<T>::grid2dirty(grid2,res2);
+ GridderConfig<T>::apply_taper(img2, res2, divide);
}
return res;
}
- pyarr_c<T> apply_taper(const pyarr_c<T> &img, bool divide) const
+ pyarr<complex<T>> grid2dirty_c(const pyarr<complex<T>> &grid) const
{
- checkArray(img, "img", {nx_dirty, ny_dirty});
- auto pin = img.data();
- auto res = makeArray<T>({nx_dirty, ny_dirty});
- auto pout = res.mutable_data();
- {
- py::gil_scoped_release release;
- if (divide)
- for (size_t i=0; i<nx_dirty; ++i)
- for (size_t j=0; j<ny_dirty; ++j)
- pout[ny_dirty*i + j] = pin[ny_dirty*i + j]/(cfu[i]*cfv[j]);
- else
- for (size_t i=0; i<nx_dirty; ++i)
- for (size_t j=0; j<ny_dirty; ++j)
- pout[ny_dirty*i + j] = pin[ny_dirty*i + j]*cfu[i]*cfv[j];
- }
- return res;
- }
- pyarr_c<complex<T>> grid2dirty_c(const pyarr_c<complex<T>> &grid) const
- {
- checkArray(grid, "grid", {nu, nv});
- auto tmp = makeArray<complex<T>>({nu, nv});
- auto ptmp = tmp.mutable_data();
- pocketfft::c2c({nu,nv},{grid.strides(0),grid.strides(1)},
- {tmp.strides(0), tmp.strides(1)}, {0,1}, pocketfft::BACKWARD,
- grid.data(), tmp.mutable_data(), T(1), nthreads);
auto res = makeArray<complex<T>>({nx_dirty, ny_dirty});
- auto pout = res.mutable_data();
+ auto grid2=make_const_mav<2>(grid);
+ auto res2=make_mav<2>(res);
{
py::gil_scoped_release release;
- for (size_t i=0; i<nx_dirty; ++i)
- for (size_t j=0; j<ny_dirty; ++j)
- {
- size_t i2 = nu-nx_dirty/2+i;
- if (i2>=nu) i2-=nu;
- size_t j2 = nv-ny_dirty/2+j;
- if (j2>=nv) j2-=nv;
- pout[ny_dirty*i + j] = ptmp[nv*i2+j2]*cfu[i]*cfv[j];
- }
+ GridderConfig<T>::grid2dirty_c(grid2,res2);
}
return res;
}
- pyarr_c<T> dirty2grid(const pyarr_c<T> &dirty) const
- {
- checkArray(dirty, "dirty", {nx_dirty, ny_dirty});
- auto pdirty = dirty.data();
- auto tmp = makeArray<T>({nu, nv});
- auto ptmp = tmp.mutable_data();
- {
- py::gil_scoped_release release;
- for (size_t i=0; i<nu*nv; ++i)
- ptmp[i] = 0.;
- for (size_t i=0; i<nx_dirty; ++i)
- for (size_t j=0; j<ny_dirty; ++j)
- {
- size_t i2 = nu-nx_dirty/2+i;
- if (i2>=nu) i2-=nu;
- size_t j2 = nv-ny_dirty/2+j;
- if (j2>=nv) j2-=nv;
- ptmp[nv*i2+j2] = pdirty[ny_dirty*i + j]*cfu[i]*cfv[j];
- }
- }
- hartley2_2D<T>(tmp, tmp);
- return tmp;
- }
- pyarr_c<complex<T>> dirty2grid_c(const pyarr_c<complex<T>> &dirty) const
+ pyarr<complex<T>> dirty2grid_c(const pyarr<complex<T>> &dirty) const
{
checkArray(dirty, "dirty", {nx_dirty, ny_dirty});
auto pdirty = dirty.data();
@@ -536,7 +388,7 @@ template<typename T> class PyGridderConfig: public GridderConfig<T>
}
return tmp;
}
- pyarr_c<complex<T>> apply_wscreen(const pyarr_c<complex<T>> &dirty_, double w, bool adjoint) const
+ pyarr<complex<T>> apply_wscreen(const pyarr<complex<T>> &dirty_, double w, bool adjoint) const
{
checkArray(dirty_, "dirty", {nx_dirty, ny_dirty});
auto dirty = dirty_.data();
@@ -600,7 +452,7 @@ Returns
np.array((nu,nv), dtype=np.complex128):
the gridded visibilities
)""";
-template<typename T> pyarr_c<complex<T>> vis2grid_c(
+template<typename T> pyarr<complex<T>> vis2grid_c(
const PyBaselines<T> &baselines, const PyGridderConfig<T> &gconf,
const pyarr<uint32_t> &idx_, const pyarr<complex<T>> &vis_,
py::object &grid_in, const py::object &wgt_)
@@ -677,10 +529,6 @@ Returns
np.array((nu,nv), dtype=np.float64):
the gridded visibilities (made real by making use of Hermitian symmetry)
)""";
-template<typename T> pyarr_c<T> vis2grid(const PyBaselines<T> &baselines,
- const PyGridderConfig<T> &gconf, const pyarr<uint32_t> &idx_,
- const pyarr<complex<T>> &vis_, py::object &grid_in, const py::object &wgt_)
- { return complex2hartley(vis2grid_c(baselines, gconf, idx_, vis_, None, wgt_), grid_in); }
constexpr auto ms2grid_c_DS = R"""(
Grids measurement set data onto a UV grid
@@ -706,7 +554,7 @@ Returns
np.array((nu,nv), dtype=np.complex128):
the gridded visibilities
)""";
-template<typename T> pyarr_c<complex<T>> ms2grid_c(
+template<typename T> pyarr<complex<T>> ms2grid_c(
const PyBaselines<T> &baselines, const PyGridderConfig<T> &gconf,
const pyarr<uint32_t> &idx_, const pyarr<complex<T>> &ms_,
py::object &grid_in, const py::object &wgt_)
@@ -764,14 +612,9 @@ template<typename T> pyarr_c<complex<T>> ms2grid_c(
return res;
}
-template<typename T> pyarr_c<T> ms2grid(const PyBaselines<T> &baselines,
- const PyGridderConfig<T> &gconf, const pyarr<uint32_t> &idx_,
- const pyarr<complex<T>> &ms_, py::object &grid_in, const py::object &wgt_)
- { return complex2hartley(ms2grid_c(baselines, gconf, idx_, ms_, None, wgt_), grid_in); }
-
-template<typename T> pyarr_c<complex<T>> grid2vis_c(
+template<typename T> pyarr<complex<T>> grid2vis_c(
const PyBaselines<T> &baselines, const PyGridderConfig<T> &gconf,
- const pyarr<uint32_t> &idx_, const pyarr_c<complex<T>> &grid_,
+ const pyarr<uint32_t> &idx_, const pyarr<complex<T>> &grid_,
const py::object &wgt_)
{
size_t nu=gconf.Nu(), nv=gconf.Nv();
@@ -845,14 +688,9 @@ Returns
np.array((nvis,), dtype=np.complex)
The degridded visibility data
)""";
-template<typename T> pyarr_c<complex<T>> grid2vis(const PyBaselines<T> &baselines,
- const PyGridderConfig<T> &gconf, const pyarr<uint32_t> &idx_,
- const pyarr_c<T> &grid_, const py::object &wgt_)
- { return grid2vis_c(baselines, gconf, idx_, hartley2complex(grid_), wgt_); }
-
-template<typename T> pyarr_c<complex<T>> grid2ms_c(const PyBaselines<T> &baselines,
+template<typename T> pyarr<complex<T>> grid2ms_c(const PyBaselines<T> &baselines,
const PyGridderConfig<T> &gconf, const pyarr<uint32_t> &idx_,
- const pyarr_c<complex<T>> &grid_, py::object &ms_in, const py::object &wgt_)
+ const pyarr<complex<T>> &grid_, py::object &ms_in, const py::object &wgt_)
{
auto nrows = baselines.Nrows();
auto nchan = baselines.Nchannels();
@@ -910,14 +748,9 @@ template<typename T> pyarr_c<complex<T>> grid2ms_c(const PyBaselines<T> &baselin
return res;
}
-template<typename T> pyarr_c<complex<T>> grid2ms(const PyBaselines<T> &baselines,
- const PyGridderConfig<T> &gconf, const pyarr<uint32_t> &idx_,
- const pyarr_c<T> &grid_, py::object &ms_in, const py::object &wgt_)
- { return grid2ms_c(baselines, gconf, idx_, hartley2complex(grid_), ms_in, wgt_); }
-
-template<typename T> pyarr_c<complex<T>> apply_holo(
+template<typename T> pyarr<complex<T>> apply_holo(
const PyBaselines<T> &baselines, const PyGridderConfig<T> &gconf,
- const pyarr<uint32_t> &idx_, const pyarr_c<complex<T>> &grid_,
+ const pyarr<uint32_t> &idx_, const pyarr<complex<T>> &grid_,
const py::object &wgt_)
{
size_t nu=gconf.Nu(), nv=gconf.Nv();
@@ -984,7 +817,7 @@ template<typename T> pyarr_c<complex<T>> apply_holo(
return res;
}
-template<typename T> pyarr_c<T> get_correlations(
+template<typename T> pyarr<T> get_correlations(
const PyBaselines<T> &baselines, const PyGridderConfig<T> &gconf,
const pyarr<uint32_t> &idx_, int du, int dv, const py::object &wgt_)
{
@@ -1077,8 +910,8 @@ np.array((nvis,), dtype=np.uint32)
the compressed indices for all entries which match the selected criteria
and are not flagged.
)""";
-template<typename T> pyarr_c<uint32_t> getIndices(const PyBaselines<T> &baselines,
- const PyGridderConfig<T> &gconf, const pyarr_c<bool> &flags_, int chbegin,
+template<typename T> pyarr<uint32_t> getIndices(const PyBaselines<T> &baselines,
+ const PyGridderConfig<T> &gconf, const pyarr<bool> &flags_, int chbegin,
int chend, T wmin, T wmax)
{
size_t nrow=baselines.Nrows(),
@@ -1134,7 +967,7 @@ template<typename T> pyarr_c<uint32_t> getIndices(const PyBaselines<T> &baseline
return res;
}
-template<typename T> pyarr_c<complex<T>> vis2dirty_wstack(const PyBaselines<T> &baselines,
+template<typename T> pyarr<complex<T>> vis2dirty_wstack(const PyBaselines<T> &baselines,
const PyGridderConfig<T> &gconf, const pyarr<uint32_t> &idx_,
const pyarr<complex<T>> &vis_)
{
@@ -1258,7 +1091,7 @@ PYBIND11_MODULE(nifty_gridder, m)
using namespace pybind11::literals;
py::class_<PyBaselines<double>> (m, "Baselines", PyBaselines_DS)
- .def(py::init<const pyarr_c<double> &, const pyarr_c<double> &>(),
+ .def(py::init<const pyarr<double> &, const pyarr<double> &>(),
"coord"_a, "freq"_a)
.def ("Nrows",&PyBaselines<double>::Nrows)
.def ("Nchannels",&PyBaselines<double>::Nchannels)
@@ -1280,11 +1113,7 @@ PYBIND11_MODULE(nifty_gridder, m)
.def("Supp", &PyGridderConfig<double>::Supp)
.def("apply_taper", &PyGridderConfig<double>::apply_taper, apply_taper_DS,
"img"_a, "divide"_a=false)
- .def("grid2dirty", &PyGridderConfig<double>::grid2dirty,
- grid2dirty_DS, "grid"_a)
.def("grid2dirty_c", &PyGridderConfig<double>::grid2dirty_c, "grid"_a)
- .def("dirty2grid", &PyGridderConfig<double>::dirty2grid,
- dirty2grid_DS, "dirty"_a)
.def("dirty2grid_c", &PyGridderConfig<double>::dirty2grid_c, "dirty"_a)
.def("apply_wscreen", &PyGridderConfig<double>::apply_wscreen,
apply_wscreen_DS, "dirty"_a, "w"_a, "adjoint"_a)
@@ -1310,14 +1139,6 @@ PYBIND11_MODULE(nifty_gridder, m)
m.def("getIndices", getIndices<double>, getIndices_DS, "baselines"_a,
"gconf"_a, "flags"_a, "chbegin"_a=-1, "chend"_a=-1,
"wmin"_a=-1e30, "wmax"_a=1e30);
- m.def("vis2grid",&vis2grid<double>, vis2grid_DS, "baselines"_a, "gconf"_a,
- "idx"_a, "vis"_a, "grid_in"_a=None, "wgt"_a=None);
- m.def("ms2grid",&ms2grid<double>, "baselines"_a, "gconf"_a, "idx"_a, "ms"_a,
- "grid_in"_a=None, "wgt"_a=None);
- m.def("grid2vis",&grid2vis<double>, grid2vis_DS, "baselines"_a, "gconf"_a,
- "idx"_a, "grid"_a, "wgt"_a=None);
- m.def("grid2ms",&grid2ms<double>, "baselines"_a, "gconf"_a, "idx"_a,
- "grid"_a, "ms_in"_a=None, "wgt"_a=None);
m.def("vis2grid_c",&vis2grid_c<double>, vis2grid_c_DS, "baselines"_a,
"gconf"_a, "idx"_a, "vis"_a, "grid_in"_a=None, "wgt"_a=None);
m.def("ms2grid_c",&ms2grid_c<double>, ms2grid_c_DS, "baselines"_a, "gconf"_a,
diff --git a/test.py b/test.py
index 801d363202a31f8714b3021acce097890a6dfcab..bac5fb4f20b576cb5ff708e524462997cc31c5db 100644
--- a/test.py
+++ b/test.py
@@ -12,7 +12,7 @@ def _init_gridder(nxdirty, nydirty, epsilon, nchan, nrow):
nydirty=nydirty,
epsilon=epsilon,
pixsize_x=pixsize,
- pixsize_y=pixsize)
+ pixsize_y=0.5*pixsize)
speedoflight, f0 = 3e8, 1e9
freq = f0 + np.arange(nchan)*(f0/nchan)
uvw = (np.random.rand(nrow, 3)-0.5)/(pixsize*f0/speedoflight)
@@ -45,11 +45,11 @@ def test_adjointness(nxdirty, nydirty, nrow, nchan, epsilon):
bl, gconf, idx = _init_gridder(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)
- grid = ng.vis2grid(bl, gconf, idx, vis)
- dirty = gconf.grid2dirty(grid)
+ grid = ng.vis2grid_c(bl, gconf, idx, vis)
+ dirty = gconf.grid2dirty_c(grid)
dirty2 = np.random.rand(*dirty.shape)-0.5
- ms2 = bl.vis2ms(ng.grid2vis(bl, gconf, idx, gconf.dirty2grid(dirty2)), idx)
- assert_allclose(np.vdot(ms, ms2).real, np.vdot(dirty, dirty2))
+ ms2 = bl.vis2ms(ng.grid2vis_c(bl, gconf, idx, gconf.dirty2grid_c(dirty2)), idx)
+ assert_allclose(np.vdot(ms, ms2).real, np.vdot(dirty.real, dirty2.real))
@pmp("nxdirty", (128,))
@@ -78,38 +78,6 @@ def test_hoisted_grid_allocation(nxdirty, nydirty, nrow, nchan, epsilon):
ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5)
vis = baselines.ms2vis(ms, idx)
- # -----------------------------
- # Test real grid case (ms2grid)
- # -----------------------------
- real_grid = np.zeros((gconf.Nu(), gconf.Nv()), dtype=np.float64)
- grid = ng.ms2grid(baselines, gconf, idx, ms, grid_in=None)
- grid2 = ng.ms2grid(baselines, gconf, idx, ms, grid_in=real_grid)
-
- assert_array_almost_equal(grid, grid2)
- assert grid.dtype == real_grid.dtype == grid2.dtype
- assert id(grid2) == id(real_grid) # Same grid object
-
- real_grid = np.zeros((gconf.Nu(), gconf.Nv()), dtype=np.float64)
- grid = ng.ms2grid(baselines, gconf, idx, ms, grid_in=None, wgt=weights)
- grid2 = ng.ms2grid(baselines, gconf, idx, ms, grid_in=real_grid,
- wgt=weights)
-
- assert_array_almost_equal(grid, grid2)
- assert grid.dtype == real_grid.dtype == grid2.dtype
- assert id(grid2) == id(real_grid) # Same grid object
-
- # ------------------------------
- # Test real grid case (vis2grid)
- # ------------------------------
- real_grid = np.zeros((gconf.Nu(), gconf.Nv()), dtype=np.float64)
- grid = ng.vis2grid(baselines, gconf, idx, vis, grid_in=None)
- grid2 = ng.vis2grid(baselines, gconf, idx, vis, grid_in=real_grid)
-
- # Almost same result
- assert_array_almost_equal(grid, grid2)
- assert grid.dtype == real_grid.dtype == grid2.dtype
- assert id(grid2) == id(real_grid) # Same grid object
-
# ----------------------------------
# Test complex grid case (ms2grid_c)
# ----------------------------------
@@ -248,7 +216,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)
- res0 = conf.grid2dirty(ng.vis2grid(bl, conf, idx, vis))
+ res0 = conf.grid2dirty_c(ng.vis2grid_c(bl, conf, idx, vis))
x, y = np.meshgrid(*[-ss/2 + np.arange(ss) for ss in [nxdirty, nydirty]],
indexing='ij')
@@ -258,7 +226,7 @@ def test_against_dft(nxdirty, nydirty, epsilon, nchan, nrow):
uvw = bl.effectiveuvw(idx)
for ii in idx:
phase = x*uvw[ii, 0] + y*uvw[ii, 1]
- res1 += (vis[ii]*np.exp(2j*np.pi*phase)).real
+ res1 += (vis[ii]*np.exp(2j*np.pi*phase))
assert_(_l2error(res0, res1) < epsilon)
@@ -275,7 +243,7 @@ def test_against_wdft(nxdirty, nydirty, nchan, nrow, fov):
nydirty=nydirty,
epsilon=epsilon,
pixsize_x=pixsize,
- pixsize_y=pixsize)
+ pixsize_y=0.5*pixsize)
speedoflight, f0 = 3e8, 1e9
freq = f0 + np.arange(nchan)*(f0/nchan)
uvw = (np.random.rand(nrow, 3)-0.5)/(pixsize*f0/speedoflight)
@@ -314,4 +282,56 @@ def test_against_wdft(nxdirty, nydirty, nchan, nrow, fov):
phase = x*uvw[ii, 0] + y*uvw[ii, 1] + uvw[ii, 2]*nm1
res1 += (vis[ii]*np.exp(2j*np.pi*phase)).real
res1 /= n
+ print(_l2error(res0, res1))
+ assert_(_l2error(res0, res1) < 1e-4)
+
+@pmp('nxdirty', [16, 64])
+@pmp('nydirty', [64])
+@pmp("nrow", (10, 100, 1000))
+@pmp("nchan", (1,))
+@pmp("fov", (50,))
+def test_against_wdft_exact(nxdirty, nydirty, nchan, nrow, fov):
+ epsilon = 1e-10
+ np.random.seed(40)
+ pixsize = fov*np.pi/180/nxdirty
+ conf = ng.GridderConfig(nxdirty=nxdirty,
+ nydirty=nydirty,
+ epsilon=epsilon,
+ pixsize_x=pixsize,
+ pixsize_y=0.2*pixsize)
+ speedoflight, f0 = 3e8, 1e9
+ freq = f0 + np.arange(nchan)*(f0/nchan)
+ uvw = (np.random.rand(nrow, 3)-0.5)/(pixsize*f0/speedoflight)
+ # use exactly one random w
+ uvw[:,2] = uvw[0,2]
+ bl = ng.Baselines(coord=uvw, freq=freq)
+ flags = np.zeros((nrow, nchan), dtype=np.bool)
+ idx = ng.getIndices(bl, conf, flags)
+ 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)
+
+ mi, ma = np.min(uvw[:, 2]), np.max(uvw[:, 2])
+ nplanes = 1
+ w = mi
+ 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).real
+
+ # Compute dft with w term
+ x, y = np.meshgrid(*[-ss/2 + np.arange(ss) for ss in [nxdirty, nydirty]],
+ indexing='ij')
+ x *= conf.Pixsize_x()
+ y *= conf.Pixsize_y()
+ res1 = np.zeros_like(res0)
+
+ eps = np.sqrt(x**2+y**2)
+ s = np.sin(eps)
+ nm1 = -s*s/(1+np.cos(eps))
+ n = nm1+1
+ for ii in idx:
+ phase = x*uvw[ii, 0] + y*uvw[ii, 1] + uvw[ii, 2]*nm1
+ res1 += (vis[ii]*np.exp(2j*np.pi*phase)).real
+ res1 /= n
+ print(_l2error(res0, res1))
assert_(_l2error(res0, res1) < 1e-4)