Merge branch 'op_twiddling' into 'NIFTy_4'

Replace InverseOperator and AdjointOperator with OperatorAdapter, and more

See merge request ift/NIFTy!237

.gitlab-ci.yml

@@ -3,32 +3,41 @@ image: debian:testing-slim
 
 stages:
   - test
+  - release
 
 variables:
   DOCKER_DRIVER: overlay
 
-before_script:
-  - apt-get update
-  - sh ci/install_basics.sh
-  - pip install --process-dependency-links -r ci/requirements.txt
-  - pip3 install --process-dependency-links -r ci/requirements.txt
-  - pip install --user .
-  - pip3 install --user .
-
-test_min:
+test_python2:
   stage: test
   script:
-    - nosetests3 -q
-    - OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 4 nosetests -q 2>/dev/null
-    - OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 4 nosetests3 -q 2>/dev/null
+    - apt-get update > /dev/null
+    - apt-get install -y git libfftw3-dev openmpi-bin libopenmpi-dev python python-pip python-dev python-nose python-numpy python-matplotlib python-future python-mpi4py python-scipy > /dev/null
+    - pip install --process-dependency-links parameterized coverage git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git > /dev/null
+    - pip install --user .
+    - OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 2 nosetests -q 2>/dev/null
     - nosetests -q --with-coverage --cover-package=nifty4 --cover-branches --cover-erase
     - >
       coverage report | grep TOTAL | awk '{ print "TOTAL: "$6; }'
 
+test_python3:
+  stage: test
+  script:
+    - apt-get update > /dev/null
+    - apt-get install -y git libfftw3-dev openmpi-bin libopenmpi-dev python3 python3-pip python3-dev python3-nose python3-numpy python3-matplotlib python3-future python3-mpi4py python3-scipy > /dev/null
+    - pip3 install --process-dependency-links parameterized git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git > /dev/null
+    - pip3 install --user .
+    - nosetests3 -q
+    - OMP_NUM_THREADS=1 mpiexec --allow-run-as-root -n 2 nosetests3 -q 2>/dev/null
+
 pages:
+  stage: release
   script:
-  - sh docs/generate.sh
-  - mv docs/build/ public/
+    - apt-get update > /dev/null
+    - apt-get install -y git libfftw3-dev python python-pip python-dev python-numpy python-future python-sphinx python-sphinx-rtd-theme python-numpydoc > /dev/null
+    - pip install --user .
+    - sh docs/generate.sh > /dev/null
+    - mv docs/build/ public/
   artifacts:
     paths:
     - public

ci/install_basics.sh

-apt-get install -y git libfftw3-dev openmpi-bin libopenmpi-dev \
-  python  python-pip  python-dev  python-nose  python-numpy  python-matplotlib  python-future  python-mpi4py python-scipy \
-  python3 python3-pip python3-dev python3-nose python3-numpy python3-matplotlib python3-future python3-mpi4py python3-scipy

ci/requirements.txt

-parameterized
-coverage
-git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
-sphinx
-sphinx_rtd_theme
-numpydoc

demos/Wiener_Filter.ipynb

@@ -651,7 +651,7 @@
     "ma = np.max(s_data)\n",
     "\n",
     "fig, axes = plt.subplots(3, 2, figsize=(15, 22.5))\n",
-    "sample = HT(curv.draw_sample()+m).to_global_data()\n",
+    "sample = HT(curv.draw_sample(from_inverse=True)+m).to_global_data()\n",
     "post_mean = (m_mean + HT(m)).to_global_data()\n",
     "\n",
     "data = [s_data, m_data, post_mean, sample, s_data - m_data, uncertainty]\n",

demos/wiener_filter_easy.py

@@ -52,8 +52,7 @@ if __name__ == "__main__":
                                     tol_abs_gradnorm=0.1)
     inverter = ift.ConjugateGradient(controller=IC)
     D = (ift.SandwichOperator(R, N.inverse) + Sh.inverse).inverse
-    # MR FIXME: we can/should provide a preconditioner here as well!
-    D = ift.InversionEnabler(D, inverter)
+    D = ift.InversionEnabler(D, inverter, approximation=Sh)
     m = D(j)
 
     # Plotting

nifty4/data_objects/distributed_do.py

@@ -236,12 +236,24 @@ class data_object(object):
     def __ipow__(self, other):
         return self._binary_helper(other, op='__ipow__')
 
-    def __eq__(self, other):
-        return self._binary_helper(other, op='__eq__')
+    def __lt__(self, other):
+        return self._binary_helper(other, op='__lt__')
+
+    def __le__(self, other):
+        return self._binary_helper(other, op='__le__')
 
     def __ne__(self, other):
         return self._binary_helper(other, op='__ne__')
 
+    def __eq__(self, other):
+        return self._binary_helper(other, op='__eq__')
+
+    def __ge__(self, other):
+        return self._binary_helper(other, op='__ge__')
+
+    def __gt__(self, other):
+        return self._binary_helper(other, op='__gt__')
+
     def __neg__(self):
         return data_object(self._shape, -self._data, self._distaxis)
 

nifty4/extra/__init__.py

 from .operator_tests import consistency_check
-from .energy_tests import check_value_gradient_consistency
+from .energy_tests import *

nifty4/extra/energy_tests.py

@@ -19,35 +19,63 @@
 import numpy as np
 from ..field import Field
 
-__all__ = ["check_value_gradient_consistency"]
+__all__ = ["check_value_gradient_consistency",
+           "check_value_gradient_curvature_consistency"]
 
 
-def check_value_gradient_consistency(E, tol=1e-6, ntries=100):
+def _get_acceptable_energy(E):
     if not np.isfinite(E.value):
         raise ValueError
+    dir = Field.from_random("normal", E.position.domain)
+    # find a step length that leads to a "reasonable" energy
+    for i in range(50):
+        try:
+            E2 = E.at(E.position+dir)
+            if np.isfinite(E2.value) and abs(E2.value) < 1e20:
+                break
+        except FloatingPointError:
+            pass
+        dir *= 0.5
+    else:
+        raise ValueError("could not find a reasonable initial step")
+    return E2
+
+
+def check_value_gradient_consistency(E, tol=1e-6, ntries=100):
     for _ in range(ntries):
-        dir = Field.from_random("normal", E.position.domain)
-        # find a step length that leads to a "reasonable" energy
+        E2 = _get_acceptable_energy(E)
+        dir = E2.position - E.position
+        Enext = E2
+        dirnorm = dir.norm()
+        dirder = E.gradient.vdot(dir)/dirnorm
         for i in range(50):
-            try:
-                E2 = E.at(E.position+dir)
-                if np.isfinite(E2.value) and abs(E2.value) < 1e20:
-                    break
-            except FloatingPointError:
-                pass
+            print(abs((E2.value-E.value)/dirnorm-dirder))
+            if abs((E2.value-E.value)/dirnorm-dirder) < tol:
+                break
             dir *= 0.5
+            dirnorm *= 0.5
+            E2 = E2.at(E.position+dir)
         else:
-            raise ValueError("could not find a reasonable initial step")
+            raise ValueError("gradient and value seem inconsistent")
+        # E = Enext
+
+
+def check_value_gradient_curvature_consistency(E, tol=1e-6, ntries=100):
+    for _ in range(ntries):
+        E2 = _get_acceptable_energy(E)
+        dir = E2.position - E.position
         Enext = E2
-        dirder = E.gradient.vdot(dir)
+        dirnorm = dir.norm()
+        dirder = E.gradient.vdot(dir)/dirnorm
+        dgrad = E.curvature(dir)/dirnorm
         for i in range(50):
-            Ediff = E2.value - E.value
-            eps = 1e-10*max(abs(E.value), abs(E2.value))
-            if abs(Ediff-dirder) < max([tol*abs(Ediff), tol*abs(dirder), eps]):
+            gdiff = E2.gradient - E.gradient
+            if abs((E2.value-E.value)/dirnorm-dirder) < tol and \
+               (abs((E2.gradient-E.gradient)/dirnorm-dgrad) < tol).all():
                 break
             dir *= 0.5
-            dirder *= 0.5
+            dirnorm *= 0.5
             E2 = E2.at(E.position+dir)
         else:
-            raise ValueError("gradient and value seem inconsistent")
-        E = Enext
+            raise ValueError("gradient, value and curvature seem inconsistent")
+        # E = Enext

nifty4/library/nonlinear_power_energy.py

@@ -70,7 +70,7 @@ class NonlinearPowerEnergy(Energy):
             if samples is None or samples == 0:
                 xi_sample_list = [xi]
             else:
-                xi_sample_list = [D.inverse_draw_sample() + xi
+                xi_sample_list = [D.draw_sample(from_inverse=True) + xi
                                   for _ in range(samples)]
         self.xi_sample_list = xi_sample_list
         self.inverter = inverter

nifty4/library/wiener_filter_curvature.py

@@ -40,4 +40,4 @@ def WienerFilterCurvature(R, N, S, inverter):
         The minimizer to use during numerical inversion
     """
     op = SandwichOperator(R, N.inverse) + S.inverse
-    return InversionEnabler(op, inverter, S)
+    return InversionEnabler(op, inverter, S.inverse)

nifty4/minimization/__init__.py

@@ -10,7 +10,7 @@ from .steepest_descent import SteepestDescent
 from .vl_bfgs import VL_BFGS
 from .l_bfgs import L_BFGS
 from .relaxed_newton import RelaxedNewton
-from .scipy_minimizer import ScipyMinimizer, NewtonCG, L_BFGS_B
+from .scipy_minimizer import ScipyMinimizer, NewtonCG, L_BFGS_B, ScipyCG
 from .energy import Energy
 from .quadratic_energy import QuadraticEnergy
 from .line_energy import LineEnergy
@@ -19,5 +19,5 @@ __all__ = ["LineSearch", "LineSearchStrongWolfe",
            "IterationController", "GradientNormController",
            "Minimizer", "ConjugateGradient", "NonlinearCG", "DescentMinimizer",
            "SteepestDescent", "VL_BFGS", "RelaxedNewton", "ScipyMinimizer",
-           "NewtonCG", "L_BFGS_B", "Energy", "QuadraticEnergy", "LineEnergy",
-           "L_BFGS"]
+           "NewtonCG", "L_BFGS_B", "ScipyCG", "Energy", "QuadraticEnergy",
+           "LineEnergy", "L_BFGS"]

nifty4/minimization/scipy_minimizer.py

@@ -24,14 +24,26 @@ from ..logger import logger
 from .iteration_controller import IterationController
 
 
+def _toNdarray(fld):
+    return fld.to_global_data().reshape(-1)
+
+
+def _toFlatNdarray(fld):
+    return fld.val.flatten()
+
+
+def _toField(arr, dom):
+    return Field.from_global_data(dom, arr.reshape(dom.shape))
+
+
 class _MinHelper(object):
     def __init__(self, energy):
         self._energy = energy
         self._domain = energy.position.domain
 
     def _update(self, x):
-        pos = Field(self._domain, x.reshape(self._domain.shape))
-        if (pos.val != self._energy.position.val).any():
+        pos = _toField(x, self._domain)
+        if (pos != self._energy.position).any():
             self._energy = self._energy.at(pos.locked_copy())
 
     def fun(self, x):
@@ -40,13 +52,12 @@ class _MinHelper(object):
 
     def jac(self, x):
         self._update(x)
-        return self._energy.gradient.val.flatten()
+        return _toFlatNdarray(self._energy.gradient)
 
     def hessp(self, x, p):
         self._update(x)
-        vec = Field(self._domain, p.reshape(self._domain.shape))
-        res = self._energy.curvature(vec)
-        return res.val.flatten()
+        res = self._energy.curvature(_toField(p, self._domain))
+        return _toFlatNdarray(res)
 
 
 class ScipyMinimizer(Minimizer):
@@ -113,3 +124,40 @@ def L_BFGS_B(ftol, gtol, maxiter, maxcor=10, disp=False, bounds=None):
     options = {"ftol": ftol, "gtol": gtol, "maxiter": maxiter,
                "maxcor": maxcor, "disp": disp}
     return ScipyMinimizer("L-BFGS-B", options, False, bounds)
+
+
+class ScipyCG(Minimizer):
+    def __init__(self, tol, maxiter):
+        super(ScipyCG, self).__init__()
+        if not dobj.is_numpy():
+            raise NotImplementedError
+        self._tol = tol
+        self._maxiter = maxiter
+
+    def __call__(self, energy, preconditioner=None):
+        from scipy.sparse.linalg import LinearOperator as scipy_linop, cg
+        from .quadratic_energy import QuadraticEnergy
+        if not isinstance(energy, QuadraticEnergy):
+            raise ValueError("need a quadratic energy for CG")
+
+        class mymatvec(object):
+            def __init__(self, op):
+                self._op = op
+
+            def __call__(self, inp):
+                return _toNdarray(self._op(_toField(inp, self._op.domain)))
+
+        op = energy._A
+        b = _toNdarray(energy._b)
+        sx = _toNdarray(energy.position)
+        sci_op = scipy_linop(shape=(op.domain.size, op.target.size),
+                             matvec=mymatvec(op))
+        prec_op = None
+        if preconditioner is not None:
+            prec_op = scipy_linop(shape=(op.domain.size, op.target.size),
+                                  matvec=mymatvec(preconditioner))
+        res, stat = cg(sci_op, b, x0=sx, tol=self._tol, M=prec_op,
+                       maxiter=self._maxiter)
+        stat = (IterationController.CONVERGED if stat >= 0 else
+                IterationController.ERROR)
+        return energy.at(_toField(res, op.domain)), stat

nifty4/operators/adjoint_operator.py

-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-# Copyright(C) 2013-2018 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
-# and financially supported by the Studienstiftung des deutschen Volkes.
-
-from .linear_operator import LinearOperator
-
-
-class AdjointOperator(LinearOperator):
-    """Adapter class representing the adjoint of a given operator."""
-
-    def __init__(self, op):
-        super(AdjointOperator, self).__init__()
-        self._op = op
-
-    @property
-    def domain(self):
-        return self._op.target
-
-    @property
-    def target(self):
-        return self._op.domain
-
-    @property
-    def capability(self):
-        return self._adjointCapability[self._op.capability]
-
-    @property
-    def adjoint(self):
-        return self._op
-
-    def apply(self, x, mode):
-        return self._op.apply(x, self._adjointMode[mode])

nifty4/operators/chain_operator.py

@@ -96,13 +96,15 @@ class ChainOperator(LinearOperator):
     def target(self):
         return self._ops[0].target
 
-    @property
-    def inverse(self):
-        return self.make([op.inverse for op in reversed(self._ops)])
-
-    @property
-    def adjoint(self):
-        return self.make([op.adjoint for op in reversed(self._ops)])
+    def _flip_modes(self, mode):
+        if mode == 0:
+            return self
+        if mode == 1 or mode == 2:
+            return self.make([op._flip_modes(mode)
+                              for op in reversed(self._ops)])
+        if mode == 3:
+            return self.make([op._flip_modes(mode) for op in self._ops])
+        raise ValueError("bad operator flipping mode")
 
     @property
     def capability(self):

nifty4/operators/diagonal_operator.py

@@ -158,35 +158,30 @@ class DiagonalOperator(EndomorphicOperator):
     def capability(self):
         return self._all_ops
 
-    @property
-    def inverse(self):
-        res = self._skeleton(())
-        res._ldiag = 1./self._ldiag
-        return res
-
-    @property
-    def adjoint(self):
-        if np.issubdtype(self._ldiag.dtype, np.floating):
+    def _flip_modes(self, mode):
+        if mode == 0:
+            return self
+        if mode == 1 and np.issubdtype(self._ldiag.dtype, np.floating):
             return self
         res = self._skeleton(())
-        res._ldiag = self._ldiag.conjugate()
-        return res
-
-    def draw_sample(self, dtype=np.float64):
-        if (np.issubdtype(self._ldiag.dtype, np.complexfloating) or
-                (self._ldiag <= 0.).any()):
-            raise ValueError("operator not positive definite")
-        res = Field.from_random(random_type="normal", domain=self._domain,
-                                dtype=dtype)
-        res.local_data[()] *= np.sqrt(self._ldiag)
+        if mode == 1:
+            res._ldiag = self._ldiag.conjugate()
+        elif mode == 2:
+            res._ldiag = 1./self._ldiag
+        elif mode == 3:
+            res._ldiag = 1./self._ldiag.conjugate()
+        else:
+            raise ValueError("bad operator flipping mode")
         return res
 
-    def inverse_draw_sample(self, dtype=np.float64):
+    def draw_sample(self, from_inverse=False, dtype=np.float64):
         if (np.issubdtype(self._ldiag.dtype, np.complexfloating) or
                 (self._ldiag <= 0.).any()):
             raise ValueError("operator not positive definite")
-
         res = Field.from_random(random_type="normal", domain=self._domain,
                                 dtype=dtype)
-        res.local_data[()] /= np.sqrt(self._ldiag)
+        if from_inverse:
+            res.local_data[()] /= np.sqrt(self._ldiag)
+        else:
+            res.local_data[()] *= np.sqrt(self._ldiag)
         return res

nifty4/operators/endomorphic_operator.py

@@ -36,31 +36,22 @@ class EndomorphicOperator(LinearOperator):
         for endomorphic operators."""
         return self.domain
 
-    def draw_sample(self, dtype=np.float64):
+    def draw_sample(self, from_inverse=False, dtype=np.float64):
         """Generate a zero-mean sample
 
         Generates a sample from a Gaussian distribution with zero mean and
         covariance given by the operator.
 
+        Parameters
+        ----------
+        from_inverse : bool (default : False)
+            if True, the sample is drawn from the inverse of the operator
+        dtype : numpy datatype (default : numpy.float64)
+            the data type to be used for the sample
+
         Returns
         -------
         Field
             A sample from the Gaussian of given covariance.
         """
         raise NotImplementedError
-
-    def inverse_draw_sample(self, dtype=np.float64):
-        """Generates a zero-mean sample
-
-        Generates a sample from a Gaussian distribution with zero mean and
-        covariance given by the inverse of the operator.
-
-        Returns
-        -------
-            A sample from the Gaussian of given covariance
-        """
-        if self.capability & self.INVERSE_TIMES:
-            x = self.draw_sample(dtype)
-            return self.inverse_times(x)
-        else:
-            raise NotImplementedError

nifty4/operators/inversion_enabler.py

@@ -20,11 +20,11 @@ from ..minimization.quadratic_energy import QuadraticEnergy
 from ..minimization.iteration_controller import IterationController
 from ..field import Field
 from ..logger import logger
-from .linear_operator import LinearOperator
+from .endomorphic_operator import EndomorphicOperator
 import numpy as np
 
 
-class InversionEnabler(LinearOperator):
+class InversionEnabler(EndomorphicOperator):
     """Class which augments the capability of another operator object via
     numerical inversion.
 
@@ -38,16 +38,18 @@ class InversionEnabler(LinearOperator):
     inverter : :class:`Minimizer`
         The minimizer to use for the iterative numerical inversion.
         Typically, this is a :class:`ConjugateGradient` object.
-    preconditioner : :class:`LinearOperator`, optional
-        if not None, this operator is used as a preconditioner during the
-        iterative inversion, to accelerate convergence.
+    approximation : :class:`LinearOperator`, optional
+        if not None, this operator should be an approximation to `op`, which
+        supports the operation modes that `op` doesn't have. It is used as a
+        preconditioner during the iterative inversion, to accelerate
+        convergence.
     """
 
-    def __init__(self, op, inverter, preconditioner=None):
+    def __init__(self, op, inverter, approximation=None):
         super(InversionEnabler, self).__init__()
         self._op = op
         self._inverter = inverter
-        self._preconditioner = preconditioner
+        self._approximation = approximation
 
     @property
     def domain(self):
@@ -66,24 +68,21 @@ class InversionEnabler(LinearOperator):
         if self._op.capability & mode:
             return self._op.apply(x, mode)
 
-        def func(x):
-            return self._op.apply(x, self._inverseMode[mode])
-
         x0 = Field.zeros(self._tgt(mode), dtype=x.dtype)
-        energy = QuadraticEnergy(A=func, b=x, position=x0)
-        r, stat = self._inverter(energy, preconditioner=self._preconditioner)
+        invmode = self._modeTable[self.INVERSE_BIT][self._ilog[mode]]
+        invop = self._op._flip_modes(self._ilog[invmode])
+        prec = self._approximation
+        if prec is not None:
+            prec = prec._flip_modes(self._ilog[mode])
+        energy = QuadraticEnergy(A=invop, b=x, position=x0)
+        r, stat = self._inverter(energy, preconditioner=prec)
         if stat != IterationController.CONVERGED:
             logger.warning("Error detected during operator inversion")
         return r.position
 
-    def draw_sample(self, dtype=np.float64):
-        try:
-            return self._op.draw_sample(dtype)
-        except:
-            return self(self._op.inverse_draw_sample(dtype))
-
-    def inverse_draw_sample(self, dtype=np.float64):
+    def draw_sample(self, from_inverse=False, dtype=np.float64):
         try:
-            return self._op.inverse_draw_sample(dtype)
+            return self._op.draw_sample(from_inverse, dtype)
         except:
-            return self.inverse_times(self._op.draw_sample(dtype))
+            samp = self._op.draw_sample(not from_inverse, dtype)
+            return self.inverse_times(samp) if from_inverse else self(samp)

nifty4/operators/linear_operator.py

@@ -48,11 +48,16 @@ class LinearOperator(NiftyMetaBase()):
     by means of a single integer number.
     """
 
+    _ilog = (-1, 0, 1, -1, 2, -1, -1, -1, 3)
     _validMode = (False, True, True, False, True, False, False, False, True)
-    _inverseMode = (0, 4, 8, 0, 1, 0, 0, 0, 2)
-    _inverseCapability = (0, 4, 8, 12, 1, 5, 9, 13, 2, 6, 10, 14, 3, 7, 11, 15)
-    _adjointMode = (0, 2, 1, 0, 8, 0, 0, 0, 4)
-    _adjointCapability = (0, 2, 1, 3, 8, 10, 9, 11, 4, 6, 5, 7, 12, 14, 13, 15)
+    _modeTable = ((1, 2, 4, 8),
+                  (2, 1, 8, 4),
+                  (4, 8, 1, 2),
+                  (8, 4, 2, 1))
+    _capTable = ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15),
+                 (0, 2, 1, 3, 8, 10, 9, 11, 4, 6, 5, 7, 12, 14, 13, 15),
+                 (0, 4, 8, 12, 1, 5, 9, 13, 2, 6, 10, 14, 3, 7, 11, 15),
+                 (0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15))
     _addInverse = (0, 5, 10, 15, 5, 5, 15, 15, 10, 15, 10, 15, 15, 15, 15, 15)
     _backwards = 6
     _all_ops = 15
@@ -61,6 +66,8 @@ class LinearOperator(NiftyMetaBase()):
     INVERSE_TIMES = 4
     ADJOINT_INVERSE_TIMES = 8