Commit 0b6d6abd authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'nightly' into better_minimizers

# Conflicts:
#	nifty/minimization/conjugate_gradient.py
#	nifty/minimization/descent_minimizer.py
parents 576aa27e 6a208374
......@@ -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());
......
......@@ -33,7 +33,7 @@ if __name__ == "__main__":
# Total side-length of the domain
L = 2.
# Grid resolution (pixels per axis)
N_pixels = 512
N_pixels = 128
signal_space = RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
harmonic_space = FFTOperator.get_default_codomain(signal_space)
......
......@@ -103,21 +103,24 @@ def clipped_exp(x):
def limited_exp(x):
return _math_helper(x, _limited_exp_helper)
def _limited_exp_helper(x):
thr = 200.
mask = x>thr
mask = x > thr
if np.count_nonzero(mask) == 0:
return np.exp(x)
result = ((1.-thr) + x)*np.exp(thr)
result[~mask] = np.exp(x[~mask])
return result
def limited_exp_deriv(x):
return _math_helper(x, _limited_exp_deriv_helper)
def _limited_exp_deriv_helper(x):
thr = 200.
mask = x>thr
mask = x > thr
if np.count_nonzero(mask) == 0:
return np.exp(x)
result = np.empty_like(x)
......
......@@ -23,7 +23,8 @@ from keepers import Loggable
from future.utils import with_metaclass
class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {}))):
class Energy(with_metaclass(NiftyMeta,
type('NewBase', (Loggable, object), {}))):
""" Provides the functional used by minimization schemes.
The Energy object is an implementation of a scalar function including its
......
......@@ -18,7 +18,6 @@
from __future__ import division
from builtins import zip
#from builtins import str
from builtins import range
import ast
......@@ -282,8 +281,8 @@ class Field(Loggable, Versionable, object):
# ---Powerspectral methods---
def power_analyze(self, spaces=None, logarithmic=None, nbin=None,
binbounds=None, keep_phase_information=False):
def power_analyze(self, spaces=None, binbounds=None,
keep_phase_information=False):
""" Computes the square root power spectrum for a subspace of `self`.
Creates a PowerSpace for the space addressed by `spaces` with the given
......@@ -297,16 +296,8 @@ class Field(Loggable, Versionable, object):
spaces : int *optional*
The subspace for which the powerspectrum shall be computed
(default : None).
logarithmic : boolean *optional*
True if the output PowerSpace should use logarithmic binning.
{default : None}
nbin : int *optional*
The number of bins the resulting PowerSpace shall have
(default : None).
if nbin==None : maximum number of bins is used
binbounds : array-like *optional*
Inner bounds of the bins (default : None).
Overrides nbin and logarithmic.
if binbounds==None : bins are inferred.
keep_phase_information : boolean, *optional*
If False, return a real-valued result containing the power spectrum
......@@ -374,8 +365,6 @@ class Field(Loggable, Versionable, object):
parts = [self._single_power_analyze(
work_field=part,
space_index=space_index,
logarithmic=logarithmic,
nbin=nbin,
binbounds=binbounds)
for part in parts]
......@@ -387,8 +376,7 @@ class Field(Loggable, Versionable, object):
return result_field
@classmethod
def _single_power_analyze(cls, work_field, space_index, logarithmic, nbin,
binbounds):
def _single_power_analyze(cls, work_field, space_index, binbounds):
if not work_field.domain[space_index].harmonic:
raise ValueError(
......@@ -407,7 +395,6 @@ class Field(Loggable, Versionable, object):
harmonic_domain = work_field.domain[space_index]
power_domain = PowerSpace(harmonic_partner=harmonic_domain,
distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds)
power_spectrum = cls._calculate_power_spectrum(
field_val=work_field.val,
......@@ -689,7 +676,7 @@ class Field(Loggable, Versionable, object):
result_list[0].val.get_axes_local_distribution_strategy(
result_list[0].domain_axes[power_space_index])
if pindex.distribution_strategy is not local_distribution_strategy:
if pindex.distribution_strategy != local_distribution_strategy:
raise AttributeError(
"The distribution_strategy of pindex does not fit the "
"slice_local distribution strategy of the synthesized field.")
......@@ -1603,7 +1590,7 @@ class Field(Loggable, Versionable, object):
new_field.dtype = np.dtype(hdf5_group.attrs['dtype'])
new_field.distribution_strategy =\
hdf5_group.attrs['distribution_strategy']
str(hdf5_group.attrs['distribution_strategy'])
return new_field
......
......@@ -26,13 +26,14 @@ class LogNormalWienerFilterCurvature(InvertibleOperatorMixin,
"""
def __init__(self, R, N, S, d, position, inverter=None,
preconditioner=None, fft4exp=None, **kwargs):
preconditioner=None, fft4exp=None, offset=None, **kwargs):
self._cache = {}
self.R = R
self.N = N
self.S = S
self.d = d
self.position = position
self.offset = offset
if preconditioner is None:
preconditioner = self.S.times
self._domain = self.S.domain
......@@ -55,6 +56,7 @@ class LogNormalWienerFilterCurvature(InvertibleOperatorMixin,
copy.N = self.N.copy()
copy.S = self.S.copy()
copy.d = self.d.copy()
copy.offset = self.offset
if 'position' in kwargs:
copy.position = kwargs['position']
else:
......@@ -92,7 +94,10 @@ class LogNormalWienerFilterCurvature(InvertibleOperatorMixin,
@property
@memo
def _expp_sspace(self):
return clipped_exp(self._fft(self.position))
xs = self._fft(self.position)
if self.offset is not None:
xs += self.offset
return clipped_exp(xs)
@property
@memo
......
......@@ -24,12 +24,14 @@ class LogNormalWienerFilterEnergy(Energy):
The prior signal covariance in harmonic space.
"""
def __init__(self, position, d, R, N, S, fft4exp=None, old_curvature=None):
def __init__(self, position, d, R, N, S, fft4exp=None, old_curvature=None,
offset=None):
super(LogNormalWienerFilterEnergy, self).__init__(position=position)
self.d = d
self.R = R
self.N = N
self.S = S
self.offset = offset
if fft4exp is None:
self._fft = create_composed_fft_operator(self.S.domain,
......@@ -43,7 +45,8 @@ class LogNormalWienerFilterEnergy(Energy):
def at(self, position):
return self.__class__(position=position, d=self.d, R=self.R, N=self.N,
S=self.S, fft4exp=self._fft,
old_curvature=self._curvature)
old_curvature=self._curvature,
offset=self.offset)
@property
@memo
......@@ -66,7 +69,8 @@ class LogNormalWienerFilterEnergy(Energy):
S=self.S,
d=self.d,
position=self.position,
fft4exp=self._fft)
fft4exp=self._fft,
offset=self.offset)
else:
self._curvature = \
self._old_curvature.copy(position=self.position)
......
......@@ -139,7 +139,7 @@ class LineSearchStrongWolfe(LineSearch):
phi_alpha1 = le_alpha1.value
if (phi_alpha1 > phi_0 + self.c1*alpha1*phiprime_0) or \
((phi_alpha1 >= phi_alpha0) and (i > 0)):
((phi_alpha1 >= phi_alpha0) and (iteration_number > 1)):
le_star = self._zoom(alpha0, alpha1, phi_0, phiprime_0,
phi_alpha0, phiprime_alpha0, phi_alpha1,
le_0)
......
......@@ -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---
......
......@@ -136,15 +136,15 @@ class DiagonalOperator(EndomorphicOperator):
def _adjoint_times(self, x, spaces):
return self._times_helper(x, spaces,
operation=lambda z: z.adjoint().__mul__)
operation=lambda z: z.conjugate().__mul__)
def _inverse_times(self, x, spaces):
return self._times_helper(
x, spaces, operation=lambda z: z.__rtruediv__)
def _adjoint_inverse_times(self, x, spaces):
return self._times_helper(x, spaces,
operation=lambda z: z.adjoint().__rtruediv__)
return self._times_helper(
x, spaces, operation=lambda z: z.conjugate().__rtruediv__)
def diagonal(self, bare=False, copy=True):
""" Returns the diagonal of the Operator.
......
......@@ -17,8 +17,6 @@ def buildIdx(nr, lmax):
new_dtype = np.float32
elif nr.dtype == np.dtype('complex128'):
new_dtype = np.float64
elif nr.dtype == np.dtype('complex256'):
new_dtype = np.float128
else:
raise TypeError("dtype of nr not supported.")
......
......@@ -42,23 +42,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.