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

don't break the interface after all

parent eb929f7b
Pipeline #80239 passed with stages
in 14 minutes and 43 seconds
......@@ -1370,14 +1370,19 @@ template<typename T> vector<idx_t> getWgtIndices(const Baselines &baselines,
template<typename T> void ms2dirty(const mav<double,2> &uvw,
const mav<double,1> &freq, const mav<complex<T>,2> &ms,
const mav<T,2> &wgt, double pixsize_x, double pixsize_y, double epsilon,
const mav<T,2> &wgt, double pixsize_x, double pixsize_y, size_t nu, size_t nv, double epsilon,
bool do_wstacking, size_t nthreads, mav<T,2> &dirty, size_t verbosity,
bool negate_v=false)
{
Baselines baselines(uvw, freq, negate_v);
// adjust for increased error when gridding in 2 or 3 dimensions
epsilon /= do_wstacking ? 3 : 2;
auto [nu, nv] = getNuNv(baselines, ms, wgt, epsilon, do_wstacking, dirty.shape(0), dirty.shape(1), pixsize_x, pixsize_y);
if (nu*nv==0)
{
auto [nu2, nv2] = getNuNv(baselines, ms, wgt, epsilon, do_wstacking, dirty.shape(0), dirty.shape(1), pixsize_x, pixsize_y);
nu = nu2;
nv = nv2;
}
GridderConfig<T> gconf(dirty.shape(0), dirty.shape(1), nu, nv, epsilon, pixsize_x, pixsize_y, baselines, nthreads);
auto idx = getWgtIndices(baselines, gconf, wgt, ms);
auto idx2 = mav<idx_t,1>(idx.data(),{idx.size()});
......@@ -1387,15 +1392,20 @@ template<typename T> void ms2dirty(const mav<double,2> &uvw,
template<typename T> void dirty2ms(const mav<double,2> &uvw,
const mav<double,1> &freq, const mav<T,2> &dirty,
const mav<T,2> &wgt, double pixsize_x, double pixsize_y, double epsilon,
bool do_wstacking, size_t nthreads, mav<complex<T>,2> &ms,
const mav<T,2> &wgt, double pixsize_x, double pixsize_y, size_t nu, size_t nv,
double epsilon, bool do_wstacking, size_t nthreads, mav<complex<T>,2> &ms,
size_t verbosity, bool negate_v=false)
{
Baselines baselines(uvw, freq, negate_v);
// adjust for increased error when gridding in 2 or 3 dimensions
epsilon /= do_wstacking ? 3 : 2;
mav<complex<T>,2> null_ms(nullptr, {0,0}, false);
auto [nu, nv] = getNuNv(baselines, null_ms, wgt, epsilon, do_wstacking, dirty.shape(0), dirty.shape(1), pixsize_x, pixsize_y);
if (nu*nv==0)
{
auto [nu2, nv2] = getNuNv(baselines, null_ms, wgt, epsilon, do_wstacking, dirty.shape(0), dirty.shape(1), pixsize_x, pixsize_y);
nu = nu2;
nv = nv2;
}
GridderConfig<T> gconf(dirty.shape(0), dirty.shape(1), nu, nv, epsilon, pixsize_x, pixsize_y, baselines, nthreads);
auto idx = getWgtIndices(baselines, gconf, wgt, null_ms);
auto idx2 = mav<idx_t,1>(idx.data(),{idx.size()});
......
......@@ -67,6 +67,7 @@ def explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, xpixsize, ypixsize,
@pmp("nxdirty", (30, 128))
@pmp("nydirty", (128, 250))
@pmp("ofactor", (0, 1.2, 1.5, 1.7, 2.0))
@pmp("nrow", (2, 27))
@pmp("nchan", (1, 5))
@pmp("epsilon", (1e-1, 1e-3, 1e-5))
......@@ -74,7 +75,7 @@ def explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, xpixsize, ypixsize,
@pmp("wstacking", (True, False))
@pmp("use_wgt", (True, False))
@pmp("nthreads", (1, 2))
def test_adjointness_ms2dirty(nxdirty, nydirty, nrow, nchan, epsilon,
def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon,
singleprec, wstacking, use_wgt, nthreads):
if singleprec and epsilon < 5e-5:
pytest.skip()
......@@ -87,14 +88,21 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, nrow, nchan, epsilon,
ms = rng.random((nrow, nchan))-0.5 + 1j*(rng.random((nrow, nchan))-0.5)
wgt = rng.uniform(0.9, 1.1, (nrow, nchan)) if use_wgt else None
dirty = rng.random((nxdirty, nydirty))-0.5
nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
if nu & 1:
nu += 1
if nv & 1:
nv += 1
if ofactor == 0:
nu = nv = 0
if singleprec:
ms = ms.astype("c8")
dirty = dirty.astype("f4")
if wgt is not None:
wgt = wgt.astype("f4")
dirty2 = ng.ms2dirty(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex,
pixsizey, epsilon, wstacking, nthreads, 0).astype("f8")
ms2 = ng.dirty2ms(uvw, freq, dirty, wgt, pixsizex, pixsizey, epsilon,
pixsizey, nu, nv, epsilon, wstacking, nthreads, 0).astype("f8")
ms2 = ng.dirty2ms(uvw, freq, dirty, wgt, pixsizex, pixsizey, nu, nv, epsilon,
wstacking, nthreads, 0).astype("c16")
ref = max(my_vdot(ms,ms).real, my_vdot(ms2,ms2).real,
my_vdot(dirty, dirty).real, my_vdot(dirty2, dirty2).real)
......@@ -104,6 +112,7 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, nrow, nchan, epsilon,
@pmp('nxdirty', [16, 64])
@pmp('nydirty', [64])
@pmp('ofactor', [0, 1.2, 1.4, 1.7, 2])
@pmp("nrow", (1, 2, 27))
@pmp("nchan", (1, 5))
@pmp("epsilon", (1e-2, 1e-3, 1e-4, 1e-7))
......@@ -112,7 +121,7 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, nrow, nchan, epsilon,
@pmp("use_wgt", (True,))
@pmp("nthreads", (1, 2))
@pmp("fov", (1., 20.))
def test_ms2dirty_against_wdft2(nxdirty, nydirty, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, fov, nthreads):
def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, fov, nthreads):
if singleprec and epsilon < 5e-5:
pytest.skip()
rng = np.random.default_rng(42)
......@@ -124,11 +133,18 @@ def test_ms2dirty_against_wdft2(nxdirty, nydirty, nrow, nchan, epsilon, singlepr
ms = rng.random((nrow, nchan))-0.5 + 1j*(rng.random((nrow, nchan))-0.5)
wgt = rng.uniform(0.9, 1.1, (nrow, 1)) if use_wgt else None
wgt = np.broadcast_to(wgt, (nrow, nchan)) if use_wgt else None
nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
if nu & 1:
nu += 1
if nv & 1:
nv += 1
if ofactor == 0:
nu = nv = 0
if singleprec:
ms = ms.astype("c8")
if wgt is not None:
wgt = wgt.astype("f4")
dirty = ng.ms2dirty(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex,
pixsizey, epsilon, wstacking, nthreads, 0).astype("f8")
pixsizey, nu, nv, epsilon, wstacking, nthreads, 0).astype("f8")
ref = explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, pixsizex, pixsizey, wstacking)
assert_allclose(_l2error(dirty, ref), 0, atol=epsilon)
......@@ -36,8 +36,8 @@ auto None = py::none();
template<typename T> py::array ms2dirty2(const py::array &uvw_,
const py::array &freq_, const py::array &ms_, const py::object &wgt_,
size_t npix_x, size_t npix_y, double pixsize_x, double pixsize_y,
double epsilon, bool do_wstacking, size_t nthreads,
size_t npix_x, size_t npix_y, double pixsize_x, double pixsize_y, size_t nu,
size_t nv, double epsilon, bool do_wstacking, size_t nthreads,
size_t verbosity)
{
auto uvw = to_mav<double,2>(uvw_, false);
......@@ -49,23 +49,23 @@ template<typename T> py::array ms2dirty2(const py::array &uvw_,
auto dirty2 = to_mav<T,2>(dirty, true);
{
py::gil_scoped_release release;
ms2dirty(uvw,freq,ms,wgt2,pixsize_x,pixsize_y,epsilon,
ms2dirty(uvw,freq,ms,wgt2,pixsize_x,pixsize_y,nu,nv,epsilon,
do_wstacking,nthreads,dirty2,verbosity);
}
return move(dirty);
}
py::array Pyms2dirty(const py::array &uvw,
const py::array &freq, const py::array &ms, const py::object &wgt,
size_t npix_x, size_t npix_y, double pixsize_x, double pixsize_y,
double epsilon, bool do_wstacking, size_t nthreads,
size_t npix_x, size_t npix_y, double pixsize_x, double pixsize_y, size_t nu,
size_t nv, double epsilon, bool do_wstacking, size_t nthreads,
size_t verbosity)
{
if (isPyarr<complex<float>>(ms))
return ms2dirty2<float>(uvw, freq, ms, wgt, npix_x, npix_y,
pixsize_x, pixsize_y, epsilon, do_wstacking, nthreads, verbosity);
pixsize_x, pixsize_y, nu, nv, epsilon, do_wstacking, nthreads, verbosity);
if (isPyarr<complex<double>>(ms))
return ms2dirty2<double>(uvw, freq, ms, wgt, npix_x, npix_y,
pixsize_x, pixsize_y, epsilon, do_wstacking, nthreads, verbosity);
pixsize_x, pixsize_y, nu, nv, epsilon, do_wstacking, nthreads, verbosity);
MR_fail("type matching failed: 'ms' has neither type 'c8' nor 'c16'");
}
constexpr auto ms2dirty_DS = R"""(
......@@ -87,6 +87,15 @@ npix_x, npix_y: int
dimensions of the dirty image
pixsize_x, pixsize_y: float
angular pixel size (in radians) of the dirty image
nu, nv: int
dimensions of the (oversampled) intermediate uv grid
These values must be >= 1.2*the dimensions of the dirty image; tupical
oversampling values lie between 1.5 and 2.
Increasing the oversampling factor decreases the kernel support width
required for the desired accuracy, so it typically reduces run-time; on the
other hand, this will increase memory consumption.
If at least one of these two values is 0, the library will automatically
pick values that result in a fast computation.
epsilon: float
accuracy at which the computation should be done. Must be larger than 2e-13.
If `ms` has type np.complex64, it must be larger than 1e-5.
......@@ -108,7 +117,7 @@ np.array((nxdirty, nydirty), dtype=float of same precision as `ms`)
template<typename T> py::array dirty2ms2(const py::array &uvw_,
const py::array &freq_, const py::array &dirty_, const py::object &wgt_,
double pixsize_x, double pixsize_y, double epsilon,
double pixsize_x, double pixsize_y, size_t nu, size_t nv, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
{
auto uvw = to_mav<double,2>(uvw_, false);
......@@ -120,22 +129,22 @@ template<typename T> py::array dirty2ms2(const py::array &uvw_,
auto ms2 = to_mav<complex<T>,2>(ms, true);
{
py::gil_scoped_release release;
dirty2ms(uvw,freq,dirty,wgt2,pixsize_x,pixsize_y,epsilon,
dirty2ms(uvw,freq,dirty,wgt2,pixsize_x,pixsize_y,nu,nv,epsilon,
do_wstacking,nthreads,ms2,verbosity);
}
return move(ms);
}
py::array Pydirty2ms(const py::array &uvw,
const py::array &freq, const py::array &dirty, const py::object &wgt,
double pixsize_x, double pixsize_y, double epsilon,
double pixsize_x, double pixsize_y, size_t nu, size_t nv, double epsilon,
bool do_wstacking, size_t nthreads, size_t verbosity)
{
if (isPyarr<float>(dirty))
return dirty2ms2<float>(uvw, freq, dirty, wgt,
pixsize_x, pixsize_y, epsilon, do_wstacking, nthreads, verbosity);
pixsize_x, pixsize_y, nu, nv, epsilon, do_wstacking, nthreads, verbosity);
if (isPyarr<double>(dirty))
return dirty2ms2<double>(uvw, freq, dirty, wgt,
pixsize_x, pixsize_y, epsilon, do_wstacking, nthreads, verbosity);
pixsize_x, pixsize_y, nu, nv, epsilon, do_wstacking, nthreads, verbosity);
MR_fail("type matching failed: 'dirty' has neither type 'f4' nor 'f8'");
}
constexpr auto dirty2ms_DS = R"""(
......@@ -155,6 +164,15 @@ wgt: np.array((nrows, nchan), same dtype as `dirty`), optional
If present, its values are multiplied to the output
pixsize_x, pixsize_y: float
angular pixel size (in radians) of the dirty image
nu, nv: int
dimensions of the (oversampled) intermediate uv grid
These values must be >= 1.2*the dimensions of the dirty image; tupical
oversampling values lie between 1.5 and 2.
Increasing the oversampling factor decreases the kernel support width
required for the desired accuracy, so it typically reduces run-time; on the
other hand, this will increase memory consumption.
If at least one of these two values is 0, the library will automatically
pick values that result in a fast computation.
epsilon: float
accuracy at which the computation should be done. Must be larger than 2e-13.
If `dirty` has type np.float32, it must be larger than 1e-5.
......@@ -180,10 +198,10 @@ void add_wgridder(py::module &msup)
auto m = msup.def_submodule("wgridder");
m.def("ms2dirty", &Pyms2dirty, ms2dirty_DS, "uvw"_a, "freq"_a, "ms"_a,
"wgt"_a=None, "npix_x"_a, "npix_y"_a, "pixsize_x"_a, "pixsize_y"_a,
"wgt"_a=None, "npix_x"_a, "npix_y"_a, "pixsize_x"_a, "pixsize_y"_a, "nu"_a, "nv"_a,
"epsilon"_a, "do_wstacking"_a=false, "nthreads"_a=1, "verbosity"_a=0);
m.def("dirty2ms", &Pydirty2ms, dirty2ms_DS, "uvw"_a, "freq"_a, "dirty"_a,
"wgt"_a=None, "pixsize_x"_a, "pixsize_y"_a, "epsilon"_a,
"wgt"_a=None, "pixsize_x"_a, "pixsize_y"_a, "nu"_a, "nv"_a, "epsilon"_a,
"do_wstacking"_a=false, "nthreads"_a=1, "verbosity"_a=0);
}
......
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