implemented WienerFilterEnergy and WienerFilterCurvature, updated wiener_filter_hamiltonian

demos/wiener_filter_hamiltonian.py

@@ -10,45 +10,20 @@ rank = comm.rank
 
 np.random.seed(42)
 
-class WienerFilterEnergy(Energy):
-    def __init__(self, position, D, j):
-        # in principle not necessary, but useful in order to make the signature
-        # explicit
-        super(WienerFilterEnergy, self).__init__(position)
-        self.D = D
-        self.j = j
-
-    def at(self, position):
-        return self.__class__(position, D=self.D, j=self.j)
-
-    @property
-    def value(self):
-        D_inv_x = self.D_inverse_x()
-        H = 0.5 * D_inv_x.dot(self.position) - self.j.dot(self.position)
-        return H.real
-
-    @property
-    def gradient(self):
-        D_inv_x = self.D_inverse_x()
-        g = D_inv_x - self.j
-        return_g = g.copy_empty(dtype=np.float)
-        return_g.val = g.val.real
-        return return_g
-
-    @property
-    def curvature(self):
-        class Dummy(object):
-            def __init__(self, x):
-                self.x = x
-            def inverse_times(self, *args, **kwargs):
-                return self.x.times(*args, **kwargs)
-        my_dummy = Dummy(self.D)
-        return my_dummy
-
-
-    @memo
-    def D_inverse_x(self):
-        return D.inverse_times(self.position)
+class AdjointFFTResponse(LinearOperator):
+    def __init__(self, FFT, R, default_spaces=None):
+        super(ResponseOperator, self).__init__(default_spaces)
+        self._domain = FFT.target
+        self.target = R.target
+        self.R = R
+        self.FFT = FFT
+
+    def _times(self, x):
+        return self.R(self.FFT.adjoint_times(x))
+
+    def _adjoint_times(self, x):
+        return self.FFT(self.R.adjoint_times(x))
+
 
 
 if __name__ == "__main__":
@@ -72,10 +47,11 @@ if __name__ == "__main__":
     ss = fft.inverse_times(sh)
 
     # model the measurement process
-    R = SmoothingOperator(s_space, sigma=0.01)
-#    R = DiagonalOperator(s_space, diagonal=1.)
-#    R._diagonal.val[200:400, 200:400] = 0
+    Instrument = SmoothingOperator(s_space, sigma=0.01)
 
+#    Instrument = DiagonalOperator(s_space, diagonal=1.)
+#    Instrument._diagonal.val[200:400, 200:400] = 0
+    R = AdjointFFTResponse(fft, Instrument)
     signal_to_noise = 1
     N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
     n = Field.from_random(domain=s_space,
@@ -84,15 +60,11 @@ if __name__ == "__main__":
                           mean=0)
 
     # create mock data
-    d = R(ss) + n
-
-    # set up reconstruction objects
-    j = R.adjoint_times(N.inverse_times(d))
-    D = PropagatorOperator(S=S, N=N, R=R)
+    d = R(sh) + n
 
     def distance_measure(energy, iteration):
-        x = energy.position
-        print (iteration, ((x-ss).norm()/ss.norm()).real)
+        x = energy.value
+        print (x, iteration)
 
 #    minimizer = SteepestDescent(convergence_tolerance=0,
 #                                iteration_limit=50,
@@ -107,13 +79,14 @@ if __name__ == "__main__":
 #                        callback=distance_measure,
 #                        max_history_length=3)
 
+    solution = energy.analytic_solution()
     m0 = Field(s_space, val=1.)
 
     energy = WienerFilterEnergy(position=m0, D=D, j=j)
 
     (energy, convergence) = minimizer(energy)
 
-    m = energy.position
+    m = fft.adjoint_times(energy.position)
 
     d_data = d.val.get_full_data().real
     if rank == 0:

nifty/energies/__init__.py

@@ -19,3 +19,4 @@
 from energy import Energy
 from line_energy import LineEnergy
 from memoization import memo
+from wiener_filter_energy import WienerFilterEnergy
\ No newline at end of file

nifty/energies/wiener_filter_energy.py

+from .energy import Energy
+from nifty.operators import WienerFilterCurvature
+
+class WienerFilterEnergy(Energy):
+    """The Energy for the Wiener filter.
+
+    It describes the situation of linear measurement with
+    Gaussian noise and Gaussain signal prior.
+
+    Parameters
+    ----------
+    d : Field,
+        the data.
+    R : Operator,
+        The response operator, describtion of the measurement process.
+    N : EndomorphicOperator,
+        The noise covariance in data space.
+    S : EndomorphicOperator,
+        The prior signal covariance in harmonic space.
+    """
+
+    def __init__(self, position, d, R, N, S):
+        super(WienerFilterEnergy, self).__init__(position)
+        self.d = d
+        self.R = R
+        self.N = N
+        self.S = S
+
+    def at(self, position):
+        return self.__class__(position, self.d, self.R, self.N, self.S)
+
+    @property
+    def value(self):
+        energy = 0.5 * self.position.dot(self.S.inverse_times(self.position))
+        energy += 0.5 * (self.d - self.R(self.position)).dot(
+            self.N.inverse_times(self.d - self.R(self.position)))
+        return energy
+
+    @property
+    def gradient(self):
+        gradient = self.S.inverse_times(self.position)
+        gradient -= self.N.inverse_times(self.d - self.R(self.position))
+        return gradient
+
+    @property
+    def curvature(self):
+        curvature = WienerFilterCurvature(R=self.R, N=self.N, S=self.S)
+        return curvature
+
+    def analytic_solution(self):
+        D_inverse = self.curvature()
+        j = self.R.adjoint_times(self.N.inverse_times(self.d))
+        new_position = D_inverse.inverse_times(j)
+        return self.at(new_position)
+

nifty/operators/__init__.py

@@ -39,3 +39,5 @@ from propagator_operator import HarmonicPropagatorOperator
 from composed_operator import ComposedOperator
 
 from response_operator import ResponseOperator
+
+from curvature_operators import WienerFilterCurvature

nifty/operators/curvature_operators/__init__.py

+from wiener_filter_curvature import WienerFilterCurvature
\ No newline at end of file

nifty/operators/curvature_operators/wiener_filter_curvature.py

+from nifty.operators import EndomorphicOperator,\
+                            InvertibleOperatorMixin
+
+
+class WienerFilterCurvature(InvertibleOperatorMixin, EndomorphicOperator):
+    def __init__(self, R, N, S, inverter=None, preconditioner=None):
+
+        self.R = R
+        self.N = N
+        self.S = S
+        if preconditioner is None:
+            preconditioner = self.S.times
+        self._domain = self.S.domain
+        super(WienerFilterCurvature, self).__init__(inverter=inverter,
+                                                 preconditioner=preconditioner)
+    @property
+    def domain(self):
+        return self._domain
+
+    @property
+    def self_adjoint(self):
+        return True
+
+    @property
+    def unitary(self):
+        return False
+
+    # ---Added properties and methods---
+
+    def _times(self, x, spaces):
+        return self.R.adjoint_times(self.N.inverse_times(self.R(x)))\
+               + self.S.inverse_times(x)