maybe CriticalPowerEnergy working, probably problem in samples or something

demos/critical_filtering.py

@@ -10,7 +10,7 @@ from mpi4py import MPI
 comm = MPI.COMM_WORLD
 rank = comm.rank
 
-np.random.seed(43)
+np.random.seed(44)
 
 
 def plot_parameters(m, t, p, p_sig,p_d):
@@ -59,7 +59,7 @@ if __name__ == "__main__":
     distribution_strategy = 'not'
 
     # Set up position space
-    dist = 1/128. *0.1
+    dist = 1/128. *10
     s_space = ift.RGSpace([128, 128], distances=[dist,dist])
     # s_space = ift.HPSpace(32)
 
@@ -74,8 +74,8 @@ if __name__ == "__main__":
                              distribution_strategy=distribution_strategy)
 
     # Choose the prior correlation structure and defining correlation operator
-    # p_spec = (lambda k: (.5 / (k + 1) ** 3))
-    p_spec = (lambda k: 1)
+    p_spec = (lambda k: (.5 / (k + 1) ** 3))
+    # p_spec = (lambda k: 1)
     S = ift.create_power_operator(h_space, power_spectrum=p_spec,
                                   distribution_strategy=distribution_strategy)
 
@@ -93,7 +93,7 @@ if __name__ == "__main__":
     # Add a harmonic transformation to the instrument
     R = AdjointFFTResponse(fft, Instrument)
 
-    noise = 1.
+    noise = .1
     ndiag = ift.Field(s_space, noise).weight(1)
     N = ift.DiagonalOperator(s_space, ndiag)
     n = ift.Field.from_random(domain=s_space,
@@ -130,7 +130,8 @@ if __name__ == "__main__":
     flat_power = ift.Field(p_space, val=1e-8)
     m0 = flat_power.power_synthesize(real_signal=True)
 
-    t0 = ift.Field(p_space, val=np.log(1./(1+p_space.kindex)**2))
+    # t0 = ift.Field(p_space, val=np.log(1./(1+p_space.kindex)**2))
+    t0 = ift.Field(p_space, val=-5)
 
     for i in range(500):
         S0 = ift.create_power_operator(h_space, power_spectrum=ift.exp(t0),
@@ -139,13 +140,13 @@ if __name__ == "__main__":
         # Initialize non-linear Wiener Filter energy
         map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
         # Solve the Wiener Filter analytically
-        # D0 = map_energy.curvature
-        # m0 = D0.inverse_times(j)
+        D0 = map_energy.curvature
+        m0 = D0.inverse_times(j)
         # Initialize power energy with updated parameters
-        power_energy = CriticalPowerEnergy(position=t0, m=sh, D=None,
-                                           smoothness_prior=1e-15, samples=3)
+        power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0,
+                                           smoothness_prior=1e-15, samples=5)
 
-        (power_energy, convergence) = minimizer2(power_energy)
+        (power_energy, convergence) = minimizer1(power_energy)
 
         # Set new power spectrum
         t0.val = power_energy.position.val.real

nifty/library/critical_filter/critical_power_energy.py

 from ...energies.energy import Energy
 from ...operators.smoothness_operator import SmoothnessOperator
 from ...operators.diagonal_operator import DiagonalOperator
+from ...operators.linear_operator import LinearOperator
+from ...operators.power_projection_operator import PowerProjection
 from . import CriticalPowerCurvature
 from ...memoization import memo
 from ...minimization import ConjugateGradient
@@ -70,6 +72,7 @@ class CriticalPowerEnergy(Energy):
                                     strength=smoothness_prior,
                                     logarithmic=logarithmic)
         self.rho = self.position.domain[0].rho
+        self.P = PowerProjection(domain=self.m.domain,target=self.position.domain)
         self._w = w if w is not None else None
         if inverter is None:
             preconditioner = DiagonalOperator(self._theta.domain,
@@ -135,15 +138,12 @@ class CriticalPowerEnergy(Energy):
                     self.logger.info("Drawing sample %i" % i)
                     posterior_sample = generate_posterior_sample(
                                                             self.m, self.D)
-                    projected_sample = posterior_sample.power_analyze(
-                     binbounds=self.position.domain[0].binbounds)
-                    w += (projected_sample) * self.rho
+                    w += self.P(abs(posterior_sample) ** 2)
+
                 w /= float(self.samples)
             else:
-                w = self.m.power_analyze(
-                     binbounds=self.position.domain[0].binbounds)
-                w *= self.rho
-            self._w = w.weight(-1)
+                w = self.P(abs(self.m)**2)
+            self._w = w
         return self._w
 
     @property

nifty/operators/power_projection_operator/__init__.py

+from .power_projection_operator import PowerProjection
\ No newline at end of file

nifty/operators/power_projection_operator/power_projection_operator.py

+from ... import Field,\
+                FieldArray
+from ..linear_operator import LinearOperator
+
+
+class PowerProjection(LinearOperator):
+    def __init__(self, domain, target, spaces=0, default_spaces=None):
+        self._domain = self._parse_domain(domain)
+        self._target = self._parse_domain(target)
+        self.pindex = self.target[spaces].pindex
+        super(PowerProjection, self).__init__(default_spaces)
+
+    def _times(self,x,spaces):
+        projected_x = self.pindex.bincount(weights=x.weight(1).val.real)
+        y = Field(self.target, val=projected_x).weight(-1)
+        return y
+
+    def _adjoint_times(self,x,spaces):
+        if spaces is None:
+            spaces = 0
+        y = Field(self.domain, val=1.)
+        axes = x.domain_axes
+        spec = x.val.get_full_data()
+
+        spec = self._spec_to_rescaler(spec, spaces, axes)
+        y.val.apply_scalar_function(lambda x: x * spec.real,
+                                            inplace=True)
+
+        return y
+
+    def _spec_to_rescaler(self, spec, power_space_index,axes):
+
+        # weight the random fields with the power spectrum
+        # therefore get the pindex from the power space
+        # take the local data from pindex. This data must be compatible to the
+        # local data of the field given the slice of the PowerSpace
+        # local_distribution_strategy = \
+        #     result.val.get_axes_local_distribution_strategy(
+        #         result.domain_axes[power_space_index])
+        #
+        # if self.pindex.distribution_strategy is not local_distribution_strategy:
+        #     self.logger.warn(
+        #         "The distribution_stragey of pindex does not fit the "
+        #         "slice_local distribution strategy of the synthesized field.")
+
+        # Now use numpy advanced indexing in order to put the entries of the
+        # power spectrum into the appropriate places of the pindex array.
+        # Do this for every 'pindex-slice' in parallel using the 'slice(None)'s
+        local_pindex = self.pindex.get_local_data(copy=False)
+
+        local_blow_up = [slice(None)]*len(spec.shape)
+        local_blow_up[axes[power_space_index][0]] = local_pindex
+        # here, the power_spectrum is distributed into the new shape
+        local_rescaler = spec[local_blow_up]
+        return local_rescaler
+    @property
+    def domain(self):
+        return self._domain
+
+    @property
+    def target(self):
+        return self._target
+
+    @property
+    def unitary(self):
+        return False
\ No newline at end of file

nifty/sugar.py

@@ -17,7 +17,7 @@
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
 import numpy as np
-
+from .operators.power_projection_operator import PowerProjection
 from . import Space,\
                   PowerSpace,\
                   Field,\
@@ -71,14 +71,14 @@ def create_power_operator(domain, power_spectrum, dtype=None,
 
     fp = Field(power_domain, val=power_spectrum, dtype=dtype,
                distribution_strategy='not')
-    f = fp.power_synthesize(mean=1, std=0, real_signal=False,
-                            distribution_strategy=distribution_strategy)
+    P = PowerProjection(domain, power_domain)
+    f = P.adjoint_times(fp)
 
     if not issubclass(fp.dtype.type, np.complexfloating):
         f = f.real
 
-    f **= 2
-    diag = Field(domain,f).weight()
+    # f **= 2
+    diag = Field(domain,f)#.weight()
     return DiagonalOperator(domain, diag)
 
 
@@ -109,9 +109,10 @@ def generate_posterior_sample(mean, covariance):
     R = covariance.R
     N = covariance.N
 
-    power = sqrt(S.diagonal().power_analyze())
-    mock_signal = power.power_synthesize(real_signal=True)
-
+    power = sqrt(S.diagonal())
+    mock_signal = Field.from_random(random_type="normal", domain=S.domain,
+                                   dtype=power.dtype)
+    mock_signal = power * mock_signal
     noise = N.diagonal().weight(-1)
 
     mock_noise = Field.from_random(random_type="normal", domain=N.domain,