diff --git a/test/test_energies/test_map.py b/test/test_energies/test_map.py index 5e56922da0a5315f6197b1dc22f5164fd05636e5..7338a0249e707b7d5266aefda425d418ef2fe217 100644 --- a/test/test_energies/test_map.py +++ b/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)