From 1c78f918188b977ab69bcc1d15ee223bc791c5e8 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Tue, 30 Jan 2018 13:57:49 +0100
Subject: [PATCH] Energy tests ready to merge
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
Sadly, it is not possible to test the curvatures other than in the linear case.
We have failed. The reason is that the curvature is approximated in a way that
it is positive definite. For details see Knollmüller and Ensslin 2017.
---
test/test_energies/test_map.py | 55 +++---------------------
test/test_energies/test_power.py | 74 +-------------------------------
2 files changed, 8 insertions(+), 121 deletions(-)
diff --git a/test/test_energies/test_map.py b/test/test_energies/test_map.py
index 2d5e18765..c4ebebb23 100644
--- a/test/test_energies/test_map.py
+++ b/test/test_energies/test_map.py
@@ -144,11 +144,11 @@ class Energy_Tests(unittest.TestCase):
S = ift.create_power_operator(hspace, power_spectrum=lambda k: 1.)
energy0 = ift.library.NonlinearWienerFilterEnergy(
position=xi0, d=d, Instrument=R, nonlinearity=f, ht=ht, power=A, N=N, S=S)
- energy1 = energy0.at(xi0)
+ energy1 = energy0.at(xi1)
a = (energy1.value - energy0.value) / eps
b = energy0.gradient.vdot(direction)
- tol = 1e-4
+ tol = 1e-5
assert_allclose(a, b, rtol=tol, atol=tol)
@@ -197,55 +197,10 @@ class Curvature_Tests(unittest.TestCase):
tol = 1e-7
assert_allclose(a.val, b.val, rtol=tol, atol=tol)
- @expand(product([ift.RGSpace(64, distances=.789),
- ift.RGSpace([32, 32], distances=.789)],
- [4, 78, 23]))
- def testLognormalMapCurvature(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)
- 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')
- sh0 = xi0 * A
- s = ht(sh0)
- Instrument = ift.ScalingOperator(10., space)
- R = Instrument * ht
- N = ift.ScalingOperator(1., space)
- d = Instrument(ift.exp(s)) + n
-
- direction = ift.Field.from_random('normal', hspace)
- direction /= np.sqrt(direction.var())
- eps = 1e-7
- sh1 = sh0 + 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=lambda k: 1.)
-
- energy0 = ift.library.LogNormalWienerFilterEnergy(
- position=sh0, d=d, R=R, N=N, S=S, inverter=inverter)
- gradient0 = energy0.gradient
- gradient1 = energy0.at(sh1).gradient
-
- a = (gradient1 - gradient0) / eps
- b = energy0.curvature(direction)
- tol = 1e-3
- assert_allclose(a.val, b.val, rtol=tol, atol=tol)
@expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
- [ift.library.Exponential, ift.library.Linear],
+ [ift.library.Linear], # Only linear case due to approximation of Hessian in the case of nontrivial nonlinearities.
[4, 78, 23]))
def testNonlinearMapCurvature(self, space, nonlinearity, seed):
np.random.seed(seed)
@@ -293,5 +248,7 @@ class Curvature_Tests(unittest.TestCase):
a = (gradient1 - gradient0) / eps
b = energy0.curvature(direction)
- tol = 1e-3
+ print(a.vdot(a))
+ print(b.vdot(b))
+ tol = 1e-7
assert_allclose(a.val, b.val, rtol=tol, atol=tol)
diff --git a/test/test_energies/test_power.py b/test/test_energies/test_power.py
index 0c45685a1..eb617ea52 100644
--- a/test/test_energies/test_power.py
+++ b/test/test_energies/test_power.py
@@ -59,15 +59,13 @@ class Energy_Tests(unittest.TestCase):
tol_abs_gradnorm=1e-5)
inverter = ift.ConjugateGradient(IC)
- S = ift.create_power_operator(hspace, power_spectrum=lambda k: 1.)
+ S = ift.create_power_operator(hspace, power_spectrum=lambda k: 1./(1+k**2))
D = ift.library.WienerFilterEnergy(position=s, d=d, R=R, N=N, S=S,
inverter=inverter).curvature
- w = ift.Field.zeros_like(tau0)
-
energy0 = ift.library.CriticalPowerEnergy(
- position=tau0, m=s, inverter=inverter, w=w, samples=10)
+ position=tau0, m=s, inverter=inverter, D=D, samples=10, smoothness_prior=1.)
energy1 = energy0.at(tau1)
a = (energy1.value - energy0.value) / eps
@@ -192,71 +190,3 @@ class Curvature_Tests(unittest.TestCase):
b = energy0.curvature(direction)
tol = 1e-5
assert_allclose(a.val, b.val, rtol=tol, atol=tol)
-
- @expand(product([ift.RGSpace(64, distances=.789),
- ift.RGSpace([32, 32], distances=.789)],
- [ift.library.Exponential, ift.library.Linear],
- [132, 42, 3]))
- def testNonlinearPowerCurvature(self, space, nonlinearity, seed):
- np.random.seed(seed)
- f = nonlinearity()
- dim = len(space.shape)
- hspace = space.get_default_codomain()
- ht = ift.HarmonicTransformOperator(hspace, space)
- binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=True)
- pspace = ift.PowerSpace(hspace, binbounds=binbounds)
- P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
- xi = ift.Field.from_random(domain=hspace, random_type='normal')
-
- def pspec(k): return 1 / (1 + k**2)**dim
- tau0 = ift.PS_field(pspace, pspec)
- A = P.adjoint_times(ift.sqrt(tau0))
- n = ift.Field.from_random(domain=space, random_type='normal')
- s = ht(xi * 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', pspace)
- direction /= np.sqrt(direction.var())
- eps = 1e-7
- tau1 = tau0 + 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=lambda k: 1.)
- D = ift.library.NonlinearWienerFilterEnergy(
- position=xi,
- d=d,
- Instrument=R,
- nonlinearity=f,
- power=A,
- N=N,
- S=S,
- ht=ht,
- inverter=inverter).curvature
-
- energy0 = ift.library.NonlinearPowerEnergy(
- position=tau0,
- d=d,
- xi=xi,
- D=D,
- Instrument=R,
- Projection=P,
- nonlinearity=f,
- ht=ht,
- N=N,
- samples=10)
-
- gradient0 = energy0.gradient
- gradient1 = energy0.at(tau1).gradient
-
- a = (gradient1 - gradient0) / eps
- b = energy0.curvature(direction)
- tol = 1e-3
- assert_allclose(a.val, b.val, rtol=tol, atol=tol)
--
GitLab