Commit 45e3eb9b authored by Martin Reinecke's avatar Martin Reinecke
Browse files

merge fix_smoothing

parents b1f4314d c639445d
......@@ -19,7 +19,6 @@ test_min:
stage: test
script:
- nosetests
- nosetests3
- nosetests -x --with-coverage --cover-package=nifty --cover-branches
- >
coverage report | grep TOTAL | awk '{ print "TOTAL: "$6; }'
......@@ -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
......@@ -60,8 +60,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
......@@ -4,17 +4,58 @@ from __future__ import division
from builtins import range
import numpy as np
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
......@@ -50,7 +91,7 @@ class DirectSmoothingOperator(SmoothingOperator):
"""
if dxmax is None:
dxmax = self.effective_smoothing_width*sigma
dxmax = self._effective_smoothing_width*sigma
x = np.asarray(x)
......@@ -139,7 +180,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]]
......@@ -165,18 +206,13 @@ 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)
return local_result
......@@ -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)
......@@ -36,17 +70,10 @@ class FFTSmoothingOperator(SmoothingOperator):
kernel = codomain.get_distance_array()
#MR FIXME: this causes calls of log(0.) which should probably be avoided
if self.log_distances:
kernel=np.log(kernel)
#kernel.apply_scalar_function(
# codomain.get_fft_smoothing_kernel_function(self.sigma),
# inplace=True)
kernel = codomain.get_fft_smoothing_kernel_function(self.sigma)(kernel)
# 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
local_kernel = kernel
......@@ -55,12 +82,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=local_transformed_x
......
# 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,13 +18,13 @@
import unittest
import numpy as np
from numpy.testing import assert_equal, assert_approx_equal,\
assert_allclose
from numpy.testing import assert_equal, assert_allclose
from nifty import Field,\
RGSpace,\
PowerSpace,\
SmoothingOperator
FFTSmoothingOperator,\
DirectSmoothingOperator
from itertools import product
from test.common import expand
......@@ -39,10 +39,9 @@ def _get_rtol(tp):
class SmoothingOperator_Tests(unittest.TestCase):
spaces = [RGSpace(128)]
@expand(product(spaces, [0., .5, 5.], [True, False]))
def test_property(self, space, sigma, log_distances):
op = SmoothingOperator.make(space, sigma=sigma,
log_distances=log_distances)
@expand(product(spaces, [0., .5, 5.]))
def test_property(self, space, sigma):
op = FFTSmoothingOperator(space, sigma=sigma)
if op.domain[0] != space:
raise TypeError
if op.unitary != False:
......@@ -51,44 +50,30 @@ class SmoothingOperator_Tests(unittest.TestCase):
raise ValueError
if op.sigma != sigma:
raise ValueError
if op.log_distances != log_distances:
raise ValueError
@expand(product(spaces, [0., .5, 5.], [True, False]))
def test_adjoint_times(self, space, sigma, log_distances):
op = SmoothingOperator.make(space, sigma=sigma,
log_distances=log_distances)
@expand(product(spaces, [0., .5, 5.]))
def test_adjoint_times(self, space, sigma):
op = FFTSmoothingOperator(space, sigma=sigma)
rand1 = Field.from_random('normal', domain=space)
rand2 = Field.from_random('normal', domain=space)
tt1 = rand1.vdot(op.times(rand2))
tt2 = rand2.vdot(op.adjoint_times(rand1))
assert_approx_equal(tt1, tt2)
assert_allclose(tt1, tt2)
@expand(product(spaces, [0., .5, 5.], [False]))
def test_times(self, space, sigma, log_distances):
op = SmoothingOperator.make(space, sigma=sigma,
log_distances=log_distances)
@expand(product(spaces, [0., .5, 5.]))
def test_times(self, space, sigma):
op = FFTSmoothingOperator(space, sigma=sigma)
rand1 = Field(space, val=0.)
rand1.val[0] = 1.
tt1 = op.times(rand1)
assert_approx_equal(1, tt1.sum())
@expand(product(spaces, [0., .5, 5.], [True, False]))
def test_inverse_adjoint_times(self, space, sigma, log_distances):
op = SmoothingOperator.make(space, sigma=sigma,
log_distances=log_distances)
rand1 = Field.from_random('normal', domain=space)
rand2 = Field.from_random('normal', domain=space)
tt1 = rand1.vdot(op.inverse_times(rand2))
tt2 = rand2.vdot(op.inverse_adjoint_times(rand1))
assert_approx_equal(tt1, tt2)
assert_allclose(1, tt1.sum())
@expand(product([128, 256], [1, 0.4], [0., 1., 3.7],
[np.float64, np.complex128]))
def test_smooth_regular1(self, sz, d, sigma, tp):
tol = _get_rtol(tp)
sp = RGSpace(sz, harmonic=True, distances=d)
smo = SmoothingOperator.make(sp, sigma=sigma)
smo = FFTSmoothingOperator(sp, sigma=sigma)
inp = Field.from_random(domain=sp, random_type='normal', std=1, mean=4,
dtype=tp)
out = smo(inp)
......@@ -99,7 +84,7 @@ class SmoothingOperator_Tests(unittest.TestCase):
def test_smooth_regular2(self, sz1, sz2, d1, d2, sigma, tp):
tol = _get_rtol(tp)
sp = RGSpace([sz1, sz2], distances=[d1, d2], harmonic=True)
smo = SmoothingOperator.make(sp, sigma=sigma)
smo = FFTSmoothingOperator(sp, sigma=sigma)
inp = Field.from_random(domain=sp, random_type='normal', std=1, mean=4,
dtype=tp)
out = smo(inp)
......@@ -111,7 +96,7 @@ class SmoothingOperator_Tests(unittest.TestCase):
tol = _get_rtol(tp)
sp = RGSpace(sz, harmonic=True)