merge

.gitlab-ci.yml

@@ -29,31 +29,23 @@ build_docker_from_cache:
     - docker build -t $CONTAINER_TEST_IMAGE .
     - docker push $CONTAINER_TEST_IMAGE
 
-test_python2_scalar:
-  stage: test
-  script:
-    - python setup.py install --user -f
-    - nosetests -q --with-coverage --cover-package=nifty4 --cover-branches --cover-erase
-    - >
-      coverage report | grep TOTAL | awk '{ print "TOTAL: "$6; }'
-
-test_python3_scalar:
-  stage: test
-  script:
-    - python3 setup.py install --user -f
-    - nosetests3 -q
-
-test_python2_mpi:
+test_python2_with_coverage:
   stage: test
   script:
     - python setup.py install --user -f
     - mpiexec -n 2 --bind-to none nosetests -q 2> /dev/null
+    - nosetests -q --with-coverage --cover-package=nifty4 --cover-erase
+    - >
+      coverage report --omit "*plotting*,*distributed_do*"
+    - >
+      coverage report --omit "*plotting*,*distributed_do*" | grep TOTAL | awk '{ print "TOTAL: "$4; }'
 
-test_python3_mpi:
+test_python3:
   stage: test
   script:
     - python3 setup.py install --user -f
     - mpiexec -n 2 --bind-to none nosetests3 -q 2> /dev/null
+    - nosetests3 -q
 
 pages:
   stage: release

Dockerfile

 FROM debian:testing-slim
 
-RUN apt-get update
-
-# Needed for gitlab tests
-RUN apt-get install -y git
-
-# Packages needed for NIFTy
-RUN apt-get install -y libfftw3-dev
-RUN apt-get install -y python python-pip python-dev python-future python-scipy
-RUN apt-get install -y python3 python3-pip python3-dev python3-future python3-scipy
-RUN pip install pyfftw
-RUN pip3 install pyfftw
-
-# Optional NIFTy dependencies
-RUN apt-get install -y openmpi-bin libopenmpi-dev python-mpi4py python3-mpi4py
-RUN pip install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
-RUN pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
-
-# Documentation build dependencies
-RUN apt-get install -y python-sphinx python-sphinx-rtd-theme python-numpydoc
-
-# Testing dependencies
-RUN apt-get install -y python-nose python-parameterized
-RUN apt-get install -y python3-nose python3-parameterized
-RUN pip install coverage
+RUN apt-get update && apt-get install -y \
+    # Needed for gitlab tests
+    git \
+    # Packages needed for NIFTy
+    libfftw3-dev \
+    python python-pip python-dev python-future python-scipy \
+    python3 python3-pip python3-dev python3-future python3-scipy \
+    # Documentation build dependencies
+    python-sphinx python-sphinx-rtd-theme python-numpydoc \
+    # Testing dependencies
+    python-nose python-parameterized \
+    python3-nose python3-parameterized \
+    # Optional NIFTy dependencies
+    openmpi-bin libopenmpi-dev python-mpi4py python3-mpi4py \
+    # Packages needed for NIFTy
+  && pip install pyfftw \
+  && pip3 install pyfftw \
+  # Optional NIFTy dependencies
+  && pip install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \
+  && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \
+  # Testing dependencies
+  && pip install coverage \
+  && rm -rf /var/lib/apt/lists/*
 
 # Create user (openmpi does not like to be run as root)
 RUN useradd -ms /bin/bash testinguser

nifty4/data_objects/distributed_do.py

@@ -59,32 +59,32 @@ class data_object(object):
         self._distaxis = distaxis
         self._data = data
 
-    def _sanity_checks(self):
-        # check whether the distaxis is consistent
-        if self._distaxis < -1 or self._distaxis >= len(self._shape):
-            raise ValueError
-        itmp = np.array(self._distaxis)
-        otmp = np.empty(ntask, dtype=np.int)
-        _comm.Allgather(itmp, otmp)
-        if np.any(otmp != self._distaxis):
-            raise ValueError
-        # check whether the global shape is consistent
-        itmp = np.array(self._shape)
-        otmp = np.empty((ntask, len(self._shape)), dtype=np.int)
-        _comm.Allgather(itmp, otmp)
-        for i in range(ntask):
-            if np.any(otmp[i, :] != self._shape):
-                raise ValueError
-        # check shape of local data
-        if self._distaxis < 0:
-            if self._data.shape != self._shape:
-                raise ValueError
-        else:
-            itmp = np.array(self._shape)
-            itmp[self._distaxis] = _shareSize(self._shape[self._distaxis],
-                                              ntask, rank)
-            if np.any(self._data.shape != itmp):
-                raise ValueError
+#     def _sanity_checks(self):
+#         # check whether the distaxis is consistent
+#         if self._distaxis < -1 or self._distaxis >= len(self._shape):
+#             raise ValueError
+#         itmp = np.array(self._distaxis)
+#         otmp = np.empty(ntask, dtype=np.int)
+#         _comm.Allgather(itmp, otmp)
+#         if np.any(otmp != self._distaxis):
+#             raise ValueError
+#         # check whether the global shape is consistent
+#         itmp = np.array(self._shape)
+#         otmp = np.empty((ntask, len(self._shape)), dtype=np.int)
+#         _comm.Allgather(itmp, otmp)
+#         for i in range(ntask):
+#             if np.any(otmp[i, :] != self._shape):
+#                 raise ValueError
+#         # check shape of local data
+#         if self._distaxis < 0:
+#             if self._data.shape != self._shape:
+#                 raise ValueError
+#         else:
+#             itmp = np.array(self._shape)
+#             itmp[self._distaxis] = _shareSize(self._shape[self._distaxis],
+#                                               ntask, rank)
+#             if np.any(self._data.shape != itmp):
+#                 raise ValueError
 
     @property
     def dtype(self):
@@ -322,6 +322,12 @@ def np_allreduce_sum(arr):
     return res
 
 
+def np_allreduce_min(arr):
+    res = np.empty_like(arr)
+    _comm.Allreduce(arr, res, MPI.MIN)
+    return res
+
+
 def distaxis(arr):
     return arr._distaxis
 

nifty4/data_objects/numpy_do.py

@@ -70,6 +70,10 @@ def np_allreduce_sum(arr):
     return arr
 
 
+def np_allreduce_min(arr):
+    return arr
+
+
 def distaxis(arr):
     return -1
 

nifty4/field.py

@@ -79,6 +79,13 @@ class Field(object):
         if locked:
             dobj.lock(self._val)
 
+    # prevent implicit conversion to bool
+    def __nonzero__(self):
+        raise TypeError("Field does not support implicit conversion to bool")
+
+    def __bool__(self):
+        raise TypeError("Field does not support implicit conversion to bool")
+
     @staticmethod
     def full(domain, val, dtype=None):
         """Creates a Field with a given domain, filled with a constant value.

nifty4/library/los_response.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.
+
 import numpy as np
 from scipy.special import erfc
 from scipy.sparse import coo_matrix

nifty4/library/nonlinearities.py

@@ -26,6 +26,9 @@ class Linear(object):
     def derivative(self, x):
         return Field.ones_like(x)
 
+    def hessian(self, x):
+        return Field.zeros_like(x)
+
 
 class Exponential(object):
     def __call__(self, x):
@@ -34,6 +37,9 @@ class Exponential(object):
     def derivative(self, x):
         return exp(x)
 
+    def hessian(self, x):
+        return exp(x)
+
 
 class Tanh(object):
     def __call__(self, x):
@@ -42,6 +48,9 @@ class Tanh(object):
     def derivative(self, x):
         return (1. - tanh(x)**2)
 
+    def hessian(self, x):
+        return - 2. * tanh(x) * (1. - tanh(x)**2)
+
 
 class PositiveTanh(object):
     def __call__(self, x):
@@ -49,3 +58,6 @@ class PositiveTanh(object):
 
     def derivative(self, x):
         return 0.5 * (1. - tanh(x)**2)
+
+    def hessian(self, x):
+        return - tanh(x) * (1. - tanh(x)**2)

nifty4/logger.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.
+
 def _logger_init():
     import logging
     from . import dobj

nifty4/minimization/steepest_descent.py

-from __future__ import division
 # 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

nifty4/multi/multi_field.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 ..field import Field
 import numpy as np
 from .multi_domain import MultiDomain
@@ -32,6 +50,13 @@ class MultiField(object):
     def dtype(self):
         return {key: val.dtype for key, val in self._val.items()}
 
+    @staticmethod
+    def from_random(random_type, domain, dtype=np.float64, **kwargs):
+        dtype = MultiField.build_dtype(dtype, domain)
+        return MultiField({key: Field.from_random(random_type, domain[key],
+                                                  dtype[key], **kwargs)
+                           for key in domain.keys()})
+
     def _check_domain(self, other):
         if other.domain != self.domain:
             raise ValueError("domains are incompatible.")

nifty4/operators/diagonal_operator.py

@@ -92,45 +92,46 @@ class DiagonalOperator(EndomorphicOperator):
             self._ldiag = self._ldiag.reshape(self._reshaper)
         else:
             self._ldiag = diagonal.local_data
+        self._update_diagmin()
+
+    def _update_diagmin(self):
         self._ldiag.flags.writeable = False
+        if not np.issubdtype(self._ldiag.dtype, np.complexfloating):
+            lmin = self._ldiag.min() if self._ldiag.size > 0 else 1.
+            self._diagmin = dobj.np_allreduce_min(np.array(lmin))[()]
 
-    def _skeleton(self, spc):
+    def _from_ldiag(self, spc, ldiag):
         res = DiagonalOperator.__new__(DiagonalOperator)
         res._domain = self._domain
         if self._spaces is None or spc is None:
             res._spaces = None
         else:
             res._spaces = tuple(set(self._spaces) | set(spc))
+        res._ldiag = ldiag
+        res._update_diagmin()
         return res
 
     def _scale(self, fct):
         if not np.isscalar(fct):
             raise TypeError("scalar value required")
-        res = self._skeleton(())
-        res._ldiag = self._ldiag*fct
-        return res
+        return self._from_ldiag((), self._ldiag*fct)
 
     def _add(self, sum):
         if not np.isscalar(sum):
             raise TypeError("scalar value required")
-        res = self._skeleton(())
-        res._ldiag = self._ldiag + sum
-        return res
+        return self._from_ldiag((), self._ldiag+sum)
 
     def _combine_prod(self, op):
         if not isinstance(op, DiagonalOperator):
             raise TypeError("DiagonalOperator required")
-        res = self._skeleton(op._spaces)
-        res._ldiag = self._ldiag*op._ldiag
-        return res
+        return self._from_ldiag(op._spaces, self._ldiag*op._ldiag)
 
     def _combine_sum(self, op, selfneg, opneg):
         if not isinstance(op, DiagonalOperator):
             raise TypeError("DiagonalOperator required")
-        res = self._skeleton(op._spaces)
-        res._ldiag = (self._ldiag * (-1 if selfneg else 1) +
-                      op._ldiag * (-1 if opneg else 1))
-        return res
+        tdiag = (self._ldiag * (-1 if selfneg else 1) +
+                 op._ldiag * (-1 if opneg else 1))
+        return self._from_ldiag(op._spaces, tdiag)
 
     def apply(self, x, mode):
         self._check_input(x, mode)
@@ -166,20 +167,20 @@ class DiagonalOperator(EndomorphicOperator):
             return self
         if trafo == ADJ and np.issubdtype(self._ldiag.dtype, np.floating):
             return self
-        res = self._skeleton(())
         if trafo == ADJ:
-            res._ldiag = self._ldiag.conjugate()
+            return self._from_ldiag((), self._ldiag.conjugate())
         elif trafo == INV:
-            res._ldiag = 1./self._ldiag
+            return self._from_ldiag((), 1./self._ldiag)
         elif trafo == ADJ | INV:
-            res._ldiag = 1./self._ldiag.conjugate()
-        else:
-            raise ValueError("invalid operator transformation")
-        return res
+            return self._from_ldiag((), 1./self._ldiag.conjugate())
+        raise ValueError("invalid operator transformation")
 
     def draw_sample(self, from_inverse=False, dtype=np.float64):
-        if (np.issubdtype(self._ldiag.dtype, np.complexfloating) or
-                (self._ldiag <= 0.).any()):
+        if np.issubdtype(self._ldiag.dtype, np.complexfloating):
+            raise ValueError("operator not positive definite")
+        if self._diagmin < 0.:
+            raise ValueError("operator not positive definite")
+        if self._diagmin == 0. and from_inverse:
             raise ValueError("operator not positive definite")
         res = Field.from_random(random_type="normal", domain=self._domain,
                                 dtype=dtype)

nifty4/operators/sandwich_operator.py

@@ -17,6 +17,7 @@
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
 from .endomorphic_operator import EndomorphicOperator
+from .scaling_operator import ScalingOperator
 import numpy as np
 
 
@@ -31,11 +32,15 @@ class SandwichOperator(EndomorphicOperator):
         the cheese part
     """
 
-    def __init__(self, bun, cheese):
+    def __init__(self, bun, cheese=None):
         super(SandwichOperator, self).__init__()
         self._bun = bun
-        self._cheese = cheese
-        self._op = bun.adjoint*cheese*bun
+        if cheese is None:
+            self._cheese = ScalingOperator(1., bun.target)
+            self._op = bun.adjoint*bun
+        else:
+            self._cheese = cheese
+            self._op = bun.adjoint*cheese*bun
 
     @property
     def domain(self):

nifty4/operators/scaling_operator.py

@@ -19,6 +19,7 @@
 from __future__ import division
 import numpy as np
 from ..field import Field
+from ..multi.multi_field import MultiField
 from ..domain_tuple import DomainTuple
 from .endomorphic_operator import EndomorphicOperator
 
@@ -61,7 +62,7 @@ class ScalingOperator(EndomorphicOperator):
         if self._factor == 1.:
             return x.copy()
         if self._factor == 0.:
-            return Field.zeros_like(x)
+            return x.zeros_like(x)
 
         if mode == self.TIMES:
             return x*self._factor
@@ -98,10 +99,11 @@ class ScalingOperator(EndomorphicOperator):
 
     def draw_sample(self, from_inverse=False, dtype=np.float64):
         fct = self._factor
-        if fct.imag != 0. or fct.real < 0.: 
+        if fct.imag != 0. or fct.real < 0.:
             raise ValueError("operator not positive definite")
-        if fct.real ==0 and from_inverse:
+        if fct.real == 0. and from_inverse:
             raise ValueError("operator not positive definite")
         fct = 1./np.sqrt(fct) if from_inverse else np.sqrt(fct)
-        return Field.from_random(
+        cls = Field if isinstance(self._domain, DomainTuple) else MultiField
+        return cls.from_random(
            random_type="normal", domain=self._domain, std=fct, dtype=dtype)

nifty4/plotting/plot.py

@@ -94,10 +94,10 @@ def _makeplot(name):
 def _limit_xy(**kwargs):
     import matplotlib.pyplot as plt
     x1, x2, y1, y2 = plt.axis()
-    x1 = _get_kw("xmin", x1, **kwargs)
-    x2 = _get_kw("xmax", x2, **kwargs)
-    y1 = _get_kw("ymin", y1, **kwargs)
-    y2 = _get_kw("ymax", y2, **kwargs)
+    x1 = kwargs.pop("xmin", x1)
+    x2 = kwargs.pop("xmax", x2)
+    y1 = kwargs.pop("ymin", y1)
+    y2 = kwargs.pop("ymax", y2)
     plt.axis((x1, x2, y1, y2))
 
 
@@ -165,12 +165,6 @@ def _register_cmaps():
     plt.register_cmap(cmap=LinearSegmentedColormap("Plus Minus", pm_cmap))
 
 
-def _get_kw(kwname, kwdefault=None, **kwargs):
-    if kwargs.get(kwname) is not None:
-        return kwargs.get(kwname)
-    return kwdefault
-
-
 def plot(f, **kwargs):
     import matplotlib.pyplot as plt
     _register_cmaps()
@@ -192,19 +186,19 @@ def plot(f, **kwargs):
                     (isinstance(dom[0], RGSpace) and len(dom[0].shape) == 1)):
                 raise ValueError("PowerSpace or 1D RGSpace required")
 
-    label = _get_kw("label", None, **kwargs)
+    label = kwargs.pop("label", None)
     if label is None:
         label = [None] * len(f)
     if not isinstance(label, list):
         label = [label]
 
-    linewidth = _get_kw("linewidth", None, **kwargs)
+    linewidth = kwargs.pop("linewidth", None)
     if linewidth is None:
         linewidth = [1.] * len(f)
     if not isinstance(linewidth, list):
         linewidth = [linewidth]
 
-    alpha = _get_kw("alpha", None, **kwargs)
+    alpha = kwargs.pop("alpha", None)
     if alpha is None:
         alpha = [None] * len(f)
     if not isinstance(alpha, list):
@@ -214,13 +208,13 @@ def plot(f, **kwargs):
     fig = plt.figure()
     ax = fig.add_subplot(1, 1, 1)
 
-    xsize = _get_kw("xsize", 6, **kwargs)
-    ysize = _get_kw("ysize", 6, **kwargs)
+    xsize = kwargs.pop("xsize", 6)
+    ysize = kwargs.pop("ysize", 6)
     fig.set_size_inches(xsize, ysize)
-    ax.set_title(_get_kw("title", "", **kwargs))
-    ax.set_xlabel(_get_kw("xlabel", "", **kwargs))
-    ax.set_ylabel(_get_kw("ylabel", "", **kwargs))
-    cmap = _get_kw("colormap", plt.rcParams['image.cmap'], **kwargs)
+    ax.set_title(kwargs.pop("title", ""))
+    ax.set_xlabel(kwargs.pop("xlabel", ""))
+    ax.set_ylabel(kwargs.pop("ylabel", ""))
+    cmap = kwargs.pop("colormap", plt.rcParams['image.cmap'])
     if isinstance(dom, RGSpace):
         if len(dom.shape) == 1:
             npoints = dom.shape[0]