Merge branch 'energy_cleanup' into 'NIFTy_5'

Energy re-design

See merge request ift/nifty-dev!25

demos/getting_started_2.py

@@ -79,7 +79,8 @@ if __name__ == '__main__':
     minimizer = ift.RelaxedNewton(ic_newton)
 
     # Minimize the Hamiltonian
-    H = ift.Hamiltonian(likelihood, ic_cg)
+    H = ift.Hamiltonian(likelihood)
+    H = H.makeInvertible(ic_cg)
     H, convergence = minimizer(H)
 
     # Plot results

demos/getting_started_3.py

@@ -64,8 +64,7 @@ if __name__ == '__main__':
     minimizer = ift.RelaxedNewton(ic_newton)
 
     # build model Hamiltonian
-    H = ift.Hamiltonian(likelihood, ic_cg,
-                        iteration_controller_sampling=ic_sampling)
+    H = ift.Hamiltonian(likelihood, ic_sampling)
 
     INITIAL_POSITION = ift.from_random('normal', H.position.domain)
     position = INITIAL_POSITION
@@ -81,7 +80,8 @@ if __name__ == '__main__':
         samples = [H.curvature.draw_sample(from_inverse=True)
                    for _ in range(N_samples)]
 
-        KL = ift.SampledKullbachLeiblerDivergence(H, samples, ic_cg)
+        KL = ift.SampledKullbachLeiblerDivergence(H, samples)
+        KL = KL.makeInvertible(ic_cg)
         KL, convergence = minimizer(KL)
         position = KL.position
 

nifty5/energies/hamiltonian.py

@@ -19,26 +19,23 @@
 from ..library.gaussian_energy import GaussianEnergy
 from ..minimization.energy import Energy
 from ..models.variable import Variable
-from ..operators import InversionEnabler, SamplingEnabler
+from ..operators.sampling_enabler import SamplingEnabler
 from ..utilities import memo
 
 
 class Hamiltonian(Energy):
-    def __init__(self, lh, iteration_controller,
-                 iteration_controller_sampling=None):
+    def __init__(self, lh, iteration_controller_sampling=None):
         """
         lh: Likelihood (energy object)
         prior:
         """
         super(Hamiltonian, self).__init__(lh.position)
         self._lh = lh
-        self._ic = iteration_controller
         self._ic_samp = iteration_controller_sampling
         self._prior = GaussianEnergy(Variable(self.position))
-        self._precond = self._prior.curvature
 
     def at(self, position):
-        return self.__class__(self._lh.at(position), self._ic, self._ic_samp)
+        return self.__class__(self._lh.at(position), self._ic_samp)
 
     @property
     @memo
@@ -55,11 +52,10 @@ class Hamiltonian(Energy):
     def curvature(self):
         prior_curv = self._prior.curvature
         if self._ic_samp is None:
-            c = self._lh.curvature + prior_curv
+            return self._lh.curvature + prior_curv
         else:
-            c = SamplingEnabler(self._lh.curvature, prior_curv.inverse,
-                                self._ic_samp, prior_curv.inverse)
-        return InversionEnabler(c, self._ic, self._precond)
+            return SamplingEnabler(self._lh.curvature, prior_curv.inverse,
+                                   self._ic_samp, prior_curv.inverse)
 
     def __str__(self):
         res = 'Likelihood:\t{:.2E}\n'.format(self._lh.value)

nifty5/energies/kl.py

+from builtins import *
 from ..minimization.energy import Energy
-from ..operators.inversion_enabler import InversionEnabler
-from ..operators.scaling_operator import ScalingOperator
-from ..utilities import memo
+from ..utilities import memo, my_sum
 
 
 class SampledKullbachLeiblerDivergence(Energy):
-    def __init__(self, h, res_samples, iteration_controller):
+    def __init__(self, h, res_samples):
         """
         h: Hamiltonian
         N: Number of samples to be used
@@ -13,41 +12,27 @@ class SampledKullbachLeiblerDivergence(Energy):
         super(SampledKullbachLeiblerDivergence, self).__init__(h.position)
         self._h = h
         self._res_samples = res_samples
-        self._iteration_controller = iteration_controller
 
-        self._energy_list = []
-        for ss in res_samples:
-            e = h.at(self.position+ss)
-            self._energy_list.append(e)
+        self._energy_list = tuple(h.at(self.position+ss)
+                                  for ss in res_samples)
 
     def at(self, position):
-        return self.__class__(self._h.at(position), self._res_samples,
-                              self._iteration_controller)
+        return self.__class__(self._h.at(position), self._res_samples)
 
     @property
     @memo
     def value(self):
-        v = self._energy_list[0].value
-        for energy in self._energy_list[1:]:
-            v += energy.value
-        return v / len(self._energy_list)
+        return (my_sum(map(lambda v: v.value, self._energy_list)) /
+                len(self._energy_list))
 
     @property
     @memo
     def gradient(self):
-        g = self._energy_list[0].gradient
-        for energy in self._energy_list[1:]:
-            g += energy.gradient
-        return g / len(self._energy_list)
+        return (my_sum(map(lambda v: v.gradient, self._energy_list)) /
+                len(self._energy_list))
 
     @property
     @memo
     def curvature(self):
-        # MR FIXME: This looks a bit strange...
-        approx = self._energy_list[-1]._prior.curvature
-        curvature_list = [e.curvature for e in self._energy_list]
-        op = curvature_list[0]
-        for curv in curvature_list[1:]:
-            op = op + curv
-        op = op * ScalingOperator(1./len(curvature_list), op.domain)
-        return InversionEnabler(op, self._iteration_controller, approx)
+       return (my_sum(map(lambda v: v.curvature, self._energy_list)) *
+               (1./len(self._energy_list)))

nifty5/library/amplitude_model.py

 import numpy as np
 
-from ..domains import PowerSpace, UnstructuredDomain
+from ..domains.power_space import PowerSpace
+from ..domains.unstructured_domain import UnstructuredDomain
 from ..field import Field
-from ..multi import MultiField
+from ..multi.multi_field import MultiField
 from ..sugar import makeOp, sqrt
 
 
@@ -26,7 +27,7 @@ def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
                         eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
                         a = ceps_a,  k0 = ceps_k0
 
-    sm, sv : slope_mean = expected exponent of powerlaw (e.g. -4),
+    sm, sv : slope_mean = expected exponent of power law (e.g. -4),
                 slope_variance (default=1)
 
     im, iv : y-intercept_mean, y-intercept_variance  of power_slope

nifty5/minimization/descent_minimizer.py

@@ -80,12 +80,9 @@ class DescentMinimizer(Minimizer):
                 return energy, controller.CONVERGED
 
             # compute a step length that reduces energy.value sufficiently
-            try:
-                new_energy, success = self.line_searcher.perform_line_search(
-                    energy=energy, pk=self.get_descent_direction(energy),
-                    f_k_minus_1=f_k_minus_1)
-            except ValueError:
-                return energy, controller.ERROR
+            new_energy, success = self.line_searcher.perform_line_search(
+                energy=energy, pk=self.get_descent_direction(energy),
+                f_k_minus_1=f_k_minus_1)
             if not success:
                 self.reset()
 

nifty5/minimization/energy.py

@@ -129,6 +129,12 @@ class Energy(NiftyMetaBase()):
         """
         return None
 
+    def makeInvertible(self, controller, preconditioner=None):
+        from .iteration_controller import IterationController
+        if not isinstance(controller, IterationController):
+            raise TypeError
+        return CurvatureInversionEnabler(self, controller, preconditioner)
+
     def __mul__(self, factor):
         from .energy_sum import EnergySum
         if isinstance(factor, (float, int)):
@@ -153,3 +159,45 @@ class Energy(NiftyMetaBase()):
     def __neg__(self):
         from .energy_sum import EnergySum
         return EnergySum.make([self], [-1.])
+
+
+class CurvatureInversionEnabler(Energy):
+    def __init__(self, ene, controller, preconditioner):
+        super(CurvatureInversionEnabler, self).__init__(ene.position)
+        self._energy = ene
+        self._controller = controller
+        self._preconditioner = preconditioner
+
+    def at(self, position):
+        if self._position.isSubsetOf(position):
+            return self
+        return CurvatureInversionEnabler(
+            self._energy.at(position), self._controller, self._preconditioner)
+
+    @property
+    def position(self):
+        return self._energy.position
+
+    @property
+    def value(self):
+        return self._energy.value
+
+    @property
+    def gradient(self):
+        return self._energy.gradient
+
+    @property
+    def curvature(self):
+        from ..operators.linear_operator import LinearOperator
+        from ..operators.inversion_enabler import InversionEnabler
+        curv = self._energy.curvature
+        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).curvature
+        return InversionEnabler(curv, self._controller, precond)
+
+    def longest_step(self, dir):
+        return self._energy.longest_step(dir)

nifty5/minimization/energy_sum.py

@@ -16,7 +16,8 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
-from ..utilities import memo
+from builtins import *
+from ..utilities import memo, my_lincomb_simple, my_lincomb
 from .energy import Energy
 
 
@@ -55,26 +56,17 @@ class EnergySum(Energy):
     @property
     @memo
     def value(self):
-        res = self._energies[0].value * self._factors[0]
-        for e, f in zip(self._energies[1:], self._factors[1:]):
-            res += e.value * f
-        return res
+        return my_lincomb_simple(map(lambda v: v.value, self._energies),
+                                 self._factors)
 
     @property
     @memo
     def gradient(self):
-        res = self._energies[0].gradient.copy() if self._factors[0] == 1. \
-            else self._energies[0].gradient * self._factors[0]
-
-        for e, f in zip(self._energies[1:], self._factors[1:]):
-            res += e.gradient if f == 1. else f*e.gradient
-        return res.lock()
+        return my_lincomb(map(lambda v: v.gradient, self._energies),
+                          self._factors).lock()
 
     @property
     @memo
     def curvature(self):
-        res = self._energies[0].curvature if self._factors[0] == 1. \
-            else self._energies[0].curvature * self._factors[0]
-        for e, f in zip(self._energies[1:], self._factors[1:]):
-            res = res + (e.curvature if f == 1. else e.curvature*f)
-        return res
+        return my_lincomb(map(lambda v: v.curvature, self._energies),
+                          self._factors)

nifty5/operators/linear_operator.py

@@ -274,9 +274,9 @@ class LinearOperator(NiftyMetaBase()):
 
     def _check_mode(self, mode):
         if not self._validMode[mode]:
-            raise ValueError("invalid operator mode specified")
+            raise NotImplementedError("invalid operator mode specified")
         if mode & self.capability == 0:
-            raise ValueError("requested operator mode is not supported")
+            raise NotImplementedError("requested operator mode is not supported")
 
     def _check_input(self, x, mode):
         self._check_mode(mode)

nifty5/utilities.py

@@ -16,15 +16,36 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
-from builtins import next, range
+from builtins import *
 import numpy as np
 from itertools import product
 import abc
 from future.utils import with_metaclass
+from functools import reduce
 
 __all__ = ["get_slice_list", "safe_cast", "parse_spaces", "infer_space",
            "memo", "NiftyMetaBase", "fft_prep", "hartley", "my_fftn_r2c",
-           "my_fftn"]
+           "my_fftn", "my_sum", "my_lincomb_simple", "my_lincomb",
+           "my_product"]
+
+
+def my_sum(terms):
+    return reduce(lambda x, y: x+y, terms)
+
+
+def my_lincomb_simple(terms, factors):
+    terms2 = map(lambda v: v[0]*v[1], zip(terms, factors))
+    return my_sum(terms2)
+
+
+def my_lincomb(terms, factors):
+    terms2 = map(lambda v: v[0] if v[1] == 1. else v[0]*v[1],
+                 zip(terms, factors))
+    return my_sum(terms2)
+
+
+def my_product(iterable):
+    return reduce(lambda x, y: x*y, iterable)
 
 
 def get_slice_list(shape, axes):