Refactored critical power files.

demos/critical_filtering.py

@@ -123,9 +123,6 @@ if __name__ == "__main__":
                        callback=convergence_measure,
                        max_history_length=3)
 
-    inverter = ConjugateGradient(convergence_level=1,
-                                 convergence_tolerance=10e-4,
-                                 preconditioner=None)
     # Setting starting position
     flat_power = Field(p_space,val=10e-8)
     m0 = flat_power.power_synthesize(real_signal=True)
@@ -144,7 +141,7 @@ if __name__ == "__main__":
         D0 = map_energy.curvature
         m0 = D0.inverse_times(j)
         # Initializing the power energy with updated parameters
-        power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=10., samples=3, inverter=inverter)
+        power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, smoothness_prior=10., samples=3)
 
         (power_energy, convergence) = minimizer1(power_energy)
 
@@ -153,6 +150,8 @@ if __name__ == "__main__":
         t0.val  = power_energy.position.val.real
 
         # Plotting current estimate
-        plot_parameters(m0,t0,log(sp), data_power)
+        print i
+        if i%50 == 0:
+            plot_parameters(m0,t0,log(sp), data_power)
 
 

nifty/library/critical_filter/critical_power_energy.py

+import numpy as np
+
 from nifty.energies.energy import Energy
 from nifty.operators.smoothness_operator import SmoothnessOperator
 from nifty.library.critical_filter import CriticalPowerCurvature
@@ -10,9 +12,10 @@ from nifty import Field, exp
 class CriticalPowerEnergy(Energy):
     """The Energy of the power spectrum according to the critical filter.
 
-    It describes the energy of the logarithmic amplitudes of the power spectrum of
-    a Gaussian random field under reconstruction uncertainty with smoothness and
-    inverse gamma prior. It is used to infer the correlation structure of a correlated signal.
+    It describes the energy of the logarithmic amplitudes of the power spectrum
+    of a Gaussian random field under reconstruction uncertainty with smoothness
+    and inverse gamma prior. It is used to infer the correlation structure of a
+    correlated signal.
 
     Parameters
     ----------
@@ -25,20 +28,21 @@ class CriticalPowerEnergy(Energy):
         If not specified, the map is assumed to be no reconstruction.
         default : None
     alpha : float
-        The spectral prior of the inverse gamma distribution. 1.0 corresponds to
-        non-informative.
+        The spectral prior of the inverse gamma distribution. 1.0 corresponds
+        to non-informative.
         default : 1.0
     q : float
-        The cutoff parameter of the inverse gamma distribution. 0.0 corresponds to
-        non-informative.
+        The cutoff parameter of the inverse gamma distribution. 0.0 corresponds
+        to non-informative.
+        default : 0.0
+    smoothness_prior : float
+        Controls the strength of the smoothness prior
         default : 0.0
-    sigma : float
-        The parameter of the smoothness prior.
-        default : ??? None? ???????
     logarithmic : boolean
         Whether smoothness acts on linear or logarithmic scale.
     samples : integer
-        Number of samples used for the estimation of the uncertainty corrections.
+        Number of samples used for the estimation of the uncertainty
+        corrections.
         default : 3
     w : Field
         The contribution from the map with or without uncertainty. It is used
@@ -49,76 +53,72 @@ class CriticalPowerEnergy(Energy):
         default : None
     """
 
-    def __init__(self, position, m, D=None, alpha =1.0, q=0., sigma=0.,
-                 logarithmic = True, samples=3, w=None, inverter=None):
-        super(CriticalPowerEnergy, self).__init__(position = position)
+    def __init__(self, position, m, D=None, alpha=1.0, q=0.,
+                 smoothness_prior=0., logarithmic=True, samples=3, w=None):
+        super(CriticalPowerEnergy, self).__init__(position=position)
         self.m = m
         self.D = D
         self.samples = samples
-        self.sigma = sigma
+        self.smoothness_prior = np.float(smoothness_prior)
         self.alpha = Field(self.position.domain, val=alpha)
         self.q = Field(self.position.domain, val=q)
-        self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma,
+        self.T = SmoothnessOperator(domain=self.position.domain[0],
+                                    strength=self.smoothness_prior,
                                     logarithmic=logarithmic)
         self.rho = self.position.domain[0].rho
-        self.inverter = inverter
-        self.w = w
-        if self.w is None:
-            self.w = self._calculate_w(self.m, self.D, self.samples)
-        self.theta = (exp(-self.position) * (self.q + self.w / 2.))
-
+        self._w = w if w is not None else None
 
     def at(self, position):
-        return self.__class__(position, self.m, D=self.D,
-                              alpha =self.alpha,
-                              q=self.q,
-                              sigma=self.sigma, w=self.w,
-                              samples=self.samples)
+        return self.__class__(position, self.m, D=self.D, alpha=self.alpha,
+                              q=self.q, smoothness_prior=self.smoothness_prior,
+                              w=self.w, samples=self.samples)
 
     @property
     def value(self):
-        energy = self._theta.vdot(Field(self.position.domain,val=1.), bare= True)
+        energy = self._theta.sum()
         energy += self.position.vdot(self._rho_prime, bare=True)
         energy += 0.5 * self.position.vdot(self._Tt)
         return energy.real
 
     @property
     def gradient(self):
-        gradient = - self._theta.weight(-1)
+        gradient = -self._theta.weight(-1)
         gradient += (self._rho_prime).weight(-1)
-        gradient +=  self._Tt
+        gradient += self._Tt
         gradient.val = gradient.val.real
         return gradient
 
     @property
     def curvature(self):
-        curvature = CriticalPowerCurvature(theta=self._theta.weight(-1), T=self.T,
-                                           inverter=self.inverter)
+        curvature = CriticalPowerCurvature(theta=self._theta.weight(-1),
+                                           T=self.T)
         return curvature
 
-    def _calculate_w(self, m, D, samples):
-        w = Field(domain=self.position.domain, val=0. , dtype=m.dtype)
-        if D is not None:
-            for i in range(samples):
-                posterior_sample = generate_posterior_sample(m, D, inverter = self.inverter)
-                projected_sample = posterior_sample.power_analyze(
-                    logarithmic=self.position.domain[0].config["logarithmic"],
-                    nbin= self.position.domain[0].config["nbin"])
-                w += (projected_sample) * self.rho
-            w /= float(samples)
-        else:
-            w = m.power_analyze(
-                    logarithmic=self.position.domain[0].config["logarithmic"],
-                        nbin=self.position.domain[0].config["nbin"])
-            w *= self.rho
-
-        return w
-
+    @property
+    def w(self):
+        if self._w is None:
+            w = Field(domain=self.position.domain, val=0., dtype=self.m.dtype)
+            if self.D is not None:
+                for i in range(self.samples):
+                    posterior_sample = generate_posterior_sample(
+                                                            self.m, self.D)
+                    projected_sample = posterior_sample.power_analyze(
+                     logarithmic=self.position.domain[0].config["logarithmic"],
+                     nbin=self.position.domain[0].config["nbin"])
+                    w += (projected_sample) * self.rho
+                w /= float(self.samples)
+            else:
+                w = self.m.power_analyze(
+                     logarithmic=self.position.domain[0].config["logarithmic"],
+                     nbin=self.position.domain[0].config["nbin"])
+                w *= self.rho
+            self._w = w
+        return self._w
 
     @property
     @memo
     def _theta(self):
-        return (exp(-self.position) * (self.q + self.w / 2.))
+        return exp(-self.position) * (self.q + self.w / 2.)
 
     @property
     @memo
@@ -129,4 +129,3 @@ class CriticalPowerEnergy(Energy):
     @memo
     def _Tt(self):
         return self.T(self.position)
-

nifty/library/wiener_filter/wiener_filter_energy.py

@@ -29,7 +29,6 @@ class WienerFilterEnergy(Energy):
         self.R = R
         self.N = N
         self.S = S
-        self.inverter = inverter
 
     def at(self, position):
         return self.__class__(position=position, d=self.d, R=self.R, N=self.N,
@@ -48,10 +47,8 @@ class WienerFilterEnergy(Energy):
     @property
     @memo
     def curvature(self):
-        return WienerFilterCurvature(R=self.R, N=self.N, S=self.S,
-                                     inverter=self.inverter)
+        return WienerFilterCurvature(R=self.R, N=self.N, S=self.S)
 
-    @property
     @memo
     def _Dx(self):
         return self.curvature(self.position)

nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py

@@ -38,7 +38,8 @@ class InvertibleOperatorMixin(object):
         (default: ConjugateGradient)
 
     preconditioner : LinearOperator
-        Preconditions the minimizaion problem
+        Preconditioner that is used by ConjugateGraduent if no minimizer was
+        given.
 
     Attributes
     ----------

nifty/operators/smoothness_operator/smoothness_operator.py

@@ -5,14 +5,17 @@ from nifty.operators.laplace_operator import LaplaceOperator
 
 
 class SmoothnessOperator(EndomorphicOperator):
-    """An operator measuring the smoothness on an irregular grid with respect to some scale.
+    """An operator measuring the smoothness on an irregular grid with respect
+    to some scale.
 
-    This operator applies the irregular LaplaceOperator and its adjoint to some Field over a
-    PowerSpace which corresponds to its smoothness and weights the result with a scale parameter sigma.
-    It is used in the smoothness prior terms of the CriticalPowerEnergy. For this purpose we
-    use free boundary conditions in the LaplaceOperator, having no curvature at both ends.
-    In addition the first entry is ignored as well, corresponding to the overall mean of the map.
-    The mean is therefore not considered in the smoothness prior.
+    This operator applies the irregular LaplaceOperator and its adjoint to some
+    Field over a PowerSpace which corresponds to its smoothness and weights the
+    result with a scale parameter sigma. It is used in the smoothness prior
+    terms of the CriticalPowerEnergy. For this purpose we use free boundary
+    conditions in the LaplaceOperator, having no curvature at both ends. In
+    addition the first entry is ignored as well, corresponding to the overall
+    mean of the map. The mean is therefore not considered in the smoothness
+    prior.
 
 
     Parameters
@@ -26,7 +29,8 @@ class SmoothnessOperator(EndomorphicOperator):
 
     # ---Overwritten properties and methods---
 
-    def __init__(self, domain, sigma, logarithmic=True, default_spaces=None):
+    def __init__(self, domain, strength=1., logarithmic=True,
+                 default_spaces=None):
 
         super(SmoothnessOperator, self).__init__(default_spaces=default_spaces)
 
@@ -37,10 +41,10 @@ class SmoothnessOperator(EndomorphicOperator):
         if not isinstance(self.domain[0], PowerSpace):
             raise TypeError("The domain must contain exactly one PowerSpace.")
 
-        if sigma <= 0:
+        if strength <= 0:
             raise ValueError("ERROR: invalid sigma.")
 
-        self._sigma = sigma
+        self._strength = strength
 
         self._laplace = LaplaceOperator(domain=self.domain,
                                         logarithmic=logarithmic)
@@ -68,11 +72,17 @@ class SmoothnessOperator(EndomorphicOperator):
         return False
 
     def _times(self, x, spaces):
-        res = self._laplace.adjoint_times(self._laplace(x, spaces), spaces)
-        return (1./self.sigma)**2*res
+        if self._strength != 0:
+            result = self._laplace.adjoint_times(self._laplace(x, spaces),
+                                                 spaces)
+            result *= self._strength**2
+        else:
+            result = x.copy_empty()
+            result.val[:] = 0
+        return result
 
     # ---Added properties and methods---
 
     @property
-    def sigma(self):
-        return self._sigma
+    def strength(self):
+        return self._strength

nifty/sugar.py

@@ -64,11 +64,13 @@ def create_power_operator(domain, power_spectrum, dtype=None,
     return DiagonalOperator(domain, diagonal=f, bare=True)
 
 
-def generate_posterior_sample(mean, covariance, inverter=None):
-    """ Generates a posterior sample from a Gaussian distribution with given mean and covariance
+def generate_posterior_sample(mean, covariance):
+    """ Generates a posterior sample from a Gaussian distribution with given
+    mean and covariance
 
-    This method generates samples by setting up the observation and reconstruction of a mock signal
-    in order to obtain residuals of the right correlation which are added to the given mean.
+    This method generates samples by setting up the observation and
+    reconstruction of a mock signal in order to obtain residuals of the right
+    correlation which are added to the given mean.
 
     Parameters
     ----------
@@ -77,9 +79,6 @@ def generate_posterior_sample(mean, covariance, inverter=None):
     covariance : WienerFilterCurvature
         The posterior correlation structure consisting of a
         response operator, noise covariance and prior signal covariance
-    inverter : ConjugateGradient *optional*
-        the conjugate gradient used to invert the curvature for the Wiener filter.
-        default : None
 
     Returns
     -------
@@ -91,8 +90,6 @@ def generate_posterior_sample(mean, covariance, inverter=None):
     S = covariance.S
     R = covariance.R
     N = covariance.N
-    if inverter is None:
-        inverter = ConjugateGradient(preconditioner=S)
 
     power = S.diagonal().power_analyze()**.5
     mock_signal = power.power_synthesize(real_signal=True)