Commit 747e2082 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

cleanup round #1

parent 5d60195f
...@@ -60,8 +60,7 @@ from .minimization.minimizer import Minimizer ...@@ -60,8 +60,7 @@ from .minimization.minimizer import Minimizer
from .minimization.conjugate_gradient import ConjugateGradient from .minimization.conjugate_gradient import ConjugateGradient
from .minimization.nonlinear_cg import NonlinearCG from .minimization.nonlinear_cg import NonlinearCG
from .minimization.descent_minimizers import ( from .minimization.descent_minimizers import (
DescentMinimizer, SteepestDescent, VL_BFGS, L_BFGS, RelaxedNewton, DescentMinimizer, SteepestDescent, VL_BFGS, L_BFGS, NewtonCG)
NewtonCG)
from .minimization.scipy_minimizer import (ScipyMinimizer, L_BFGS_B, ScipyCG) from .minimization.scipy_minimizer import (ScipyMinimizer, L_BFGS_B, ScipyCG)
from .minimization.energy import Energy from .minimization.energy import Energy
from .minimization.quadratic_energy import QuadraticEnergy from .minimization.quadratic_energy import QuadraticEnergy
......
...@@ -75,7 +75,7 @@ class ConjugateGradient(Minimizer): ...@@ -75,7 +75,7 @@ class ConjugateGradient(Minimizer):
return energy, controller.CONVERGED return energy, controller.CONVERGED
while True: while True:
q = energy.metric(d) q = energy.apply_metric(d)
ddotq = d.vdot(q).real ddotq = d.vdot(q).real
if ddotq == 0.: if ddotq == 0.:
logger.error("Error: ConjugateGradient: ddotq==0.") logger.error("Error: ConjugateGradient: ddotq==0.")
......
...@@ -135,24 +135,6 @@ class SteepestDescent(DescentMinimizer): ...@@ -135,24 +135,6 @@ class SteepestDescent(DescentMinimizer):
return -energy.gradient return -energy.gradient
class RelaxedNewton(DescentMinimizer):
""" Calculates the descent direction according to a Newton scheme.
The descent direction is determined by weighting the gradient at the
current parameter position with the inverse local metric.
"""
def __init__(self, controller, line_searcher=None):
if line_searcher is None:
line_searcher = LineSearchStrongWolfe(
preferred_initial_step_size=1.)
super(RelaxedNewton, self).__init__(controller=controller,
line_searcher=line_searcher)
def get_descent_direction(self, energy):
return -energy.metric.inverse_times(energy.gradient)
class NewtonCG(DescentMinimizer): class NewtonCG(DescentMinimizer):
""" Calculates the descent direction according to a Newton-CG scheme. """ Calculates the descent direction according to a Newton-CG scheme.
...@@ -180,7 +162,7 @@ class NewtonCG(DescentMinimizer): ...@@ -180,7 +162,7 @@ class NewtonCG(DescentMinimizer):
while True: while True:
if abs(ri).sum() <= termcond: if abs(ri).sum() <= termcond:
return xsupi return xsupi
Ap = energy.metric(psupi) Ap = energy.apply_metric(psupi)
# check curvature # check curvature
curv = psupi.vdot(Ap) curv = psupi.vdot(Ap)
if 0 <= curv <= 3*float64eps: if 0 <= curv <= 3*float64eps:
......
...@@ -100,12 +100,17 @@ class Energy(NiftyMetaBase()): ...@@ -100,12 +100,17 @@ class Energy(NiftyMetaBase()):
self._gradnorm = self.gradient.norm() self._gradnorm = self.gradient.norm()
return self._gradnorm return self._gradnorm
@property def apply_metric(self, x):
def metric(self):
""" """
LinearOperator : implicitly defined metric. Parameters
A positive semi-definite operator or function describing the ----------
metric of the potential at the given `position`. x: Field/MultiField
Argument for the metric operator
Returns
-------
Field/MultiField:
Output of the metric operator
""" """
raise NotImplementedError raise NotImplementedError
......
...@@ -8,18 +8,23 @@ from ..operators.scaling_operator import ScalingOperator ...@@ -8,18 +8,23 @@ from ..operators.scaling_operator import ScalingOperator
class EnergyAdapter(Energy): class EnergyAdapter(Energy):
def __init__(self, position, op, controller=None, preconditioner=None, def __init__(self, position, op, constants=[]):
constants=[]):
super(EnergyAdapter, self).__init__(position) super(EnergyAdapter, self).__init__(position)
self._op = op self._op = op
self._val = self._grad = self._metric = None
self._controller = controller
self._preconditioner = preconditioner
self._constants = constants self._constants = constants
if len(self._constants) == 0:
tmp = self._op(Linearization.make_var(self._position))
else:
ops = [ScalingOperator(0. if key in self._constants else 1., dom)
for key, dom in self._position.domain.items()]
bdop = BlockDiagonalOperator(self._position.domain, tuple(ops))
tmp = self._op(Linearization(self._position, bdop))
self._val = tmp.val.local_data[()]
self._grad = tmp.gradient
self._metric = tmp._metric
def at(self, position): def at(self, position):
return EnergyAdapter(position, self._op, self._controller, return EnergyAdapter(position, self._op, self._constants)
self._preconditioner, self._constants)
def _fill_all(self): def _fill_all(self):
if len(self._constants) == 0: if len(self._constants) == 0:
...@@ -31,35 +36,15 @@ class EnergyAdapter(Energy): ...@@ -31,35 +36,15 @@ class EnergyAdapter(Energy):
tmp = self._op(Linearization(self._position, bdop)) tmp = self._op(Linearization(self._position, bdop))
self._val = tmp.val.local_data[()] self._val = tmp.val.local_data[()]
self._grad = tmp.gradient self._grad = tmp.gradient
if self._controller is not None: self._metric = tmp._metric
from ..operators.linear_operator import LinearOperator
from ..operators.inversion_enabler import InversionEnabler
if self._preconditioner is None:
precond = None
elif isinstance(self._preconditioner, LinearOperator):
precond = self._preconditioner
elif isinstance(self._preconditioner, Energy):
precond = self._preconditioner.at(self._position).metric
self._metric = InversionEnabler(tmp._metric, self._controller,
precond)
else:
self._metric = tmp._metric
@property @property
def value(self): def value(self):
if self._val is None:
self._val = self._op(self._position).local_data[()]
return self._val return self._val
@property @property
def gradient(self): def gradient(self):
if self._grad is None:
self._fill_all()
return self._grad return self._grad
@property def apply_metric(self, x):
def metric(self): return self._metric(x)
if self._metric is None:
self._fill_all()
return self._metric
...@@ -74,6 +74,5 @@ class QuadraticEnergy(Energy): ...@@ -74,6 +74,5 @@ class QuadraticEnergy(Energy):
def gradient(self): def gradient(self):
return self._grad return self._grad
@property def apply_metric(self, x):
def metric(self): return self._A(x)
return self._A
...@@ -93,7 +93,7 @@ class _MinHelper(object): ...@@ -93,7 +93,7 @@ class _MinHelper(object):
def hessp(self, x, p): def hessp(self, x, p):
self._update(x) self._update(x)
res = self._energy.metric(_toField(p, self._energy.position)) res = self._energy.apply_metric(_toField(p, self._energy.position))
return _toArray_rw(res) return _toArray_rw(res)
......
...@@ -38,7 +38,6 @@ minimizers = ['ift.VL_BFGS(IC)', ...@@ -38,7 +38,6 @@ minimizers = ['ift.VL_BFGS(IC)',
'ift.L_BFGS(IC)', 'ift.L_BFGS(IC)',
'ift.NewtonCG(IC)'] 'ift.NewtonCG(IC)']
newton_minimizers = ['ift.RelaxedNewton(IC)']
quadratic_only_minimizers = ['ift.ConjugateGradient(IC)', quadratic_only_minimizers = ['ift.ConjugateGradient(IC)',
'ift.ScipyCG(tol=1e-5, maxiter=300)'] 'ift.ScipyCG(tol=1e-5, maxiter=300)']
slow_minimizers = ['ift.SteepestDescent(IC)'] slow_minimizers = ['ift.SteepestDescent(IC)']
...@@ -46,8 +45,8 @@ slow_minimizers = ['ift.SteepestDescent(IC)'] ...@@ -46,8 +45,8 @@ slow_minimizers = ['ift.SteepestDescent(IC)']
class Test_Minimizers(unittest.TestCase): class Test_Minimizers(unittest.TestCase):
@expand(product(minimizers + newton_minimizers + @expand(product(minimizers + quadratic_only_minimizers + slow_minimizers,
quadratic_only_minimizers + slow_minimizers, spaces)) spaces))
def test_quadratic_minimization(self, minimizer, space): def test_quadratic_minimization(self, minimizer, space):
np.random.seed(42) np.random.seed(42)
starting_point = ift.Field.from_random('normal', domain=space)*10 starting_point = ift.Field.from_random('normal', domain=space)*10
...@@ -70,7 +69,7 @@ class Test_Minimizers(unittest.TestCase): ...@@ -70,7 +69,7 @@ class Test_Minimizers(unittest.TestCase):
1./covariance_diagonal.local_data, 1./covariance_diagonal.local_data,
rtol=1e-3, atol=1e-3) rtol=1e-3, atol=1e-3)
@expand(product(minimizers+newton_minimizers)) @expand(product(minimizers))
def test_rosenbrock(self, minimizer): def test_rosenbrock(self, minimizer):
try: try:
from scipy.optimize import rosen, rosen_der, rosen_hess_prod from scipy.optimize import rosen, rosen_der, rosen_hess_prod
...@@ -94,24 +93,11 @@ class Test_Minimizers(unittest.TestCase): ...@@ -94,24 +93,11 @@ class Test_Minimizers(unittest.TestCase):
out = ift.Field.from_global_data(space, rosen_der(inp)) out = ift.Field.from_global_data(space, rosen_der(inp))
return out return out
@property def apply_metric(self, x):
def metric(self): inp = x.to_global_data_rw()
class RBCurv(ift.EndomorphicOperator): pos = self._position.to_global_data_rw()
def __init__(self, loc): return ift.Field.from_global_data(
self._loc = loc.to_global_data_rw() space, rosen_hess_prod(pos, inp))
self._capability = self.TIMES
self._domain = space
def apply(self, x, mode):
self._check_input(x, mode)
inp = x.to_global_data_rw()
out = ift.Field.from_global_data(
space, rosen_hess_prod(self._loc.copy(), inp))
return out
t1 = ift.GradientNormController(tol_abs_gradnorm=1e-5,
iteration_limit=1000)
return ift.InversionEnabler(RBCurv(self._position), t1)
try: try:
minimizer = eval(minimizer) minimizer = eval(minimizer)
...@@ -145,12 +131,11 @@ class Test_Minimizers(unittest.TestCase): ...@@ -145,12 +131,11 @@ class Test_Minimizers(unittest.TestCase):
return ift.Field.full(self.position.domain, return ift.Field.full(self.position.domain,
2*x*np.exp(-(x**2))) 2*x*np.exp(-(x**2)))
@property def apply_metric(self, x):
def metric(self): p = self.position.to_global_data()[0]
x = self.position.to_global_data()[0] v = (2 - 4*p*p)*np.exp(-p**2)
v = (2 - 4*x*x)*np.exp(-x**2)
return ift.DiagonalOperator( return ift.DiagonalOperator(
ift.Field.full(self.position.domain, v)) ift.Field.full(self.position.domain, v))(x)
try: try:
minimizer = eval(minimizer) minimizer = eval(minimizer)
...@@ -164,7 +149,7 @@ class Test_Minimizers(unittest.TestCase): ...@@ -164,7 +149,7 @@ class Test_Minimizers(unittest.TestCase):
assert_allclose(energy.position.local_data, 0., assert_allclose(energy.position.local_data, 0.,
atol=1e-3) atol=1e-3)
@expand(product(minimizers+newton_minimizers+slow_minimizers)) @expand(product(minimizers+slow_minimizers))
def test_cosh(self, minimizer): def test_cosh(self, minimizer):
space = ift.UnstructuredDomain((1,)) space = ift.UnstructuredDomain((1,))
starting_point = ift.Field.full(space, 3.) starting_point = ift.Field.full(space, 3.)
...@@ -183,12 +168,11 @@ class Test_Minimizers(unittest.TestCase): ...@@ -183,12 +168,11 @@ class Test_Minimizers(unittest.TestCase):
x = self.position.to_global_data()[0] x = self.position.to_global_data()[0]
return ift.Field.full(self.position.domain, np.sinh(x)) return ift.Field.full(self.position.domain, np.sinh(x))
@property def apply_metric(self, x):
def metric(self): p = self.position.to_global_data()[0]
x = self.position.to_global_data()[0] v = np.cosh(p)
v = np.cosh(x)
return ift.DiagonalOperator( return ift.DiagonalOperator(
ift.Field.full(self.position.domain, v)) ift.Field.full(self.position.domain, v))(x)
try: try:
minimizer = eval(minimizer) minimizer = eval(minimizer)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment