Commit 80086a53 authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'nightly' into byebye_zerocenter

# Conflicts:
#	test/test_spaces/test_rg_space.py
parents 10e3efd2 fb6d3eb1
......@@ -68,7 +68,7 @@ if __name__ == "__main__":
h_space = fft.target[0]
# Set up power space
p_space = PowerSpace(h_space, logarithmic=True,
p_space = PowerSpace(h_space, binbounds=PowerSpace.useful_binbounds(h_space,logarithmic=True),
distribution_strategy=distribution_strategy)
# Choose the prior correlation structure and defining correlation operator
......@@ -91,7 +91,8 @@ if __name__ == "__main__":
R = AdjointFFTResponse(fft, Instrument)
noise = 1.
N = DiagonalOperator(s_space, diagonal=noise, bare=True)
ndiag = Field(s_space, noise).weight(1)
N = DiagonalOperator(s_space, ndiag)
n = Field.from_random(domain=s_space,
random_type='normal',
std=sqrt(noise),
......
......@@ -39,7 +39,8 @@ if __name__ == "__main__":
R_harmonic = ComposedOperator([fft, R], default_spaces=[0, 0])
# Setting up the noise covariance and drawing a random noise realization
N = DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True)
ndiag = Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
N = DiagonalOperator(data_domain, ndiag)
noise = Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(exp(mock_signal)) + noise #|\label{code:wf_mock_data}|
......
......@@ -93,8 +93,8 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=(0, 1, 0, 1))
# Setting up the noise covariance and drawing a random noise realization
N = ift.DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise,
bare=True)
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
N = ift.DiagonalOperator(data_domain, ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise),
mean=0)
......@@ -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());
......
......@@ -41,7 +41,8 @@ if __name__ == "__main__":
R_harmonic = ift.ComposedOperator([fft, R], default_spaces=[0, 0])
# Setting up the noise covariance and drawing a random noise realization
N = ift.DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True)
ndiag = ift.Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
N = ift.DiagonalOperator(data_domain, ndiag)
noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(mock_signal) + noise #|\label{code:wf_mock_data}|
......@@ -57,7 +58,7 @@ if __name__ == "__main__":
proby = Proby(signal_space, probe_count=800)
proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z)))) #|\label{code:wf_variance_fft_wrap}|
sm = ift.SmoothingOperator.make(signal_space, sigma=0.03)
sm = ift.FFTSmoothingOperator(signal_space, sigma=0.03)
variance = ift.sqrt(sm(proby.diagonal.weight(-1))) #|\label{code:wf_variance_weighting}|
repo = Repository('repo_800.h5')
......
......@@ -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)
......@@ -54,9 +54,8 @@ if __name__ == "__main__":
data_domain = R.target[0]
R_harmonic = ComposedOperator([fft, R], default_spaces=[0, 0])
N = DiagonalOperator(data_domain,
diagonal=mock_signal.var()/signal_to_noise,
bare=True)
ndiag = Field(data_domain, mock_signal.var()/signal_to_noise).weight(1)
N = DiagonalOperator(data_domain, ndiag)
noise = Field.from_random(domain=data_domain,
random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise),
......
import d2o
from nifty import *
import plotly.offline as pl
......@@ -8,7 +11,7 @@ from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
np.random.seed(42)
d2o.random.seed(42)
class AdjointFFTResponse(LinearOperator):
......@@ -40,7 +43,8 @@ class AdjointFFTResponse(LinearOperator):
if __name__ == "__main__":
distribution_strategy = 'not'
nifty_configuration['default_distribution_strategy'] = 'fftw'
nifty_configuration['harmonic_rg_base'] = 'real'
# Set up position space
s_space = RGSpace([128, 128])
......@@ -51,18 +55,16 @@ if __name__ == "__main__":
h_space = fft.target[0]
# Setting up power space
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
p_space = PowerSpace(h_space)
# Choosing the prior correlation structure and defining
# correlation operator
p_spec = (lambda k: (42 / (k + 1) ** 3))
S = create_power_operator(h_space, power_spectrum=p_spec,
distribution_strategy=distribution_strategy)
S = create_power_operator(h_space, power_spectrum=p_spec)
# Drawing a sample sh from the prior distribution in harmonic space
sp = Field(p_space, val=p_spec,
distribution_strategy=distribution_strategy)
sp = Field(p_space, val=p_spec)
sh = sp.power_synthesize(real_signal=True)
ss = fft.adjoint_times(sh)
......@@ -74,7 +76,8 @@ if __name__ == "__main__":
# Adding a harmonic transformation to the instrument
R = AdjointFFTResponse(fft, Instrument)
signal_to_noise = 1.
N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
ndiag = Field(s_space, val=ss.var()/signal_to_noise).weight()
N = DiagonalOperator(s_space, diagonal=ndiag)
n = Field.from_random(domain=s_space,
random_type='normal',
std=ss.std()/np.sqrt(signal_to_noise),
......@@ -94,9 +97,17 @@ if __name__ == "__main__":
# iteration_limit=50,
# callback=convergence_measure)
minimizer = RelaxedNewton(convergence_tolerance=0,
iteration_limit=1,
callback=convergence_measure)
controller = GradientNormController(iteration_limit=50,
callback=convergence_measure)
minimizer = VL_BFGS(controller=controller)
controller = GradientNormController(iteration_limit=1,
callback=convergence_measure)
minimizer = RelaxedNewton(controller=controller)
#
# minimizer = VL_BFGS(convergence_tolerance=0,
# iteration_limit=50,
......@@ -104,9 +115,6 @@ if __name__ == "__main__":
# max_history_length=3)
#
inverter = ConjugateGradient(convergence_level=3,
convergence_tolerance=1e-5,
preconditioner=None)
# Setting starting position
m0 = Field(h_space, val=.0)
......@@ -115,15 +123,22 @@ if __name__ == "__main__":
D0 = energy.curvature
# Solving the problem analytically
m0 = D0.inverse_times(j)
sample_variance = Field(sh.domain, val=0.)
sample_mean = Field(sh.domain, val=0.)
# sampling the uncertainty map
n_samples = 10
for i in range(n_samples):
sample = fft(sugar.generate_posterior_sample(0., D0))
sample_variance += sample**2
sample_mean += sample
variance = (sample_variance - sample_mean**2)/n_samples
# m0 = D0.inverse_times(j)
m = minimizer(energy)[0].position
plotter = plotting.RG2DPlotter()
plotter(ss, path='signal.html')
plotter(fft.inverse_times(m), path='m.html')
# sample_variance = Field(sh.domain, val=0.)
# sample_mean = Field(sh.domain, val=0.)
# # sampling the uncertainty map
# n_samples = 10
# for i in range(n_samples):
# sample = fft(sugar.generate_posterior_sample(0., D0))
# sample_variance += sample**2
# sample_mean += sample
# variance = (sample_variance - sample_mean**2)/n_samples
......@@ -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)
......
......@@ -17,5 +17,6 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from .energy import Energy
from .quadratic_energy import QuadraticEnergy
from .line_energy import LineEnergy
from .memoization import memo
......@@ -17,12 +17,14 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from ..nifty_meta import NiftyMeta
from .memoization import memo
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
......@@ -41,7 +43,7 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {})))
value : np.float
The value of the energy functional at given `position`.
gradient : Field
The gradient at given `position` in parameter direction.
The gradient at given `position`.
curvature : LinearOperator, callable
A positive semi-definite operator or function describing the curvature
of the potential at the given `position`.
......@@ -64,12 +66,18 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {})))
"""
def __init__(self, position):
def __init__(self, position, gradient=None, curvature=None):
super(Energy, self).__init__()
self._cache = {}
self._position = position.copy()
def at(self, position):
self._cache = {}
if gradient is not None:
self._cache['gradient'] = gradient
if curvature is not None:
self._cache['curvature'] = curvature
def at(self, position, gradient=None, curvature=None):
""" Initializes and returns a new Energy object at the new position.
Parameters
......@@ -83,10 +91,12 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {})))
Energy object at new position.
"""
return self.__class__(position)
return self.__class__(position,
gradient=gradient,
curvature=curvature)
@property
@memo
def position(self):
"""
The Field location in parameter space where value, gradient and
......@@ -97,6 +107,7 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {})))
return self._position
@property
@memo
def value(self):
"""
The value of the energy functional at given `position`.
......@@ -106,15 +117,37 @@ class Energy(with_metaclass(NiftyMeta, type('NewBase', (Loggable, object), {})))
raise NotImplementedError
@property
@memo
def gradient(self):
"""
The gradient at given `position` in parameter direction.
The gradient at given `position`.
"""
raise NotImplementedError
@property
@memo
def gradient_norm(self):
"""
The length of the gradient at given `position`.
"""
return self.gradient.norm()
@property
@memo
def gradient_infnorm(self):
"""
The infinity norm of the gradient at given `position`.
"""
return abs(self.gradient).max()
@property
@memo
def curvature(self):
"""
A positive semi-definite operator or function describing the curvature
......
......@@ -18,8 +18,13 @@
from __future__ import print_function
from keepers import Loggable
from ..nifty_meta import NiftyMeta
from future.utils import with_metaclass
class LineEnergy(object):
class LineEnergy((with_metaclass(NiftyMeta,
type('NewBase', (Loggable, object), {})))):
""" Evaluates an underlying Energy along a certain line direction.
Given an Energy class and a line direction, its position is parametrized by
......@@ -74,7 +79,7 @@ class LineEnergy(object):
self._line_position = float(line_position)
self._line_direction = line_direction
if self._line_position==float(offset):
if self._line_position == float(offset):
self.energy = energy
else:
pos = energy.position \
......@@ -116,6 +121,6 @@ class LineEnergy(object):
def directional_derivative(self):
res = self.energy.gradient.vdot(self.line_direction)
if abs(res.imag) / max(abs(res.real), 1.) > 1e-12:
print ("directional derivative has non-negligible "
"imaginary part:", res)
self.logger.warn("directional derivative has non-negligible "
"imaginary part:", res)
return res.real
......@@ -18,9 +18,11 @@
def memo(f):
name = id(f)
name = f.__name__
def wrapped_f(self):
if not hasattr(self, "_cache"):
self._cache = {}
try:
return self._cache[name]
except KeyError:
......
from .energy import Energy
from .memoization import memo
class QuadraticEnergy(Energy):
"""The Energy for a quadratic form.
The most important aspect of this energy is that its curvature must be
position-independent.
"""
def __init__(self, position, A, b, gradient=None, curvature=None):
super(QuadraticEnergy, self).__init__(position=position,
gradient=gradient,
curvature=curvature)
self._A = A
self._b = b
if gradient is not None:
self._Ax = gradient + self._b
else:
self._Ax = self._A(self.position)
def at(self, position, gradient=None, curvature=None):
return self.__class__(position=position, A=self._A, b=self._b,
gradient=gradient, curvature=curvature)
@property
@memo
def value(self):
return 0.5*self.position.vdot(self._Ax) - self._b.vdot(self.position)
@property
@memo
def gradient(self):
return self._Ax - self._b
@property
def curvature(self):
return self._A
......@@ -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.")
......@@ -827,7 +814,7 @@ class Field(Loggable, Versionable, object):
try:
return int(reduce(lambda x, y: x * y, dim_tuple))
except TypeError:
return 0
return 1
@property
def dof(self):
......@@ -866,9 +853,9 @@ class Field(Loggable, Versionable, object):
def imag(self):
""" The imaginary part of the field (data is not copied).
"""
real_part = self.val.imag
result = self.copy_empty(dtype=real_part.dtype)
result.set_val(new_val=real_part, copy=False)
imag_part = self.val.imag
result = self.copy_empty(dtype=imag_part.dtype)
result.set_val(new_val=imag_part, copy=False)
return result
# ---Special unary/binary operations---
......@@ -1085,7 +1072,7 @@ class Field(Loggable, Versionable, object):
new_field.set_val(new_val=new_val, copy=False)
return new_field
def vdot(self, x=None, spaces=None, bare=False):
def vdot(self, x=None, spaces=None):
""" Computes the volume-factor-aware dot product of 'self' with x.
Parameters
......@@ -1097,9 +1084,6 @@ class Field(Loggable, Versionable, object):
If the domain of `self` and `x` are not the same, `spaces` specfies
the mapping.
bare : boolean
If true, no volume factors will be included in the computation.
Returns
-------
out : float, complex
......@@ -1110,10 +1094,7 @@ class Field(Loggable, Versionable, object):
"the NIFTy field class")
# Compute the dot respecting the fact of discrete/continuous spaces
if bare:
y = self
else:
y = self.weight(power=1)
y = self.weight(power=1)
if spaces is None:
x_val = x.get_val(copy=False)
......@@ -1603,7 +1584,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
......
from ...operators.endomorphic_operator import EndomorphicOperator
from ...operators.invertible_operator_mixin import InvertibleOperatorMixin
from ...operators.diagonal_operator import DiagonalOperator
from ...minimization import ConjugateGradient
class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
......@@ -21,34 +22,16 @@ class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
# ---Overwritten properties and methods---
def __init__(self, theta, T, inverter=None, preconditioner=None, **kwargs):
def __init__(self, theta, T, inverter=None, **kwargs):
self.theta = DiagonalOperator(theta.domain, diagonal=theta)
self.theta = DiagonalOperator(theta.domain, diagonal=theta, copy=False)
self.T = T
if preconditioner is None:
preconditioner = self.theta.inverse_times
if inverter is None:
inverter = ConjugateGradient(
preconditioner=self.theta.inverse_times)
self._domain = self.theta.domain
super(CriticalPowerCurvature, self).__init__(
inverter=inverter,
preconditioner=preconditioner,
**kwargs)
def _add_attributes_to_copy(self, copy, **kwargs):
copy._domain = self._domain
if 'theta' in kwargs:
theta = kwargs['theta']
copy.theta = DiagonalOperator(theta.domain, diagonal=theta)
else:
copy.theta = self.theta.copy()
if 'T' in kwargs:
copy.T = kwargs['T']
else:
copy.T = self.T