Skip to content

GitLab

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

tweaks

parent 18c0cba6
Hide whitespace changes
Inline Side-by-side
gridder_cxx.h
View file @ c7d0c24a
......@@ -30,6 +30,7 @@
#include "pocketfft_hdronly.h"
#include <array>
#define NOINLINE __attribute__((noinline))
namespace gridder {
......@@ -279,7 +280,7 @@ template<typename T> class EC_Kernel_with_correction: public EC_Kernel<T>
T tmp=0;
for (int i=0; i<p; ++i)
tmp += wgt[i]*psi[i]*cos(pi*supp*v*x[i]);
return T(1)/(supp*tmp);
return T(1./(supp*tmp));
}
};
/* Compute correction factors for the ES gridding kernel
......@@ -1003,7 +1004,6 @@ template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &base
wmin = min(wmin,wval[ipart]);
wmax = max(wmax,wval[ipart]);
}
cout << "data w range: " << wmin << " to " << wmax << endl;
double x0 = -0.5*nx_dirty*psx,
y0 = -0.5*ny_dirty*psy;
......@@ -1017,7 +1017,6 @@ cout << "data w range: " << wmin << " to " << wmax << endl;
wmin -= (0.5*w_supp-1)*dw;
wmax += (0.5*w_supp-1)*dw;
nplanes += w_supp-2;
cout << "nplanes: " << nplanes << endl;
vector<size_t> nvis_plane(nplanes,0);
vector<int> minplane(nvis);
for (size_t ipart=0; ipart<nvis; ++ipart)
......@@ -1040,8 +1039,6 @@ cout << "nplanes: " << nplanes << endl;
auto tdirty=tdirty_.getMav();
for (size_t iw=0; iw<nplanes; ++iw)
{
cout << "working on w plane #" << iw << " containing " << nvis_plane[iw]
<< " visibilities" << endl;
if (nvis_plane[iw]==0) continue;
double wcur = wmin+iw*dw;
size_t cnt=0;
......@@ -1062,14 +1059,13 @@ myassert(cnt==nvis_plane[iw],"must not happen 2");
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.apply_wscreen(tdirty, tdirty, wcur, false);
#pragma omp parallel for num_threads(nthreads)
for (size_t i=0; i<nx_dirty; ++i)
for (size_t j=0; j<ny_dirty; ++j)
dirty(i,j) += tdirty(i,j);
}
}
// correct for w gridding
cout << "applying correction for gridding in w direction" << endl;
apply_wcorr(gconf, dirty, kernel, dw);
}
template<typename T> void vis2dirty_wstack(const Baselines<T> &baselines,
......@@ -1101,7 +1097,6 @@ template<typename T, typename Serv> void dirty2x_wstack(const Baselines<T> &base
wmin = min(wmin,wval[ipart]);
wmax = max(wmax,wval[ipart]);
}
//cout << "data w range: " << wmin << " to " << wmax << endl;
double x0 = -0.5*nx_dirty*psx,
y0 = -0.5*ny_dirty*psy;
......@@ -1115,7 +1110,6 @@ template<typename T, typename Serv> void dirty2x_wstack(const Baselines<T> &base
wmin -= (0.5*w_supp-1)*dw;
wmax += (0.5*w_supp-1)*dw;
nplanes += w_supp-2;
cout << "nplanes: " << nplanes << endl;
vector<size_t> nvis_plane(nplanes,0);
vector<int> minplane(nvis);
for (size_t ipart=0; ipart<nvis; ++ipart)
......@@ -1135,7 +1129,6 @@ cout << "nplanes: " << nplanes << endl;
for (size_t j=0; j<ny_dirty; ++j)
tdirty(i,j) = dirty(i,j);
// correct for w gridding
//cout << "applying correction for gridding in w direction" << endl;
apply_wcorr(gconf, tdirty, kernel, dw);
tmpStorage<complex<T>,2> grid_({nu,nv});
......@@ -1146,8 +1139,6 @@ cout << "nplanes: " << nplanes << endl;
vector<uint32_t> idx_loc_;
for (size_t iw=0; iw<nplanes; ++iw)
{
cout << "working on w plane #" << iw << " containing " << nvis_plane[iw]
<< " visibilities" << endl;
if (nvis_plane[iw]==0) continue;
double wcur = wmin+iw*dw;
gconf.apply_wscreen(tdirty, tdirty2, wcur, true);
......@@ -1222,8 +1213,9 @@ template<typename T> vector<uint32_t> getWgtIndices(const Baselines<T> &baseline
template<typename T> void full_gridding(const const_mav<T,2> &uvw,
const const_mav<T,1> &freq, const const_mav<complex<T>,2> &ms,
const const_mav<T,2> &wgt, T pixsize_x, T pixsize_y, double epsilon,
const mav<complex<T>,2> &dirty)
size_t nthreads, const mav<complex<T>,2> &dirty)
{
set_nthreads(nthreads);
Baselines<T> baselines(uvw, freq);
GridderConfig<T> gconf(dirty.shape(0), dirty.shape(1), epsilon, pixsize_x, pixsize_y);
auto idx = getWgtIndices(baselines, gconf, wgt);
......@@ -1235,8 +1227,9 @@ template<typename T> void full_gridding(const const_mav<T,2> &uvw,
template<typename T> void full_degridding(const const_mav<T,2> &uvw,
const const_mav<T,1> &freq, const const_mav<complex<T>,2> &dirty,
const const_mav<T,2> &wgt, T pixsize_x, T pixsize_y, double epsilon,
const mav<complex<T>,2> &ms)
size_t nthreads, const mav<complex<T>,2> &ms)
{
set_nthreads(nthreads);
Baselines<T> baselines(uvw, freq);
GridderConfig<T> gconf(dirty.shape(0), dirty.shape(1), epsilon, pixsize_x, pixsize_y);
auto idx = getWgtIndices(baselines, gconf, wgt);
......
nifty_gridder.cc
View file @ c7d0c24a
......@@ -706,7 +706,7 @@ template<typename T> pyarr<complex<T>> Pydirty2vis_wstack(const PyBaselines<T> &
template<typename T> pyarr<complex<T>> Pyfull_gridding(const pyarr<T> &uvw,
const pyarr<T> &freq, const pyarr<complex<T>> &ms, const py::object &wgt_,
size_t npix_x, size_t npix_y, T pixsize_x, T pixsize_y, double epsilon)
size_t npix_x, size_t npix_y, T pixsize_x, T pixsize_y, double epsilon, size_t nthreads)
{
auto uvw2 = make_const_mav<2>(uvw);
auto freq2 = make_const_mav<1>(freq);
......@@ -715,13 +715,13 @@ template<typename T> pyarr<complex<T>> Pyfull_gridding(const pyarr<T> &uvw,
auto wgt2 = make_const_mav<2>(wgt);
auto dirty=makeArray<complex<T>>({npix_x,npix_y});
auto dirty2 =make_mav<2>(dirty);
full_gridding(uvw2,freq2,ms2,wgt2,pixsize_x,pixsize_y,epsilon,dirty2);
full_gridding(uvw2,freq2,ms2,wgt2,pixsize_x,pixsize_y,epsilon,nthreads,dirty2);
return dirty;
}
template<typename T> pyarr<complex<T>> Pyfull_degridding(const pyarr<T> &uvw,
const pyarr<T> &freq, const pyarr<complex<T>> &dirty, const py::object &wgt_,
T pixsize_x, T pixsize_y, double epsilon)
T pixsize_x, T pixsize_y, double epsilon, size_t nthreads)
{
auto uvw2 = make_const_mav<2>(uvw);
auto freq2 = make_const_mav<1>(freq);
......@@ -730,7 +730,7 @@ template<typename T> pyarr<complex<T>> Pyfull_degridding(const pyarr<T> &uvw,
auto wgt2 = make_const_mav<2>(wgt);
auto ms=makeArray<complex<T>>({uvw2.shape(0),freq2.shape(0)});
auto ms2 =make_mav<2>(ms);
full_degridding(uvw2,freq2,dirty2,wgt2,pixsize_x,pixsize_y,epsilon,ms2);
full_degridding(uvw2,freq2,dirty2,wgt2,pixsize_x,pixsize_y,epsilon,nthreads,ms2);
return ms;
}
......@@ -809,11 +809,11 @@ PYBIND11_MODULE(nifty_gridder, m)
"idx"_a, "dirty"_a);
m.def("full_gridding",&Pyfull_gridding<double>,"uvw"_a,"freq"_a,"vis"_a,
"wgt"_a,"npix_x"_a,"npix_y"_a,"pixsize_x"_a,"pixsize_y"_a,"epsilon"_a);
"wgt"_a,"npix_x"_a,"npix_y"_a,"pixsize_x"_a,"pixsize_y"_a,"epsilon"_a,"nthreads"_a=1);
m.def("full_gridding_f",&Pyfull_gridding<float>,"uvw"_a,"freq"_a,"vis"_a,
"wgt"_a,"npix_x"_a,"npix_y"_a,"pixsize_x"_a,"pixsize_y"_a,"epsilon"_a);
"wgt"_a,"npix_x"_a,"npix_y"_a,"pixsize_x"_a,"pixsize_y"_a,"epsilon"_a,"nthreads"_a=1);
m.def("full_degridding",&Pyfull_degridding<double>,"uvw"_a,"freq"_a,"dirty"_a,
"wgt"_a,"pixsize_x"_a,"pixsize_y"_a, "epsilon"_a);
"wgt"_a,"pixsize_x"_a,"pixsize_y"_a, "epsilon"_a,"nthreads"_a=1);
m.def("full_degridding_f",&Pyfull_degridding<float>,"uvw"_a,"freq"_a,"dirty"_a,
"wgt"_a,"pixsize_x"_a,"pixsize_y"_a, "epsilon"_a);
"wgt"_a,"pixsize_x"_a,"pixsize_y"_a, "epsilon"_a, "nthreads"_a=1);
}
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