started critical filtering

demos/nonlinear_wiener_filter.py

@@ -54,7 +54,7 @@ if __name__ == "__main__":
     distribution_strategy = 'not'
 
     # Set up position space
-    s_space = RGSpace([128,128])
+    s_space = RGSpace([256,256])
     # s_space = HPSpace(32)
 
     # Define harmonic transformation and associated harmonic space
@@ -173,3 +173,9 @@ if __name__ == "__main__":
     if rank == 0:
         pl.plot([go.Heatmap(z=m_data)], filename='map.html')
 
+    f_m_data = function(m).val.get_full_data().real
+    if rank == 0:
+        pl.plot([go.Heatmap(z=f_m_data)], filename='f_map.html')
+    f_ss_data = function(ss).val.get_full_data().real
+    if rank == 0:
+        pl.plot([go.Heatmap(z=f_ss_data)], filename='f_ss.html')

nifty/library/energy_library/__init__.py

 from wiener_filter_energy import WienerFilterEnergy
-from nonlinear_wiener_filter_energy import NonlinearWienerFilterEnergy
\ No newline at end of file
+from nonlinear_wiener_filter_energy import NonlinearWienerFilterEnergy
+from critical_power_energy import CriticalPowerEnergy
\ No newline at end of file

nifty/library/energy_library/critical_power_energy.py

+from nifty.energies.energy import Energy
+from nifty.library.operator_library import CriticalPowerCurvature
+from nifty.sugar import generate_posterior_sample
+from nifty import Field
+
+class CriticalPowerEnergy(Energy):
+    """The Energy for the Gaussian lognormal case.
+
+    It describes the situation of linear measurement  of a
+    lognormal signal with Gaussian noise and Gaussain signal prior.
+
+    Parameters
+    ----------
+    d : Field,
+        the data.
+    R : Operator,
+        The nonlinear 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, m, D, alpha, beta, w = None, samples=3):
+        super(CriticalPowerEnergy, self).__init__(position = position)
+        self.al = 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.real
+
+    @property
+    def gradient(self):
+        gradient = self.S.inverse_times(self.position)
+        gradient -= self.R.derived_adjoint_times(
+                    self.N.inverse_times(self.d - self.R(self.position)), self.position)
+        return gradient
+
+    @property
+    def curvature(self):
+        curvature =CriticalPowerCurvature(R=self.R,
+                                                   N=self.N,
+                                                   S=self.S,
+                                                   position=self.position)
+        return curvature
+
+    def _calculate_w(self, m, D, samples):
+        w = Field(domain=self.position.domain, val=0)
+        for i in range(samples):
+            posterior_sample = generate_posterior_sample(m, D)
+            projected_sample =posterior_sample.power_analyze()**2
+            w += projected_sample
+        return w / float(samples)
+
+
+

nifty/library/energy_library/nonlinear_wiener_filter_energy.py

 from nifty.energies.energy import Energy
 from nifty.library.operator_library import NonlinearWienerFilterCurvature
-
+from nifty import Field
 class NonlinearWienerFilterEnergy(Energy):
     """The Energy for the Gaussian lognormal case.
 
@@ -51,4 +51,3 @@ class NonlinearWienerFilterEnergy(Energy):
                                                    S=self.S,
                                                    position=self.position)
         return curvature
-

nifty/library/operator_library/__init__.py

 from wiener_filter_curvature import WienerFilterCurvature
-from nonlinear_wiener_filter_curvature import NonlinearWienerFilterCurvature
\ No newline at end of file
+from nonlinear_wiener_filter_curvature import NonlinearWienerFilterCurvature
+from critical_power_curvature import CriticalPowerCurvature
\ No newline at end of file

nifty/library/operator_library/critical_power_curvature.py

+from nifty.operators import EndomorphicOperator,\
+                            InvertibleOperatorMixin
+
+
+class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
+    def __init__(self, R, N, S, position, inverter=None, preconditioner=None):
+
+        self.R = R
+        self.N = N
+        self.S = S
+        self.position = position
+        if preconditioner is None:
+            preconditioner = self.S.times
+        self._domain = self.S.domain
+        super(CriticalPowerCurvature, 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.derived_adjoint_times(
+                self.N.inverse_times(self.R.derived_times(
+                    x, self.position)), self.position)\
+               + self.S.inverse_times(x)

nifty/sugar.py

@@ -19,7 +19,8 @@
 from nifty import PowerSpace,\
                   Field,\
                   DiagonalOperator,\
-                  FFTOperator
+                  FFTOperator,\
+                  sqrt
 
 __all__ = ['create_power_operator']
 
@@ -41,4 +42,23 @@ def create_power_operator(domain, power_spectrum, dtype=None,
 
     power_operator = DiagonalOperator(domain, diagonal=f, bare=True)
 
-    return power_operator
\ No newline at end of file
+    return power_operator
+
+
+def generate_posterior_sample(mean, covariance):
+    S = covariance.S
+    R = covariance.R
+    N = covariance.N
+    power = S.diagonal()
+    noise = N.diagonal().val
+    mock_signal = Field.from_random(random_type="normal", domain=S.domain,
+                                    std = sqrt(power), dtype = power.dtype)
+    mock_noise = Field.from_random(random_type="normal", domain=N.domain,
+                                   std = sqrt(noise), dtype = noise.dtype)
+    mock_data = R(mock_signal) + mock_noise
+
+    mock_j = R.adjoint_times(N.inverse_times(mock_data))
+    mock_m = covariance.inverse_times(mock_j)
+    sample = mock_signal - mock_m + mean
+    return sample
+