Merge branch 'noise_energy' into 'NIFTy_4'

Noise energy

See merge request ift/NIFTy!226

nifty4/library/noise_energy.py

@@ -16,52 +16,27 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
-from .. import Field, exp
-from ..operators.diagonal_operator import DiagonalOperator
+from ..field import Field, exp
 from ..minimization.energy import Energy
-
-# TODO Take only residual_sample_list as argument
+from ..operators.diagonal_operator import DiagonalOperator
 
 
 class NoiseEnergy(Energy):
-    def __init__(self, position, d, xi, D, t, ht, Instrument,
-                 nonlinearity, alpha, q, Distributor, samples=3,
-                 xi_sample_list=None, inverter=None):
-        super(NoiseEnergy, self).__init__(position=position)
-        self.xi = xi
-        self.D = D
-        self.d = d
-        self.N = DiagonalOperator(diagonal=exp(self.position))
-        self.t = t
-        self.samples = samples
-        self.ht = ht
-        self.Instrument = Instrument
-        self.nonlinearity = nonlinearity
+    def __init__(self, position, alpha, q, res_sample_list):
+        super(NoiseEnergy, self).__init__(position)
 
+        self.N = DiagonalOperator(diagonal=exp(self.position))
         self.alpha = alpha
         self.q = q
-        self.Distributor = Distributor
-        self.power = self.Distributor(exp(0.5 * self.t))
-        if xi_sample_list is None:
-            if samples is None or samples == 0:
-                xi_sample_list = [xi]
-            else:
-                xi_sample_list = [D.draw_sample() + xi
-                                  for _ in range(samples)]
-        self.xi_sample_list = xi_sample_list
-        self.inverter = inverter
-
-        A = Distributor(exp(.5*self.t))
+        alpha_field = Field(self.position.domain, val=alpha)
+        q_field = Field(self.position.domain, val=q)
+        self.res_sample_list = res_sample_list
 
         self._gradient = None
-        for sample in self.xi_sample_list:
-            map_s = self.ht(A * sample)
-
-            residual = self.d - \
-                self.Instrument(self.nonlinearity(map_s))
-            lh = .5 * residual.vdot(self.N.inverse_times(residual))
-            grad = -.5 * self.N.inverse_times(residual.conjugate()*residual)
 
+        for s in self.res_sample_list:
+            lh = .5 * s.vdot(self.N.inverse_times(s))
+            grad = -.5 * self.N.inverse_times(s.conjugate()*s)
             if self._gradient is None:
                 self._value = lh
                 self._gradient = grad.copy()
@@ -69,20 +44,19 @@ class NoiseEnergy(Energy):
                 self._value += lh
                 self._gradient += grad
 
-        self._value *= 1. / len(self.xi_sample_list)
+        expmpos = exp(-position)
+        self._value *= 1./len(self.res_sample_list)
         self._value += .5 * self.position.sum()
-        self._value += (self.alpha - 1.).vdot(self.position) + \
-            self.q.vdot(exp(-self.position))
+        self._value += (alpha_field-1.).vdot(self.position) + \
+            q_field.vdot(expmpos)
 
-        self._gradient *= 1. / len(self.xi_sample_list)
-        self._gradient += (self.alpha-0.5) - self.q*(exp(-self.position))
+        self._gradient *= 1./len(self.res_sample_list)
+        self._gradient += (alpha_field-0.5) - q_field*expmpos
+        self._gradient.lock()
 
     def at(self, position):
-        return self.__class__(
-            position, self.d, self.xi, self.D, self.t, self.ht,
-            self.Instrument, self.nonlinearity, self.alpha, self.q,
-            self.Distributor, xi_sample_list=self.xi_sample_list,
-            samples=self.samples, inverter=self.inverter)
+        return self.__class__(position, self.alpha, self.q,
+                              self.res_sample_list)
 
     @property
     def value(self):

test/test_energies/test_noise.py

@@ -73,7 +73,7 @@ class Noise_Energy_Tests(unittest.TestCase):
         inverter = ift.ConjugateGradient(IC)
 
         S = ift.create_power_operator(hspace, power_spectrum=_flat_PS)
-        D = ift.library.NonlinearWienerFilterEnergy(
+        C = ift.library.NonlinearWienerFilterEnergy(
             position=xi,
             d=d,
             Instrument=R,
@@ -84,10 +84,10 @@ class Noise_Energy_Tests(unittest.TestCase):
             S=S,
             inverter=inverter).curvature
 
-        energy0 = ift.library.NoiseEnergy(
-            position=eta0, d=d, xi=xi, D=D, t=tau, Instrument=R,
-            alpha=alpha, q=q, Distributor=Dist, nonlinearity=f,
-            ht=ht, samples=10)
+        res_sample_list = [d - R(f(ht(C.draw_sample() + xi)))
+                           for _ in range(10)]
+
+        energy0 = ift.library.NoiseEnergy(eta0, alpha, q, res_sample_list)
         energy1 = energy0.at(eta1)
 
         a = (energy1.value - energy0.value) / eps