Commit 943a3944 authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'better_minimizers' into 'nightly'

Better minimizers

See merge request !191
parents bf3848ee edf74c43
......@@ -54,8 +54,8 @@ if __name__ == "__main__":
data_domain = R.target[0]
R_harmonic = ComposedOperator([fft, R], default_spaces=[0, 0])
ndiag = Field(data_domain,mock_signal.var()/signal_to_noise).weight(1)
N = DiagonalOperator(data_domain,ndiag)
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)
......@@ -95,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,
......@@ -105,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)
......@@ -116,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
......@@ -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,6 +17,7 @@
# 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
......@@ -42,7 +43,7 @@ class Energy(with_metaclass(NiftyMeta,
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`.
......@@ -65,12 +66,18 @@ class Energy(with_metaclass(NiftyMeta,
"""
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
......@@ -84,10 +91,12 @@ class Energy(with_metaclass(NiftyMeta,
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
......@@ -98,6 +107,7 @@ class Energy(with_metaclass(NiftyMeta,
return self._position
@property
@memo
def value(self):
"""
The value of the energy functional at given `position`.
......@@ -107,15 +117,37 @@ class Energy(with_metaclass(NiftyMeta,
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
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
copy = super(CriticalPowerCurvature,
self)._add_attributes_to_copy(copy, **kwargs)
return copy
inverter=inverter, **kwargs)
def _times(self, x, spaces):
return self.T(x) + self.theta(x)
......
from ...energies.energy import Energy
from ...operators.smoothness_operator import SmoothnessOperator
from ...operators.diagonal_operator import DiagonalOperator
from . import CriticalPowerCurvature
from ...energies.memoization import memo
from ...minimization import ConjugateGradient
from ...sugar import generate_posterior_sample
from ... import Field, exp
......@@ -55,8 +57,10 @@ class CriticalPowerEnergy(Energy):
def __init__(self, position, m, D=None, alpha=1.0, q=0.,
smoothness_prior=0., logarithmic=True, samples=3, w=None,
old_curvature=None):
super(CriticalPowerEnergy, self).__init__(position=position)
inverter=None, gradient=None, curvature=None):
super(CriticalPowerEnergy, self).__init__(position=position,
gradient=gradient,
curvature=curvature)
self.m = m
self.D = D
self.samples = samples
......@@ -67,8 +71,16 @@ class CriticalPowerEnergy(Energy):
logarithmic=logarithmic)
self.rho = self.position.domain[0].rho
self._w = w if w is not None else None
self._old_curvature = old_curvature
self._curvature = None
if inverter is None:
preconditioner = DiagonalOperator(self._theta.domain,
diagonal=self._theta.weight(-1),
copy=False)
inverter = ConjugateGradient(preconditioner=preconditioner)
self._inverter = inverter
@property
def inverter(self):
return self._inverter
# ---Mandatory properties and methods---
......@@ -77,7 +89,7 @@ class CriticalPowerEnergy(Energy):
q=self.q, smoothness_prior=self.smoothness_prior,
logarithmic=self.logarithmic,
w=self.w, samples=self.samples,
old_curvature=self._curvature)
inverter=self.inverter)
@property
@memo
......@@ -97,15 +109,10 @@ class CriticalPowerEnergy(Energy):
return gradient
@property
@memo
def curvature(self):
if self._curvature is None:
if self._old_curvature is None:
self._curvature = CriticalPowerCurvature(
theta=self._theta.weight(-1), T=self.T)
else:
self._curvature = self._old_curvature.copy(
theta=self._theta.weight(-1), T=self.T)
return self._curvature
return CriticalPowerCurvature(theta=self._theta.weight(-1), T=self.T,
inverter=self.inverter)
# ---Added properties and methods---
......
......@@ -3,6 +3,7 @@ from ...operators import EndomorphicOperator,\
from ...energies.memoization import memo
from ...basic_arithmetics import clipped_exp
from ...sugar import create_composed_fft_operator
from ...minimization import ConjugateGradient
class LogNormalWienerFilterCurvature(InvertibleOperatorMixin,
......@@ -25,17 +26,16 @@ class LogNormalWienerFilterCurvature(InvertibleOperatorMixin,
"""
def __init__(self, R, N, S, d, position, inverter=None,
preconditioner=None, fft4exp=None, offset=None, **kwargs):
self._cache = {}
def __init__(self, R, N, S, d, position, inverter=None, fft4exp=None,
offset=None, **kwargs):
self.R = R
self.N = N
self.S = S
self.d = d
if inverter is None:
inverter = ConjugateGradient(preconditioner=self.S.times)
self.position = position
self.offset = offset
if preconditioner is None:
preconditioner = self.S.times
self._domain = self.S.domain
if fft4exp is None:
......@@ -45,27 +45,7 @@ class LogNormalWienerFilterCurvature(InvertibleOperatorMixin,
self._fft = fft4exp
super(LogNormalWienerFilterCurvature, self).__init__(
inverter=inverter,
preconditioner=preconditioner,
**kwargs)
def _add_attributes_to_copy(self, copy, **kwargs):
copy._cache = {}
copy._domain = self._domain
copy.R = self.R.copy()
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:
copy.position = self.position.copy()
copy._fft = self._fft
copy = super(LogNormalWienerFilterCurvature,
self)._add_attributes_to_copy(copy, **kwargs)
return copy
inverter=inverter, **kwargs)
@property
def domain(self):
......
from ...energies.energy import Energy
from ...energies.memoization import memo
from ...minimization import ConjugateGradient
from . import LogNormalWienerFilterCurvature
from ...sugar import create_composed_fft_operator
......@@ -24,29 +25,35 @@ class LogNormalWienerFilterEnergy(Energy):
The prior signal covariance in harmonic space.
"""
def __init__(self, position, d, R, N, S, fft4exp=None, old_curvature=None,
offset=None):
super(LogNormalWienerFilterEnergy, self).__init__(position=position)
def __init__(self, position, d, R, N, S, offset=None, fft4exp=None,
inverter=None, gradient=None, curvature=None):
super(LogNormalWienerFilterEnergy, self).__init__(position=position,
gradient=gradient,
curvature=curvature)
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,
all_to='position')
else:
self._fft = fft4exp
self._old_curvature = old_curvature
self._curvature = None
if inverter is None:
inverter = ConjugateGradient(preconditioner=self.S.times)
self._inverter = inverter
@property
def inverter(self):
return self._inverter
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,
offset=self.offset)
def at(self, position, gradient=None, curvature=None):
return self.__class__(position=position,
d=self.d, R=self.R, N=self.N, S=self.S,
offset=self.offset, fft4exp=self._fft,
gradient=gradient, curvature=curvature)
@property
@memo
......@@ -60,21 +67,16 @@ class LogNormalWienerFilterEnergy(Energy):
return self._Sp + self._exppRNRexppd
@property
@memo
def curvature(self):
if self._curvature is None:
if self._old_curvature is None:
self._curvature = LogNormalWienerFilterCurvature(
R=self.R,
N=self.N,
S=self.S,
d=self.d,
position=self.position,
fft4exp=self._fft,
offset=self.offset)
else:
self._curvature = \
self._old_curvature.copy(position=self.position)
return self._curvature
return LogNormalWienerFilterCurvature(R=self.R,
N=self.N,
S=self.S,
d=self.d,
position=self.position,
fft4exp=self._fft,
offset=self.offset,
inverter=self.inverter)
@property
def _expp(self):
......
from ...operators import EndomorphicOperator,\
InvertibleOperatorMixin
from ...minimization import ConjugateGradient
class WienerFilterCurvature(InvertibleOperatorMixin, EndomorphicOperator):
......@@ -22,27 +23,16 @@ class WienerFilterCurvature(InvertibleOperatorMixin, EndomorphicOperator):
"""
def __init__(self, R, N, S, inverter=None, preconditioner=None, **kwargs):
def __init__(self, R, N, S, inverter=None, **kwargs):
self.R = R
self.N = N
self.S = S
if preconditioner is None:
preconditioner = self.S.times