diff --git a/nifty/operators/smoothing_operator/smooth_util.pyx b/nifty/operators/smoothing_operator/smooth_util.pyx index db1676543516362c864f911b0497769747de42c2..f1e714562c9c907774315bb2ee87cd39d93b1d18 100644 --- a/nifty/operators/smoothing_operator/smooth_util.pyx +++ b/nifty/operators/smoothing_operator/smooth_util.pyx @@ -1,17 +1,38 @@ +#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 -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)) +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)) +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): +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] @@ -19,16 +40,29 @@ cpdef np.ndarray[np.float64_t, ndim=1] apply_kernel_along_array(np.ndarray[np.fl 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.int i + 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): - l = max(i-int(2*smooth_length)-1,0) - u = min(i+int(2*smooth_length)+2,len(p_smooth)) - p_smooth[i-startindex] = gaussianKernel(power[l:u], distances[l:u], distances[i], smooth_length) + 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.float32_t smooth_length): +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] @@ -36,20 +70,33 @@ cpdef np.ndarray[np.float32_t, ndim=1] apply_kernel_along_array_f(np.ndarray[np. 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.float32_t) - cdef np.int i + 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): - l = max(i-int(2*smooth_length)-1,0) - u = min(i+int(2*smooth_length)+2,len(p_smooth)) - p_smooth[i-startindex] = gaussianKernel_f(power[l:u], distances[l:u], distances[i], smooth_length) + 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): +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 @@ -62,7 +109,8 @@ cpdef np.ndarray[np.float64_t, ndim=1] apply_along_axis(np.int_t axis, np.ndarra 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) + cdef np.ndarray[np.int_t, ndim=1] outshape = \ + np.asarray(shape).take(indlist) i.put(indlist, ind) @@ -73,7 +121,12 @@ cpdef np.ndarray[np.float64_t, ndim=1] apply_along_axis(np.int_t axis, np.ndarra 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) + res = apply_kernel_along_array(slicedArr, + startindex, + endindex, + distances, + smooth_length, + smoothing_width) outshape = np.asarray(getShape(arr)) outshape[axis] = endindex - startindex @@ -82,7 +135,7 @@ cpdef np.ndarray[np.float64_t, ndim=1] apply_along_axis(np.int_t axis, np.ndarra cdef np.int_t n, k = 1 while k < Ntot: # increment the index - ind[-1] += 1 + ind[nd-1] += 1 n = -1 while (ind[n] >= holdshape[n]) and (n > (1-nd)): ind[n-1] += 1 @@ -90,14 +143,26 @@ cpdef np.ndarray[np.float64_t, ndim=1] apply_along_axis(np.int_t axis, np.ndarra n -= 1 i.put(indlist, ind) slicedArr = arr[tuple(i.tolist())] - res = apply_kernel_along_array(slicedArr, startindex, endindex, distances, smooth_length) + 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.float32_t smooth_length): +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 @@ -110,7 +175,8 @@ cpdef np.ndarray[np.float32_t, ndim=1] apply_along_axis_f(np.int_t axis, np.ndar 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) + cdef np.ndarray[np.int_t, ndim=1] outshape = \ + np.asarray(shape).take(indlist) i.put(indlist, ind) @@ -121,16 +187,21 @@ cpdef np.ndarray[np.float32_t, ndim=1] apply_along_axis_f(np.int_t axis, np.ndar 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) + 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_t) + 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[-1] += 1 + ind[nd-1] += 1 n = -1 while (ind[n] >= holdshape[n]) and (n > (1-nd)): ind[n-1] += 1 @@ -138,7 +209,12 @@ cpdef np.ndarray[np.float32_t, ndim=1] apply_along_axis_f(np.int_t axis, np.ndar n -= 1 i.put(indlist, ind) slicedArr = arr[tuple(i.tolist())] - res = apply_kernel_along_array(slicedArr, startindex, endindex, distances, smooth_length) + res = apply_kernel_along_array(slicedArr, + startindex, + endindex, + distances, + smooth_length, + smoothing_width) outarr[tuple(i.tolist())] = res k += 1 diff --git a/nifty/operators/smoothing_operator/smoothing_operator.py b/nifty/operators/smoothing_operator/smoothing_operator.py index b7f3451cd730f5d923c1f812506b186c52d0f8f6..f8cac64e28491d94b03def89b2c78fc77f87f5e4 100644 --- a/nifty/operators/smoothing_operator/smoothing_operator.py +++ b/nifty/operators/smoothing_operator/smoothing_operator.py @@ -31,7 +31,7 @@ class SmoothingOperator(EndomorphicOperator): self.sigma = sigma self.log_distances = log_distances - self._direct_smoothing_width = 2. + self._direct_smoothing_width = 3. def _inverse_times(self, x, spaces, types): return self._smoothing_helper(x, spaces, inverse=True) @@ -147,14 +147,18 @@ class SmoothingOperator(EndomorphicOperator): # we rely on the knowledge, that `spaces` is a tuple with length 1. affected_axes = x.domain_axes[spaces[0]] - axes_local_distribution_strategy = \ - x.val.get_axes_local_distribution_strategy(axes=affected_axes) + 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=axes_local_distribution_strategy) + distribution_strategy='not') + distance_array = distance_array.get_local_data(copy=False) if self.log_distances: - distance_array.apply_scalar_function(np.log, inplace=True) + np.log(distance_array, out=distance_array) # collect the local data + ghost cells local_data_Q = False @@ -164,88 +168,100 @@ class SmoothingOperator(EndomorphicOperator): elif x.distribution_strategy in STRATEGIES['slicing']: # infer the local start/end based on the slicing information of # x's d2o. Only gets non-trivial for axis==0. - if 0 not in affected_axes: + if 0 != affected_axis: local_data_Q = True else: - # we rely on the fact, that the content of x.domain_axes is - # sorted - true_starts = [x.val.distributor.local_start] - true_starts += [0] * (len(affected_axes) - 1) - true_ends = [x.val.distributor.local_end] - true_ends += [x.shape[i] for i in affected_axes[1:]] - - augmented_start = max(0, - true_starts[0] - - self._direct_smoothing_width * self.sigma) - augmented_end = min(x.shape[affected_axes[0]], - true_ends[0] + - self._direct_smoothing_width * self.sigma) - augmented_slice = slice(augmented_start, augmented_end) + start_index = x.val.distributor.local_start + start_distance = distance_array[start_index] + augmented_start_distance = \ + (start_distance - self._direct_smoothing_width*self.sigma) + augmented_start_index = \ + np.searchsorted(distance_array, augmented_start_distance) + true_start = start_index - augmented_start_index + end_index = x.val.distributor.local_end + end_distance = distance_array[end_index-1] + augmented_end_distance = \ + (end_distance + self._direct_smoothing_width*self.sigma) + augmented_end_index = \ + np.searchsorted(distance_array, augmented_end_distance) + true_end = true_start + x.val.distributor.local_length + augmented_slice = slice(augmented_start_index, + augmented_end_index) + augmented_data = x.val.get_data(augmented_slice, local_keys=True, copy=False) augmented_data = augmented_data.get_local_data(copy=False) - augmented_distance_array = distance_array.get_data( - augmented_slice, - local_keys=True, - copy=False) - augmented_distance_array = \ - augmented_distance_array.get_local_data(copy=False) + augmented_distance_array = distance_array[augmented_slice] else: - raise ValueError(about._errors.cstring( - "ERROR: Direct smoothing not implemented for given" - "distribution strategy.")) + raise ValueError("Direct smoothing not implemented for given" + "distribution strategy.") if local_data_Q: # if the needed data resides on the nodes already, the necessary # are the same; no matter what the distribution strategy was. augmented_data = x.val.get_local_data(copy=False) - augmented_distance_array = \ - distance_array.get_local_data(copy=False) - true_starts = [0] * len(affected_axes) - true_ends = [x.shape[i] for i in affected_axes] + augmented_distance_array = distance_array + true_start = 0 + true_end = x.shape[affected_axis] # perform the convolution along the affected axes - local_result = augmented_data - for index in range(len(affected_axes)): - data_axis = affected_axes[index] - distances_axis = index - true_start = true_starts[index] - true_end = true_ends[index] - - local_result = self._direct_smoothing_single_axis( - local_result, - data_axis, - augmented_distance_array, - distances_axis, - true_start, - true_end, - inverse) - + # currently only one axis is supported + data_axis = affected_axes[0] + local_result = self._direct_smoothing_single_axis( + augmented_data, + data_axis, + augmented_distance_array, + true_start, + true_end, + inverse) result = x.copy_empty() result.val.set_local_data(local_result, copy=False) return result def _direct_smoothing_single_axis(self, data, data_axis, distances, - distances_axis, true_start, true_end, - inverse): + true_start, true_end, inverse): if inverse: - true_sigma = 1 / self.sigma + true_sigma = 1. / self.sigma else: true_sigma = self.sigma - if (data.dtype == np.dtype('float32')): - smoothed_data = su.apply_along_axis_f(data_axis, data, - startindex=true_start, - endindex=true_end, - distances=distances, - smooth_length=true_sigma) + if data.dtype is np.dtype('float32'): + distances = distances.astype(np.float32, copy=False) + smoothed_data = su.apply_along_axis_f( + data_axis, data, + startindex=true_start, + endindex=true_end, + distances=distances, + smooth_length=true_sigma, + smoothing_width=self._direct_smoothing_width) + elif data.dtype is np.dtype('float64'): + distances = distances.astype(np.float64, copy=False) + smoothed_data = su.apply_along_axis( + data_axis, data, + startindex=true_start, + endindex=true_end, + distances=distances, + smooth_length=true_sigma, + smoothing_width=self._direct_smoothing_width) + + elif np.issubdtype(data.dtype, np.complexfloating): + real = self._direct_smoothing_single_axis(data.real, + data_axis, + distances, + true_start, + true_end, inverse) + imag = self._direct_smoothing_single_axis(data.imag, + data_axis, + distances, + true_start, + true_end, inverse) + + return real + 1j*imag + else: - smoothed_data = su.apply_along_axis(data_axis, data, - startindex=true_start, - endindex=true_end, - distances=distances, - smooth_length=true_sigma) + raise TypeError("Dtype %s not supported" % str(data.dtype)) + return smoothed_data diff --git a/nifty/spaces/gl_space/gl_space.py b/nifty/spaces/gl_space/gl_space.py index 81cfeff6b47c2379b62cf6b74d8ab961a0485b40..381fff54daeb1082ad045eb1da6a55c170279230 100644 --- a/nifty/spaces/gl_space/gl_space.py +++ b/nifty/spaces/gl_space/gl_space.py @@ -3,7 +3,8 @@ from __future__ import division import itertools import numpy as np -from d2o import arange, STRATEGIES as DISTRIBUTION_STRATEGIES +import d2o +from d2o import STRATEGIES as DISTRIBUTION_STRATEGIES from nifty.spaces.space import Space from nifty.config import nifty_configuration as gc,\ @@ -154,8 +155,8 @@ class GLSpace(Space): return result_x def get_distance_array(self, distribution_strategy): - dists = arange(start=0, stop=self.shape[0], - distribution_strategy=distribution_strategy) + dists = d2o.arange(start=0, stop=self.shape[0], + distribution_strategy=distribution_strategy) dists = dists.apply_scalar_function( lambda x: self._distance_array_helper(divmod(x, self.nlon)), @@ -172,7 +173,7 @@ class GLSpace(Space): return np.arctan(numerator / denominator) - def get_smoothing_kernel_function(self, sigma): + def get_fft_smoothing_kernel_function(self, sigma): if sigma is None: sigma = np.sqrt(2) * np.pi diff --git a/nifty/spaces/hp_space/hp_space.py b/nifty/spaces/hp_space/hp_space.py index c8d1fffcca87418ba662bb04d4ad0558eb4449cb..ae74476f29f5abe4d4ba30decbd0ec303d68ad39 100644 --- a/nifty/spaces/hp_space/hp_space.py +++ b/nifty/spaces/hp_space/hp_space.py @@ -35,6 +35,8 @@ from __future__ import division import numpy as np +import d2o + from nifty.spaces.space import Space from nifty.config import nifty_configuration as gc, \ dependency_injector as gdi @@ -169,7 +171,7 @@ class HPSpace(Space): ------- dists: distributed_data_object """ - dists = arange( + dists = d2o.arange( start=0, stop=self.shape[0], distribution_strategy=distribution_strategy ) @@ -183,7 +185,7 @@ class HPSpace(Space): return dists - def get_smoothing_kernel_function(self, sigma): + def get_fft_smoothing_kernel_function(self, sigma): if sigma is None: sigma = np.sqrt(2) * np.pi diff --git a/nifty/spaces/lm_space/lm_space.py b/nifty/spaces/lm_space/lm_space.py index 0a34c939eccd6a895ca68c81cd85eb69385bd9d7..4a33a354e892ad7d8b74c2566c559ad5206909be 100644 --- a/nifty/spaces/lm_space/lm_space.py +++ b/nifty/spaces/lm_space/lm_space.py @@ -156,7 +156,7 @@ class LMSpace(Space): return dists - def get_smoothing_kernel_function(self, sigma): + def get_fft_smoothing_kernel_function(self, sigma): if sigma is None: sigma = np.sqrt(2) * np.pi / (self.lmax + 1) diff --git a/nifty/spaces/power_space/power_space.py b/nifty/spaces/power_space/power_space.py index 54ed22824e48cca65df6affed858f7fdf91313fa..6da672751e918d646e63913a3cff4712a18b5af6 100644 --- a/nifty/spaces/power_space/power_space.py +++ b/nifty/spaces/power_space/power_space.py @@ -2,6 +2,8 @@ import numpy as np +import d2o + from power_index_factory import PowerIndexFactory from nifty.spaces.space import Space @@ -106,13 +108,14 @@ class PowerSpace(Space): return result_x def get_distance_array(self, distribution_strategy): - raise NotImplementedError( - "There is no get_distance_array implementation for " - "PowerSpace.") + result = d2o.distributed_data_object( + self.kindex, + distribution_strategy=distribution_strategy) + return result - def get_smoothing_kernel_function(self, sigma): + def get_fft_smoothing_kernel_function(self, sigma): raise NotImplementedError( - "There is no smoothing function for PowerSpace.") + "There is no fft smoothing function for PowerSpace.") # ---Added properties and methods--- diff --git a/nifty/spaces/rg_space/rg_space.py b/nifty/spaces/rg_space/rg_space.py index 1bd9bfae949965e1ff9e663532f89e4433b0527d..9c99617dfead8e429d80d78ec34b3044a1e64c66 100644 --- a/nifty/spaces/rg_space/rg_space.py +++ b/nifty/spaces/rg_space/rg_space.py @@ -282,7 +282,7 @@ class RGSpace(Space): dists = np.sqrt(dists) return dists - def get_smoothing_kernel_function(self, sigma): + def get_fft_smoothing_kernel_function(self, sigma): if sigma is None: sigma = np.sqrt(2) * np.max(self.distances) diff --git a/nifty/spaces/space/space.py b/nifty/spaces/space/space.py index 802ae2f0c7ba997f670c8d95420edabd8bf528f5..eaa9286c2ec72d5f0b823ce4fcb4f1ce9c63b04d 100644 --- a/nifty/spaces/space/space.py +++ b/nifty/spaces/space/space.py @@ -278,7 +278,7 @@ class Space(object, Loggable): raise NotImplementedError( "There is no generic distance structure for Space base class.") - def get_smoothing_kernel_function(self, sigma): + def get_fft_smoothing_kernel_function(self, sigma): raise NotImplementedError( "There is no generic co-smoothing kernel for Space base class.") diff --git a/setup.py b/setup.py index cb20468e86ffaacbc906c8ebf5875a369d8b308c..83e64ae9b8d2f7e35b8fb03df51d7da2b79fd188 100644 --- a/setup.py +++ b/setup.py @@ -39,8 +39,9 @@ setup(name="ift_nifty", zip_safe=False, ext_modules=cythonize([ "nifty/operators/fft_operator/transformations/" - "lm_transformation_factory.pyx", - "nifty/spaces/lm_space/lm_helper.pyx" + "lm_transformation_factory.pyx", + "nifty/spaces/lm_space/lm_helper.pyx", + "nifty/operators/smoothing_operator/smooth_util.pyx" ]), include_dirs=[numpy.get_include()], dependency_links=[