Commit 6d217b7b authored by Theo Steininger's avatar Theo Steininger
Browse files

Decomposed the SmoothingOperator class.

parent 93067284
Pipeline #12421 passed with stage
in 6 minutes and 2 seconds
......@@ -24,7 +24,7 @@ from diagonal_operator import DiagonalOperator
from endomorphic_operator import EndomorphicOperator
from smoothing_operator import SmoothingOperator
from smoothing_operator import *
from fft_operator import *
......
......@@ -75,7 +75,8 @@ class LinearOperator(Loggable, object):
def __init__(self, default_spaces=None):
self.default_spaces = default_spaces
def _parse_domain(self, domain):
@staticmethod
def _parse_domain(domain):
return utilities.parse_domain(domain)
@abc.abstractproperty
......
......@@ -6,6 +6,7 @@ from nifty.operators.smoothing_operator import SmoothingOperator
from nifty.operators.composed_operator import ComposedOperator
from nifty.operators.diagonal_operator import DiagonalOperator
class ResponseOperator(LinearOperator):
""" NIFTy ResponseOperator (example)
......
......@@ -16,4 +16,4 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
from smoothing_operator import SmoothingOperator
from .smoothing_operator import SmoothingOperator
# -*- coding: utf8 -*-
import numpy as np
from d2o import STRATEGIES
from .smoothing_operator import SmoothingOperator
class DirectSmoothingOperator(SmoothingOperator):
def __init__(self, domain, sigma, log_distances=False,
default_spaces=None):
super(DirectSmoothingOperator, self).__init__(domain,
sigma,
log_distances,
default_spaces)
self.effective_smoothing_width = 3.01
def _precompute(self, x, sigma, dxmax=None):
""" Does precomputations for Gaussian smoothing on a 1D irregular grid.
Parameters
----------
x: 1D floating point array or list containing the individual grid
positions. Points must be given in ascending order.
sigma: The sigma of the Gaussian with which the function living on x
should be smoothed, in the same units as x.
dxmax: (optional) The maximum distance up to which smoothing is
performed, in the same units as x. Default is 3.01*sigma.
Returns
-------
ibegin: integer array of the same size as x
ibegin[i] is the minimum grid index to consider when computing the
smoothed value at grid index i
nval: integer array of the same size as x
nval[i] is the number of indices to consider when computing the
smoothed value at grid index i.
wgt: list with the same number of entries as x
wgt[i] is an array with nval[i] entries containing the
normalized smoothing weights.
"""
if dxmax is None:
dxmax = self.effective_smoothing_width*sigma
x = np.asarray(x)
ibegin = np.searchsorted(x, x-dxmax)
nval = np.searchsorted(x, x+dxmax) - ibegin
wgt = []
expfac = 1. / (2. * sigma*sigma)
for i in range(x.size):
t = x[ibegin[i]:ibegin[i]+nval[i]]-x[i]
t = np.exp(-t*t*expfac)
t *= 1./np.sum(t)
wgt.append(t)
return ibegin, nval, wgt
def _apply_kernel_along_array(self, power, startindex, endindex,
distances, smooth_length, smoothing_width,
ibegin, nval, wgt):
if smooth_length == 0.0:
return power[startindex:endindex]
p_smooth = np.zeros(endindex-startindex, dtype=power.dtype)
for i in xrange(startindex, endindex):
imin = max(startindex, ibegin[i])
imax = min(endindex, ibegin[i]+nval[i])
p_smooth[imin:imax] += (power[i] *
wgt[i][imin-ibegin[i]:imax-imin+ibegin[i]])
return p_smooth
def _apply_along_axis(self, axis, arr, startindex, endindex, distances,
smooth_length, smoothing_width):
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))
ibegin, nval, wgt = self._precompute(
distances, smooth_length, smooth_length*smoothing_width)
ind = np.zeros(nd-1, dtype=np.int)
i = np.zeros(nd, dtype=object)
shape = arr.shape
indlist = np.asarray(range(nd))
indlist = np.delete(indlist, axis)
i[axis] = slice(None, None)
outshape = np.asarray(shape).take(indlist)
i.put(indlist, ind)
Ntot = np.product(outshape)
holdshape = outshape
slicedArr = arr[tuple(i.tolist())]
res = self._apply_kernel_along_array(
slicedArr, startindex, endindex, distances,
smooth_length, smoothing_width, ibegin, nval, wgt)
outshape = np.asarray(arr.shape)
outshape[axis] = endindex - startindex
outarr = np.zeros(outshape, dtype=arr.dtype)
outarr[tuple(i.tolist())] = res
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 = self._apply_kernel_along_array(
slicedArr, startindex, endindex, distances,
smooth_length, smoothing_width, ibegin, nval, wgt)
outarr[tuple(i.tolist())] = res
k += 1
return outarr
def _smooth(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
local_data_Q = False
if x.distribution_strategy == 'not':
local_data_Q = True
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 != affected_axis:
local_data_Q = True
else:
start_index = x.val.distributor.local_start
start_distance = distance_array[start_index]
augmented_start_distance = \
(start_distance -
self.effective_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.effective_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[augmented_slice]
else:
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
true_start = 0
true_end = x.shape[affected_axis]
# perform the convolution along the affected axes
# currently only one axis is supported
data_axis = affected_axes[0]
if inverse:
true_sigma = 1. / self.sigma
else:
true_sigma = self.sigma
local_result = self._apply_along_axis(
data_axis,
augmented_data,
startindex=true_start,
endindex=true_end,
distances=augmented_distance_array,
smooth_length=true_sigma,
smoothing_width=self.effective_smoothing_width)
result = x.copy_empty()
result.val.set_local_data(local_result, copy=False)
return result
# -*- coding: utf-8 -*-
import numpy as np
from nifty.operators.fft_operator import FFTOperator
from .smoothing_operator import SmoothingOperator
class FFTSmoothingOperator(SmoothingOperator):
def _smooth(self, x, spaces, inverse):
Transformator = FFTOperator(x.domain[spaces[0]])
# transform to the (global-)default codomain and perform all remaining
# steps therein
transformed_x = Transformator(x, spaces=spaces)
codomain = transformed_x.domain[spaces[0]]
coaxes = transformed_x.domain_axes[spaces[0]]
# create the kernel using the knowledge of codomain about itself
axes_local_distribution_strategy = \
transformed_x.val.get_axes_local_distribution_strategy(axes=coaxes)
kernel = codomain.get_distance_array(
distribution_strategy=axes_local_distribution_strategy)
if self.log_distances:
kernel.apply_scalar_function(np.log, inplace=True)
kernel.apply_scalar_function(
codomain.get_fft_smoothing_kernel_function(self.sigma),
inplace=True)
# now, apply the kernel to transformed_x
# this is done node-locally utilizing numpys reshaping in order to
# apply the kernel to the correct axes
local_transformed_x = transformed_x.val.get_local_data(copy=False)
local_kernel = kernel.get_local_data(copy=False)
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)
# apply the kernel
if inverse:
local_transformed_x /= local_kernel
else:
local_transformed_x *= local_kernel
transformed_x.val.set_local_data(local_transformed_x, copy=False)
smoothed_x = Transformator.adjoint_times(transformed_x, spaces=spaces)
result = x.copy_empty()
result.set_val(smoothed_x, copy=False)
return result
......@@ -16,12 +16,12 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
import abc
import numpy as np
import nifty.nifty_utilities as utilities
from nifty.operators.endomorphic_operator import EndomorphicOperator
from nifty.operators.fft_operator import FFTOperator
from d2o import STRATEGIES
from nifty.spaces import RGSpace, GLSpace, HPSpace, PowerSpace
class SmoothingOperator(EndomorphicOperator):
......@@ -75,29 +75,81 @@ class SmoothingOperator(EndomorphicOperator):
"""
_fft_smoothing_spaces = [RGSpace,
GLSpace,
HPSpace]
_direct_smoothing_spaces = [PowerSpace]
def __new__(cls, domain, *args, **kwargs):
if cls is SmoothingOperator:
domain = cls._parse_domain(domain)
if len(domain) != 1:
raise ValueError("SmoothingOperator only accepts exactly one "
"space as input domain.")
if np.any([isinstance(domain[0], sp)
for sp in cls._fft_smoothing_spaces]):
from .fft_smoothing_operator import FFTSmoothingOperator
return super(SmoothingOperator, cls).__new__(
FFTSmoothingOperator, domain, *args, **kwargs)
elif np.any([isinstance(domain[0], sp)
for sp in cls._direct_smoothing_spaces]):
from .direct_smoothing_operator import DirectSmoothingOperator
return super(SmoothingOperator, cls).__new__(
DirectSmoothingOperator, domain, *args, **kwargs)
else:
raise NotImplementedError("For the given Space smoothing "
" is not available.")
else:
print 'new 4'
return super(SmoothingOperator, cls).__new__(cls,
domain,
*args,
**kwargs)
# ---Overwritten properties and methods---
def __init__(self, domain, sigma, log_distances=False,
default_spaces=None):
super(SmoothingOperator, self).__init__(default_spaces)
# # the _parse_domain is already done in the __new__ method
# self._domain = self._parse_domain(domain)
# if len(self.domain) != 1:
# raise ValueError("SmoothingOperator only accepts exactly one "
# "space as input domain.")
self._domain = self._parse_domain(domain)
if len(self.domain) != 1:
raise ValueError(
'ERROR: SmoothOperator accepts only exactly one '
'space as input domain.'
)
self.sigma = sigma
self.log_distances = log_distances
self._direct_smoothing_width = 3.
def _inverse_times(self, x, spaces):
return self._smoothing_helper(x, spaces, inverse=True)
if self.sigma == 0:
return x.copy()
# the domain of the smoothing operator contains exactly one space.
# Hence, if spaces is None, but we passed LinearOperator's
# _check_input_compatibility, we know that x is also solely defined
# on that space
if spaces is None:
spaces = (0,)
return self._smooth(x, spaces, inverse=True)
def _times(self, x, spaces):
return self._smoothing_helper(x, spaces, inverse=False)
if self.sigma == 0:
return x.copy()
# the domain of the smoothing operator contains exactly one space.
# Hence, if spaces is None, but we passed LinearOperator's
# _check_input_compatibility, we know that x is also solely defined
# on that space
if spaces is None:
spaces = (0,)
return self._smooth(x, spaces, inverse=False)
# ---Mandatory properties and methods---
@property
......@@ -130,285 +182,6 @@ class SmoothingOperator(EndomorphicOperator):
def log_distances(self, log_distances):
self._log_distances = bool(log_distances)
@staticmethod
def _precompute(x, sigma, dxmax=None):
""" Performs the precomputations for Gaussian smoothing on a 1D irregular grid.
Parameters
----------
x: 1D floating point array or list containing the individual grid
positions. Points must be given in ascending order.
sigma: The sigma of the Gaussian with which the function living on x
should be smoothed, in the same units as x.
dxmax: (optional) The maximum distance up to which smoothing is
performed, in the same units as x. Default is 3.01*sigma.
Returns
-------
ibegin: integer array of the same size as x
ibegin[i] is the minimum grid index to consider when computing the
smoothed value at grid index i
nval: integer array of the same size as x
nval[i] is the number of indices to consider when computing the
smoothed value at grid index i.
wgt: list with the same number of entries as x
wgt[i] is an array with nval[i] entries containing the
normalized smoothing weights.
"""
if dxmax is None:
dxmax = 3.01*sigma
x = np.asarray(x)
ibegin = np.searchsorted(x, x-dxmax)
nval = np.searchsorted(x, x+dxmax) - ibegin
wgt = []
expfac = 1. / (2. * sigma*sigma)
for i in range(x.size):
t = x[ibegin[i]:ibegin[i]+nval[i]]-x[i]
t = np.exp(-t*t*expfac)
t *= 1./np.sum(t)
wgt.append(t)
return ibegin, nval, wgt
@staticmethod
def _apply_kernel_along_array(
power, startindex, endindex, distances, smooth_length,
smoothing_width, ibegin, nval, wgt):
if smooth_length == 0.0:
return power[startindex:endindex]
p_smooth = np.zeros(endindex-startindex, dtype=power.dtype)
for i in xrange(startindex, endindex):
imin = max(startindex, ibegin[i])
imax = min(endindex, ibegin[i]+nval[i])
p_smooth[imin:imax] += (power[i] *
wgt[i][imin-ibegin[i]:imax-imin+ibegin[i]])
return p_smooth
@staticmethod
def _getShape(a):
return tuple(a.shape)
@staticmethod
def _apply_along_axis(axis, arr, startindex, endindex, distances,
smooth_length, smoothing_width):
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))
ibegin, nval, wgt = SmoothingOperator._precompute(
distances, smooth_length, smooth_length*smoothing_width)
ind = np.zeros(nd-1, dtype=np.int)
i = np.zeros(nd, dtype=object)
shape = SmoothingOperator._getShape(arr)
indlist = np.asarray(range(nd))
indlist = np.delete(indlist, axis)
i[axis] = slice(None, None)
outshape = np.asarray(shape).take(indlist)
i.put(indlist, ind)
Ntot = np.product(outshape)
holdshape = outshape
slicedArr = arr[tuple(i.tolist())]
res = SmoothingOperator._apply_kernel_along_array(
slicedArr, startindex, endindex, distances,
smooth_length, smoothing_width, ibegin, nval, wgt)
outshape = np.asarray(SmoothingOperator._getShape(arr))
outshape[axis] = endindex - startindex
outarr = np.zeros(outshape, dtype=arr.dtype)
outarr[tuple(i.tolist())] = res
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 = SmoothingOperator._apply_kernel_along_array(
slicedArr, startindex, endindex, distances,
smooth_length, smoothing_width, ibegin, nval, wgt)
outarr[tuple(i.tolist())] = res
k += 1
return outarr
def _smoothing_helper(self, x, spaces, inverse):
if self.sigma == 0:
return x.copy()
# the domain of the smoothing operator contains exactly one space.
# Hence, if spaces is None, but we passed LinearOperator's
# _check_input_compatibility, we know that x is also solely defined
# on that space
if spaces is None:
spaces = (0,)
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
# steps therein
transformed_x = Transformator(x, spaces=spaces)
codomain = transformed_x.domain[spaces[0]]
coaxes = transformed_x.domain_axes[spaces[0]]
# create the kernel using the knowledge of codomain about itself
axes_local_distribution_strategy = \
transformed_x.val.get_axes_local_distribution_strategy(axes=coaxes)
kernel = codomain.get_distance_array(
distribution_strategy=axes_local_distribution_strategy)
if self.log_distances:
kernel.apply_scalar_function(np.log, inplace=True)
kernel.apply_scalar_function(
codomain.get_fft_smoothing_kernel_function(self.sigma),
inplace=True)