Commit c710bb41 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

tweak InversionEnabler, move generate_posterior_sample() into WienerFilterCurvature

parent e7273eea
Pipeline #22437 passed with stage
in 4 minutes and 46 seconds
import nifty2go as ift
from nifty2go.library import NonlinearWienerFilterEnergy, NonlinearPowerEnergy
from nifty2go.library.nonlinearities import *
from nifty2go.library.nonlinearities import Exponential
import numpy as np
np.random.seed(42)
......@@ -32,9 +31,10 @@ if __name__ == "__main__":
p_space = ift.PowerSpace(h_space,
binbounds=ift.PowerSpace.useful_binbounds(
h_space, logarithmic=True))
s_spec = ift.Field(p_space,val=1.)
# Choosing the prior correlation structure and defining correlation operator
p = ift.Field(p_space,val = p_spec(p_space.k_lengths))
s_spec = ift.Field(p_space, val=1.)
# Choosing the prior correlation structure and defining
# correlation operator
p = ift.Field(p_space, val=p_spec(p_space.k_lengths))
log_p = ift.log(p)
S = ift.create_power_operator(h_space, power_spectrum=s_spec)
......@@ -42,40 +42,40 @@ if __name__ == "__main__":
sp = ift.Field(p_space, val=s_spec)
sh = ift.power_synthesize(sp)
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
mask = np.ones(s_space.shape)
mask[6000:8000] = 0.
mask = ift.Field(s_space,val=mask)
mask = ift.Field(s_space, val=mask)
MaskOperator = ift.DiagonalOperator(mask)
InstrumentResponse = ift.ResponseOperator(s_space, sigma=[0.0])
MeasurementOperator = ift.ComposedOperator([MaskOperator, InstrumentResponse])
MeasurementOperator = ift.ComposedOperator([MaskOperator,
InstrumentResponse])
d_space = MeasurementOperator.target
noise_covariance = ift.Field(d_space, val=noise_level**2).weight()
N = ift.DiagonalOperator(noise_covariance)
n = ift.Field.from_random(domain=d_space,
random_type='normal',
std=noise_level,
mean=0)
Projection = ift.PowerProjectionOperator(domain= h_space, power_space=p_space)
random_type='normal',
std=noise_level,
mean=0)
Projection = ift.PowerProjectionOperator(domain=h_space,
power_space=p_space)
power = Projection.adjoint_times(ift.exp(0.5 * log_p))
# Creating the mock data
true_sky = nonlinearity(FFT.adjoint_times(power * sh))
d = MeasurementOperator(true_sky) + n
flat_power = ift.Field(p_space,val=10e-8)
flat_power = ift.Field(p_space, val=1e-7)
m0 = ift.power_synthesize(flat_power)
t0 = ift.Field(p_space,val=-4.)
t0 = ift.Field(p_space, val=-4.)
power0 = Projection.adjoint_times(ift.exp(0.5 * t0))
IC1 = ift.GradientNormController(verbose=True, iteration_limit=100,
tol_abs_gradnorm=1e-3)
minimizer = ift.RelaxedNewton(IC1)
LS = ift.LineSearchStrongWolfe(c2=0.02)
minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
ICI = ift.GradientNormController(verbose=False, name="ICI",
iteration_limit=500,
......@@ -85,22 +85,22 @@ if __name__ == "__main__":
for i in range(20):
power0 = Projection.adjoint_times(ift.exp(0.5*t0))
map0_energy = NonlinearWienerFilterEnergy(m0, d, MeasurementOperator,
nonlinearity, FFT, power0, N, S, inverter=inverter)
map0_energy = ift.library.NonlinearWienerFilterEnergy(
m0, d, MeasurementOperator, nonlinearity, FFT, power0, N, S,
inverter=inverter)
# Minimization with chosen minimizer
map0_energy, convergence = minimizer(map0_energy)
m0 = map0_energy.position
# Updating parameters for correlation structure reconstruction
D0 = map0_energy.curvature
# Initializing the power energy with updated parameters
power0_energy = NonlinearPowerEnergy(position=t0, d=d, N=N, m=m0,
D=D0, FFT=FFT, Instrument=MeasurementOperator,
nonlinearity=nonlinearity,Projection=Projection,
sigma=1., samples=2, inverter=inverter)
power0_energy = ift.library.NonlinearPowerEnergy(
position=t0, d=d, N=N, m=m0, D=D0, FFT=FFT,
Instrument=MeasurementOperator, nonlinearity=nonlinearity,
Projection=Projection, sigma=1., samples=2, inverter=inverter)
(power0_energy, convergence) = minimizer(power0_energy)
......@@ -108,9 +108,11 @@ if __name__ == "__main__":
t0_ = t0.copy()
t0 = power0_energy.position
# break degeneracy between amplitude and excitation by setting excitation monopole to 1
m0, t0 = adjust_zero_mode(m0,t0)
# break degeneracy between amplitude and excitation by setting
# excitation monopole to 1
m0, t0 = adjust_zero_mode(m0, t0)
ift.plotting.plot(true_sky)
ift.plotting.plot(nonlinearity(FFT.adjoint_times(power0*m0)),title='reconstructed_sky')
ift.plotting.plot(nonlinearity(FFT.adjoint_times(power0*m0)),
title='reconstructed_sky')
ift.plotting.plot(MeasurementOperator.adjoint_times(d))
......@@ -93,9 +93,8 @@ if __name__ == "__main__":
j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N,
R=R_harmonic)
wiener_curvature = ift.InversionEnabler(wiener_curvature, inverter)
wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R_harmonic, inverter=inverter)
m_k = wiener_curvature.inverse_times(j)
m = fft(m_k)
......
......@@ -49,9 +49,8 @@ if __name__ == "__main__":
j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N,
R=R_harmonic)
wiener_curvature = ift.InversionEnabler(wiener_curvature, inverter)
wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R_harmonic, inverter=inverter)
m_k = wiener_curvature.inverse_times(j)
m = fft(m_k)
......
......@@ -5,9 +5,10 @@ import nifty2go as ift
# Note that the constructor of PropagatorOperator takes as arguments the
# response R and noise covariance N operating on signal space and signal
# covariance operating on harmonic space.
class PropagatorOperator(ift.EndomorphicOperator):
def __init__(self, R, N, Sh):
super(PropagatorOperator, self).__init__()
class PropagatorOperator(ift.InversionEnabler, ift.EndomorphicOperator):
def __init__(self, R, N, Sh, inverter):
ift.InversionEnabler.__init__(self, inverter)
ift.EndomorphicOperator.__init__(self)
self.R = R
self.N = N
......@@ -84,7 +85,6 @@ if __name__ == "__main__":
j = R.adjoint_times(N.inverse_times(d))
IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=IC)
D = ift.InversionEnabler(PropagatorOperator(Sh=Sh, N=N, R=R),
inverter=inverter)
D = PropagatorOperator(Sh=Sh, N=N, R=R, inverter=inverter)
m = D(j)
......@@ -69,10 +69,9 @@ if __name__ == "__main__":
j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(
verbose=True, tol_abs_gradnorm=1e-4/nu.K/(nu.m**(0.5*dimensionality)))
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N,
R=R_harmonic)
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.InversionEnabler(wiener_curvature, inverter)
wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R_harmonic, inverter=inverter)
m = wiener_curvature.inverse_times(j)
m_s = fft(m)
......
......@@ -71,7 +71,7 @@ if __name__ == "__main__":
n_samples = 50
for i in range(n_samples):
sample = fft(ift.sugar.generate_posterior_sample(m, D))
sample = fft(D.generate_posterior_sample(m))
sample_variance += sample**2
sample_mean += sample
sample_mean /= n_samples
......
from ..operators.endomorphic_operator import EndomorphicOperator
from ..operators.diagonal_operator import DiagonalOperator
from ..operators import EndomorphicOperator, InversionEnabler, DiagonalOperator
class CriticalPowerCurvature(EndomorphicOperator):
class CriticalPowerCurvature(InversionEnabler, EndomorphicOperator):
"""The curvature of the CriticalPowerEnergy.
This operator implements the second derivative of the
......@@ -17,15 +16,12 @@ class CriticalPowerCurvature(EndomorphicOperator):
The smoothness prior contribution to the curvature.
"""
def __init__(self, theta, T):
super(CriticalPowerCurvature, self).__init__()
def __init__(self, theta, T, inverter):
EndomorphicOperator.__init__(self)
self._theta = DiagonalOperator(theta)
InversionEnabler.__init__(self, inverter, self._theta.inverse_times)
self._T = T
@property
def preconditioner(self):
return self._theta.inverse_times
def _times(self, x):
return self._T(x) + self._theta(x)
......
from ..minimization.energy import Energy
from ..operators.smoothness_operator import SmoothnessOperator
from ..operators.power_projection_operator import PowerProjectionOperator
from ..operators.inversion_enabler import InversionEnabler
from .critical_power_curvature import CriticalPowerCurvature
from ..utilities import memo
from .. import Field, exp
from ..sugar import generate_posterior_sample
class CriticalPowerEnergy(Energy):
......@@ -67,14 +65,12 @@ class CriticalPowerEnergy(Energy):
self._inverter = inverter
if w is None:
# self.logger.info("Initializing w")
P = PowerProjectionOperator(domain=self.m.domain,
power_space=self.position.domain[0])
if self.D is not None:
w = Field.zeros(self.position.domain, dtype=self.m.dtype)
for i in range(self.samples):
# self.logger.info("Drawing sample %i" % i)
sample = generate_posterior_sample(self.m, self.D)
sample = self.D.generate_posterior_sample(self.m)
w += P(abs(sample)**2)
w *= 1./self.samples
......@@ -82,15 +78,15 @@ class CriticalPowerEnergy(Energy):
w = P(abs(self.m)**2)
self._w = w
theta = exp(-self.position) * (self.q + self._w*0.5)
self._theta = exp(-self.position) * (self.q + self._w*0.5)
Tt = self.T(self.position)
energy = theta.integrate()
energy = self._theta.integrate()
energy += self.position.integrate()*(self.alpha-0.5)
energy += 0.5*self.position.vdot(Tt)
self._value = energy.real
gradient = -theta
gradient = -self._theta
gradient += self.alpha-0.5
gradient += Tt
self._gradient = gradient.real
......@@ -99,7 +95,7 @@ class CriticalPowerEnergy(Energy):
return self.__class__(position, self.m, D=self.D, alpha=self.alpha,
q=self.q, smoothness_prior=self.smoothness_prior,
logarithmic=self.logarithmic,
samples=self.samples, w=self.w,
samples=self.samples, w=self._w,
inverter=self._inverter)
@property
......@@ -113,8 +109,8 @@ class CriticalPowerEnergy(Energy):
@property
@memo
def curvature(self):
curv = CriticalPowerCurvature(theta=self._theta, T=self.T)
return InversionEnabler(curv, inverter=self._inverter)
return CriticalPowerCurvature(theta=self._theta, T=self.T,
inverter=self._inverter)
@property
def logarithmic(self):
......
from ..operators import EndomorphicOperator
from ..operators import EndomorphicOperator, InversionEnabler
from ..utilities import memo
from ..field import exp
class LogNormalWienerFilterCurvature(EndomorphicOperator):
class LogNormalWienerFilterCurvature(InversionEnabler, EndomorphicOperator):
"""The curvature of the LogNormalWienerFilterEnergy.
This operator implements the second derivative of the
......@@ -21,8 +21,9 @@ class LogNormalWienerFilterCurvature(EndomorphicOperator):
The prior signal covariance
"""
def __init__(self, R, N, S, position, fft4exp):
super(LogNormalWienerFilterCurvature, self).__init__()
def __init__(self, R, N, S, position, fft4exp, inverter):
InversionEnabler.__init__(self, inverter)
EndomorphicOperator.__init__(self)
self.R = R
self.N = N
self.S = S
......
......@@ -2,7 +2,6 @@ from ..minimization.energy import Energy
from ..utilities import memo
from .log_normal_wiener_filter_curvature import LogNormalWienerFilterCurvature
from ..sugar import create_composed_fft_operator
from ..operators.inversion_enabler import InversionEnabler
class LogNormalWienerFilterEnergy(Energy):
......@@ -58,11 +57,9 @@ class LogNormalWienerFilterEnergy(Energy):
@property
@memo
def curvature(self):
return InversionEnabler(
LogNormalWienerFilterCurvature(R=self.R, N=self.N, S=self.S,
position=self.position,
fft4exp=self._fft),
inverter=self._inverter)
return LogNormalWienerFilterCurvature(
R=self.R, N=self.N, S=self.S, position=self.position,
fft4exp=self._fft, inverter=self._inverter)
@property
@memo
......@@ -72,7 +69,7 @@ class LogNormalWienerFilterEnergy(Energy):
@property
@memo
def _Rexppd(self):
expp = self._fft.adjoint_times(self.curvature.op._expp_sspace)
expp = self._fft.adjoint_times(self.curvature._expp_sspace)
return self.R(expp) - self.d
@property
......@@ -84,5 +81,5 @@ class LogNormalWienerFilterEnergy(Energy):
@memo
def _exppRNRexppd(self):
return self._fft.adjoint_times(
self.curvature.op._expp_sspace *
self.curvature._expp_sspace *
self._fft(self.R.adjoint_times(self._NRexppd)))
from .. import Field, exp
from ..operators.diagonal_operator import DiagonalOperator
from ..sugar import generate_posterior_sample
from ..minimization.energy import Energy
from ..utilities import memo
......@@ -31,7 +30,7 @@ class NoiseEnergy(Energy):
sample_list.append(self.m)
else:
for i in range(samples):
sample = generate_posterior_sample(m, D)
sample = D.generate_posterior_sample(m)
sample = FFT(Field(FFT.domain, val=(
FFT.adjoint_times(sample).val)))
sample_list.append(sample)
......
from ..operators.endomorphic_operator import EndomorphicOperator
from ..operators import EndomorphicOperator, InversionEnabler
from .response_operators import LinearizedPowerResponse
class NonlinearPowerCurvature(EndomorphicOperator):
class NonlinearPowerCurvature(InversionEnabler, EndomorphicOperator):
def __init__(self, position, FFT, Instrument, nonlinearity,
Projection, N, T, sample_list):
super(NonlinearPowerCurvature, self).__init__()
Projection, N, T, sample_list, inverter):
InversionEnabler.__init__(self, inverter)
EndomorphicOperator.__init__(self)
self.N = N
self.FFT = FFT
self.Instrument = Instrument
......
from .. import exp
from ..utilities import memo
from ..operators.smoothness_operator import SmoothnessOperator
from ..sugar import generate_posterior_sample
from .nonlinear_power_curvature import NonlinearPowerCurvature
from .response_operators import LinearizedPowerResponse
from ..minimization.energy import Energy
from ..operators.inversion_enabler import InversionEnabler
class NonlinearPowerEnergy(Energy):
......@@ -56,7 +54,7 @@ class NonlinearPowerEnergy(Energy):
if samples is None or samples == 0:
sample_list = [m]
else:
sample_list = [generate_posterior_sample(m, D)
sample_list = [D.generate_posterior_sample(m)
for _ in range(samples)]
self.sample_list = sample_list
self.inverter = inverter
......@@ -105,7 +103,7 @@ class NonlinearPowerEnergy(Energy):
@property
@memo
def curvature(self):
curvature = NonlinearPowerCurvature(
return NonlinearPowerCurvature(
self.position, self.FFT, self.Instrument, self.nonlinearity,
self.Projection, self.N, self.T, self.sample_list)
return InversionEnabler(curvature, inverter=self.inverter)
self.Projection, self.N, self.T, self.sample_list,
inverter=self.inverter)
from .wiener_filter_curvature import WienerFilterCurvature
from ..utilities import memo
from ..minimization.energy import Energy
from ..operators.inversion_enabler import InversionEnabler
from .response_operators import LinearizedSignalResponse
......@@ -45,6 +44,5 @@ class NonlinearWienerFilterEnergy(Energy):
@property
@memo
def curvature(self):
curvature = WienerFilterCurvature(R=self.LinearizedResponse,
N=self.N, S=self.S)
return InversionEnabler(curvature, inverter=self.inverter)
return WienerFilterCurvature(R=self.LinearizedResponse, N=self.N,
S=self.S, inverter=self.inverter)
......@@ -42,7 +42,7 @@ class LinearizedPowerResponse(LinearOperator):
self.Instrument = Instrument
self.FFT = FFT
self.Projection = Projection
self.power = exp(0.5 * t)
self.power = exp(0.5*t)
self.m = m
position = FFT.adjoint_times(
self.Projection.adjoint_times(self.power) * self.m)
......
from ..operators import EndomorphicOperator
from ..operators import EndomorphicOperator, InversionEnabler
from ..field import Field, sqrt
from ..sugar import power_analyze, power_synthesize
class WienerFilterCurvature(EndomorphicOperator):
class WienerFilterCurvature(InversionEnabler, EndomorphicOperator):
"""The curvature of the WienerFilterEnergy.
This operator implements the second derivative of the
......@@ -19,16 +21,13 @@ class WienerFilterCurvature(EndomorphicOperator):
The prior signal covariance
"""
def __init__(self, R, N, S):
super(WienerFilterCurvature, self).__init__()
def __init__(self, R, N, S, inverter):
EndomorphicOperator.__init__(self)
InversionEnabler.__init__(self, inverter, S.times)
self.R = R
self.N = N
self.S = S
@property
def preconditioner(self):
return self.S.times
@property
def domain(self):
return self.S.domain
......@@ -45,3 +44,40 @@ class WienerFilterCurvature(EndomorphicOperator):
res = self.R.adjoint_times(self.N.inverse_times(self.R(x)))
res += self.S.inverse_times(x)
return res
def generate_posterior_sample(self, mean):
""" Generates a posterior sample from a Gaussian distribution with
given mean and covariance.
This method generates samples by setting up the observation and
reconstruction of a mock signal in order to obtain residuals of the
right correlation which are added to the given mean.
Parameters
----------
mean : Field
the mean of the posterior Gaussian distribution
Returns
-------
sample : Field
Returns the a sample from the Gaussian of given mean and
covariance.
"""
power = sqrt(power_analyze(self.S.diagonal()))
mock_signal = power_synthesize(power, real_signal=True)
noise = self.N.diagonal().weight(-1)
mock_noise = Field.from_random(random_type="normal",
domain=self.N.domain,
dtype=noise.dtype.type)
mock_noise *= sqrt(noise)
mock_data = self.R(mock_signal) + mock_noise
mock_j = self.R.adjoint_times(self.N.inverse_times(mock_data))
mock_m = self.inverse_times(mock_j)
sample = mock_signal - mock_m + mean
return sample
from ..minimization.energy import Energy
from ..operators.inversion_enabler import InversionEnabler
from .wiener_filter_curvature import WienerFilterCurvature
......@@ -28,8 +27,7 @@ class WienerFilterEnergy(Energy):
self.R = R
self.N = N
self.S = S
self._curvature = InversionEnabler(WienerFilterCurvature(R, N, S),
inverter=inverter)
self._curvature = WienerFilterCurvature(R, N, S, inverter=inverter)
self._inverter = inverter
if _j is None:
_j = self.R.adjoint_times(self.N.inverse_times(d))
......
......@@ -32,9 +32,6 @@ class ConjugateGradient(Minimizer):
----------
controller : IterationController
Object that decides when to terminate the minimization.
preconditioner : Operator *optional*
This operator can be provided which transforms the variables of the
system to improve the conditioning (default: None).
References
----------
......@@ -42,8 +39,7 @@ class ConjugateGradient(Minimizer):
2006, Springer-Verlag New York
"""
def __init__(self, controller, preconditioner=None):
self._preconditioner = preconditioner
def __init__(self, controller):
self._controller = controller
def __call__(self, energy, preconditioner=None):
......@@ -54,6 +50,9 @@ class ConjugateGradient(Minimizer):
energy : Energy object at the starting point of the iteration.
Its curvature operator must be independent of position, otherwise
linear conjugate gradient minimization will fail.
preconditioner : Operator *optional*
This operator can be provided which transforms the variables of the
system to improve the conditioning (default: None).
Returns
-------
......@@ -62,9 +61,6 @@ class ConjugateGradient(Minimizer):
status : integer
Can be controller.CONVERGED or controller.ERROR
"""
if preconditioner is None:
preconditioner = self._preconditioner
controller = self._controller
status = controller.start(energy)
if status != controller.CONTINUE:
......
......@@ -18,57 +18,29 @@
from ..minimization.quadratic_energy import QuadraticEnergy
from ..field import Field
from .linear_operator import LinearOperator
class InversionEnabler(LinearOperator):
class InversionEnabler(object):
def __init__(self, op, inverter, preconditioner=None):
self._op = op
self._inverter = inverter
if preconditioner is None and hasattr(op, "preconditioner"):
self._preconditioner = op.preconditioner
else:
self._preconditioner = preconditioner
def __init__(self, inverter, preconditioner=None):
super(InversionEnabler, self).__init__()
self._inverter = inverter