Fixed bug in LMSpace.get_distance_array and HealpixPlotter. Worked on demos.

demos/wiener_filter_easy.py

-import numpy as np
-
-from nifty import RGSpace, PowerSpace, Field, FFTOperator, ComposedOperator,\
-                  SmoothingOperator, DiagonalOperator, create_power_operator
-from nifty.library import WienerFilterCurvature
-
-#import plotly.offline as pl
-#import plotly.graph_objs as go
-
-from mpi4py import MPI
-comm = MPI.COMM_WORLD
-rank = comm.rank
-
-
-if __name__ == "__main__":
-
-    distribution_strategy = 'fftw'
-
-    # Setting up physical constants
-    # total length of Interval or Volume the field lives on, e.g. in meters
-    L = 2.
-    # typical distance over which the field is correlated (in same unit as L)
-    correlation_length = 0.1
-    # variance of field in position space sqrt(<|s_x|^2>) (in unit of s)
-    field_variance = 2.
-    # smoothing length of response (in same unit as L)
-    response_sigma = 0.1
-
-    # defining resolution (pixels per dimension)
-    N_pixels = 512
-
-    # Setting up derived constants
-    k_0 = 1./correlation_length
-    # note that field_variance**2 = a*k_0/4. for this analytic form of power
-    # spectrum
-    a = field_variance**2/k_0*4.
-    pow_spec = (lambda k: a / (1 + k/k_0) ** 4)
-    pixel_length = L/N_pixels
-
-    # Setting up the geometry
-    s_space = RGSpace([N_pixels, N_pixels], distances=pixel_length)
-    fft = FFTOperator(s_space, domain_dtype=np.float, target_dtype=np.complex)
-    h_space = fft.target[0]
-    inverse_fft = FFTOperator(h_space, target=s_space,
-                              domain_dtype=np.complex, target_dtype=np.float)
-    p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
-
-    # Creating the mock data
-
-    S = create_power_operator(h_space, power_spectrum=pow_spec,
-                              distribution_strategy=distribution_strategy)
-
-    sp = Field(p_space, val=pow_spec,
-               distribution_strategy=distribution_strategy)
-    sh = sp.power_synthesize(real_signal=True)
-    ss = fft.inverse_times(sh)
-
-    R = SmoothingOperator(s_space, sigma=response_sigma)
-    R_harmonic = ComposedOperator([inverse_fft, R], default_spaces=[0, 0])
-
-    signal_to_noise = 1
-    N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
-    n = Field.from_random(domain=s_space,
-                          random_type='normal',
-                          std=ss.std()/np.sqrt(signal_to_noise),
-                          mean=0)
-
-    d = R(ss) + n
-
-    # Wiener filter
-
-    j = R_harmonic.adjoint_times(N.inverse_times(d))
-    wiener_curvature = WienerFilterCurvature(S=S, N=N, R=R_harmonic)
-
-    m = wiener_curvature.inverse_times(j)
-    m_s = inverse_fft(m)
-

demos/wiener_filter_via_curvature.py

+import numpy as np
+
+from nifty import RGSpace, PowerSpace, Field, FFTOperator, ComposedOperator,\
+                  DiagonalOperator, ResponseOperator, plotting,\
+                  create_power_operator
+from nifty.library import WienerFilterCurvature
+
+
+if __name__ == "__main__":
+
+    distribution_strategy = 'not'
+
+    # Setting up variable parameters
+
+    # Typical distance over which the field is correlated
+    correlation_length = 0.01
+    # Variance of field in position space sqrt(<|s_x|^2>)
+    field_variance = 2.
+    # smoothing length of response (in same unit as L)
+    response_sigma = 0.1
+    # The signal to noise ratio
+    signal_to_noise = 0.7
+
+    # note that field_variance**2 = a*k_0/4. for this analytic form of power
+    # spectrum
+    def power_spectrum(k):
+        a = 4 * correlation_length * field_variance**2
+        return a / (1 + k * correlation_length) ** 4
+
+    # Setting up the geometry
+
+    # Total side-length of the domain
+    L = 2.
+    # Grid resolution (pixels per axis)
+    N_pixels = 512
+
+    signal_space = RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
+    harmonic_space = FFTOperator.get_default_codomain(signal_space)
+    fft = FFTOperator(harmonic_space, target=signal_space,
+                      domain_dtype=np.complex, target_dtype=np.float)
+    power_space = PowerSpace(harmonic_space,
+                             distribution_strategy=distribution_strategy)
+
+    # Creating the mock data
+    S = create_power_operator(harmonic_space, power_spectrum=power_spectrum,
+                              distribution_strategy=distribution_strategy)
+
+    mock_power = Field(power_space, val=power_spectrum,
+                       distribution_strategy=distribution_strategy)
+    np.random.seed(43)
+    mock_harmonic = mock_power.power_synthesize(real_signal=True)
+    mock_signal = fft(mock_harmonic)
+
+    R = ResponseOperator(signal_space, sigma=(response_sigma,))
+    data_domain = R.target[0]
+    R_harmonic = ComposedOperator([fft, R], default_spaces=[0, 0])
+
+    N = DiagonalOperator(data_domain,
+                         diagonal=mock_signal.var()/signal_to_noise,
+                         bare=True)
+    noise = Field.from_random(domain=data_domain,
+                              random_type='normal',
+                              std=mock_signal.std()/np.sqrt(signal_to_noise),
+                              mean=0)
+    data = R(mock_signal) + noise
+
+    # Wiener filter
+
+    j = R_harmonic.adjoint_times(N.inverse_times(data))
+    wiener_curvature = WienerFilterCurvature(S=S, N=N, R=R_harmonic)
+
+    m = wiener_curvature.inverse_times(j)
+    m_s = fft(m)
+
+    plotter = plotting.RG2DPlotter()
+    plotter.title = 'mock_signal.html';
+    plotter(mock_signal)
+    plotter.title = 'data.html'
+    plotter(Field(signal_space,
+                  val=data.val.get_full_data().reshape(signal_space.shape)))
+    plotter.title = 'map.html'; plotter(m_s)
\ No newline at end of file

nifty/operators/smoothing_operator/fft_smoothing_operator.py

@@ -21,7 +21,6 @@ class FFTSmoothingOperator(SmoothingOperator):
         # transform to the (global-)default codomain and perform all remaining
         # steps therein
         transformator = self._get_transformator(x.dtype)
-
         transformed_x = transformator(x, spaces=spaces)
         codomain = transformed_x.domain[spaces[0]]
         coaxes = transformed_x.domain_axes[spaces[0]]

nifty/plotting/plotter/healpix_plotter.py

@@ -21,3 +21,6 @@ class HealpixPlotter(PlotterBase):
 
     def _initialize_figure(self):
         return Figure2D(plots=None)
+
+    def _parse_data(self, data, field, spaces):
+        return data

nifty/spaces/lm_space/lm_space.py

@@ -22,7 +22,7 @@ import numpy as np
 
 from nifty.spaces.space import Space
 
-from d2o import arange
+from d2o import arange, distributed_data_object
 
 
 class LMSpace(Space):
@@ -127,6 +127,24 @@ class LMSpace(Space):
             return x.copy()
 
     def get_distance_array(self, distribution_strategy):
+        if distribution_strategy == 'not':  # short cut
+            lmax = self.lmax
+            ldist = np.empty((self.dim,), dtype=np.float64)
+            ldist[0:lmax+1] = np.arange(lmax+1, dtype=np.float64)
+            tmp = np.empty((2*lmax+2), dtype=np.float64)
+            tmp[0::2] = np.arange(lmax+1)
+            tmp[1::2] = np.arange(lmax+1)
+            idx = lmax+1
+            for l in range(1, lmax+1):
+                ldist[idx:idx+2*(lmax+1-l)] = tmp[2*l:]
+                idx += 2*(lmax+1-l)
+            dists = distributed_data_object(
+                            global_shape=self.shape,
+                            dtype=np.float,
+                            distribution_strategy=distribution_strategy)
+            dists.set_local_data(ldist)
+            return dists
+
         dists = arange(start=0, stop=self.shape[0],
                        distribution_strategy=distribution_strategy)