more tests

test/test_field.py

@@ -19,8 +19,7 @@
 import unittest
 
 import numpy as np
-from numpy.testing import assert_,\
-                          assert_almost_equal,\
+from numpy.testing import assert_equal,\
                           assert_allclose
 from nose.plugins.skip import SkipTest
 
@@ -32,6 +31,8 @@ from nifty2go import Field,\
                   PowerSpace,\
                   DomainTuple
 
+from nifty2go.sugar import create_power_operator
+
 from test.common import expand
 
 
@@ -50,7 +51,7 @@ class Test_Interface(unittest.TestCase):
         attribute = attribute_desired_type[0]
         desired_type = attribute_desired_type[1]
         f = Field(domain=domain)
-        assert_(isinstance(getattr(f, attribute), desired_type))
+        assert_equal(isinstance(getattr(f, attribute), desired_type), True)
 
 
 class Test_Functionality(unittest.TestCase):
@@ -85,6 +86,41 @@ class Test_Functionality(unittest.TestCase):
         assert_allclose(ps1.val/samples, fp1.val, rtol=0.2)
         assert_allclose(ps2.val/samples, fp2.val, rtol=0.2)
 
+
+    @expand(product([RGSpace((8,), harmonic=True),
+                     RGSpace((8, 8), harmonic=True, distances=0.123)],
+                    [RGSpace((8,), harmonic=True),
+                     LMSpace(12)]))
+    def test_DiagonalOperator_power_analyze(self, space1, space2):
+        np.random.seed(11)
+
+        p1 = PowerSpace(space1)
+        spec1 = lambda k: 42/(1+k)**2
+        fp1 = Field(p1, val=spec1(p1.kindex))
+
+        p2 = PowerSpace(space2)
+        spec2 = lambda k: 42/(1+k)**3
+        fp2 = Field(p2, val=spec2(p2.kindex))
+
+        S_1 = create_power_operator(space1, lambda x: np.sqrt(spec1(x)))
+        S_2 = create_power_operator(space2, lambda x: np.sqrt(spec2(x)))
+        S_1.set_diagonal(S_1.diagonal().weight(-1))
+        S_2.set_diagonal(S_2.diagonal().weight(-1))
+
+        samples = 500
+        ps1 = 0.
+        ps2 = 0.
+
+        for ii in range(samples):
+            rand_k = Field.from_random('normal', domain=(space1, space2))
+            sk = S_1.times(S_2.times(rand_k, spaces=1), spaces=0)
+            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/samples, fp1.val, rtol=0.2)
+        assert_allclose(ps2.val/samples, fp2.val, rtol=0.2)
+
     def test_vdot(self):
         s=RGSpace((10,))
         f1=Field.from_random("normal",domain=s,dtype=np.complex128)