Merge branch 'mmm5' of gitlab.mpcdf.mpg.de:ift/NIFTy into mmm5

.gitlab-ci.yml

@@ -14,7 +14,7 @@ before_script:
   - >
     apt-get install -y build-essential python python-pip python-dev git
     autoconf gsl-bin libgsl-dev wget python-numpy cython
-  - pip install -r ci/requirements_base.txt
+  - pip install --upgrade -r ci/requirements_base.txt
   - chmod +x ci/*.sh
 
 test_min:
@@ -26,7 +26,8 @@ test_min:
 test_mpi:
   stage: test
   script:
-    - apt-get install -y openmpi-bin libopenmpi-dev python-mpi4py
+    - apt-get install -y openmpi-bin libopenmpi-dev
+    - pip install mpi4py
     - ci/install_pyHealpix.sh
     - python setup.py build_ext --inplace
     - nosetests -vv
@@ -38,8 +39,10 @@ test_mpi_fftw:
     - >
       apt-get install -y libatlas-base-dev libfftw3-bin libfftw3-dev
       libfftw3-double3 libfftw3-long3 libfftw3-mpi-dev libfftw3-mpi3
-      libfftw3-quad3 libfftw3-single3 python-mpi4py python-pyfftw
+      libfftw3-quad3 libfftw3-single3
+    - pip install mpi4py
     - ci/install_pyHealpix.sh
+    - ci/install_pyfftw.sh
     - python setup.py build_ext --inplace
     - nosetests -vv
 
@@ -53,8 +56,11 @@ test_mpi_fftw_hdf5:
       libfftw3-quad3 libfftw3-single3
     - >
       apt-get install -y libhdf5-10 libhdf5-dev libhdf5-openmpi-10
-      libhdf5-openmpi-dev hdf5-tools python-mpi4py python-pyfftw python-h5py
+      libhdf5-openmpi-dev hdf5-tools
+    - pip install mpi4py
     - ci/install_pyHealpix.sh
+    - ci/install_h5py.sh
+    - ci/install_pyfftw.sh
     - python setup.py build_ext --inplace
     - nosetests -vv --with-coverage --cover-package=nifty --cover-branches
     - >

Dockerfile

@@ -4,19 +4,23 @@ FROM ubuntu:latest
 RUN \
     apt-get update && \
     apt-get install -y build-essential python python-pip python-dev git \
-    gfortran autoconf gsl-bin libgsl-dev python-matplotlib openmpi-bin \
+    autoconf gsl-bin libgsl-dev openmpi-bin wget \
     libopenmpi-dev libatlas-base-dev libfftw3-bin libfftw3-dev \
     libfftw3-double3 libfftw3-long3 libfftw3-mpi-dev libfftw3-mpi3 \
     libfftw3-quad3 libfftw3-single3 libhdf5-10 libhdf5-dev \
-    libhdf5-openmpi-10 libhdf5-openmpi-dev hdf5-tools python-h5py python-pyfftw
+    libhdf5-openmpi-10 libhdf5-openmpi-dev hdf5-tools
 
 # python dependencies
 ADD ci/requirements.txt /tmp/requirements.txt
-RUN pip install -r /tmp/requirements.txt
+RUN pip install --upgrade -r /tmp/requirements.txt
 
-# install pyHealpix
+# install pyHealpix, pyfftw and h5py
 ADD ci/install_pyHealpix.sh /tmp/install_pyHealpix.sh
 RUN cd /tmp && chmod +x install_pyHealpix.sh && ./install_pyHealpix.sh
+ADD ci/install_pyfftw.sh /tmp/install_pyfftw.sh
+RUN cd /tmp && chmod +x install_pyfftw.sh && ./install_pyfftw.sh
+ADD ci/install_h5py.sh /tmp/install_h5py.sh
+RUN cd /tmp && chmod +x install_h5py.sh && ./install_h5py.sh
 
 # copy sources and install nifty
 COPY . /tmp/NIFTy

ci/install_h5py.sh

+#!/bin/bash
+
+wget https://api.github.com/repos/h5py/h5py/tags -O - | grep tarball_url | grep -v rc | head -n 1 | cut -d '"' -f 4 | wget -i - -O h5py.tar.gz
+tar xzf h5py.tar.gz
+cd h5py-h5py*
+export CC=mpicc
+export HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/openmpi
+python setup.py configure --mpi
+python setup.py build
+python setup.py install
+cd ..
+rm -r h5py-h5py*
+rm h5py.tar.gz

ci/install_pyfftw.sh

+#!/bin/bash
+
+git clone -b mpi https://github.com/fredRos/pyFFTW.git
+cd pyFFTW/
+CC=mpicc python setup.py build_ext install
+cd ..
+rm -r pyFFTW

ci/requirements.txt

+numpy
+cython
+mpi4py
+matplotlib
+ipython==5.3.0
 nose
 parameterized
 coverage

ci/requirements_base.txt

+numpy
 nose
 parameterized
 coverage

nifty/basic_arithmetics.py

@@ -31,10 +31,8 @@ def _math_helper(x, function):
         result_val = x.val.apply_scalar_function(function)
         result = x.copy_empty(dtype=result_val.dtype)
         result.val = result_val
-
     elif isinstance(x, distributed_data_object):
         result = x.apply_scalar_function(function, inplace=False)
-
     else:
         result = function(np.asarray(x))
 

nifty/config/nifty_config.py

@@ -17,7 +17,6 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 import os
-from distutils.version import LooseVersion as lv
 
 import numpy as np
 import keepers
@@ -36,22 +35,6 @@ variable_fft_module = keepers.Variable(
                                lambda z: (('pyfftw' in dependency_injector)
                                           if z == 'fftw' else True))
 
-def _pyHealpix_validator(use_pyHealpix):
-    if not isinstance(use_pyHealpix, bool):
-        return False
-    if not use_pyHealpix:
-        return True
-    if 'pyHealpix' not in dependency_injector:
-        return False
-    pyHealpix = dependency_injector['pyHealpix']
-    return True
-
-
-variable_use_pyHealpix = keepers.Variable(
-                          'use_pyHealpix',
-                          [True, False],
-                          _pyHealpix_validator,
-                          genus='boolean')
 
 def _dtype_validator(dtype):
     try:
@@ -69,7 +52,7 @@ variable_default_field_dtype = keepers.Variable(
 
 variable_default_distribution_strategy = keepers.Variable(
                               'default_distribution_strategy',
-                              ['fftw', 'equal'],
+                              ['fftw', 'equal', 'not'],
                               lambda z: (('pyfftw' in dependency_injector)
                                          if z == 'fftw' else True),
                               genus='str')
@@ -77,7 +60,6 @@ variable_default_distribution_strategy = keepers.Variable(
 nifty_configuration = keepers.get_Configuration(
                  name='NIFTy',
                  variables=[variable_fft_module,
-                            variable_use_pyHealpix,
                             variable_default_field_dtype,
                             variable_default_distribution_strategy],
                  file_name='NIFTy.conf',

nifty/minimization/vl_bfgs.py

@@ -77,8 +77,6 @@ class InformationStore(object):
         self._sy_store = {}
         self._yy_store = {}
 
-#        self.dot_matrix = {}
-
     @property
     def history_length(self):
         return min(self.k, self.max_history_length)
@@ -193,48 +191,6 @@ class InformationStore(object):
         self.last_x = x.copy()
         self.last_gradient = gradient.copy()
 
-#
-#        k = self.k
-#        m = self.actual_history_length
-#        big_m = self.history_length
-#
-#        # compute dot products
-#        for i in xrange(k-1, k-m-1, -1):
-#            # new_s with s
-#            key = (big_m+m, big_m+1+i)
-#            self.dot_matrix[key] = new_s.dot(self.s[i])
-#
-#            # new_s with y
-#            key = (big_m+m, i+1)
-#            self.dot_matrix[key] = new_s.dot(self.y[i])
-#
-#            # new_y with s
-#            if i != k-1:
-#                key = (big_m+1+i, k)
-#                self.dot_matrix[key] = new_y.dot(self.s[i])
-#
-#            # new_y with y
-#            # actually key = (i+1, k) but the convention is that the first
-#            # index is larger than the second one
-#            key = (k, i+1)
-#            self.dot_matrix[key] = new_y.dot(self.y[i])
-#
-#            # gradient with s
-#            key = (big_m+1+i, 0)
-#            self.dot_matrix[key] = gradient.dot(self.s[i])
-#
-#            # gradient with y
-#            key = (i+1, 0)
-#            self.dot_matrix[key] = gradient.dot(self.y[i])
-#
-#        # gradient with gradient
-#        key = (0, 0)
-#        self.dot_matrix[key] = gradient.dot(gradient)
-#
-#        self.last_x = x
-#        self.last_gradient = gradient
-#
-
 
 class LimitedList(object):
     def __init__(self, history_length):

nifty/operators/fft_operator/fft_operator.py

@@ -55,33 +55,25 @@ class FFTOperator(LinearOperator):
     def __init__(self, domain=(), target=None, module=None,
                  domain_dtype=None, target_dtype=None):
 
-        self._domain = self._parse_domain(domain)
-
         # Initialize domain and target
+
+        self._domain = self._parse_domain(domain)
         if len(self.domain) != 1:
-            raise ValueError(
-                    'ERROR: TransformationOperator accepts only exactly one '
-                    'space as input domain.')
+            raise ValueError("TransformationOperator accepts only exactly one "
+                             "space as input domain.")
 
         if target is None:
             target = (self.get_default_codomain(self.domain[0]), )
         self._target = self._parse_domain(target)
+        if len(self.target) != 1:
+            raise ValueError("TransformationOperator accepts only exactly one "
+                             "space as output target.")
 
         # Create transformation instances
-        try:
-            forward_class = self.transformation_dictionary[
+        forward_class = self.transformation_dictionary[
                 (self.domain[0].__class__, self.target[0].__class__)]
-        except KeyError:
-            raise ValueError(
-                "No forward transformation for domain-target pair "
-                "found.")
-        try:
-            backward_class = self.transformation_dictionary[
+        backward_class = self.transformation_dictionary[
                 (self.target[0].__class__, self.domain[0].__class__)]
-        except KeyError:
-            raise ValueError(
-                "No backward transformation for domain-target pair "
-                "found.")
 
         self._forward_transformation = TransformationCache.create(
             forward_class, self.domain[0], self.target[0], module=module)

nifty/operators/fft_operator/transformations/gllmtransformation.py

@@ -21,7 +21,7 @@ import numpy as np
 from nifty.config import dependency_injector as gdi
 from nifty import GLSpace, LMSpace
 from slicing_transformation import SlicingTransformation
-import lm_transformation_factory as ltf
+import lm_transformation_factory
 
 pyHealpix = gdi.get('pyHealpix')
 
@@ -58,30 +58,21 @@ class GLLMTransformation(SlicingTransformation):
         """
 
         if not isinstance(domain, GLSpace):
-            raise TypeError(
-                "domain needs to be a GLSpace")
+            raise TypeError("domain needs to be a GLSpace")
 
         nlat = domain.nlat
         lmax = nlat - 1
-        mmax = nlat - 1
-        if domain.dtype == np.dtype('float32'):
-            return_dtype = np.float32
-        else:
-            return_dtype = np.float64
 
-        result = LMSpace(lmax=lmax, dtype=return_dtype)
-        cls.check_codomain(domain, result)
+        result = LMSpace(lmax=lmax, dtype=domain.dtype)
         return result
 
-    @staticmethod
-    def check_codomain(domain, codomain):
+    @classmethod
+    def check_codomain(cls, domain, codomain):
         if not isinstance(domain, GLSpace):
-            raise TypeError(
-                "domain is not a GLSpace")
+            raise TypeError("domain is not a GLSpace")
 
         if not isinstance(codomain, LMSpace):
-            raise TypeError(
-                "codomain must be a LMSpace.")
+            raise TypeError("codomain must be a LMSpace.")
 
         nlat = domain.nlat
         nlon = domain.nlon
@@ -89,18 +80,15 @@ class GLLMTransformation(SlicingTransformation):
         mmax = codomain.mmax
 
         if lmax != mmax:
-            raise ValueError(
-                'ERROR: codomain has lmax != mmax.')
+            cls.Logger.warn("Unrecommended: codomain has lmax != mmax.")
 
         if lmax != nlat - 1:
-            raise ValueError(
-                'ERROR: codomain has lmax != nlat - 1.')
+            cls.Logger.warn("Unrecommended: codomain has lmax != nlat - 1.")
 
-        if nlon != 2 * nlat - 1:
-            raise ValueError(
-                'ERROR: domain has nlon != 2 * nlat - 1.')
+        if nlon != 2*nlat - 1:
+            cls.Logger.warn("Unrecommended: domain has nlon != 2*nlat - 1.")
 
-        return None
+        super(GLLMTransformation, cls).check_codomain(domain, codomain)
 
     def _transformation_of_slice(self, inp, **kwargs):
         nlat = self.domain.nlat
@@ -108,20 +96,20 @@ class GLLMTransformation(SlicingTransformation):
         lmax = self.codomain.lmax
         mmax = self.codomain.mmax
 
-        sjob=pyHealpix.sharpjob_d()
-        sjob.set_Gauss_geometry(nlat,nlon)
-        sjob.set_triangular_alm_info(lmax,mmax)
+        sjob = pyHealpix.sharpjob_d()
+        sjob.set_Gauss_geometry(nlat, nlon)
+        sjob.set_triangular_alm_info(lmax, mmax)
         if issubclass(inp.dtype.type, np.complexfloating):
-
             [resultReal, resultImag] = [sjob.map2alm(x)
                                         for x in (inp.real, inp.imag)]
 
-            [resultReal, resultImag] = [ltf.buildIdx(x, lmax=lmax)
+            [resultReal,
+             resultImag] = [lm_transformation_factory.buildIdx(x, lmax=lmax)
                                         for x in [resultReal, resultImag]]
 
             result = self._combine_complex_result(resultReal, resultImag)
         else:
             result = sjob.map2alm(inp)
-            result = ltf.buildIdx(result, lmax=lmax)
+            result = lm_transformation_factory.buildIdx(result, lmax=lmax)
 
         return result

nifty/operators/fft_operator/transformations/hplmtransformation.py

@@ -22,7 +22,7 @@ from nifty.config import dependency_injector as gdi
 from nifty import HPSpace, LMSpace
 from slicing_transformation import SlicingTransformation
 
-import lm_transformation_factory as ltf
+import lm_transformation_factory
 
 pyHealpix = gdi.get('pyHealpix')
 
@@ -59,49 +59,50 @@ class HPLMTransformation(SlicingTransformation):
         """
 
         if not isinstance(domain, HPSpace):
-            raise TypeError(
-                "domain needs to be a HPSpace")
+            raise TypeError("domain needs to be a HPSpace")
 
-        lmax = 2 * domain.nside
+        lmax = 2*domain.nside
 
-        result = LMSpace(lmax=lmax, dtype=np.dtype('float64'))
-        cls.check_codomain(domain, result)
+        result = LMSpace(lmax=lmax, dtype=domain.dtype)
         return result
 
-    @staticmethod
-    def check_codomain(domain, codomain):
+    @classmethod
+    def check_codomain(cls, domain, codomain):
         if not isinstance(domain, HPSpace):
-            raise TypeError(
-                'ERROR: domain is not a HPSpace')
+            raise TypeError("domain is not a HPSpace")
 
         if not isinstance(codomain, LMSpace):
-            raise TypeError(
-                'ERROR: codomain must be a LMSpace.')
+            raise TypeError("codomain must be a LMSpace.")
 
-        nside = domain.nside
         lmax = codomain.lmax
+        nside = domain.nside
+
+        if lmax != 2*nside:
+            cls.Logger.warn("Unrecommended: lmax != 2*nside.")
 
-        return None
+        super(HPLMTransformation, cls).check_codomain(domain, codomain)
 
     def _transformation_of_slice(self, inp, **kwargs):
+        nside = self.domain.nside
         lmax = self.codomain.lmax
         mmax = lmax
         nside= self.domain.nside
 
-        sjob=pyHealpix.sharpjob_d()
+        sjob = pyHealpix.sharpjob_d()
         sjob.set_Healpix_geometry(nside)
-        sjob.set_triangular_alm_info(lmax,mmax)
+        sjob.set_triangular_alm_info(lmax, mmax)
         if issubclass(inp.dtype.type, np.complexfloating):
             [resultReal, resultImag] = [sjob.map2alm(x)
                                         for x in (inp.real, inp.imag)]
 
-            [resultReal, resultImag] = [ltf.buildIdx(x, lmax=lmax)
-                                        for x in [resultReal, resultImag]]
+            [resultReal,
+             resultImag] = [lm_transformation_factory.buildIdx(x, lmax=lmax)
+                            for x in [resultReal, resultImag]]
 
             result = self._combine_complex_result(resultReal, resultImag)
 
         else:
             result = sjob.map2alm(inp)
-            result = ltf.buildIdx(result, lmax=lmax)
+            result = lm_transformation_factory.buildIdx(result, lmax=lmax)
 
         return result

nifty/operators/fft_operator/transformations/lm_transformation_factory.py

 import numpy as np
 
-def buildLm(inp, **kwargs):
-    if inp.dtype == np.dtype('float32'):
-        return _buildLm_f(inp, **kwargs)
-    else:
-        return _buildLm(inp, **kwargs)
 
-def buildIdx(inp, **kwargs):
-    if inp.dtype == np.dtype('complex64'):
-        return _buildIdx_f(inp, **kwargs)
-    else:
-        return _buildIdx(inp, **kwargs)
+def buildLm(nr, lmax):
+    new_dtype = np.result_type(nr.dtype, np.complex64)
 
-def _buildIdx_f(nr, lmax):
-    size = (lmax+1)*(lmax+1)
+    size = (len(nr)-lmax-1)/2+lmax+1
+    res = np.zeros([size], dtype=new_dtype)
+    res[0:lmax+1] = nr[0:lmax+1]
+    res[lmax+1:] = np.sqrt(0.5)*(nr[lmax+1::2] + 1j*nr[lmax+2::2])
+    return res
 
-    final=np.zeros([size], dtype=np.float32)
-    final[0:lmax+1] = nr[0:lmax+1].real
-    final[lmax+1::2] = np.sqrt(2)*nr[lmax+1:].real
-    final[lmax+2::2] = np.sqrt(2)*nr[lmax+1:].imag
-    return final
 
-def _buildIdx(nr, lmax):
-    size = (lmax+1)*(lmax+1)
+def buildIdx(nr, lmax):
+    if nr.dtype == np.dtype('complex64'):
+        new_dtype = np.float32
+    elif nr.dtype == np.dtype('complex128'):
+        new_dtype = np.float64
+    elif nr.dtype == np.dtype('complex256'):
+        new_dtype = np.float128
+    else:
+        raise TypeError("dtype of nr not supported.")
 
-    final=np.zeros([size], dtype=np.float64)
+    size = (lmax+1)*(lmax+1)
+    final = np.zeros([size], dtype=new_dtype)
     final[0:lmax+1] = nr[0:lmax+1].real
     final[lmax+1::2] = np.sqrt(2)*nr[lmax+1:].real
     final[lmax+2::2] = np.sqrt(2)*nr[lmax+1:].imag
     return final
-
-def _buildLm_f(nr, lmax):
-    size = (len(nr)-lmax-1)/2+lmax+1
-
-    res=np.zeros([size], dtype=np.complex64)
-    res[0:lmax+1] = nr[0:lmax+1]
-    res[lmax+1:] = np.sqrt(0.5)*(nr[lmax+1::2] + 1j*nr[lmax+2::2])
-    return res
-
-def _buildLm(nr, lmax):
-    size = (len(nr)-lmax-1)/2+lmax+1
-
-    res=np.zeros([size], dtype=np.complex128)
-    res[0:lmax+1] = nr[0:lmax+1]
-    res[lmax+1:] = np.sqrt(0.5)*(nr[lmax+1::2] + 1j*nr[lmax+2::2])
-    return res

nifty/operators/fft_operator/transformations/lmgltransformation.py

@@ -21,7 +21,7 @@ from nifty.config import dependency_injector as gdi
 from nifty import GLSpace, LMSpace
 
 from slicing_transformation import SlicingTransformation
-import lm_transformation_factory as ltf
+import lm_transformation_factory
 
 pyHealpix = gdi.get('pyHealpix')
 
@@ -64,30 +64,21 @@ class LMGLTransformation(SlicingTransformation):
                    `arXiv:1303.4945 <http://www.arxiv.org/abs/1303.4945>`_
         """
         if not isinstance(domain, LMSpace):
-            raise TypeError(
-                'ERROR: domain needs to be a LMSpace')
-
-        if domain.dtype is np.dtype('float32'):
-            new_dtype = np.float32
-        else:
-            new_dtype = np.float64
+            raise TypeError("domain needs to be a LMSpace")
 
         nlat = domain.lmax + 1
-        nlon = domain.lmax * 2 + 1
+        nlon = domain.lmax*2 + 1
 
-        result = GLSpace(nlat=nlat, nlon=nlon, dtype=new_dtype)
-        cls.check_codomain(domain, result)
+        result = GLSpace(nlat=nlat, nlon=nlon, dtype=domain.dtype)
         return result
 
-    @staticmethod
-    def check_codomain(domain, codomain):
+    @classmethod
+    def check_codomain(cls, domain, codomain):
         if not isinstance(domain, LMSpace):
-            raise TypeError(
-                'ERROR: domain is not a LMSpace')
+            raise TypeError("domain is not a LMSpace")
 
         if not isinstance(codomain, GLSpace):
-            raise TypeError(
-                'ERROR: codomain must be a GLSpace.')
+            raise TypeError("codomain must be a GLSpace.")
 
         nlat = codomain.nlat
         nlon = codomain.nlon
@@ -95,18 +86,15 @@ class LMGLTransformation(SlicingTransformation):
         mmax = domain.mmax
 
         if lmax != mmax:
-            raise ValueError(
-                'ERROR: domain has lmax != mmax.')
+            cls.Logger.warn("Unrecommended: codomain has lmax != mmax.")
 
         if nlat != lmax + 1:
-            raise ValueError(
-                'ERROR: codomain has nlat != lmax + 1.')
+            cls.Logger.warn("Unrecommended: codomain has nlat != lmax + 1.")
 
-        if nlon != 2 * lmax + 1:
-            raise ValueError(
-                'ERROR: domain has nlon != 2 * lmax + 1.')
+        if nlon != 2*lmax + 1:
+            cls.Logger.warn("Unrecommended: domain has nlon != 2*lmax + 1.")
 
-        return None
+        super(LMGLTransformation, cls).check_codomain(domain, codomain)
 
     def _transformation_of_slice(self, inp, **kwargs):
         nlat = self.codomain.nlat
@@ -114,11 +102,12 @@ class LMGLTransformation(SlicingTransformation):
         lmax = self.domain.lmax
         mmax = self.domain.mmax
 
-        sjob=pyHealpix.sharpjob_d()
-        sjob.set_Gauss_geometry(nlat,nlon)
-        sjob.set_triangular_alm_info(lmax,mmax)
+        sjob = pyHealpix.sharpjob_d()
+        sjob.set_Gauss_geometry(nlat, nlon)
+        sjob.set_triangular_alm_info(lmax, mmax)
         if issubclass(inp.dtype.type, np.complexfloating):
-            [resultReal, resultImag] = [ltf.buildLm(x, lmax=lmax)
+            [resultReal,
+             resultImag] = [lm_transformation_factory.buildLm(x, lmax=lmax)
                                         for x in (inp.real, inp.imag)]