simplifications

demos/critical_filtering.py

@@ -65,9 +65,7 @@ if __name__ == "__main__":
     noiseless_data = MeasurementOperator(true_sky)
     noise_amplitude = noiseless_data.val.std()*noise_level
     N = ift.ScalingOperator(noise_amplitude**2, d_space)
-    n = ift.Field.from_random(
-        domain=d_space, random_type='normal',
-        std=noise_amplitude, mean=0)
+    n = N.draw_sample()
     # Creating the mock data
     d = noiseless_data + n
 

demos/nonlinear_critical_filter.py

@@ -66,9 +66,7 @@ if __name__ == "__main__":
     noiseless_data = MeasurementOperator(true_sky)
     noise_amplitude = noiseless_data.val.std()*noise_level
     N = ift.ScalingOperator(noise_amplitude**2, d_space)
-    n = ift.Field.from_random(
-        domain=d_space, random_type='normal',
-        std=noise_amplitude, mean=0)
+    n = N.draw_sample()
     # Creating the mock data
     d = noiseless_data + n
 

demos/paper_demos/cartesian_wiener_filter.py

@@ -83,9 +83,7 @@ if __name__ == "__main__":
     noise_amplitude = noiseless_data.val.std()/signal_to_noise
     # Setting up the noise covariance and drawing a random noise realization
     N = ift.ScalingOperator(noise_amplitude**2, data_domain)
-    noise = ift.Field.from_random(
-        domain=data_domain, random_type='normal',
-        std=noise_amplitude, mean=0)
+    noise = N.draw_sample()
     data = noiseless_data + noise
 
     # Wiener filter

demos/paper_demos/wiener_filter.py

@@ -46,9 +46,7 @@ if __name__ == "__main__":
     noise_amplitude = noiseless_data.val.std()/signal_to_noise
     # Setting up the noise covariance and drawing a random noise realization
     N = ift.ScalingOperator(noise_amplitude**2, data_domain)
-    noise = ift.Field.from_random(
-        domain=data_domain, random_type='normal',
-        std=noise_amplitude, mean=0)
+    noise = N.draw_sample()
     data = noiseless_data + noise
 
     # Wiener filter

demos/wiener_filter_via_curvature.py

@@ -73,9 +73,7 @@ if __name__ == "__main__":
     noiseless_data = R(mock_signal)
     noise_amplitude = noiseless_data.val.std()/signal_to_noise
     N = ift.ScalingOperator(noise_amplitude**2, data_domain)
-    noise = ift.Field.from_random(
-        domain=data_domain, random_type='normal',
-        std=noise_amplitude, mean=0)
+    noise = N.draw_sample()
     data = noiseless_data + noise
 
     j = R.adjoint_times(N.inverse_times(data))

demos/wiener_filter_via_hamiltonian.py

@@ -43,10 +43,7 @@ if __name__ == "__main__":
     signal_to_noise = 1.
     noise_amplitude = noiseless_data.val.std()/signal_to_noise
     N = ift.ScalingOperator(noise_amplitude**2, s_space)
-    n = ift.Field.from_random(domain=s_space,
-                              random_type='normal',
-                              std=noise_amplitude,
-                              mean=0)
+    n = N.draw_sample()
 
     # Create mock data
     d = noiseless_data + n

nifty4/operators/diagonal_operator.py

@@ -134,11 +134,11 @@ class DiagonalOperator(EndomorphicOperator):
                                 self._spaces)
 
     def draw_sample(self):
-        if self._spaces is not None:
-            raise ValueError("Cannot draw (yet) from this operator")
+        if np.issubdtype(self._ldiag.dtype, np.complexfloating):
+            raise ValueError("cannot draw sample from complex-valued operator")
 
         res = Field.from_random(random_type="normal",
                                 domain=self._domain,
                                 dtype=self._diagonal.dtype)
-        res *= sqrt(self._diagonal)
+        res.val[()] *= np.sqrt(self._ldiag)
         return res

nifty4/sugar.py

@@ -63,8 +63,9 @@ def power_analyze(field, spaces=None, binbounds=None,
     ----------
     field : Field
         The field to be analyzed
-    spaces : None or int or tuple of int , optional
-        The set of subdomains for which the powerspectrum shall be computed.
+    spaces : None or int or tuple of int, optional
+        The indices of subdomains for which the power spectrum shall be
+        computed.
         If None, all subdomains will be converted.
         (default : None).
     binbounds : None or array-like, optional
@@ -85,7 +86,7 @@ def power_analyze(field, spaces=None, binbounds=None,
     -------
     Field
         The output object. Its domain is a PowerSpace and it contains
-        the power spectrum of 'field'.
+        the power spectrum of `field`.
     """
 
     for sp in field.domain:
@@ -133,7 +134,7 @@ def power_synthesize_nonrandom(field, spaces=None):
 
 
 def power_synthesize(field, spaces=None, real_power=True, real_signal=True):
-    """Returns a sampled field with `field`**2 as its power spectrum.
+    """Returns a sampled field with `field`\**2 as its power spectrum.
 
     This method draws a Gaussian random field in the harmonic partner
     domain of this field's domains, using this field as power spectrum.
@@ -231,12 +232,7 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None):
         An operator that implements the given power spectrum.
     """
     domain = DomainTuple.make(domain)
-    if space is None:
-        if len(domain) != 1:
-            raise ValueError("space keyword must be set")
-        else:
-            space = 0
-    space = int(space)
+    space = utilities.infer_space(domain, space)
     return DiagonalOperator(
         create_power_field(domain[space], power_spectrum, dtype),
         domain=domain, spaces=space)