Fixed bug in _spec_to_rescaler and added stochastic tests for power_analyze/power_synthesize

nifty/field.py

@@ -672,8 +672,11 @@ class Field(Loggable, Versionable, object):
         # Do this for every 'pindex-slice' in parallel using the 'slice(None)'s
         local_pindex = pindex.get_local_data(copy=False)
 
-        local_blow_up = [slice(None)]*len(self.shape)
-        local_blow_up[self.domain_axes[power_space_index][0]] = local_pindex
+        local_blow_up = [slice(None)]*len(spec.shape)
+        # it is important to count from behind, since spec potentially grows
+        # with every iteration
+        index = self.domain_axes[power_space_index][0]-len(self.shape)
+        local_blow_up[index] = local_pindex
         # here, the power_spectrum is distributed into the new shape
         local_rescaler = spec[local_blow_up]
         return local_rescaler

test/test_field.py

@@ -20,12 +20,15 @@ import unittest
 
 import numpy as np
 from numpy.testing import assert_,\
-                          assert_almost_equal
+                          assert_almost_equal,\
+                          assert_allclose
 
 from itertools import product
 
 from nifty import Field,\
-                  RGSpace
+                  RGSpace,\
+                  LMSpace,\
+                  PowerSpace
 
 from d2o import distributed_data_object
 
@@ -80,3 +83,39 @@ class Test_Functionality(unittest.TestCase):
         assert_almost_equal(a1.get_full_data(), a2.get_full_data())
         assert_almost_equal(h1.get_full_data(), h3.get_full_data())
         assert_almost_equal(a1.get_full_data(), a3.get_full_data())
+
+    @expand(product([RGSpace((8,), harmonic=True,
+                             zerocenter=False),
+                     RGSpace((8, 8), harmonic=True, distances=0.123,
+                             zerocenter=True)],
+                    [RGSpace((8,), harmonic=True,
+                             zerocenter=False),
+                     LMSpace(12)]))
+    def test_power_synthesize_analyze(self, space1, space2):
+        p1 = PowerSpace(space1)
+        spec1 = lambda k: 42/(1+k)**2
+        fp1 = Field(p1, val=spec1)
+
+        p2 = PowerSpace(space2)
+        spec2 = lambda k: 42/(1+k)**3
+        fp2 = Field(p2, val=spec2)
+
+        outer = np.outer(fp1.val.get_full_data(), fp2.val.get_full_data())
+        fp = Field((p1, p2), val=outer)
+
+        samples = 1000
+        ps1 = 0.
+        ps2 = 0.
+        for ii in xrange(samples):
+            sk = fp.power_synthesize(spaces=(0, 1), real_signal=True)
+
+            sp = sk.power_analyze(spaces=(0, 1), keep_phase_information=False)
+            ps1 += sp.sum(spaces=1)/fp2.sum()
+            ps2 += sp.sum(spaces=0)/fp1.sum()
+
+        assert_allclose(ps1.val.get_full_data()/samples,
+                        fp1.val.get_full_data(),
+                        rtol=0.1)
+        assert_allclose(ps2.val.get_full_data()/samples,
+                        fp2.val.get_full_data(),
+                        rtol=0.1)