Commit cc7f44e4 authored by Martin Reinecke's avatar Martin Reinecke

more critical filter fallout

parent d43c72e0
Pipeline #25127 passed with stages
in 5 minutes and 38 seconds
......@@ -44,7 +44,7 @@ if __name__ == "__main__":
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
mask = np.ones(s_space.shape)
mask[6000:8000] = 0.
mask[30:70,30:70] = 0.
mask = ift.Field.from_global_data(s_space, mask)
MaskOperator = ift.DiagonalOperator(mask)
......@@ -111,7 +111,9 @@ if __name__ == "__main__":
# excitation monopole to 1
m0, t0 = adjust_zero_mode(m0, t0)
ift.plot(true_sky)
ift.plot(nonlinearity(HT(power0*m0)),
title='reconstructed_sky')
ift.plot(MeasurementOperator.adjoint_times(d))
plotdict = {"colormap": "Planck-like"}
ift.plot(true_sky, title="True sky", name="true_sky.png", **plotdict)
ift.plot(nonlinearity(HT(power0*m0)), title="Reconstructed sky",
name="reconstructed_sky.png", **plotdict)
ift.plot(MeasurementOperator.adjoint_times(d), title="Data",
name="data.png", **plotdict)
......@@ -26,8 +26,7 @@ def SmoothnessOperator(domain, strength=1., logarithmic=True, space=None):
This operator applies the irregular LaplaceOperator and its adjoint to some
Field over a PowerSpace which corresponds to its smoothness and weights the
result with a scale parameter sigma. It is used in the smoothness prior
terms of the CriticalPowerEnergy. For this purpose we use free boundary
result with a scale parameter sigma. For this purpose we use free boundary
conditions in the LaplaceOperator, having no curvature at both ends. In
addition the first entry is ignored as well, corresponding to the overall
mean of the map. The mean is therefore not considered in the smoothness
......
......@@ -29,56 +29,6 @@ def _flat_PS(k):
class Energy_Tests(unittest.TestCase):
@expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[132, 42, 3]))
def testLinearPower(self, space, seed):
np.random.seed(seed)
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)
Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
xi = ift.Field.from_random(domain=hspace, random_type='normal')
def pspec(k): return 64 / (1 + k**2)**dim
pspec = ift.PS_field(pspace, pspec)
tau0 = ift.log(pspec)
A = Dist(ift.sqrt(pspec))
n = ift.Field.from_random(domain=space, random_type='normal', std=.01)
N = ift.DiagonalOperator(n**2)
s = xi * A
Instrument = ift.ScalingOperator(1., space)
R = Instrument * ht
d = R(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. / (1 + k**2))
D = ift.library.WienerFilterEnergy(position=s, d=d, R=R, N=N, S=S,
inverter=inverter).curvature
energy0 = ift.library.CriticalPowerEnergy(
position=tau0, m=s, inverter=inverter, D=D, samples=10,
smoothness_prior=1.)
energy1 = energy0.at(tau1)
a = (energy1.value - energy0.value) / eps
b = energy0.gradient.vdot(direction)
tol = 1e-4
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],
......@@ -143,61 +93,3 @@ class Energy_Tests(unittest.TestCase):
b = energy0.gradient.vdot(direction)
tol = 1e-4
assert_allclose(a, b, rtol=tol, atol=tol)
class Curvature_Tests(unittest.TestCase):
# Note: It is only possible to test the linear power 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)],
[132, 42, 3]))
def testLinearPowerCurvature(self, space, seed):
np.random.seed(seed)
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)
Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
xi = ift.Field.from_random(domain=hspace, random_type='normal')
def pspec(k): return 64 / (1 + k**2)**dim
pspec = ift.PS_field(pspace, pspec)
tau0 = ift.log(pspec)
A = Dist(ift.sqrt(pspec))
n = ift.Field.from_random(domain=space, random_type='normal', std=.01)
N = ift.DiagonalOperator(n**2)
s = xi * A
diag = ift.Field.ones(space)
Instrument = ift.DiagonalOperator(diag)
R = Instrument * ht
d = R(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=_flat_PS)
D = ift.library.WienerFilterEnergy(position=s, d=d, R=R, N=N, S=S,
inverter=inverter).curvature
energy0 = ift.library.CriticalPowerEnergy(
position=tau0, m=s, inverter=inverter, samples=10, D=D, alpha=2.)
gradient0 = energy0.gradient
gradient1 = energy0.at(tau1).gradient
a = (gradient1 - gradient0) / eps
b = energy0.curvature(direction)
tol = 1e-5
assert_allclose(a.to_global_data(), b.to_global_data(),
rtol=tol, atol=tol)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment