diff --git a/nifty5/extra.py b/nifty5/extra.py index 1fccf4a48e4b7dbeeca9e8402f51d7b46c07336a..fb83ba415de4047428c3ed53648d78cf02edaa37 100644 --- a/nifty5/extra.py +++ b/nifty5/extra.py @@ -22,8 +22,7 @@ from .linearization import Linearization from .operators.linear_operator import LinearOperator from .sugar import from_random -__all__ = ["consistency_check", "check_value_gradient_consistency", - "check_value_gradient_metric_consistency"] +__all__ = ["consistency_check", "check_jacobian_consistency"] def _assert_allclose(f1, f2, atol, rtol): @@ -92,12 +91,12 @@ def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64, The data type of the random vectors in the operator's target. Default is `np.float64`. atol : float - Absolute tolerance for the check. If rtol is specified, - then satisfying any tolerance will let the check pass. + Absolute tolerance for the check. If rtol is specified, + then satisfying any tolerance will let the check pass. Default: 0. rtol : float - Relative tolerance for the check. If atol is specified, - then satisfying any tolerance will let the check pass. + Relative tolerance for the check. If atol is specified, + then satisfying any tolerance will let the check pass. Default: 0. """ if not isinstance(op, LinearOperator): @@ -134,25 +133,37 @@ def _get_acceptable_location(op, loc, lin): return loc2, lin2 -def _check_consistency(op, loc, tol, ntries, do_metric): +def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100): + """ + Checks the Jacobian of an operator against its finite difference + approximation. + + Computes the Jacobian with finite differences and compares it to the + implemented Jacobian. + + Parameters + ---------- + op : Operator + Operator which shall be checked. + loc : Field or MultiField + An Field or MultiField instance which has the same domain + as op. The location at which the gradient is checked + tol : float + Tolerance for the check. + """ for _ in range(ntries): - lin = op(Linearization.make_var(loc, do_metric)) + lin = op(Linearization.make_var(loc)) loc2, lin2 = _get_acceptable_location(op, loc, lin) dir = loc2-loc locnext = loc2 dirnorm = dir.norm() for i in range(50): locmid = loc + 0.5*dir - linmid = op(Linearization.make_var(locmid, do_metric)) + linmid = op(Linearization.make_var(locmid)) dirder = linmid.jac(dir) numgrad = (lin2.val-lin.val) xtol = tol * dirder.norm() / np.sqrt(dirder.size) - cond = (abs(numgrad-dirder) <= xtol).all() - if do_metric: - dgrad = linmid.metric(dir) - dgrad2 = (lin2.gradient-lin.gradient) - cond = cond and (abs(dgrad-dgrad2) <= xtol).all() - if cond: + if (abs(numgrad-dirder) <= xtol).all(): break dir = dir*0.5 dirnorm *= 0.5 @@ -161,33 +172,3 @@ def _check_consistency(op, loc, tol, ntries, do_metric): raise ValueError("gradient and value seem inconsistent") loc = locnext - -def check_value_gradient_consistency(op, loc, tol=1e-8, ntries=100): - """ - Checks the gradient (jacobian) of an operator against its value. - - Computes the gradient (jacobian) with finite differences and compares - it to the implemented gradient (jacobian). - - Parameters - ---------- - op : Operator - Operator which shall be checked. - loc : Field or MultiField - An Field or MultiField instance which has the same domain - as op. The location at which the gradient is checked - atol : float - Absolute tolerance for the check. If rtol is specified, - then satisfying any tolerance will let the check pass. - Default: 0. - rtol : float - Relative tolerance for the check. If atol is specified, - then satisfying any tolerance will let the check pass. - Default: 0 - """ - _check_consistency(op, loc, tol, ntries, False) - - -def check_value_gradient_metric_consistency(op, loc, tol=1e-8, ntries=100): - """FIXME""" - _check_consistency(op, loc, tol, ntries, True) diff --git a/test/test_energy_gradients.py b/test/test_energy_gradients.py index 129e40292cb1d4ff9aad3edf62c8cc945497e9f1..cd4d1b50e43cc90a5f40b6e19d84a605a3f5b5a7 100644 --- a/test/test_energy_gradients.py +++ b/test/test_energy_gradients.py @@ -44,7 +44,7 @@ def field(request): def test_gaussian(field): energy = ift.GaussianEnergy(domain=field.domain) - ift.extra.check_value_gradient_consistency(energy, field) + ift.extra.check_jacobian_consistency(energy, field) def test_inverse_gamma(field): @@ -53,7 +53,7 @@ def test_inverse_gamma(field): d = np.random.normal(10, size=space.shape)**2 d = ift.Field.from_global_data(space, d) energy = ift.InverseGammaLikelihood(d) - ift.extra.check_value_gradient_consistency(energy, field, tol=1e-7) + ift.extra.check_jacobian_consistency(energy, field, tol=1e-7) def testPoissonian(field): @@ -62,7 +62,7 @@ def testPoissonian(field): d = np.random.poisson(120, size=space.shape) d = ift.Field.from_global_data(space, d) energy = ift.PoissonianEnergy(d) - ift.extra.check_value_gradient_consistency(energy, field, tol=1e-7) + ift.extra.check_jacobian_consistency(energy, field, tol=1e-7) def test_hamiltonian_and_KL(field): @@ -70,11 +70,11 @@ def test_hamiltonian_and_KL(field): space = field.domain lh = ift.GaussianEnergy(domain=space) hamiltonian = ift.StandardHamiltonian(lh) - ift.extra.check_value_gradient_consistency(hamiltonian, field) + ift.extra.check_jacobian_consistency(hamiltonian, field) S = ift.ScalingOperator(1., space) samps = [S.draw_sample() for i in range(3)] kl = ift.AveragedEnergy(hamiltonian, samps) - ift.extra.check_value_gradient_consistency(kl, field) + ift.extra.check_jacobian_consistency(kl, field) def test_bernoulli(field): @@ -83,4 +83,4 @@ def test_bernoulli(field): d = np.random.binomial(1, 0.1, size=space.shape) d = ift.Field.from_global_data(space, d) energy = ift.BernoulliEnergy(d) - ift.extra.check_value_gradient_consistency(energy, field, tol=1e-6) + ift.extra.check_jacobian_consistency(energy, field, tol=1e-6) diff --git a/test/test_gaussian_metric.py b/test/test_gaussian_energy.py similarity index 85% rename from test/test_gaussian_metric.py rename to test/test_gaussian_energy.py index aa9f4beb76ef8ee90fb8dde44daacf268202cdb5..e77604ed3297fa34f1ecaeec825460ce04b9fd66 100644 --- a/test/test_gaussian_metric.py +++ b/test/test_gaussian_energy.py @@ -45,8 +45,6 @@ def test_gaussian_energy(space, nonlinearity, noise, seed): pspace = ift.PowerSpace(hspace, binbounds=binbounds) Dist = ift.PowerDistributor(target=hspace, power_space=pspace) xi0 = ift.Field.from_random(domain=hspace, random_type='normal') - # FIXME Needed? - xi0_var = ift.Linearization.make_var(xi0) def pspec(k): return 1/(1 + k**2)**dim @@ -55,8 +53,6 @@ def test_gaussian_energy(space, nonlinearity, noise, seed): A = Dist(ift.sqrt(pspec)) N = ift.ScalingOperator(noise, space) n = N.draw_sample() - # FIXME Needed? - s = ht(ift.makeOp(A)(xi0_var)) R = ift.ScalingOperator(10., space) def d_model(): @@ -73,9 +69,5 @@ def test_gaussian_energy(space, nonlinearity, noise, seed): N = None energy = ift.GaussianEnergy(d, N)(d_model()) - if nonlinearity == "": - ift.extra.check_value_gradient_metric_consistency( - energy, xi0, ntries=10) - else: - ift.extra.check_value_gradient_consistency( - energy, xi0, ntries=10, tol=5e-8) + ift.extra.check_jacobian_consistency( + energy, xi0, ntries=10, tol=5e-8) diff --git a/test/test_operators/test_gradients.py b/test/test_operators/test_jacobian.py similarity index 82% rename from test/test_operators/test_gradients.py rename to test/test_operators/test_jacobian.py index c988c09ea83848fab7e8f33de03c3dbfae17c880..4b825b064a26dcb8c313461bc51b4879bfc75303 100644 --- a/test/test_operators/test_gradients.py +++ b/test/test_operators/test_jacobian.py @@ -47,7 +47,7 @@ def _make_linearization(type, space, seed): def testBasics(space, seed): var = _make_linearization("Variable", space, seed) model = ift.ScalingOperator(6., var.target) - ift.extra.check_value_gradient_consistency(model, var.val) + ift.extra.check_jacobian_consistency(model, var.val) @pmp('type1', ['Variable', 'Constant']) @@ -65,32 +65,32 @@ def testBinary(type1, type2, space, seed): select_s2 = ift.ducktape(None, dom, "s2") model = select_s1*select_s2 pos = ift.from_random("normal", dom) - ift.extra.check_value_gradient_consistency(model, pos, ntries=20) + ift.extra.check_jacobian_consistency(model, pos, ntries=20) model = select_s1 + select_s2 pos = ift.from_random("normal", dom) - ift.extra.check_value_gradient_consistency(model, pos, ntries=20) + ift.extra.check_jacobian_consistency(model, pos, ntries=20) model = select_s1.scale(3.) pos = ift.from_random("normal", dom1) - ift.extra.check_value_gradient_consistency(model, pos, ntries=20) + ift.extra.check_jacobian_consistency(model, pos, ntries=20) model = ift.ScalingOperator(2.456, space)(select_s1*select_s2) pos = ift.from_random("normal", dom) - ift.extra.check_value_gradient_consistency(model, pos, ntries=20) + ift.extra.check_jacobian_consistency(model, pos, ntries=20) model = ift.sigmoid(2.456*(select_s1*select_s2)) pos = ift.from_random("normal", dom) - ift.extra.check_value_gradient_consistency(model, pos, ntries=20) + ift.extra.check_jacobian_consistency(model, pos, ntries=20) pos = ift.from_random("normal", dom) model = ift.OuterProduct(pos['s1'], ift.makeDomain(space)) - ift.extra.check_value_gradient_consistency(model, pos['s2'], ntries=20) + ift.extra.check_jacobian_consistency(model, pos['s2'], ntries=20) model = select_s1**2 pos = ift.from_random("normal", dom1) - ift.extra.check_value_gradient_consistency(model, pos, ntries=20) + ift.extra.check_jacobian_consistency(model, pos, ntries=20) model = select_s1.clip(-1, 1) pos = ift.from_random("normal", dom1) - ift.extra.check_value_gradient_consistency(model, pos, ntries=20) + ift.extra.check_jacobian_consistency(model, pos, ntries=20) if isinstance(space, ift.RGSpace): model = ift.FFTOperator(space)(select_s1*select_s2) pos = ift.from_random("normal", dom) - ift.extra.check_value_gradient_consistency(model, pos, ntries=20) + ift.extra.check_jacobian_consistency(model, pos, ntries=20) def testModelLibrary(space, seed): @@ -102,18 +102,18 @@ def testModelLibrary(space, seed): assert_(isinstance(model, ift.Operator)) S = ift.ScalingOperator(1., model.domain) pos = S.draw_sample() - ift.extra.check_value_gradient_consistency(model, pos, ntries=20) + ift.extra.check_jacobian_consistency(model, pos, ntries=20) model2 = ift.CorrelatedField(space, model) S = ift.ScalingOperator(1., model2.domain) pos = S.draw_sample() - ift.extra.check_value_gradient_consistency(model2, pos, ntries=20) + ift.extra.check_jacobian_consistency(model2, pos, ntries=20) domtup = ift.DomainTuple.make((space, space)) model3 = ift.MfCorrelatedField(domtup, [model, model]) S = ift.ScalingOperator(1., model3.domain) pos = S.draw_sample() - ift.extra.check_value_gradient_consistency(model3, pos, ntries=20) + ift.extra.check_jacobian_consistency(model3, pos, ntries=20) def testPointModel(space, seed): @@ -123,7 +123,7 @@ def testPointModel(space, seed): q = 0.73 model = ift.InverseGammaOperator(space, alpha, q) # FIXME All those cdfs and ppfs are not very accurate - ift.extra.check_value_gradient_consistency(model, pos, tol=1e-2, ntries=20) + ift.extra.check_jacobian_consistency(model, pos, tol=1e-2, ntries=20) @pmp('target', [ @@ -148,7 +148,7 @@ def testDynamicModel(target, causal, minimum_phase, seed): S = ift.ScalingOperator(1., model.domain) pos = S.draw_sample() # FIXME I dont know why smaller tol fails for 3D example - ift.extra.check_value_gradient_consistency(model, pos, tol=1e-5, ntries=20) + ift.extra.check_jacobian_consistency(model, pos, tol=1e-5, ntries=20) if len(target.shape) > 1: dct = { 'target': target, @@ -169,5 +169,5 @@ def testDynamicModel(target, causal, minimum_phase, seed): S = ift.ScalingOperator(1., model.domain) pos = S.draw_sample() # FIXME I dont know why smaller tol fails for 3D example - ift.extra.check_value_gradient_consistency( + ift.extra.check_jacobian_consistency( model, pos, tol=1e-5, ntries=20)