Commit 4d421738 authored by Theo Steininger's avatar Theo Steininger
Browse files

Refactored minimizers. Removed LinearOperator.copy

parent 715b01f3
......@@ -66,12 +66,20 @@ 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):
if gradient is not None:
key = id(self.gradient)
self._cache[key] = gradient
if curvature is not None:
key = id(self.curvature)
self._cache[key] = curvature
def at(self, position, gradient=None, curvature=None):
""" Initializes and returns a new Energy object at the new position.
Parameters
......@@ -85,10 +93,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
......@@ -99,6 +109,7 @@ class Energy(with_metaclass(NiftyMeta,
return self._position
@property
@memo
def value(self):
"""
The value of the energy functional at given `position`.
......@@ -108,6 +119,7 @@ class Energy(with_metaclass(NiftyMeta,
raise NotImplementedError
@property
@memo
def gradient(self):
"""
The gradient at given `position`.
......@@ -137,6 +149,7 @@ class Energy(with_metaclass(NiftyMeta,
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
......@@ -8,21 +8,20 @@ class QuadraticEnergy(Energy):
position-independent.
"""
def __init__(self, position, A, b, grad=None):
super(QuadraticEnergy, self).__init__(position=position)
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 grad is not None:
self._Ax = grad + self._b
if gradient is not None:
self._Ax = gradient + self._b
else:
self._Ax = self._A(self.position)
def at(self, position):
return self.__class__(position=position, A=self._A, b=self._b)
def at_with_grad(self, position, grad):
def at(self, position, gradient=None, curvature=None):
return self.__class__(position=position, A=self._A, b=self._b,
grad=grad)
gradient=gradient, curvature=curvature)
@property
@memo
......
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,17 @@ class LogNormalWienerFilterCurvature(InvertibleOperatorMixin,
"""
def __init__(self, R, N, S, d, position, inverter=None,
preconditioner=None, fft4exp=None, offset=None, **kwargs):
def __init__(self, R, N, S, d, position, inverter=None, fft4exp=None,
offset=None, **kwargs):
self._cache = {}
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 +46,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
if inverter is None:
inverter = ConjugateGradient(preconditioner=self.S.times)
self._domain = self.S.domain
super(WienerFilterCurvature, self).__init__(
inverter=inverter,
preconditioner=preconditioner,
**kwargs)
def _add_attributes_to_copy(self, copy, **kwargs):
copy._domain = self._domain
copy.R = self.R.copy()
copy.N = self.N.copy()
copy.S = self.S.copy()
copy = super(WienerFilterCurvature, self)._add_attributes_to_copy(
copy, **kwargs)
return copy
super(WienerFilterCurvature, self).__init__(inverter=inverter,
**kwargs)
@property
def domain(self):
......
from ...energies.energy import Energy
from ...energies.memoization import memo
from ...minimization import ConjugateGradient
from . import WienerFilterCurvature
......@@ -23,17 +25,26 @@ class WienerFilterEnergy(Energy):
The prior signal covariance in harmonic space.
"""
def __init__(self, position, d, R, N, S, old_curvature=None):
super(WienerFilterEnergy, self).__init__(position=position)
def __init__(self, position, d, R, N, S, inverter=None,
gradient=None, curvature=None):
super(WienerFilterEnergy, self).__init__(position=position,
gradient=gradient,
curvature=curvature)
self.d = d
self.R = R
self.N = N
self.S = S
self._curvature = old_curvature
if inverter is None:
inverter = ConjugateGradient(preconditioner=self.S.times)
self._inverter = inverter
@property
def inverter(self):
return self._inverter
def at(self, position):
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, old_curvature=self.curvature)
S=self.S, inverter=self.inverter)
@property
@memo
......@@ -46,12 +57,10 @@ class WienerFilterEnergy(Energy):
return self._Dx - self._j
@property
@memo
def curvature(self):
if self._curvature is None:
self._curvature = WienerFilterCurvature(R=self.R,
N=self.N,
S=self.S)
return self._curvature
return WienerFilterCurvature(R=self.R, N=self.N, S=self.S,
inverter=self.inverter)
@property
@memo
......
......@@ -17,8 +17,7 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from .line_searching import *
from .iteration_controller import IterationController
from .default_iteration_controller import DefaultIterationController
from .iteration_controlling import *
from .minimizer import Minimizer
from .conjugate_gradient import ConjugateGradient
from .nonlinear_cg import NonlinearCG
......
......@@ -17,9 +17,11 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from __future__ import division
import numpy as np
from .minimizer import Minimizer
from .iteration_controlling import GradientNormController
class ConjugateGradient(Minimizer):
......@@ -43,7 +45,9 @@ class ConjugateGradient(Minimizer):
"""
def __init__(self, controller, preconditioner=None):
def __init__(self,
controller=GradientNormController(iteration_limit=100),
preconditioner=None):
self._preconditioner = preconditioner
self._controller = controller
......@@ -66,50 +70,44 @@ class ConjugateGradient(Minimizer):
"""
controller = self._controller
status = controller.start(energy)
if status != controller.CONTINUE:
return energy, status
status = controller.reset(energy)
r = -energy.gradient
if self._preconditioner is not None:
d = self._preconditioner(r)
else:
d = r.copy()
previous_gamma = (r.vdot(d)).real
if previous_gamma == 0:
return energy, controller.CONVERGED
previous_gamma = np.inf
d = r.copy_empty()
d.val[:] = 0.
while True:
q = energy.curvature(d)
ddotq = d.vdot(q).real
if ddotq==0.:
self.logger.error("Alpha became infinite! Stopping.")
return energy, controller.ERROR
alpha = previous_gamma/ddotq
if alpha < 0:
self.logger.warn("Positive definiteness of A violated!")
return energy, controller.ERROR
r -= q * alpha
energy = energy.at_with_grad(energy.position+d*alpha,-r)
status = self._controller.check(energy)
if status != controller.CONTINUE:
return energy, status
if self._preconditioner is not None:
s = self._preconditioner(r)
else:
s = r
gamma = r.vdot(s).real
if gamma < 0:
self.logger.warn(
"Positive definiteness of preconditioner violated!")
if gamma == 0:
self.logger.info("Gamma == 0. Stopping.")
return energy, controller.CONVERGED
d = s + d * max(0, gamma/previous_gamma)
previous_gamma = gamma
status = self._controller.check(energy)
if status != controller.CONTINUE:
return energy, status
q = energy.curvature(d)
ddotq = d.vdot(q).real
if ddotq == 0.:
self.logger.error("Alpha became infinite! Stopping.")
return energy, controller.ERROR
alpha = previous_gamma/ddotq
if alpha < 0:
self.logger.error(
"Positive definiteness of A violated! Stopping.")
return energy, controller.ERROR