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

documentation

parent d0018a6c
......@@ -348,6 +348,28 @@ template<typename T> pyarr_c<complex<T>> hartley2complex
return res;
}
template<typename T> void hartley2_2D(const pyarr_c<T> &in, pyarr_c<T> &out)
{
size_t nu=in.shape(0), nv=in.shape(1);
pocketfft::r2r_hartley({nu, nv},
{in.strides(0), in.strides(1)},
{out.strides(0), out.strides(1)}, {0,1},
in.data(), out.mutable_data(), T(1), 0);
auto ptmp = out.mutable_data();
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);
}
}
/* Compute correction factors for the ES gridding kernel
This implementation follows eqs. (3.8) to (3.10) of Barnett et al. 2018 */
vector<double> correction_factors (size_t n, size_t nval, size_t w)
......@@ -380,6 +402,7 @@ template<typename T> struct UV
UV operator* (T fct) const
{ return UV(u*fct, v*fct); }
};
template<typename T> struct UVW
{
T u, v, w;
......@@ -564,23 +587,7 @@ template<typename T> class GridderConfig
{
auto tmp = makearray<T>({nu, nv});
auto ptmp = tmp.mutable_data();
pocketfft::r2r_hartley({nu, nv},
{grid.strides(0), grid.strides(1)},
{tmp.strides(0), tmp.strides(1)}, {0,1},
grid.data(),
ptmp, T(1), 0);
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);
}
hartley2_2D<T>(grid, tmp);
auto res = makearray<T>({nx_dirty, ny_dirty});
auto pout = res.mutable_data();
for (size_t i=0; i<nx_dirty; ++i)
......@@ -610,22 +617,7 @@ template<typename T> class GridderConfig
if (j2>=nv) j2-=nv;
ptmp[nv*i2+j2] = pdirty[ny_dirty*i + j]*cfu[i]*cfv[j];
}
pocketfft::r2r_hartley({nu, nv},
{tmp.strides(0), tmp.strides(1)},
{tmp.strides(0), tmp.strides(1)}, {0,1},
tmp.data(), tmp.mutable_data(), T(1), 0);
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);
}
hartley2_2D<T>(tmp, tmp);
return tmp;
}
};
......@@ -981,6 +973,119 @@ np.array((nrows, nchannels), dtype=np.complex)
the updated MS
)""";
const char *GridderConfig_DS = R"""(
Class storing information related to the gridding/degridding process.
Parameters
==========
nxdirty: int
x resolution of the dirty image; must be even
nydirty: int
y resolution of the dirty image; must be even
epsilon: float
required accuracy for the gridding/degridding step
Must be >= 2e-13.
urange: float
vrange: float
)""";
const char *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
)""";
const char *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
)""";
const char *getIndices_DS = R"""(
Selects a subset of entries from a `Baselines` object.
Parameters
==========
baselines: Baselines
the Baselines object
gconf: GridderConf
the GridderConf object to be used with the returned indices.
(used to optimize the ordering of the indices)
chbegin: int
first channel to use (-1: start with the first available channel)
chend: int
one-past last channel to use (-1: one past the last available channel)
wmin: float
only select entries with w>=wmin
wmax: float
only select entries with w<wmax
Returns
=======
np.array((nvis,), dtype=np.uint32)
the compressed indices for all entries which match the selected criteria
and are not flagged.
)""";
const char *vis2grid_DS = R"""(
Grids visibilities onto a UV grid
Parameters
==========
baselines: Baselines
the Baselines object
gconf: GridderConf
the GridderConf object to be used
(used to optimize the ordering of the indices)
idx: np.array((nvis,), dtype=np.uint32)
the indices for the entries to be gridded
vis: np.array((nvis,), dtype=np.complex)
The visibility data for the index array
Returns
=======
np.array((nu,nv), dtype=np.float64):
the gridded visibilities (made real by making use of Hermitian symmetry)
)""";
const char *grid2vis_DS = R"""(
Degrids visibilities from a UV grid
Parameters
==========
baselines: Baselines
the Baselines object
gconf: GridderConf
the GridderConf object to be used
(used to optimize the ordering of the indices)
idx: np.array((nvis,), dtype=np.uint32)
the indices for the entries to be degridded
grid: np.array((nu,nv), dtype=np.float64):
the gridded visibilities (made real by making use of Hermitian symmetry)
vis: np.array((nvis,), dtype=np.complex)
The visibility data for the index array
Returns
=======
np.array((nvis,), dtype=np.complex)
The degridded visibility data
)""";
} // unnamed namespace
PYBIND11_MODULE(nifty_gridder, m)
......@@ -995,17 +1100,17 @@ PYBIND11_MODULE(nifty_gridder, m)
.def ("vis2ms",&Baselines<double>::vis2ms<complex<double>>, BL_vis2ms_DS, "vis"_a, "idx"_a)
.def ("add_vis_to_ms",&Baselines<double>::add_vis_to_ms<complex<double>>, BL_add_vis_to_ms_DS,
"vis"_a, "idx"_a, "ms"_a.noconvert());
py::class_<GridderConfig<double>> (m, "GridderConfig")
py::class_<GridderConfig<double>> (m, "GridderConfig", GridderConfig_DS)
.def(py::init<size_t, size_t, double, double, double>(),"nxdirty"_a,
"nydirty"_a, "epsilon"_a, "urange"_a, "vrange"_a)
.def("Nu", &GridderConfig<double>::Nu)
.def("Nv", &GridderConfig<double>::Nv)
.def("grid2dirty", &GridderConfig<double>::grid2dirty, "grid"_a)
.def("dirty2grid", &GridderConfig<double>::dirty2grid, "dirty"_a);
m.def("getIndices", getIndices<double>, "baselines"_a, "gconf"_a, "chbegin"_a=-1,
"chend"_a=-1, "wmin"_a=-1e30, "wmax"_a=1e30);
m.def("vis2grid",&vis2grid<double>, "baselines"_a, "gconf"_a, "idx"_a, "vis"_a);
m.def("grid2vis",&grid2vis<double>, "baselines"_a, "gconf"_a, "idx"_a, "grid"_a);
.def("grid2dirty", &GridderConfig<double>::grid2dirty, grid2dirty_DS, "grid"_a)
.def("dirty2grid", &GridderConfig<double>::dirty2grid, dirty2grid_DS, "dirty"_a);
m.def("getIndices", getIndices<double>, getIndices_DS, "baselines"_a,
"gconf"_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);
m.def("grid2vis",&grid2vis<double>, grid2vis_DS, "baselines"_a, "gconf"_a, "idx"_a, "grid"_a);
py::class_<Baselines<float>> (m, "Baselines_f", Baselines_DS)
.def(py::init<const pyarr_c<float> &, const pyarr_c<float> &,
......
Supports Markdown
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