Merge branch 'PowerGrid' into 'NIFTy_5'

Power grid

See merge request !327

Dockerfile

@@ -13,7 +13,7 @@ RUN apt-get update && apt-get install -y \
     python3-mpi4py python3-matplotlib \
   # more optional NIFTy dependencies
   && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \
-  && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git \
+  && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/nifty_gridder.git@better_params \
   && pip3 install git+https://gitlab.mpcdf.mpg.de/mtr/pypocketfft.git \
   && pip3 install jupyter \
   && rm -rf /var/lib/apt/lists/*

demos/bench_gridder.py

@@ -9,7 +9,7 @@ np.random.seed(40)
 
 N0s, a0s, b0s, c0s = [], [], [], []
 
-for ii in range(10, 23):
+for ii in range(10, 26):
     nu = 1024
     nv = 1024
     N = int(2**ii)
@@ -27,17 +27,15 @@ for ii in range(10, 23):
     img = ift.from_global_data(uvspace, img)
 
     t0 = time()
-    GM = ift.GridderMaker(uvspace, eps=1e-7)
-    idx = GM.getReordering(uv)
-    uv = uv[idx]
-    vis = vis[idx]
+    GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv)
     vis = ift.from_global_data(visspace, vis)
-    op = GM.getFull(uv).adjoint
+    op = GM.getFull().adjoint
     t1 = time()
     op(img).to_global_data()
     t2 = time()
     op.adjoint(vis).to_global_data()
     t3 = time()
+    print(t2-t1, t3-t2)
     N0s.append(N)
     a0s.append(t1 - t0)
     b0s.append(t2 - t1)

demos/getting_started_3.py

@@ -109,7 +109,8 @@ if __name__ == '__main__':
     minimizer = ift.NewtonCG(ic_newton)
 
     # Set up likelihood and information Hamiltonian
-    likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
+    likelihood = ift.GaussianEnergy(mean=data,
+                                    inverse_covariance=N.inverse)(signal_response)
     H = ift.StandardHamiltonian(likelihood, ic_sampling)
 
     initial_mean = ift.MultiField.full(H.domain, 0.)

nifty5/domains/log_rg_space.py

@@ -80,8 +80,8 @@ class LogRGSpace(StructuredDomain):
         return np.array(self._t_0)
 
     def __repr__(self):
-        return ("LogRGSpace(shape={}, harmonic={})".format(
-            self.shape, self.harmonic))
+        return ("LogRGSpace(shape={}, bindistances={}, t_0={}, harmonic={})".format(
+            self.shape, self.bindistances, self.t_0, self.harmonic))
 
     def get_default_codomain(self):
         """Returns a :class:`LogRGSpace` object representing the (position or

nifty5/library/gridder.py

@@ -15,108 +15,96 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import numpy as np
-
 from ..domain_tuple import DomainTuple
 from ..domains.rg_space import RGSpace
 from ..domains.unstructured_domain import UnstructuredDomain
-from ..fft import hartley
 from ..operators.linear_operator import LinearOperator
 from ..sugar import from_global_data, makeDomain
+import numpy as np
 
 
 class GridderMaker(object):
-    def __init__(self, domain, eps=2e-13):
-        from nifty_gridder import get_w
-        domain = makeDomain(domain)
-        if (len(domain) != 1 or not isinstance(domain[0], RGSpace) or
-                not len(domain.shape) == 2):
-            raise ValueError("need domain with exactly one 2D RGSpace")
-        nu, nv = domain.shape
-        if nu % 2 != 0 or nv % 2 != 0:
-            raise ValueError("dimensions must be even")
-        nu2, nv2 = 2*nu, 2*nv
-        w = get_w(eps)
-        nsafe = (w+1)//2
-        nu2 = max([nu2, 2*nsafe])
-        nv2 = max([nv2, 2*nsafe])
-
-        oversampled_domain = RGSpace(
-            [nu2, nv2], distances=[1, 1], harmonic=False)
-
-        self._eps = eps
-        self._rest = _RestOperator(domain, oversampled_domain, eps)
-
-    def getReordering(self, uv):
-        from nifty_gridder import peanoindex
-        nu2, nv2 = self._rest._domain.shape
-        return peanoindex(uv, nu2, nv2)
-
-    def getGridder(self, uv):
-        return RadioGridder(self._rest.domain, self._eps, uv)
+    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")
+        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, uv):
-        return self.getRest() @ self.getGridder(uv)
+    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, domain, oversampled_domain, eps):
-        from nifty_gridder import correction_factors
-        self._domain = makeDomain(oversampled_domain)
-        self._target = domain
-        nu, nv = domain.shape
-        nu2, nv2 = oversampled_domain.shape
 
-        fu = correction_factors(nu2, nu//2+1, eps)
-        fv = correction_factors(nv2, nv//2+1, eps)
-        # compute deconvolution operator
-        rng = np.arange(nu)
-        k = np.minimum(rng, nu-rng)
-        self._deconv_u = np.roll(fu[k], -nu//2).reshape((-1, 1))
-        rng = np.arange(nv)
-        k = np.minimum(rng, nv-rng)
-        self._deconv_v = np.roll(fv[k], -nv//2).reshape((1, -1))
+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)
-        nu, nv = self._target.shape
         res = x.to_global_data()
         if mode == self.TIMES:
-            res = hartley(res)
-            res = np.roll(res, (nu//2, nv//2), axis=(0, 1))
-            res = res[:nu, :nv]
-            res *= self._deconv_u
-            res *= self._deconv_v
+            res = self._gconf.grid2dirty(res)
         else:
-            res = res*self._deconv_u
-            res *= self._deconv_v
-            nu2, nv2 = self._domain.shape
-            res = np.pad(res, ((0, nu2-nu), (0, nv2-nv)), mode='constant',
-                         constant_values=0)
-            res = np.roll(res, (-nu//2, -nv//2), axis=(0, 1))
-            res = hartley(res)
+            res = self._gconf.dirty2grid(res)
         return from_global_data(self._tgt(mode), res)
 
 
 class RadioGridder(LinearOperator):
-    def __init__(self, target, eps, uv):
+    def __init__(self, grid_domain, bl, gconf, idx):
         self._domain = DomainTuple.make(
-            UnstructuredDomain((uv.shape[0],)))
-        self._target = DomainTuple.make(target)
+            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
-        self._eps = float(eps)
-        self._uv = uv  # FIXME: should we write-protect this?
 
     def apply(self, x, mode):
-        from nifty_gridder import to_grid, from_grid
+        import nifty_gridder
         self._check_input(x, mode)
         if mode == self.TIMES:
-            nu2, nv2 = self._target.shape
-            res = to_grid(self._uv, x.to_global_data(), nu2, nv2, self._eps)
+            x = self._bl.ms2vis(x.to_global_data().reshape((-1, 1)), self._idx)
+            res = nifty_gridder.vis2grid(self._bl, self._gconf, self._idx, x)
         else:
-            res = from_grid(self._uv, x.to_global_data(), self._eps)
+            res = nifty_gridder.grid2vis(self._bl, self._gconf, self._idx,
+                                         x.to_global_data())
+            res = self._bl.vis2ms(res, self._idx).reshape((-1,))
         return from_global_data(self._tgt(mode), res)

nifty5/operators/adder.py

@@ -28,12 +28,15 @@ class Adder(Operator):
     field : Field or MultiField
         The field by which the input is shifted.
     """
-    def __init__(self, field):
+    def __init__(self, field, neg=False):
         if not isinstance(field, (Field, MultiField)):
             raise TypeError
         self._field = field
         self._domain = self._target = field.domain
+        self._neg = bool(neg)
 
     def apply(self, x):
         self._check_input(x)
+        if self._neg:
+            return x - self._field
         return x + self._field

nifty5/operators/energy_operators.py

@@ -110,8 +110,8 @@ class GaussianEnergy(EnergyOperator):
     ----------
     mean : Field
         Mean of the Gaussian. Default is 0.
-    covariance : LinearOperator
-        Covariance of the Gaussian. Default is the identity operator.
+    inverse_covariance : LinearOperator
+        Inverse covariance of the Gaussian. Default is the identity operator.
     domain : Domain, DomainTuple, tuple of Domain or MultiDomain
         Operator domain. By default it is inferred from `mean` or
         `covariance` if specified
@@ -121,28 +121,27 @@ class GaussianEnergy(EnergyOperator):
     At least one of the arguments has to be provided.
     """
 
-    def __init__(self, mean=None, covariance=None, domain=None):
+    def __init__(self, mean=None, inverse_covariance=None, domain=None):
         if mean is not None and not isinstance(mean, (Field, MultiField)):
             raise TypeError
-        if covariance is not None and not isinstance(covariance,
-                                                     LinearOperator):
+        if inverse_covariance is not None and not isinstance(inverse_covariance, LinearOperator):
             raise TypeError
 
         self._domain = None
         if mean is not None:
             self._checkEquivalence(mean.domain)
-        if covariance is not None:
-            self._checkEquivalence(covariance.domain)
+        if inverse_covariance is not None:
+            self._checkEquivalence(inverse_covariance.domain)
         if domain is not None:
             self._checkEquivalence(domain)
         if self._domain is None:
             raise ValueError("no domain given")
         self._mean = mean
-        if covariance is None:
+        if inverse_covariance is None:
             self._op = SquaredNormOperator(self._domain).scale(0.5)
         else:
-            self._op = QuadraticFormOperator(covariance.inverse)
-        self._icov = None if covariance is None else covariance.inverse
+            self._op = QuadraticFormOperator(inverse_covariance)
+        self._icov = None if inverse_covariance is None else inverse_covariance
 
     def _checkEquivalence(self, newdom):
         newdom = makeDomain(newdom)

test/test_operators/test_nft.py

@@ -39,23 +39,53 @@ def test_gridding(nu, nv, N, eps):
     vis = np.random.randn(N) + 1j*np.random.randn(N)
 
     # Nifty
-    GM = ift.GridderMaker(ift.RGSpace((nu, nv)), eps=eps)
-    # re-order for performance
-    idx = GM.getReordering(uv)
-    uv, vis = uv[idx], vis[idx]
+    dom = ift.RGSpace((nu, nv), distances=(0.2, 1.12))
+    dstx, dsty = dom.distances
+    uv[:, 0] = uv[:, 0]/dstx
+    uv[:, 1] = uv[:, 1]/dsty
+    GM = ift.GridderMaker(dom, uv=uv, eps=eps)
     vis2 = ift.from_global_data(ift.UnstructuredDomain(vis.shape), vis)
 
-    Op = GM.getFull(uv)
+    Op = GM.getFull()
     pynu = Op(vis2).to_global_data()
     # DFT
     x, y = np.meshgrid(
         *[-ss/2 + np.arange(ss) for ss in [nu, nv]], indexing='ij')
     dft = pynu*0.
     for i in range(N):
-        dft += (vis[i]*np.exp(2j*np.pi*(x*uv[i, 0] + y*uv[i, 1]))).real
+        dft += (
+            vis[i]*np.exp(2j*np.pi*(x*uv[i, 0]*dstx + y*uv[i, 1]*dsty))).real
     assert_(_l2error(dft, pynu) < eps)
 
 
+def test_cartesian():
+    nx, ny = 2, 6
+    dstx, dsty = 0.3, 0.2
+    dom = ift.RGSpace((nx, ny), (dstx, dsty))
+
+    kx = np.fft.fftfreq(nx, dstx)
+    ky = np.fft.fftfreq(ny, dsty)
+    uu, vv = np.meshgrid(kx, ky)
+    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
+
+    fld = ift.from_random('normal', dom)
+    arr = fld.to_global_data()
+
+    fld2 = ift.from_global_data(dom, np.roll(arr, (nx//2, ny//2), axis=(0, 1)))
+    res = op(fld2).to_global_data().reshape(nx, ny)
+
+    fft = ift.FFTOperator(dom.get_default_codomain(), target=dom).adjoint
+    vol = ift.full(dom, 1.).integrate()
+    res1 = fft(fld).to_global_data()
+
+    # FIXME: we don't understand the conjugate() yet
+    np.testing.assert_allclose(res, res1.conjugate()*vol)
+
+
 @pmp('eps', [1e-2, 1e-6, 2e-13])
 @pmp('nu', [12, 128])
 @pmp('nv', [4, 12, 128])
@@ -63,14 +93,11 @@ def test_gridding(nu, nv, N, eps):
 def test_build(nu, nv, N, eps):
     dom = ift.RGSpace([nu, nv])
     uv = np.random.rand(N, 2) - 0.5
-    GM = ift.GridderMaker(dom, eps=eps)
-    # re-order for performance
-    idx = GM.getReordering(uv)
-    uv = uv[idx]
-    R0 = GM.getGridder(uv)
+    GM = ift.GridderMaker(dom, uv=uv, eps=eps)
+    R0 = GM.getGridder()
     R1 = GM.getRest()
     R = R1@R0
-    RF = GM.getFull(uv)
+    RF = GM.getFull()
 
     # Consistency checks
     flt = np.float64