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

nifty/operators/fft_operator/fft_operator.py

@@ -59,14 +59,14 @@ class FFTOperator(LinearOperator):
             forward_class = self.transformation_dictionary[
                 (self.domain[0].__class__, self.target[0].__class__)]
         except KeyError:
-            raise TypeError(
+            raise ValueError(
                 "No forward transformation for domain-target pair "
                 "found.")
         try:
             backward_class = self.transformation_dictionary[
                 (self.target[0].__class__, self.domain[0].__class__)]
         except KeyError:
-            raise TypeError(
+            raise ValueError(
                 "No backward transformation for domain-target pair "
                 "found.")
 
@@ -156,13 +156,13 @@ class FFTOperator(LinearOperator):
         try:
             codomain_class = cls.default_codomain_dictionary[domain_class]
         except KeyError:
-            raise TypeError("unknown domain")
+            raise ValueError("Unknown domain")
 
         try:
             transform_class = cls.transformation_dictionary[(domain_class,
                                                              codomain_class)]
         except KeyError:
-            raise TypeError(
+            raise ValueError(
                 "No transformation for domain-codomain pair found.")
 
         return transform_class.get_codomain(domain)

nifty/operators/fft_operator/transformations/rg_transforms.py

@@ -164,24 +164,26 @@ class FFTW(Transform):
             self.centering_mask_dict[temp_id] = centering_mask
         return self.centering_mask_dict[temp_id]
 
-    def _get_transform_info(self, domain, codomain, local_shape,
+    def _get_transform_info(self, domain, codomain, axes, local_shape,
                             local_offset_Q, is_local, transform_shape=None,
                             **kwargs):
         # generate a id-tuple which identifies the domain-codomain setting
         temp_id = (domain.__hash__() ^
                    (101 * codomain.__hash__()) ^
-                   (211 * transform_shape.__hash__()))
+                   (211 * transform_shape.__hash__()) ^
+                   (131 * is_local.__hash__())
+                   )
 
         # generate the plan_and_info object if not already there
         if temp_id not in self.info_dict:
             if is_local:
                 self.info_dict[temp_id] = FFTWLocalTransformInfo(
-                    domain, codomain, local_shape,
+                    domain, codomain, axes, local_shape,
                     local_offset_Q, self, **kwargs
                 )
             else:
                 self.info_dict[temp_id] = FFTWMPITransfromInfo(
-                    domain, codomain, local_shape,
+                    domain, codomain, axes, local_shape,
                     local_offset_Q, self, transform_shape, **kwargs
                 )
 
@@ -248,17 +250,16 @@ class FFTW(Transform):
         # val must be numpy array or d2o with slicing distributor
         ###
 
-        local_offset_Q = False
         try:
             local_val = val.get_local_data(copy=False)
-            if axes is None or 0 in axes:
-                local_offset_Q = val.distributor.local_shape[0] % 2
         except(AttributeError):
             local_val = val
+
         current_info = self._get_transform_info(self.domain,
                                                 self.codomain,
+                                                axes,
                                                 local_shape=local_val.shape,
-                                                local_offset_Q=local_offset_Q,
+                                                local_offset_Q=False,
                                                 is_local=True,
                                                 **kwargs)
 
@@ -309,14 +310,10 @@ class FFTW(Transform):
 
     def _mpi_transform(self, val, axes, **kwargs):
 
-        if axes is None or 0 in axes:
-            local_offset_list = np.cumsum(
-                np.concatenate([[0, ], val.distributor.all_local_slices[:, 2]])
-            )
-            local_offset_Q = bool(
-                local_offset_list[val.distributor.comm.rank] % 2)
-        else:
-            local_offset_Q = False
+        local_offset_list = np.cumsum(
+            np.concatenate([[0, ], val.distributor.all_local_slices[:, 2]])
+        )
+        local_offset_Q = bool(local_offset_list[val.distributor.comm.rank] % 2)
 
         return_val = val.copy_empty(global_shape=val.shape,
                                     dtype=self.codomain.dtype)
@@ -362,6 +359,7 @@ class FFTW(Transform):
                 current_info = self._get_transform_info(
                     self.domain,
                     self.codomain,
+                    axes,
                     local_shape=val.local_shape,
                     local_offset_Q=local_offset_Q,
                     is_local=False,
@@ -446,20 +444,22 @@ class FFTW(Transform):
 
 
 class FFTWTransformInfo(object):
-    def __init__(self, domain, codomain, local_shape,
+    def __init__(self, domain, codomain, axes, local_shape,
                  local_offset_Q, fftw_context, **kwargs):
         if pyfftw is None:
             raise ImportError("The module pyfftw is needed but not available.")
 
-        self.cmask_domain = fftw_context.get_centering_mask(
-            domain.zerocenter,
-            local_shape,
-            local_offset_Q)
+        shape = (local_shape if axes is None else
+                 [y for x, y in enumerate(local_shape) if x in axes])
+
+        self.cmask_domain = fftw_context.get_centering_mask(domain.zerocenter,
+                                                            shape,
+                                                            local_offset_Q)
 
         self.cmask_codomain = fftw_context.get_centering_mask(
-            codomain.zerocenter,
-            local_shape,
-            local_offset_Q)
+                                                        codomain.zerocenter,
+                                                        shape,
+                                                        local_offset_Q)
 
         # If both domain and codomain are zero-centered the result,
         # will get a global minus. Store the sign to correct it.
@@ -493,10 +493,11 @@ class FFTWTransformInfo(object):
 
 
 class FFTWLocalTransformInfo(FFTWTransformInfo):
-    def __init__(self, domain, codomain, local_shape,
+    def __init__(self, domain, codomain, axes, local_shape,
                  local_offset_Q, fftw_context, **kwargs):
         super(FFTWLocalTransformInfo, self).__init__(domain,
                                                      codomain,
+                                                     axes,
                                                      local_shape,
                                                      local_offset_Q,
                                                      fftw_context,
@@ -512,10 +513,11 @@ class FFTWLocalTransformInfo(FFTWTransformInfo):
 
 
 class FFTWMPITransfromInfo(FFTWTransformInfo):
-    def __init__(self, domain, codomain, local_shape,
+    def __init__(self, domain, codomain, axes, local_shape,
                  local_offset_Q, fftw_context, transform_shape, **kwargs):
         super(FFTWMPITransfromInfo, self).__init__(domain,
                                                    codomain,
+                                                   axes,
                                                    local_shape,
                                                    local_offset_Q,
                                                    fftw_context,

nifty/operators/fft_operator/transformations/rgrgtransformation.py

@@ -7,29 +7,29 @@ from nifty import RGSpace, nifty_configuration
 
 class RGRGTransformation(Transformation):
     def __init__(self, domain, codomain=None, module=None):
-        super(RGRGTransformation, self).__init__(domain, codomain, module)
+        super(RGRGTransformation, self).__init__(domain, codomain)
 
         if module is None:
             if nifty_configuration['fft_module'] == 'pyfftw':
-                self._transform = FFTW(domain, codomain)
+                self._transform = FFTW(self.domain, self.codomain)
             elif (nifty_configuration['fft_module'] == 'gfft' or
                   nifty_configuration['fft_module'] == 'gfft_dummy'):
                 self._transform = \
-                    GFFT(domain,
-                         codomain,
+                    GFFT(self.domain,
+                         self.codomain,
                          gdi.get(nifty_configuration['fft_module']))
             else:
                 raise ValueError('ERROR: unknow default FFT module:' +
                                  nifty_configuration['fft_module'])
         else:
             if module == 'pyfftw':
-                self._transform = FFTW(domain, codomain)
+                self._transform = FFTW(self.domain, self.codomain)
             elif module == 'gfft':
                 self._transform = \
-                    GFFT(domain, codomain, gdi.get('gfft'))
+                    GFFT(self.domain, self.codomain, gdi.get('gfft'))
             elif module == 'gfft_dummy':
                 self._transform = \
-                    GFFT(domain, codomain, gdi.get('gfft_dummy'))
+                    GFFT(self.domain, self.codomain, gdi.get('gfft_dummy'))
             else:
                 raise ValueError('ERROR: unknow FFT module:' + module)
 

nifty/operators/fft_operator/transformations/transformation.py

@@ -11,7 +11,7 @@ class Transformation(object, Loggable):
     """
     __metaclass__ = abc.ABCMeta
 
-    def __init__(self, domain, codomain=None, module=None):
+    def __init__(self, domain, codomain):
         if codomain is None:
             self.domain = domain
             self.codomain = self.get_codomain(domain)

nifty/operators/smoothing_operator/smooth_util.pyx

+#cython: nonecheck=False
+#cython: boundscheck=True
+#cython: wraparound=False
+#cython: cdivision=True
+
+import numpy as np
+cimport numpy as np
+#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+
+cdef int SIGMA_RANGE = 3
+
+cpdef np.float32_t gaussianKernel_f(np.ndarray[np.float32_t, ndim=1] mpower,
+                                    np.ndarray[np.float32_t, ndim=1] mk,
+                                    np.float32_t mu,
+                                    np.float32_t smooth_length):
+    cdef np.ndarray[np.float32_t, ndim=1] C = \
+        np.exp(-(mk - mu) ** 2 / (2. * smooth_length ** 2))
+    return np.sum(C * mpower) / np.sum(C)
+
+cpdef np.float64_t gaussianKernel(np.ndarray[np.float64_t, ndim=1] mpower,
+                                  np.ndarray[np.float64_t, ndim=1] mk,
+                                  np.float64_t mu,
+                                  np.float64_t smooth_length):
+    cdef np.ndarray[np.float64_t, ndim=1] C = \
+        np.exp(-(mk - mu) ** 2 / (2. * smooth_length ** 2))
+    return np.sum(C * mpower) / np.sum(C)
+
+
+cpdef np.ndarray[np.float64_t, ndim=1] apply_kernel_along_array(
+                                    np.ndarray[np.float64_t, ndim=1] power,
+                                    np.int_t startindex,
+                                    np.int_t endindex,
+                                    np.ndarray[np.float64_t, ndim=1] distances,
+                                    np.float64_t smooth_length,
+                                    np.float64_t smoothing_width):
+
+    if smooth_length == 0.0:
+        return power[startindex:endindex]
+
+    if (smooth_length is None) or (smooth_length < 0):
+        smooth_length = distances[1]-distances[0]
+
+    cdef np.ndarray[np.float64_t, ndim=1] p_smooth = \
+        np.empty(endindex-startindex, dtype=np.float64)
+    cdef np.int64_t i, l, u
+    for i in xrange(startindex, endindex):
+        current_location = distances[i]
+        l = np.searchsorted(distances,
+                            current_location - smoothing_width*smooth_length)
+        u = np.searchsorted(distances,
+                            current_location + smoothing_width*smooth_length)
+        p_smooth[i-startindex] = gaussianKernel(power[l:u],
+                                                distances[l:u],
+                                                distances[i],
+                                                smooth_length)
+
+    return p_smooth
+
+cpdef np.ndarray[np.float32_t, ndim=1] apply_kernel_along_array_f(
+                                    np.ndarray[np.float32_t, ndim=1] power,
+                                    np.int_t startindex,
+                                    np.int_t endindex,
+                                    np.ndarray[np.float32_t, ndim=1] distances,
+                                    np.float64_t smooth_length,
+                                    np.float64_t smoothing_width):
+
+    if smooth_length == 0.0:
+        return power[startindex:endindex]
+
+    if (smooth_length is None) or (smooth_length < 0):
+        smooth_length = distances[1]-distances[0]
+
+    cdef np.ndarray[np.float32_t, ndim=1] p_smooth = \
+        np.empty(endindex-startindex, dtype=np.float64)
+    cdef np.int64_t i, l, u
+    for i in xrange(startindex, endindex):
+        current_location = distances[i]
+        l = np.searchsorted(distances,
+                            current_location - smoothing_width*smooth_length)
+        u = np.searchsorted(distances,
+                            current_location + smoothing_width*smooth_length)
+        p_smooth[i-startindex] = gaussianKernel(power[l:u],
+                                                distances[l:u],
+                                                distances[i],
+                                                smooth_length)
+
+    return p_smooth
+
+def getShape(a):
+    return tuple(a.shape)
+
+cpdef np.ndarray[np.float64_t, ndim=1] apply_along_axis(
+                                    np.int_t axis,
+                                    np.ndarray arr,
+                                    np.int_t startindex,
+                                    np.int_t endindex,
+                                    np.ndarray[np.float64_t, ndim=1] distances,
+                                    np.float64_t smooth_length,
+                                    np.float64_t smoothing_width):
+    cdef np.int_t nd = arr.ndim
+    if axis < 0:
+        axis += nd
+    if (axis >= nd):
+        raise ValueError("axis must be less than arr.ndim; axis=%d, rank=%d."
+            % (axis, nd))
+    cdef np.ndarray[np.int_t, ndim=1] ind = np.zeros(nd-1, dtype=np.int)
+    cdef np.ndarray i = np.zeros(nd, dtype=object)
+    cdef tuple shape = getShape(arr)
+    cdef np.ndarray[np.int_t, ndim=1] indlist = np.asarray(range(nd))
+    indlist = np.delete(indlist,axis)
+    i[axis] = slice(None, None)
+    cdef np.ndarray[np.int_t, ndim=1] outshape = \
+        np.asarray(shape).take(indlist)
+
+    i.put(indlist, ind)
+
+    cdef np.ndarray[np.float64_t, ndim=1] slicedArr
+    cdef np.ndarray[np.float64_t, ndim=1] res
+
+    cdef np.int_t Ntot = np.product(outshape)
+    cdef np.ndarray[np.int_t, ndim=1] holdshape = outshape
+    slicedArr = arr[tuple(i.tolist())]
+
+    res = apply_kernel_along_array(slicedArr,
+                                   startindex,
+                                   endindex,
+                                   distances,
+                                   smooth_length,
+                                   smoothing_width)
+
+    outshape = np.asarray(getShape(arr))
+    outshape[axis] = endindex - startindex
+    outarr = np.zeros(outshape, dtype=np.float64)
+    outarr[tuple(i.tolist())] = res
+    cdef np.int_t n, k = 1
+    while k < Ntot:
+        # increment the index
+        ind[nd-1] += 1
+        n = -1
+        while (ind[n] >= holdshape[n]) and (n > (1-nd)):
+            ind[n-1] += 1
+            ind[n] = 0
+            n -= 1
+        i.put(indlist, ind)
+        slicedArr = arr[tuple(i.tolist())]
+        res = apply_kernel_along_array(slicedArr,
+                                       startindex,
+                                       endindex,
+                                       distances,
+                                       smooth_length,
+                                       smoothing_width)
+        outarr[tuple(i.tolist())] = res
+        k += 1
+
+    return outarr
+
+
+cpdef np.ndarray[np.float32_t, ndim=1] apply_along_axis_f(
+                                    np.int_t axis,
+                                    np.ndarray arr,
+                                    np.int_t startindex,
+                                    np.int_t endindex,
+                                    np.ndarray[np.float32_t, ndim=1] distances,
+                                    np.float64_t smooth_length,
+                                    np.float64_t smoothing_width):
+    cdef np.int_t nd = arr.ndim
+    if axis < 0:
+        axis += nd
+    if (axis >= nd):
+        raise ValueError("axis must be less than arr.ndim; axis=%d, rank=%d."
+            % (axis, nd))
+    cdef np.ndarray[np.int_t, ndim=1] ind = np.zeros(nd-1, dtype=np.int)
+    cdef np.ndarray i = np.zeros(nd, dtype=object)
+    cdef tuple shape = getShape(arr)
+    cdef np.ndarray[np.int_t, ndim=1] indlist = np.asarray(range(nd))
+    indlist = np.delete(indlist,axis)
+    i[axis] = slice(None, None)
+    cdef np.ndarray[np.int_t, ndim=1] outshape = \
+        np.asarray(shape).take(indlist)
+
+    i.put(indlist, ind)
+
+    cdef np.ndarray[np.float32_t, ndim=1] slicedArr
+    cdef np.ndarray[np.float32_t, ndim=1] res
+
+    cdef np.int_t Ntot = np.product(outshape)
+    cdef np.ndarray[np.int_t, ndim=1] holdshape = outshape
+    slicedArr = arr[tuple(i.tolist())]
+
+    res = apply_kernel_along_array_f(slicedArr,
+                                     startindex,
+                                     endindex,
+                                     distances,
+                                     smooth_length,
+                                     smoothing_width)
+
+    outshape = np.asarray(getShape(arr))
+    outshape[axis] = endindex - startindex
+    outarr = np.zeros(outshape, dtype=np.float32)
+    outarr[tuple(i.tolist())] = res
+    cdef np.int_t n, k = 1
+    while k < Ntot:
+        # increment the index
+        ind[nd-1] += 1
+        n = -1
+        while (ind[n] >= holdshape[n]) and (n > (1-nd)):
+            ind[n-1] += 1
+            ind[n] = 0
+            n -= 1
+        i.put(indlist, ind)
+        slicedArr = arr[tuple(i.tolist())]
+        res = apply_kernel_along_array(slicedArr,
+                                       startindex,
+                                       endindex,
+                                       distances,
+                                       smooth_length,
+                                       smoothing_width)
+        outarr[tuple(i.tolist())] = res
+        k += 1
+
+    return outarr
+

nifty/operators/smoothing_operator/smoothing_operator.py

 import numpy as np
 
 import nifty.nifty_utilities as utilities
-from nifty import RGSpace, LMSpace
 from nifty.operators.endomorphic_operator import EndomorphicOperator
 from nifty.operators.fft_operator import FFTOperator
+import smooth_util as su
+from d2o import STRATEGIES
 
 
 class SmoothingOperator(EndomorphicOperator):
-
     # ---Overwritten properties and methods---
-    def __init__(self, domain=(), field_type=(), sigma=0):
+    def __init__(self, domain=(), field_type=(), sigma=0,
+                 log_distances=False):
 
         self._domain = self._parse_domain(domain)
         self._field_type = self._parse_field_type(field_type)
@@ -17,8 +18,8 @@ class SmoothingOperator(EndomorphicOperator):
         if len(self.domain) != 1:
             raise ValueError(
 
-                    'ERROR: SmoothOperator accepts only exactly one '
-                    'space as input domain.'
+                'ERROR: SmoothOperator accepts only exactly one '
+                'space as input domain.'
             )
 
         if self.field_type != ():
@@ -27,13 +28,16 @@ class SmoothingOperator(EndomorphicOperator):
                 'empty tuple.'
             )
 
-        self._sigma = sigma
+        self.sigma = sigma
+        self.log_distances = log_distances
+
+        self._direct_smoothing_width = 3.
 
     def _inverse_times(self, x, spaces, types):
-        return self._smooth_helper(x, spaces, types, inverse=True)
+        return self._smoothing_helper(x, spaces, inverse=True)
 
     def _times(self, x, spaces, types):
-        return self._smooth_helper(x, spaces, types)
+        return self._smoothing_helper(x, spaces, inverse=False)
 
     # ---Mandatory properties and methods---
     @property
@@ -50,18 +54,31 @@ class SmoothingOperator(EndomorphicOperator):
 
     @property
     def symmetric(self):
-        return True
+        return False
 
     @property
     def unitary(self):
         return False
 
     # ---Added properties and methods---
+
     @property
     def sigma(self):
         return self._sigma
 
-    def _smooth_helper(self, x, spaces, types, inverse=False):
+    @sigma.setter
+    def sigma(self, sigma):
+        self._sigma = np.float(sigma)
+
+    @property
+    def log_distances(self):
+        return self._log_distances
+
+    @log_distances.setter
+    def log_distances(self, log_distances):
+        self._log_distances = bool(log_distances)
+
+    def _smoothing_helper(self, x, spaces, inverse):
         if self.sigma == 0:
             return x.copy()
 
@@ -74,6 +91,13 @@ class SmoothingOperator(EndomorphicOperator):
         else:
             spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
 
+        try:
+            result = self._fft_smoothing(x, spaces, inverse)
+        except ValueError:
+            result = self._direct_smoothing(x, spaces, inverse)
+        return result
+
+    def _fft_smoothing(self, x, spaces, inverse):
         Transformator = FFTOperator(x.domain[spaces[0]])
 
         # transform to the (global-)default codomain and perform all remaining
@@ -87,9 +111,13 @@ class SmoothingOperator(EndomorphicOperator):
             transformed_x.val.get_axes_local_distribution_strategy(axes=coaxes)
 
         kernel = codomain.get_distance_array(
-                        distribution_strategy=axes_local_distribution_strategy)
+            distribution_strategy=axes_local_distribution_strategy)
+
+        if self.log_distances:
+            kernel.apply_scalar_function(np.log, inplace=True)
+
         kernel.apply_scalar_function(
-            codomain.get_smoothing_kernel_function(self.sigma),
+            codomain.get_fft_smoothing_kernel_function(self.sigma),
             inplace=True)
 
         # now, apply the kernel to transformed_x
@@ -98,7 +126,7 @@ class SmoothingOperator(EndomorphicOperator):
         local_transformed_x = transformed_x.val.get_local_data(copy=False)
         local_kernel = kernel.get_local_data(copy=False)
 
-        reshaper = [local_transformed_x.shape[i] if i in coaxes else 1
+        reshaper = [transformed_x.shape[i] if i in coaxes else 1
                     for i in xrange(len(transformed_x.shape))]
         local_kernel = np.reshape(local_kernel, reshaper)
 
@@ -113,3 +141,127 @@ class SmoothingOperator(EndomorphicOperator):
         result = Transformator.inverse_times(transformed_x, spaces=spaces)
 
         return result
+
+    def _direct_smoothing(self, x, spaces, inverse):
+        # infer affected axes
+        # we rely on the knowledge, that `spaces` is a tuple with length 1.
+        affected_axes = x.domain_axes[spaces[0]]
+
+        if len(affected_axes) > 1:
+            raise ValueError("By this implementation only one-dimensional "
+                             "spaces can be smoothed directly.")
+
+        affected_axis = affected_axes[0]
+
+        distance_array = x.domain[spaces[0]].get_distance_array(
+            distribution_strategy='not')
+        distance_array = distance_array.get_local_data(copy=False)
+
+        if self.log_distances:
+            np.log(distance_array, out=distance_array)
+
+        # collect the local data + ghost cells