diff --git a/demos/wiener_filter_easy.py b/demos/wiener_filter_easy.py
deleted file mode 100644
index 1693f568bfd7fab56fb0d7c92b3fa514ac3c3638..0000000000000000000000000000000000000000
--- a/demos/wiener_filter_easy.py
+++ /dev/null
@@ -1,77 +0,0 @@
-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)
-
diff --git a/demos/wiener_filter_via_curvature.py b/demos/wiener_filter_via_curvature.py
new file mode 100644
index 0000000000000000000000000000000000000000..0224333321b1bf1621bc0c10ae3a759d9d0fb0f6
--- /dev/null
+++ b/demos/wiener_filter_via_curvature.py
@@ -0,0 +1,81 @@
+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
diff --git a/demos/wiener_filter_advanced.py b/demos/wiener_filter_via_hamiltonian.py
similarity index 100%
rename from demos/wiener_filter_advanced.py
rename to demos/wiener_filter_via_hamiltonian.py
diff --git a/nifty/operators/smoothing_operator/fft_smoothing_operator.py b/nifty/operators/smoothing_operator/fft_smoothing_operator.py
index f9e825db33757fb058ae909b94934cdeab4e39f4..9cc5d156ede00525faaa4b185d94f917171adb56 100644
--- a/nifty/operators/smoothing_operator/fft_smoothing_operator.py
+++ b/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]]
diff --git a/nifty/plotting/plotter/healpix_plotter.py b/nifty/plotting/plotter/healpix_plotter.py
index 33eb6863971ab9afcefe08781ee7c282186a8c31..506b9c08351ff61ab7d4894205a4da7b538b4817 100644
--- a/nifty/plotting/plotter/healpix_plotter.py
+++ b/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
diff --git a/nifty/spaces/lm_space/lm_space.py b/nifty/spaces/lm_space/lm_space.py
index 42022f056978952ab90e011e9a911b4e7a2be3d2..d6ec2da05927508397d577d642544e52a82f8c43 100644
--- a/nifty/spaces/lm_space/lm_space.py
+++ b/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)