Commit 9fd1fd9a authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge remote-tracking branch 'remotes/origin/fix_smoothing' into nightly

# Conflicts:
#	test/test_operators/test_smoothing_operator.py
parents d6d8935d c639445d
......@@ -130,7 +130,7 @@ if __name__ == "__main__":
plotter.plot.zmin = 0.
plotter.plot.zmax = 3.
sm = ift.SmoothingOperator.make(plot_space, sigma=0.03)
sm = ift.FFTSmoothingOperator(plot_space, sigma=0.03)
plotter(ift.log(ift.sqrt(sm(ift.Field(plot_space, val=variance.val.real)))), path='uncertainty.html')
plotter.plot.zmin = np.real(mock_signal.min());
......
......@@ -72,12 +72,6 @@ class ComposedOperator(LinearOperator):
>>> f = Field.from_random('normal', domain=(x1,x2))
>>> FFT.times(f)
See Also
--------
EndomorphicOperator, ProjectionOperator,
DiagonalOperator, SmoothingOperator, ResponseOperator,
PropagatorOperator, ComposedOperator
"""
# ---Overwritten properties and methods---
......
......@@ -41,23 +41,6 @@ class InvertibleOperatorMixin(object):
Preconditioner that is used by ConjugateGradient if no minimizer was
given.
Attributes
----------
Raises
------
Notes
-----
Examples
--------
The PropagatorOperator inherits from InvertibleOperatorMixin.
See Also
--------
PropagatorOperator
"""
def __init__(self, inverter=None, preconditioner=None,
......
......@@ -66,12 +66,6 @@ class LinearOperator(with_metaclass(
the LinearOperator. A LinearOperator must have the attributes domain,
target and unitary to be properly defined.
See Also
--------
EndomorphicOperator, ProjectionOperator,
DiagonalOperator, SmoothingOperator, ResponseOperator,
PropagatorOperator, ComposedOperator
"""
def __init__(self, default_spaces=None):
......
......@@ -3,7 +3,7 @@ import numpy as np
from ... import Field,\
FieldArray
from ..linear_operator import LinearOperator
from ..smoothing_operator import SmoothingOperator
from ..smoothing_operator import FFTSmoothingOperator
from ..composed_operator import ComposedOperator
from ..diagonal_operator import DiagonalOperator
......@@ -81,8 +81,8 @@ class ResponseOperator(LinearOperator):
"exposure do not match")
for ii in range(len(kernel_smoothing)):
kernel_smoothing[ii] = SmoothingOperator.make(self._domain[ii],
sigma=sigma[ii])
kernel_smoothing[ii] = FFTSmoothingOperator(self._domain[ii],
sigma=sigma[ii])
kernel_exposure[ii] = DiagonalOperator(self._domain[ii],
diagonal=exposure[ii])
......
......@@ -16,4 +16,5 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from .smoothing_operator import SmoothingOperator
from .fft_smoothing_operator import FFTSmoothingOperator
from .direct_smoothing_operator import DirectSmoothingOperator
......@@ -6,17 +6,58 @@ import numpy as np
from d2o import STRATEGIES
from .smoothing_operator import SmoothingOperator
from ..endomorphic_operator import EndomorphicOperator
class DirectSmoothingOperator(SmoothingOperator):
class DirectSmoothingOperator(EndomorphicOperator):
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
super(DirectSmoothingOperator, self).__init__(default_spaces)
self._domain = self._parse_domain(domain)
if len(self._domain) != 1:
raise ValueError("DirectSmoothingOperator only accepts exactly one"
" space as input domain.")
self._sigma = sigma
self._log_distances = log_distances
self._effective_smoothing_width = 3.01
def _times(self, x, spaces):
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)
# ---Mandatory properties and methods---
@property
def domain(self):
return self._domain
@property
def self_adjoint(self):
return True
@property
def unitary(self):
return False
# ---Added properties and methods---
@property
def sigma(self):
return self._sigma
@property
def log_distances(self):
return self._log_distances
def _add_attributes_to_copy(self, copy, **kwargs):
copy.effective_smoothing_width = self.effective_smoothing_width
......@@ -52,7 +93,7 @@ class DirectSmoothingOperator(SmoothingOperator):
"""
if dxmax is None:
dxmax = self.effective_smoothing_width*sigma
dxmax = self._effective_smoothing_width*sigma
x = np.asarray(x)
......@@ -141,7 +182,7 @@ class DirectSmoothingOperator(SmoothingOperator):
return outarr
def _smooth(self, x, spaces, inverse):
def _smooth(self, x, spaces):
# infer affected axes
# we rely on the knowledge, that `spaces` is a tuple with length 1.
affected_axes = x.domain_axes[spaces[0]]
......@@ -175,14 +216,14 @@ class DirectSmoothingOperator(SmoothingOperator):
start_distance = distance_array[start_index]
augmented_start_distance = \
(start_distance -
self.effective_smoothing_width*self.sigma)
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)
(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
......@@ -212,19 +253,14 @@ class DirectSmoothingOperator(SmoothingOperator):
# 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)
smooth_length=self.sigma,
smoothing_width=self._effective_smoothing_width)
result = x.copy_empty()
result.val.set_local_data(local_result, copy=False)
......
......@@ -3,22 +3,56 @@
from builtins import range
import numpy as np
from nifty.operators.fft_operator import FFTOperator
from ..endomorphic_operator import EndomorphicOperator
from ..fft_operator import FFTOperator
from .smoothing_operator import SmoothingOperator
class FFTSmoothingOperator(EndomorphicOperator):
class FFTSmoothingOperator(SmoothingOperator):
def __init__(self, domain, sigma, log_distances=False,
def __init__(self, domain, sigma,
default_spaces=None):
super(FFTSmoothingOperator, self).__init__(
domain=domain,
sigma=sigma,
log_distances=log_distances,
default_spaces=default_spaces)
super(FFTSmoothingOperator, self).__init__(default_spaces)
self._domain = self._parse_domain(domain)
if len(self._domain) != 1:
raise ValueError("SmoothingOperator only accepts exactly one "
"space as input domain.")
self._sigma = sigma
self._transformator_cache = {}
def _times(self, x, spaces):
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)
# ---Mandatory properties and methods---
@property
def domain(self):
return self._domain
@property
def self_adjoint(self):
return True
@property
def unitary(self):
return False
# ---Added properties and methods---
@property
def sigma(self):
return self._sigma
def _add_attributes_to_copy(self, copy, **kwargs):
copy._transformator_cache = self._transformator_cache
......@@ -26,7 +60,7 @@ class FFTSmoothingOperator(SmoothingOperator):
copy, **kwargs)
return copy
def _smooth(self, x, spaces, inverse):
def _smooth(self, x, spaces):
# transform to the (global-)default codomain and perform all remaining
# steps therein
transformator = self._get_transformator(x.dtype)
......@@ -41,16 +75,12 @@ class FFTSmoothingOperator(SmoothingOperator):
kernel = codomain.get_distance_array(
distribution_strategy=axes_local_distribution_strategy)
#MR FIXME: this causes calls of log(0.) which should probably be avoided
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
# this is done node-locally utilizing numpy's 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)
......@@ -59,12 +89,7 @@ class FFTSmoothingOperator(SmoothingOperator):
for i in range(len(transformed_x.shape))]
local_kernel = np.reshape(local_kernel, reshaper)
# apply the kernel
if inverse:
#MR FIXME: danger of having division by zero or overflows
local_transformed_x /= local_kernel
else:
local_transformed_x *= local_kernel
local_transformed_x *= local_kernel
transformed_x.val.set_local_data(local_transformed_x, copy=False)
......
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import abc
import numpy as np
from nifty.operators.endomorphic_operator import EndomorphicOperator
from nifty.spaces import RGSpace, GLSpace, HPSpace, PowerSpace
class SmoothingOperator(EndomorphicOperator):
""" NIFTY class for smoothing operators.
The NIFTy SmoothingOperator smooths Fields, with a given kernel length.
Fields which are not living over a PowerSpace are smoothed
via a gaussian convolution. Fields living over the PowerSpace are directly
smoothed.
Parameters
----------
domain : DomainObject, i.e. Space or FieldType
The Space on which the operator acts. The SmoothingOperator
can only live on one space or FieldType
sigma : float
Sets the length of the Gaussian convolution kernel
log_distances : boolean *optional*
States whether the convolution happens on the logarithmic grid or not
(default: False).
default_spaces : tuple of ints *optional*
Defines on which space(s) of a given field the Operator acts by
default (default: None).
Attributes
----------
domain : DomainObject, i.e. Space or FieldType
The domain on which the Operator's input Field lives.
target : tuple of DomainObjects, i.e. Spaces and FieldTypes
The domain in which the outcome of the operator lives. As the Operator
is endomorphic this is the same as its domain.
unitary : boolean
Indicates whether the Operator is unitary or not.
self_adjoint : boolean
Indicates whether the operator is self_adjoint or not.
sigma : float
Sets the length of the Gaussian convolution kernel
log_distances : boolean
States whether the convolution happens on the logarithmic grid or not.
Raises
------
ValueError
Raised if
* the given domain inherits more than one space. The
SmoothingOperator acts only on one Space.
Notes
-----
Examples
--------
>>> x = RGSpace(5)
>>> S = SmoothingOperator(x, sigma=1.)
>>> f = Field(x, val=[1,2,3,4,5])
>>> S.times(f).val
<distributed_data_object>
array([ 3., 3., 3., 3., 3.])
See Also
--------
DiagonalOperator, SmoothingOperator,
PropagatorOperator, ProjectionOperator,
ComposedOperator
"""
@staticmethod
def make(domain, sigma, log_distances=False, default_spaces=None):
_fft_smoothing_spaces = [RGSpace,
GLSpace,
HPSpace]
_direct_smoothing_spaces = [PowerSpace]
domain = SmoothingOperator._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 _fft_smoothing_spaces]):
from .fft_smoothing_operator import FFTSmoothingOperator
return FFTSmoothingOperator (domain, sigma, log_distances,\
default_spaces)
elif np.any([isinstance(domain[0], sp)
for sp in _direct_smoothing_spaces]):
from .direct_smoothing_operator import DirectSmoothingOperator
return DirectSmoothingOperator (domain, sigma, log_distances,\
default_spaces)
else:
raise NotImplementedError("For the given Space smoothing "
" is not available.")
# ---Overwritten properties and methods---
def __init__(self, domain, sigma, log_distances=False,
default_spaces=None):
super(SmoothingOperator, self).__init__(default_spaces)
self._domain = self._parse_domain(domain)
if len(self._domain) != 1:
raise ValueError("SmoothingOperator only accepts exactly one "
"space as input domain.")
self._sigma = sigma
self._log_distances = log_distances
def _add_attributes_to_copy(self, copy, **kwargs):
copy._domain = self._domain
copy._sigma = self._sigma
copy._log_distances = self._log_distances
copy = super(SmoothingOperator, self)._add_attributes_to_copy(copy,
**kwargs)
return copy
def _inverse_times(self, x, spaces):
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):
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
def domain(self):
return self._domain
@property
def self_adjoint(self):
return True
@property
def unitary(self):
return False
# ---Added properties and methods---
@property
def sigma(self):
return self._sigma
@property
def log_distances(self):
return self._log_distances
@abc.abstractmethod
def _smooth(self, x, spaces, inverse):
raise NotImplementedError
......@@ -18,7 +18,7 @@
import unittest
import numpy as np
from numpy.testing import assert_approx_equal, assert_allclose
from numpy.testing import assert_allclose
from nifty import Field,\
RGSpace,\
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment