Commit 1c53be26 authored by theos's avatar theos

Connected the Cython code for explicit smoothing with the SmoothingOperator.

Fixed several bugs.
parent ad0a6d12
#cython: nonecheck=False
#cython: boundscheck=True
#cython: wraparound=False
#cython: cdivision=True
import numpy as np import numpy as np
cimport numpy as np cimport numpy as np
#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #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 int SIGMA_RANGE = 3
cdef np.ndarray[np.float32_t, ndim=1] C = np.exp(-(mk - mu) ** 2 / (2. * smooth_length ** 2))
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) 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): cpdef np.float64_t gaussianKernel(np.ndarray[np.float64_t, ndim=1] mpower,
cdef np.ndarray[np.float64_t, ndim=1] C = np.exp(-(mk - mu) ** 2 / (2. * smooth_length ** 2)) 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) 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: if smooth_length == 0.0:
return power[startindex:endindex] return power[startindex:endindex]
...@@ -19,16 +40,29 @@ cpdef np.ndarray[np.float64_t, ndim=1] apply_kernel_along_array(np.ndarray[np.fl ...@@ -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): if (smooth_length is None) or (smooth_length < 0):
smooth_length = distances[1]-distances[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.ndarray[np.float64_t, ndim=1] p_smooth = \
cdef np.int i np.empty(endindex-startindex, dtype=np.float64)
cdef np.int64_t i, l, u
for i in xrange(startindex, endindex): for i in xrange(startindex, endindex):
l = max(i-int(2*smooth_length)-1,0) current_location = distances[i]
u = min(i+int(2*smooth_length)+2,len(p_smooth)) l = np.searchsorted(distances,
p_smooth[i-startindex] = gaussianKernel(power[l:u], distances[l:u], distances[i], smooth_length) 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 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: if smooth_length == 0.0:
return power[startindex:endindex] return power[startindex:endindex]
...@@ -36,20 +70,33 @@ cpdef np.ndarray[np.float32_t, ndim=1] apply_kernel_along_array_f(np.ndarray[np. ...@@ -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): if (smooth_length is None) or (smooth_length < 0):
smooth_length = distances[1]-distances[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.ndarray[np.float32_t, ndim=1] p_smooth = \
cdef np.int i np.empty(endindex-startindex, dtype=np.float64)
cdef np.int64_t i, l, u
for i in xrange(startindex, endindex): for i in xrange(startindex, endindex):
l = max(i-int(2*smooth_length)-1,0) current_location = distances[i]
u = min(i+int(2*smooth_length)+2,len(p_smooth)) l = np.searchsorted(distances,
p_smooth[i-startindex] = gaussianKernel_f(power[l:u], distances[l:u], distances[i], smooth_length) 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 return p_smooth
def getShape(a): def getShape(a):
return tuple(a.shape) 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 cdef np.int_t nd = arr.ndim
if axis < 0: if axis < 0:
axis += nd axis += nd
...@@ -62,7 +109,8 @@ cpdef np.ndarray[np.float64_t, ndim=1] apply_along_axis(np.int_t axis, np.ndarra ...@@ -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)) cdef np.ndarray[np.int_t, ndim=1] indlist = np.asarray(range(nd))
indlist = np.delete(indlist,axis) indlist = np.delete(indlist,axis)
i[axis] = slice(None, None) 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) 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 ...@@ -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 cdef np.ndarray[np.int_t, ndim=1] holdshape = outshape
slicedArr = arr[tuple(i.tolist())] 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 = np.asarray(getShape(arr))
outshape[axis] = endindex - startindex 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 ...@@ -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 cdef np.int_t n, k = 1
while k < Ntot: while k < Ntot:
# increment the index # increment the index
ind[-1] += 1 ind[nd-1] += 1
n = -1 n = -1
while (ind[n] >= holdshape[n]) and (n > (1-nd)): while (ind[n] >= holdshape[n]) and (n > (1-nd)):
ind[n-1] += 1 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 ...@@ -90,14 +143,26 @@ cpdef np.ndarray[np.float64_t, ndim=1] apply_along_axis(np.int_t axis, np.ndarra
n -= 1 n -= 1
i.put(indlist, ind) i.put(indlist, ind)
slicedArr = arr[tuple(i.tolist())] 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 outarr[tuple(i.tolist())] = res
k += 1 k += 1
return outarr 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 cdef np.int_t nd = arr.ndim
if axis < 0: if axis < 0:
axis += nd axis += nd
...@@ -110,7 +175,8 @@ cpdef np.ndarray[np.float32_t, ndim=1] apply_along_axis_f(np.int_t axis, np.ndar ...@@ -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)) cdef np.ndarray[np.int_t, ndim=1] indlist = np.asarray(range(nd))
indlist = np.delete(indlist,axis) indlist = np.delete(indlist,axis)
i[axis] = slice(None, None) 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) 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 ...@@ -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 cdef np.ndarray[np.int_t, ndim=1] holdshape = outshape
slicedArr = arr[tuple(i.tolist())] 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 = np.asarray(getShape(arr))
outshape[axis] = endindex - startindex outshape[axis] = endindex - startindex
outarr = np.zeros(outshape, dtype=np.float32_t) outarr = np.zeros(outshape, dtype=np.float32)
outarr[tuple(i.tolist())] = res outarr[tuple(i.tolist())] = res
cdef np.int_t n, k = 1 cdef np.int_t n, k = 1
while k < Ntot: while k < Ntot:
# increment the index # increment the index
ind[-1] += 1 ind[nd-1] += 1
n = -1 n = -1
while (ind[n] >= holdshape[n]) and (n > (1-nd)): while (ind[n] >= holdshape[n]) and (n > (1-nd)):
ind[n-1] += 1 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 ...@@ -138,7 +209,12 @@ cpdef np.ndarray[np.float32_t, ndim=1] apply_along_axis_f(np.int_t axis, np.ndar
n -= 1 n -= 1
i.put(indlist, ind) i.put(indlist, ind)
slicedArr = arr[tuple(i.tolist())] 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 outarr[tuple(i.tolist())] = res
k += 1 k += 1
......
...@@ -31,7 +31,7 @@ class SmoothingOperator(EndomorphicOperator): ...@@ -31,7 +31,7 @@ class SmoothingOperator(EndomorphicOperator):
self.sigma = sigma self.sigma = sigma
self.log_distances = log_distances self.log_distances = log_distances
self._direct_smoothing_width = 2. self._direct_smoothing_width = 3.
def _inverse_times(self, x, spaces, types): def _inverse_times(self, x, spaces, types):
return self._smoothing_helper(x, spaces, inverse=True) return self._smoothing_helper(x, spaces, inverse=True)
...@@ -147,14 +147,18 @@ class SmoothingOperator(EndomorphicOperator): ...@@ -147,14 +147,18 @@ class SmoothingOperator(EndomorphicOperator):
# we rely on the knowledge, that `spaces` is a tuple with length 1. # we rely on the knowledge, that `spaces` is a tuple with length 1.
affected_axes = x.domain_axes[spaces[0]] affected_axes = x.domain_axes[spaces[0]]
axes_local_distribution_strategy = \ if len(affected_axes) > 1:
x.val.get_axes_local_distribution_strategy(axes=affected_axes) 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( 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: 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 # collect the local data + ghost cells
local_data_Q = False local_data_Q = False
...@@ -164,88 +168,100 @@ class SmoothingOperator(EndomorphicOperator): ...@@ -164,88 +168,100 @@ class SmoothingOperator(EndomorphicOperator):
elif x.distribution_strategy in STRATEGIES['slicing']: elif x.distribution_strategy in STRATEGIES['slicing']:
# infer the local start/end based on the slicing information of # infer the local start/end based on the slicing information of
# x's d2o. Only gets non-trivial for axis==0. # x's d2o. Only gets non-trivial for axis==0.
if 0 not in affected_axes: if 0 != affected_axis:
local_data_Q = True local_data_Q = True
else: else:
# we rely on the fact, that the content of x.domain_axes is start_index = x.val.distributor.local_start
# sorted start_distance = distance_array[start_index]
true_starts = [x.val.distributor.local_start] augmented_start_distance = \
true_starts += [0] * (len(affected_axes) - 1) (start_distance - self._direct_smoothing_width*self.sigma)
true_ends = [x.val.distributor.local_end] augmented_start_index = \
true_ends += [x.shape[i] for i in affected_axes[1:]] np.searchsorted(distance_array, augmented_start_distance)
true_start = start_index - augmented_start_index
augmented_start = max(0, end_index = x.val.distributor.local_end
true_starts[0] - end_distance = distance_array[end_index-1]
self._direct_smoothing_width * self.sigma) augmented_end_distance = \
augmented_end = min(x.shape[affected_axes[0]], (end_distance + self._direct_smoothing_width*self.sigma)
true_ends[0] + augmented_end_index = \
self._direct_smoothing_width * self.sigma) np.searchsorted(distance_array, augmented_end_distance)
augmented_slice = slice(augmented_start, augmented_end) 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, augmented_data = x.val.get_data(augmented_slice,
local_keys=True, local_keys=True,
copy=False) copy=False)
augmented_data = augmented_data.get_local_data(copy=False) augmented_data = augmented_data.get_local_data(copy=False)
augmented_distance_array = distance_array.get_data( augmented_distance_array = distance_array[augmented_slice]
augmented_slice,
local_keys=True,
copy=False)
augmented_distance_array = \
augmented_distance_array.get_local_data(copy=False)
else: else:
raise ValueError(about._errors.cstring( raise ValueError("Direct smoothing not implemented for given"
"ERROR: Direct smoothing not implemented for given" "distribution strategy.")
"distribution strategy."))
if local_data_Q: if local_data_Q:
# if the needed data resides on the nodes already, the necessary # if the needed data resides on the nodes already, the necessary
# are the same; no matter what the distribution strategy was. # are the same; no matter what the distribution strategy was.
augmented_data = x.val.get_local_data(copy=False) augmented_data = x.val.get_local_data(copy=False)
augmented_distance_array = \ augmented_distance_array = distance_array
distance_array.get_local_data(copy=False) true_start = 0
true_starts = [0] * len(affected_axes) true_end = x.shape[affected_axis]
true_ends = [x.shape[i] for i in affected_axes]
# perform the convolution along the affected axes # perform the convolution along the affected axes
local_result = augmented_data # currently only one axis is supported
for index in range(len(affected_axes)): data_axis = affected_axes[0]
data_axis = affected_axes[index] local_result = self._direct_smoothing_single_axis(
distances_axis = index augmented_data,
true_start = true_starts[index] data_axis,
true_end = true_ends[index] augmented_distance_array,
true_start,
local_result = self._direct_smoothing_single_axis( true_end,
local_result, inverse)
data_axis,
augmented_distance_array,
distances_axis,
true_start,
true_end,
inverse)
result = x.copy_empty() result = x.copy_empty()
result.val.set_local_data(local_result, copy=False) result.val.set_local_data(local_result, copy=False)
return result return result
def _direct_smoothing_single_axis(self, data, data_axis, distances, def _direct_smoothing_single_axis(self, data, data_axis, distances,
distances_axis, true_start, true_end, true_start, true_end, inverse):
inverse):
if inverse: if inverse:
true_sigma = 1 / self.sigma true_sigma = 1. / self.sigma
else: else:
true_sigma = self.sigma true_sigma = self.sigma
if (data.dtype == np.dtype('float32')): if data.dtype is np.dtype('float32'):
smoothed_data = su.apply_along_axis_f(data_axis, data, distances = distances.astype(np.float32, copy=False)
startindex=true_start, smoothed_data = su.apply_along_axis_f(
endindex=true_end, data_axis, data,
distances=distances, startindex=true_start,
smooth_length=true_sigma) 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: else:
smoothed_data = su.apply_along_axis(data_axis, data, raise TypeError("Dtype %s not supported" % str(data.dtype))
startindex=true_start,
endindex=true_end,
distances=distances,
smooth_length=true_sigma)
return smoothed_data return smoothed_data
...@@ -3,7 +3,8 @@ from __future__ import division ...@@ -3,7 +3,8 @@ from __future__ import division
import itertools import itertools
import numpy as np 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.spaces.space import Space
from nifty.config import nifty_configuration as gc,\ from nifty.config import nifty_configuration as gc,\
...@@ -154,8 +155,8 @@ class GLSpace(Space): ...@@ -154,8 +155,8 @@ class GLSpace(Space):
return result_x return result_x
def get_distance_array(self, distribution_strategy): def get_distance_array(self, distribution_strategy):
dists = arange(start=0, stop=self.shape[0], dists = d2o.arange(start=0, stop=self.shape[0],
distribution_strategy=distribution_strategy) distribution_strategy=distribution_strategy)
dists = dists.apply_scalar_function( dists = dists.apply_scalar_function(
lambda x: self._distance_array_helper(divmod(x, self.nlon)), lambda x: self._distance_array_helper(divmod(x, self.nlon)),
...@@ -172,7 +173,7 @@ class GLSpace(Space): ...@@ -172,7 +173,7 @@ class GLSpace(Space):
return np.arctan(numerator / denominator) return np.arctan(numerator / denominator)
def get_smoothing_kernel_function(self, sigma): def get_fft_smoothing_kernel_function(self, sigma):
if sigma is None: if sigma is None:
sigma = np.sqrt(2) * np.pi sigma = np.sqrt(2) * np.pi
......
...@@ -35,6 +35,8 @@ from __future__ import division ...@@ -35,6 +35,8 @@ from __future__ import division
import numpy as np import numpy as np
import d2o
from nifty.spaces.space import Space from nifty.spaces.space import Space
from nifty.config import nifty_configuration as gc, \ from nifty.config import nifty_configuration as gc, \
dependency_injector as gdi dependency_injector as gdi
...@@ -169,7 +171,7 @@ class HPSpace(Space): ...@@ -169,7 +171,7 @@ class HPSpace(Space):
------- -------
dists: distributed_data_object dists: distributed_data_object
""" """
dists = arange( dists = d2o.arange(
start=0, stop=self.shape[0], start=0, stop=self.shape[0],
distribution_strategy=distribution_strategy distribution_strategy=distribution_strategy
) )
...@@ -183,7 +185,7 @@ class HPSpace(Space): ...@@ -183,7 +185,7 @@ class HPSpace(Space):
return dists return dists
def get_smoothing_kernel_function(self, sigma): def get_fft_smoothing_kernel_function(self, sigma):
if sigma is None: if sigma is None:
sigma = np.sqrt(2) * np.pi sigma = np.sqrt(2) * np.pi
......