diff --git a/nifty4/extra/__init__.py b/nifty4/extra/__init__.py index be34ded633e43270be20b89d005029453de7ba7f..d0299cd62e06a2f68d7932a1290c48b97f04ea04 100644 --- a/nifty4/extra/__init__.py +++ b/nifty4/extra/__init__.py @@ -1,2 +1,2 @@ from .operator_tests import consistency_check -from .energy_tests import check_value_gradient_consistency +from .energy_tests import * diff --git a/nifty4/extra/energy_tests.py b/nifty4/extra/energy_tests.py index f4c74270e38a10a01e2db5e3739bb068f7788977..38acc18b8fc8c18d95d1917068655183542cfbb3 100644 --- a/nifty4/extra/energy_tests.py +++ b/nifty4/extra/energy_tests.py @@ -19,35 +19,63 @@ import numpy as np from ..field import Field -__all__ = ["check_value_gradient_consistency"] +__all__ = ["check_value_gradient_consistency", + "check_value_gradient_curvature_consistency"] -def check_value_gradient_consistency(E, tol=1e-6, ntries=100): +def _get_acceptable_energy(E): if not np.isfinite(E.value): raise ValueError + dir = Field.from_random("normal", E.position.domain) + # find a step length that leads to a "reasonable" energy + for i in range(50): + try: + E2 = E.at(E.position+dir) + if np.isfinite(E2.value) and abs(E2.value) < 1e20: + break + except FloatingPointError: + pass + dir *= 0.5 + else: + raise ValueError("could not find a reasonable initial step") + return E2 + + +def check_value_gradient_consistency(E, tol=1e-6, ntries=100): for _ in range(ntries): - dir = Field.from_random("normal", E.position.domain) - # find a step length that leads to a "reasonable" energy + E2 = _get_acceptable_energy(E) + dir = E2.position - E.position + Enext = E2 + dirnorm = dir.norm() + dirder = E.gradient.vdot(dir)/dirnorm for i in range(50): - try: - E2 = E.at(E.position+dir) - if np.isfinite(E2.value) and abs(E2.value) < 1e20: - break - except FloatingPointError: - pass + print(abs((E2.value-E.value)/dirnorm-dirder)) + if abs((E2.value-E.value)/dirnorm-dirder) < tol: + break dir *= 0.5 + dirnorm *= 0.5 + E2 = E2.at(E.position+dir) else: - raise ValueError("could not find a reasonable initial step") + raise ValueError("gradient and value seem inconsistent") + # E = Enext + + +def check_value_gradient_curvature_consistency(E, tol=1e-6, ntries=100): + for _ in range(ntries): + E2 = _get_acceptable_energy(E) + dir = E2.position - E.position Enext = E2 - dirder = E.gradient.vdot(dir) + dirnorm = dir.norm() + dirder = E.gradient.vdot(dir)/dirnorm + dgrad = E.curvature(dir)/dirnorm for i in range(50): - Ediff = E2.value - E.value - eps = 1e-10*max(abs(E.value), abs(E2.value)) - if abs(Ediff-dirder) < max([tol*abs(Ediff), tol*abs(dirder), eps]): + gdiff = E2.gradient - E.gradient + if abs((E2.value-E.value)/dirnorm-dirder) < tol and \ + (abs((E2.gradient-E.gradient)/dirnorm-dgrad) < tol).all(): break dir *= 0.5 - dirder *= 0.5 + dirnorm *= 0.5 E2 = E2.at(E.position+dir) else: - raise ValueError("gradient and value seem inconsistent") - E = Enext + raise ValueError("gradient, value and curvature seem inconsistent") + # E = Enext diff --git a/test/test_energies/test_map.py b/test/test_energies/test_map.py index 8e0c671dae9c306dbbdaddff5d33d4e9b0514a0b..311feb1884b5dc2e0e0937bd0c170de27d7f837b 100644 --- a/test/test_energies/test_map.py +++ b/test/test_energies/test_map.py @@ -52,11 +52,6 @@ class Energy_Tests(unittest.TestCase): N = ift.ScalingOperator(1., space) d = R(s0) + n - direction = ift.Field.from_random('normal', hspace) - direction /= np.sqrt(direction.var()) - eps = 1e-7 - s1 = s0 + eps * direction - IC = ift.GradientNormController( iteration_limit=100, tol_abs_gradnorm=1e-5) @@ -65,7 +60,8 @@ class Energy_Tests(unittest.TestCase): S = ift.create_power_operator(hspace, power_spectrum=_flat_PS) energy = ift.library.WienerFilterEnergy( position=s0, d=d, R=R, N=N, S=S, inverter=inverter) - ift.extra.check_value_gradient_consistency(energy, tol=1e-8, ntries=10) + ift.extra.check_value_gradient_curvature_consistency( + energy, tol=1e-4, ntries=10) @expand(product([ift.RGSpace(64, distances=.789), ift.RGSpace([32, 32], distances=.789)], @@ -92,120 +88,13 @@ class Energy_Tests(unittest.TestCase): N = ift.ScalingOperator(1., space) d = R(f(s)) + n - direction = ift.Field.from_random('normal', hspace) - direction /= np.sqrt(direction.var()) - eps = 1e-7 - xi1 = xi0 + eps * direction - S = ift.create_power_operator(hspace, power_spectrum=_flat_PS) energy = ift.library.NonlinearWienerFilterEnergy( position=xi0, d=d, Instrument=R, nonlinearity=f, ht=ht, power=A, N=N, S=S) - ift.extra.check_value_gradient_consistency(energy, tol=1e-8, ntries=10) - - -class Curvature_Tests(unittest.TestCase): - # Note: It is only possible to test linear curvatures since the non-linear - # curvatures are not the exact second derivative but only a part of it. One - # term is neglected which would render the second derivative non-positive - # definite. - @expand(product([ift.RGSpace(64, distances=.789), - ift.RGSpace([32, 32], distances=.789)], - [4, 78, 23])) - def testLinearMapCurvature(self, space, seed): - np.random.seed(seed) - dim = len(space.shape) - 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) - Dist = ift.PowerDistributor(target=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 = Dist(ift.sqrt(pspec)) - n = ift.Field.from_random(domain=space, random_type='normal') - s0 = xi0 * A - Instrument = ift.ScalingOperator(10., space) - R = Instrument * ht - N = ift.ScalingOperator(1., space) - d = R(s0) + n - - direction = ift.Field.from_random('normal', hspace) - direction /= np.sqrt(direction.var()) - eps = 1e-7 - s1 = s0 + eps * direction - - IC = ift.GradientNormController( - iteration_limit=100, - tol_abs_gradnorm=1e-5) - inverter = ift.ConjugateGradient(IC) - - S = ift.create_power_operator(hspace, power_spectrum=_flat_PS) - energy0 = ift.library.WienerFilterEnergy( - position=s0, d=d, R=R, N=N, S=S, inverter=inverter) - gradient0 = energy0.gradient - gradient1 = energy0.at(s1).gradient - - a = (gradient1 - gradient0) / eps - b = energy0.curvature(direction) - tol = 1e-7 - assert_allclose(a.to_global_data(), b.to_global_data(), - rtol=tol, atol=tol) - - @expand(product([ift.RGSpace(64, distances=.789), - ift.RGSpace([32, 32], distances=.789)], - # Only linear case due to approximation of Hessian in the - # case of nontrivial nonlinearities. - [ift.library.Linear], - [4, 78, 23])) - def testNonlinearMapCurvature(self, space, nonlinearity, seed): - np.random.seed(seed) - f = nonlinearity() - dim = len(space.shape) - 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) - Dist = ift.PowerDistributor(target=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 = Dist(ift.sqrt(pspec)) - n = ift.Field.from_random(domain=space, random_type='normal') - s = ht(xi0 * A) - R = ift.ScalingOperator(10., space) - N = ift.ScalingOperator(1., space) - d = R(f(s)) + n - - direction = ift.Field.from_random('normal', hspace) - direction /= np.sqrt(direction.var()) - eps = 1e-7 - xi1 = xi0 + eps * direction - - S = ift.create_power_operator(hspace, power_spectrum=_flat_PS) - - IC = ift.GradientNormController( - iteration_limit=500, - tol_abs_gradnorm=1e-7) - inverter = ift.ConjugateGradient(IC) - energy0 = ift.library.NonlinearWienerFilterEnergy( - position=xi0, - d=d, - Instrument=R, - nonlinearity=f, - ht=ht, - power=A, - N=N, - S=S, - inverter=inverter) - gradient0 = energy0.gradient - gradient1 = energy0.at(xi1).gradient - - a = (gradient1 - gradient0) / eps - b = energy0.curvature(direction) - tol = 1e-7 - assert_allclose(a.to_global_data(), b.to_global_data(), - rtol=tol, atol=tol) + if isinstance(nonlinearity, ift.library.Linear): + ift.extra.check_value_gradient_curvature_consistency( + energy, tol=1e-4, ntries=10) + else: + ift.extra.check_value_gradient_consistency( + energy, tol=1e-4, ntries=10) diff --git a/test/test_energies/test_noise.py b/test/test_energies/test_noise.py index 324bf750e447ae3d9924b58d0690f271aaeaab4d..9a4253452b357ee5bb546f2cde6befa77cdbb66e 100644 --- a/test/test_energies/test_noise.py +++ b/test/test_energies/test_noise.py @@ -62,11 +62,6 @@ class Noise_Energy_Tests(unittest.TestCase): alpha = ift.Field.full(d.domain, 2.) q = ift.Field.full(d.domain, 1e-5) - direction = ift.Field.from_random('normal', d.domain) - direction /= np.sqrt(direction.var()) - eps = 1e-8 - eta1 = eta0 + eps * direction - IC = ift.GradientNormController( iteration_limit=100, tol_abs_gradnorm=1e-5) @@ -88,4 +83,4 @@ class Noise_Energy_Tests(unittest.TestCase): for _ in range(10)] energy = ift.library.NoiseEnergy(eta0, alpha, q, res_sample_list) - ift.extra.check_value_gradient_consistency(energy, tol=1e-8, ntries=10) + ift.extra.check_value_gradient_consistency(energy, tol=1e-6, ntries=10) diff --git a/test/test_energies/test_power.py b/test/test_energies/test_power.py index 5c05ac456badb0d5d097e9690ad98bb6f796a571..1062f65fc851bb70cdbc830e453089474538513e 100644 --- a/test/test_energies/test_power.py +++ b/test/test_energies/test_power.py @@ -54,11 +54,6 @@ class Energy_Tests(unittest.TestCase): N = ift.DiagonalOperator(diag) d = R(f(s)) + n - direction = ift.Field.from_random('normal', pspace) - direction /= np.sqrt(direction.var()) - eps = 1e-7 - tau1 = tau0 + eps * direction - IC = ift.GradientNormController( iteration_limit=100, tol_abs_gradnorm=1e-5) @@ -87,4 +82,4 @@ class Energy_Tests(unittest.TestCase): ht=ht, N=N, samples=10) - ift.extra.check_value_gradient_consistency(energy, tol=1e-8, ntries=10) + ift.extra.check_value_gradient_consistency(energy, tol=1e-5, ntries=10)