Even more tests. Now all energy classes are covered.

test/test_energies/test_map.py

@@ -25,17 +25,16 @@ from numpy.testing import assert_allclose
 
 
 # TODO Add also other space types
+# TODO Set tolerances to reasonable values
 
 
 class Map_Energy_Tests(unittest.TestCase):
     @expand(product([ift.RGSpace(64, distances=.789),
-                     ift.RGSpace([32, 32], distances=.789)],
-                    [ift.library.Exponential, ift.library.Linear]))
-    def testNonlinearMap(self, space, nonlinearity):
-        f = nonlinearity()
+                     ift.RGSpace([32, 32], distances=.789)]))
+    def testLinearMap(self, space):
         dim = len(space.shape)
-        fft = ift.FFTOperator(space)
-        hspace = fft.target[0]
+        hspace = space.get_default_codomain()
+        ht = ift.HarmonicTransformOperator(hspace, target=space)
         binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
         pspace = ift.PowerSpace(hspace, binbounds=binbounds)
         P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
@@ -45,32 +44,40 @@ class Map_Energy_Tests(unittest.TestCase):
         pspec = ift.PS_field(pspace, pspec)
         A = P.adjoint_times(ift.sqrt(pspec))
         n = ift.Field.from_random(domain=space, random_type='normal')
-        s = fft.inverse_times(xi0 * A)
+        s0 = xi0 * A
         diag = ift.Field.ones(space) * 10
-        R = ift.DiagonalOperator(diag)
+        Instrument = ift.DiagonalOperator(diag)
+        R = Instrument * ht
         diag = ift.Field.ones(space)
         N = ift.DiagonalOperator(diag)
-        d = R(f(s)) + n
+        d = R(s0) + n
 
         direction = ift.Field.from_random('normal', hspace)
         direction /= np.sqrt(direction.var())
         eps = 1e-10
-        xi1 = xi0 + eps * direction
+        s1 = s0 + eps * direction
+
+        IC = ift.GradientNormController(
+            name='IC',
+            verbose=False,
+            iteration_limit=100,
+            tol_abs_gradnorm=1e-5)
+        inverter = ift.ConjugateGradient(IC)
 
         S = ift.create_power_operator(hspace, power_spectrum=lambda k: 1.)
-        energy0 = ift.library.NonlinearWienerFilterEnergy(
-            position=xi0, d=d, Instrument=R, nonlinearity=f, FFT=fft, power=A, N=N, S=S)
-        energy1 = ift.library.NonlinearWienerFilterEnergy(
-            position=xi1, d=d, Instrument=R, nonlinearity=f, FFT=fft, power=A, N=N, S=S)
+        energy0 = ift.library.WienerFilterEnergy(
+            position=s0, d=d, R=R, N=N, S=S, inverter=inverter)
+        energy1 = ift.library.WienerFilterEnergy(
+            position=s1, d=d, R=R, N=N, S=S, inverter=inverter)
 
         a = (energy1.value - energy0.value) / eps
         b = energy0.gradient.vdot(direction)
-        tol = 1e-2
+        tol = 1e-3
         assert_allclose(a, b, rtol=tol, atol=tol)
 
     @expand(product([ift.RGSpace(64, distances=.789),
                      ift.RGSpace([32, 32], distances=.789)]))
-    def testLinearMap(self, space):
+    def testLognormalMap(self, space):
         dim = len(space.shape)
         hspace = space.get_default_codomain()
         ht = ift.HarmonicTransformOperator(hspace, target=space)
@@ -83,18 +90,19 @@ class Map_Energy_Tests(unittest.TestCase):
         pspec = ift.PS_field(pspace, pspec)
         A = P.adjoint_times(ift.sqrt(pspec))
         n = ift.Field.from_random(domain=space, random_type='normal')
-        s0 = xi0 * A
+        sh0 = xi0 * A
+        s = ht(sh0)
         diag = ift.Field.ones(space) * 10
         Instrument = ift.DiagonalOperator(diag)
         R = Instrument * ht
         diag = ift.Field.ones(space)
         N = ift.DiagonalOperator(diag)
-        d = R(s0) + n
+        d = Instrument(ift.exp(s)) + n
 
         direction = ift.Field.from_random('normal', hspace)
         direction /= np.sqrt(direction.var())
         eps = 1e-10
-        s1 = s0 + eps * direction
+        sh1 = sh0 + eps * direction
 
         IC = ift.GradientNormController(
             name='IC',
@@ -104,12 +112,52 @@ class Map_Energy_Tests(unittest.TestCase):
         inverter = ift.ConjugateGradient(IC)
 
         S = ift.create_power_operator(hspace, power_spectrum=lambda k: 1.)
-        energy0 = ift.library.WienerFilterEnergy(
-            position=s0, d=d, R=R, N=N, S=S, inverter=inverter)
-        energy1 = ift.library.WienerFilterEnergy(
-            position=s1, d=d, R=R, N=N, S=S, inverter=inverter)
+        energy0 = ift.library.LogNormalWienerFilterEnergy(
+            position=sh0, d=d, R=R, N=N, S=S, inverter=inverter)
+        energy1 = ift.library.LogNormalWienerFilterEnergy(
+            position=sh1, d=d, R=R, N=N, S=S, inverter=inverter)
 
         a = (energy1.value - energy0.value) / eps
         b = energy0.gradient.vdot(direction)
-        tol = 1e-3
+        tol = 1e-2
+        assert_allclose(a, b, rtol=tol, atol=tol)
+
+    @expand(product([ift.RGSpace(64, distances=.789),
+                     ift.RGSpace([32, 32], distances=.789)],
+                    [ift.library.Exponential, ift.library.Linear]))
+    def testNonlinearMap(self, space, nonlinearity):
+        f = nonlinearity()
+        dim = len(space.shape)
+        fft = ift.FFTOperator(space)
+        hspace = fft.target[0]
+        binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
+        pspace = ift.PowerSpace(hspace, binbounds=binbounds)
+        P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
+        xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
+
+        def pspec(k): return 1 / (1 + k**2)**dim
+        pspec = ift.PS_field(pspace, pspec)
+        A = P.adjoint_times(ift.sqrt(pspec))
+        n = ift.Field.from_random(domain=space, random_type='normal')
+        s = fft.inverse_times(xi0 * A)
+        diag = ift.Field.ones(space) * 10
+        R = ift.DiagonalOperator(diag)
+        diag = ift.Field.ones(space)
+        N = ift.DiagonalOperator(diag)
+        d = R(f(s)) + n
+
+        direction = ift.Field.from_random('normal', hspace)
+        direction /= np.sqrt(direction.var())
+        eps = 1e-10
+        xi1 = xi0 + eps * direction
+
+        S = ift.create_power_operator(hspace, power_spectrum=lambda k: 1.)
+        energy0 = ift.library.NonlinearWienerFilterEnergy(
+            position=xi0, d=d, Instrument=R, nonlinearity=f, FFT=fft, power=A, N=N, S=S)
+        energy1 = ift.library.NonlinearWienerFilterEnergy(
+            position=xi1, d=d, Instrument=R, nonlinearity=f, FFT=fft, power=A, N=N, S=S)
+
+        a = (energy1.value - energy0.value) / eps
+        b = energy0.gradient.vdot(direction)
+        tol = 1e-2
         assert_allclose(a, b, rtol=tol, atol=tol)