Merge branch 'critical_nightly' into 'nightly'

Add PowerProjection operator; rewrite CriticalPowerEnergy

See merge request !204

demos/critical_filtering.py

@@ -10,20 +10,21 @@ from mpi4py import MPI
 comm = MPI.COMM_WORLD
 rank = comm.rank
 
-np.random.seed(42)
+np.random.seed(44)
 
 
-def plot_parameters(m, t, p, p_d):
+def plot_parameters(m, t, p, p_sig,p_d):
 
     x = np.log(t.domain[0].kindex)
     m = fft.adjoint_times(m)
     m = m.val.get_full_data().real
     t = t.val.get_full_data().real
     p = p.val.get_full_data().real
+    pd_sig = p_sig.val.get_full_data()
     p_d = p_d.val.get_full_data().real
     pl.plot([go.Heatmap(z=m)], filename='map.html', auto_open=False)
     pl.plot([go.Scatter(x=x, y=t), go.Scatter(x=x, y=p),
-             go.Scatter(x=x, y=p_d)], filename="t.html", auto_open=False)
+             go.Scatter(x=x, y=p_d),go.Scatter(x=x, y=pd_sig)], filename="t.html", auto_open=False)
 
 
 class AdjointFFTResponse(ift.LinearOperator):
@@ -58,7 +59,8 @@ if __name__ == "__main__":
     distribution_strategy = 'not'
 
     # Set up position space
-    s_space = ift.RGSpace([128, 128])
+    dist = 1/128. *0.1
+    s_space = ift.RGSpace([128, 128], distances=[dist,dist])
     # s_space = ift.HPSpace(32)
 
     # Define harmonic transformation and associated harmonic space
@@ -73,6 +75,7 @@ if __name__ == "__main__":
 
     # Choose the prior correlation structure and defining correlation operator
     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)
 
@@ -90,7 +93,7 @@ if __name__ == "__main__":
     # Add a harmonic transformation to the instrument
     R = AdjointFFTResponse(fft, Instrument)
 
-    noise = 1.
+    noise = 10.
     ndiag = ift.Field(s_space, noise).weight(1)
     N = ift.DiagonalOperator(s_space, ndiag)
     n = ift.Field.from_random(domain=s_space,
@@ -123,12 +126,12 @@ if __name__ == "__main__":
     IC3 = ift.GradientNormController(iteration_limit=100,
                                      tol_abs_gradnorm=0.1)
     minimizer3 = ift.SteepestDescent(IC3)
-
     # Set starting position
     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),
@@ -141,14 +144,15 @@ if __name__ == "__main__":
         m0 = D0.inverse_times(j)
         # Initialize power energy with updated parameters
         power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0,
-                                           smoothness_prior=10., samples=3)
+                                           smoothness_prior=1., 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
 
         # Plot current estimate
         print(i)
-        if i % 5 == 0:
-            plot_parameters(m0, t0, ift.log(sp), data_power)
+        if i % 1 == 0:
+            plot_parameters(sh, t0, ift.log(sp), ift.log(sh.power_analyze(binbounds=p_space.binbounds)),data_power)
+            print ift.log(sh.power_analyze(binbounds=p_space.binbounds)).val - t0.val

demos/wiener_filter_via_hamiltonian.py

@@ -11,7 +11,7 @@ from mpi4py import MPI
 comm = MPI.COMM_WORLD
 rank = comm.rank
 
-d2o.random.seed(42)
+# d2o.random.seed(42)
 
 
 class AdjointFFTResponse(LinearOperator):
@@ -71,7 +71,7 @@ if __name__ == "__main__":
     # Choosing the measurement instrument
     # Instrument = SmoothingOperator(s_space, sigma=0.05)
     Instrument = DiagonalOperator(s_space, diagonal=1.)
-#    Instrument._diagonal.val[200:400, 200:400] = 0
+    Instrument._diagonal.val[64:80, 32:80] = 0.
 
     # Adding a harmonic transformation to the instrument
     R = AdjointFFTResponse(fft, Instrument)
@@ -123,8 +123,10 @@ if __name__ == "__main__":
     D0 = energy.curvature
 
     # Solving the problem analytically
-#    m0 = D0.inverse_times(j)
-
+    m0 = D0.inverse_times(j)
+    energy, convergence = minimizer(energy)
+    m = energy.position
+    D = energy.curvature
     m = minimizer(energy)[0].position
 
     plotter = plotting.RG2DPlotter()
@@ -132,13 +134,15 @@ if __name__ == "__main__":
     plotter(fft.inverse_times(m), path='m.html')
 
 
-#    sample_variance = Field(sh.domain, val=0.)
-#    sample_mean = Field(sh.domain, val=0.)
+    sample_variance = Field(s_space, val=0.)
+    sample_mean = Field(s_space, val=0.)
 
-#    # sampling the uncertainty map
-#    n_samples = 10
-#    for i in range(n_samples):
-#        sample = fft(sugar.generate_posterior_sample(0., D0))
-#        sample_variance += sample**2
-#        sample_mean += sample
-#    variance = (sample_variance - sample_mean**2)/n_samples
+    # sampling the uncertainty map
+    n_samples = 50
+    for i in range(n_samples):
+        sample = fft.adjoint_times(sugar.generate_posterior_sample(m, D))
+        sample_variance += sample**2
+        sample_mean += sample
+    sample_mean /= n_samples
+    sample_variance /= n_samples
+    variance = (sample_variance - sample_mean**2)

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.power_projection_operator import PowerProjection
 from . import CriticalPowerCurvature
 from ...memoization import memo
 from ...minimization import ConjugateGradient
@@ -70,13 +71,17 @@ 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,
-                                              diagonal=self._theta.weight(-1),
+                                              diagonal=self._theta,
                                               copy=False)
             inverter = ConjugateGradient(preconditioner=preconditioner)
         self._inverter = inverter
+        self.one = Field(self.position.domain, val=1.)
+        self.constants = self.one/2. + self.alpha - 1
 
     @property
     def inverter(self):
@@ -94,16 +99,16 @@ class CriticalPowerEnergy(Energy):
     @property
     @memo
     def value(self):
-        energy = self._theta.sum()
-        energy += self.position.weight(-1).vdot(self._rho_prime)
+        energy = self.one.vdot(self._theta)
+        energy += self.position.vdot(self.constants)
         energy += 0.5 * self.position.vdot(self._Tt)
         return energy.real
 
     @property
     @memo
     def gradient(self):
-        gradient = -self._theta.weight(-1)
-        gradient += (self._rho_prime).weight(-1)
+        gradient = -self._theta
+        gradient += (self.constants)
         gradient += self._Tt
         gradient.val = gradient.val.real
         return gradient
@@ -111,7 +116,7 @@ class CriticalPowerEnergy(Energy):
     @property
     @memo
     def curvature(self):
-        return CriticalPowerCurvature(theta=self._theta.weight(-1), T=self.T,
+        return CriticalPowerCurvature(theta=self._theta, T=self.T,
                                       inverter=self.inverter)
 
     # ---Added properties and methods---
@@ -134,14 +139,11 @@ 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
+                w = self.P(abs(self.m)**2)
             self._w = w
         return self._w
 
@@ -150,11 +152,6 @@ class CriticalPowerEnergy(Energy):
     def _theta(self):
         return exp(-self.position) * (self.q + self.w / 2.)
 
-    @property
-    @memo
-    def _rho_prime(self):
-        return self.alpha - 1. + self.rho / 2.
-
     @property
     @memo
     def _Tt(self):

nifty/operators/__init__.py

@@ -39,3 +39,5 @@ from .response_operator import ResponseOperator
 from .laplace_operator import LaplaceOperator
 
 from .smoothness_operator import SmoothnessOperator
+
+from .power_projection_operator import PowerProjection

nifty/operators/power_projection_operator/__init__.py

+from .power_projection_operator import PowerProjection

nifty/operators/power_projection_operator/power_projection_operator.py

+from ... import Field
+from ..linear_operator import LinearOperator
+from ...spaces.power_space import PowerSpace
+
+
+class PowerProjection(LinearOperator):
+    def __init__(self, domain, target, default_spaces=None):
+        self._domain = self._parse_domain(domain)
+        self._target = self._parse_domain(target)
+        if len(self._domain) != 1 or len(self._target) != 1:
+            raise ValueError("Operator only works over one space")
+        if not self._domain[0].harmonic:
+            raise ValueError("domain must be a harmonic space")
+        if not isinstance(self._target[0], PowerSpace):
+            raise ValueError("target must be a PowerSpace")
+        self.pindex = self.target[0].pindex
+        super(PowerProjection, self).__init__(default_spaces)
+
+    def _times(self, x, spaces):
+        if spaces is None:
+            spaces = 0
+        projected_x = self.pindex.bincount(
+            weights=x.weight(1, spaces=spaces).val.real,
+            axis=x.domain_axes[spaces])
+        tgt_domain = list(x.domain)
+        tgt_domain[spaces] = self._target[0]
+        y = Field(tgt_domain, val=projected_x).weight(-1, spaces=spaces)
+        return y
+
+    def _adjoint_times(self, x, spaces):
+        if spaces is None:
+            spaces = 0
+        tgt_domain = list(x.domain)
+        tgt_domain[spaces] = self._domain[0]
+        y = Field(tgt_domain, val=1.)
+        axes = x.domain_axes[spaces]
+        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]] = 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

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().weight(-1))
+    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,