Commit 4ba6126c authored by Martin Reinecke's avatar Martin Reinecke

switch to ducc0.wgridder

parent d9a04481
Pipeline #77116 passed with stages
in 23 minutes and 42 seconds
......@@ -13,7 +13,6 @@ RUN apt-get update && apt-get install -y \
python3-mpi4py python3-matplotlib \
# more optional NIFTy dependencies
&& pip3 install ducc0 \
&& pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git \
&& pip3 install jupyter \
&& rm -rf /var/lib/apt/lists/*
......
......@@ -50,9 +50,7 @@ Installation
Optional dependencies:
- [DUCC0](https://gitlab.mpcdf.mpg.de/mtr/ducc) for faster FFTs, spherical
harmonic transforms, and non-uniform Fourier transforms
- [nifty_gridder](https://gitlab.mpcdf.mpg.de/ift/nifty_gridder) (for radio
interferometry responses)
harmonic transforms, and radio interferometry gridding support
- [mpi4py](https://mpi4py.scipy.org) (for MPI-parallel execution)
- [matplotlib](https://matplotlib.org/) (for field plotting)
......@@ -86,10 +84,6 @@ If this library is present, NIFTy will detect it automatically and prefer
DUCC's FFT is compiled with optimizations for the host CPU and can provide
significantly faster transforms.
Support for the radio interferometry gridder is added via:
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
MPI support is added via:
sudo apt-get install python3-mpi4py
......
......@@ -45,9 +45,8 @@ def main():
img = ift.makeField(uvspace, img)
t0 = time()
GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv)
op = ift.Gridder(uvspace, eps=1e-7, uv=uv).adjoint
vis = ift.makeField(visspace, vis)
op = GM.getFull().adjoint
t1 = time()
op(img).val
t2 = time()
......@@ -59,15 +58,6 @@ def main():
b0s.append(t2 - t1)
c0s.append(t3 - t2)
print('Measure rest operator')
sc = ift.StatCalculator()
op = GM.getRest().adjoint
for _ in range(10):
t0 = time()
res = op(img)
sc.add(time() - t0)
print('FFT shape', res.shape)
plt.scatter(N0s, a0s, label='Gridder mr')
plt.legend()
# no idea why this is necessary, but if it is omitted, the range is wrong
......@@ -80,9 +70,6 @@ def main():
plt.scatter(N0s, b0s, color='k', marker='^', label='Gridder mr times')
plt.scatter(N0s, c0s, color='k', label='Gridder mr adjoint times')
plt.axhline(sc.mean, label='FFT')
plt.axhline(sc.mean + np.sqrt(sc.var))
plt.axhline(sc.mean - np.sqrt(sc.var))
plt.legend()
plt.ylabel('time [s]')
plt.title('Apply')
......
......@@ -23,10 +23,6 @@ If this library is present, NIFTy will detect it automatically and prefer
DUCC's FFT is compiled with optimizations for the host CPU and can provide
significantly faster transforms.
Support for the radio interferometry gridder is added via::
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git
MPI support is added via::
sudo apt-get install python3-mpi4py
......
......@@ -86,7 +86,7 @@ from .library.light_cone_operator import LightConeOperator
from .library.wiener_filter_curvature import WienerFilterCurvature
from .library.adjust_variances import (make_adjust_variances_hamiltonian,
do_adjust_variances)
from .library.gridder import GridderMaker
from .library.gridder import Gridder
from .library.correlated_fields import CorrelatedFieldMaker
from . import extra
......
......@@ -24,88 +24,42 @@ from ..operators.linear_operator import LinearOperator
from ..sugar import makeDomain, makeField
class GridderMaker(object):
def __init__(self, dirty_domain, uv, eps=2e-13):
import nifty_gridder
dirty_domain = makeDomain(dirty_domain)
if (len(dirty_domain) != 1 or not isinstance(dirty_domain[0], RGSpace)
or not len(dirty_domain.shape) == 2):
raise ValueError("need dirty_domain with exactly one 2D RGSpace")
class Gridder(LinearOperator):
def __init__(self, target, uv, eps=2e-10, nthreads=1):
self._capability = self.TIMES | self.ADJOINT_TIMES
self._target = makeDomain(target)
if (len(self._target) != 1 or not isinstance(self._target[0], RGSpace)
or not len(self._target.shape) == 2):
raise ValueError("need target with exactly one 2D RGSpace")
if uv.ndim != 2:
raise ValueError("uv must be a 2D array")
if uv.shape[1] != 2:
raise ValueError("second dimension of uv must have length 2")
dstx, dsty = dirty_domain[0].distances
# wasteful hack to adjust to shape required by nifty_gridder
uvw = np.empty((uv.shape[0], 3), dtype=np.float64)
uvw[:, 0:2] = uv
uvw[:, 2] = 0.
# Scale uv such that 0<uv<=1 which is assumed by nifty_gridder
uvw[:, 0] = uvw[:, 0]*dstx
uvw[:, 1] = uvw[:, 1]*dsty
speedOfLight = 299792458.
bl = nifty_gridder.Baselines(uvw, np.array([speedOfLight]))
nxdirty, nydirty = dirty_domain.shape
gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.)
nu, nv = gconf.Nu(), gconf.Nv()
self._idx = nifty_gridder.getIndices(
bl, gconf, np.zeros((uv.shape[0], 1), dtype=np.bool))
self._bl = bl
du, dv = 1./(nu*dstx), 1./(nv*dsty)
grid_domain = RGSpace([nu, nv], distances=[du, dv], harmonic=True)
self._rest = _RestOperator(dirty_domain, grid_domain, gconf)
self._gridder = RadioGridder(grid_domain, bl, gconf, self._idx)
def getGridder(self):
return self._gridder
def getRest(self):
return self._rest
def getFull(self):
return self.getRest() @ self._gridder
def ms2vis(self, x):
return self._bl.ms2vis(x, self._idx)
class _RestOperator(LinearOperator):
def __init__(self, dirty_domain, grid_domain, gconf):
self._domain = makeDomain(grid_domain)
self._target = makeDomain(dirty_domain)
self._gconf = gconf
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self, x, mode):
self._check_input(x, mode)
res = x.val
if mode == self.TIMES:
res = self._gconf.grid2dirty(res)
else:
res = self._gconf.dirty2grid(res)
return makeField(self._tgt(mode), res)
class RadioGridder(LinearOperator):
def __init__(self, grid_domain, bl, gconf, idx):
self._domain = DomainTuple.make(
UnstructuredDomain((bl.Nrows())))
self._target = DomainTuple.make(grid_domain)
self._bl = bl
self._gconf = gconf
self._idx = idx
self._capability = self.TIMES | self.ADJOINT_TIMES
UnstructuredDomain((uv.shape[0])))
# wasteful hack to adjust to shape required by ducc0.wgridder
self._uvw = np.empty((uv.shape[0], 3), dtype=np.float64)
self._uvw[:, 0:2] = uv
self._uvw[:, 2] = 0.
self._eps = float(eps)
self._nthreads = int(nthreads)
def apply(self, x, mode):
import nifty_gridder
self._check_input(x, mode)
speedOfLight = 299792458.
freq = np.array([speedOfLight])
x = x.val
nxdirty, nydirty = self._target[0].shape
nu, nv = max(2*nxdirty, 16), max(2*nydirty, 16)
dstx, dsty = self._target[0].distances
from ducc0.wgridder import ms2dirty, dirty2ms
if mode == self.TIMES:
x = self._bl.ms2vis(x.val.reshape((-1, 1)), self._idx)
res = nifty_gridder.vis2grid(self._bl, self._gconf, self._idx, x)
res = ms2dirty(self._uvw, freq, x.reshape((-1,1)), None, nxdirty,
nydirty, dstx, dsty, nu, nv,
self._eps, False, self._nthreads, 0)
else:
res = nifty_gridder.grid2vis(self._bl, self._gconf, self._idx,
x.val)
res = self._bl.vis2ms(res, self._idx).reshape((-1,))
res = dirty2ms(self._uvw, freq, x, None, dstx, dsty, nu, nv,
self._eps, False, self._nthreads, 0)
res = res.reshape((-1,))
return makeField(self._tgt(mode), res)
......@@ -46,10 +46,9 @@ def test_gridding(nu, nv, N, eps):
dstx, dsty = dom.distances
uv[:, 0] = uv[:, 0]/dstx
uv[:, 1] = uv[:, 1]/dsty
GM = ift.GridderMaker(dom, uv=uv, eps=eps)
Op = ift.Gridder(dom, uv=uv, eps=eps)
vis2 = ift.makeField(ift.UnstructuredDomain(vis.shape), vis)
Op = GM.getFull()
pynu = Op(vis2).val
# DFT
x, y = np.meshgrid(
......@@ -72,8 +71,7 @@ def test_cartesian():
tmp = np.vstack([uu[None, :], vv[None, :]])
uv = np.transpose(tmp, (2, 1, 0)).reshape(-1, 2)
GM = ift.GridderMaker(dom, uv=uv)
op = GM.getFull().adjoint
op = ift.Gridder(dom, uv=uv).adjoint
fld = ift.from_random(dom, 'normal')
arr = fld.val
......@@ -95,16 +93,9 @@ def test_cartesian():
def test_build(nu, nv, N, eps):
dom = ift.RGSpace([nu, nv])
uv = ift.random.current_rng().random((N, 2)) - 0.5
GM = ift.GridderMaker(dom, uv=uv, eps=eps)
R0 = GM.getGridder()
R1 = GM.getRest()
R = R1@R0
RF = GM.getFull()
RF = ift.Gridder(dom, uv=uv, eps=eps)
# Consistency checks
flt = np.float64
cmplx = np.complex128
ift.extra.check_linear_operator(R0, cmplx, flt, only_r_linear=True)
ift.extra.check_linear_operator(R1, flt, flt)
ift.extra.check_linear_operator(R, cmplx, flt, only_r_linear=True)
ift.extra.check_linear_operator(RF, cmplx, flt, only_r_linear=True)
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