Commit 3ccb329c authored by Martin Reinecke's avatar Martin Reinecke
Browse files

cosmetics

parent efdfcfff
Pipeline #22317 passed with stage
in 4 minutes and 46 seconds
......@@ -8,4 +8,4 @@ from .noise_curvature import NoiseCurvature
from .noise_energy import NoiseEnergy
from .nonlinear_power_curvature import NonlinearPowerCurvature
from .nonlinear_power_energy import NonlinearPowerEnergy
from .nonlinear_signal_energy import NonlinearWienerFilterEnergy
from .nonlinear_wiener_filter_energy import NonlinearWienerFilterEnergy
from .. import Field
from ..operators.endomorphic_operator import EndomorphicOperator
from .response_operators import LinearizedPowerResponse
......@@ -7,6 +6,7 @@ class NonlinearPowerCurvature(EndomorphicOperator):
def __init__(self, position, FFT, Instrument, nonlinearity,
Projection, N, T, sample_list):
super(NonlinearPowerCurvature, self).__init__()
self.N = N
self.FFT = FFT
self.Instrument = Instrument
......@@ -16,10 +16,6 @@ class NonlinearPowerCurvature(EndomorphicOperator):
self.Projection = Projection
self.nonlinearity = nonlinearity
# if preconditioner is None:
# preconditioner = self.theta.inverse_times
super(NonlinearPowerCurvature, self).__init__()
@property
def domain(self):
return self.position.domain
......@@ -43,7 +39,8 @@ class NonlinearPowerCurvature(EndomorphicOperator):
return result + self.T(x)
def _sample_times(self, x, sample):
LinearizedResponse = LinearizedPowerResponse(self.Instrument, self.nonlinearity,
self.FFT, self.Projection, self.position, sample)
LinearizedResponse = LinearizedPowerResponse(
self.Instrument, self.nonlinearity, self.FFT, self.Projection,
self.position, sample)
return LinearizedResponse.adjoint_times(
self.N.inverse_times(LinearizedResponse(x)))
from .. import Field, exp
from .. import exp
from ..utilities import memo
from ..operators.smoothness_operator import SmoothnessOperator
from ..sugar import generate_posterior_sample
......@@ -11,9 +11,10 @@ from ..operators.inversion_enabler import InversionEnabler
class NonlinearPowerEnergy(Energy):
"""The Energy of the power spectrum according to the critical filter.
It describes the energy of the logarithmic amplitudes of the power spectrum of
a Gaussian random field under reconstruction uncertainty with smoothness and
inverse gamma prior. It is used to infer the correlation structure of a correlated signal.
It describes the energy of the logarithmic amplitudes of the power spectrum
of a Gaussian random field under reconstruction uncertainty with smoothness
and inverse gamma prior. It is used to infer the correlation structure of a
correlated signal.
Parameters
----------
......@@ -29,47 +30,50 @@ class NonlinearPowerEnergy(Energy):
The parameter of the smoothness prior.
default : ??? None? ???????
samples : integer
Number of samples used for the estimation of the uncertainty corrections.
Number of samples used for the estimation of the uncertainty
corrections.
default : 3
"""
def __init__(self, position, d, N, m, D, FFT, Instrument, nonlinearity, Projection,
sigma=0., samples=3, sample_list=None, inverter=None):
super(NonlinearPowerEnergy, self).__init__(position=position.copy())
dummy = self.position.norm()
def __init__(self, position, d, N, m, D, FFT, Instrument, nonlinearity,
Projection, sigma=0., samples=3, sample_list=None,
inverter=None):
super(NonlinearPowerEnergy, self).__init__(position)
self.m = m
self.D = D
self.d = d
self.N = N
self.sigma = sigma
self.T = SmoothnessOperator(domain=self.position.domain[0], strength=self.sigma,
logarithmic=True)
self.T = SmoothnessOperator(domain=self.position.domain[0],
strength=self.sigma, logarithmic=True)
self.samples = samples
self.FFT = FFT
self.Instrument = Instrument
self.nonlinearity = nonlinearity
self.Projection = Projection
self.LinearizedResponse = LinearizedPowerResponse(self.Instrument, self.nonlinearity,
self.FFT, self.Projection, self.position, self.m)
self.LinearizedResponse = LinearizedPowerResponse(
self.Instrument, self.nonlinearity, self.FFT, self.Projection,
self.position, self.m)
self.power = self.Projection.adjoint_times(
exp(0.5 * self.position))
self.power = self.Projection.adjoint_times(exp(0.5*self.position))
if sample_list is None:
sample_list = []
if samples is None:
sample_list.append(self.m)
if samples is None or samples == 0:
sample_list = [self.m]
else:
for i in range(samples):
sample = generate_posterior_sample(m, D)
sample = FFT(Field(FFT.domain, val=(
FFT.adjoint_times(sample).val)))
sample_list.append(sample)
sample_list = [generate_posterior_sample(m, D)
for _ in range(samples)]
# sample_list = []
# for i in range(samples):
# sample = generate_posterior_sample(m, D)
# # MR FIXME: really needed? For RGSpaces this is a no-op
# sample = FFT(FFT.adjoint_times(sample))
# sample_list.append(sample)
self.sample_list = sample_list
self.inverter = inverter
def at(self, position):
return self.__class__(position, self.d, self.N, self.m,
self.D, self.FFT, self.Instrument, self.nonlinearity,
return self.__class__(position, self.d, self.N, self.m, self.D,
self.FFT, self.Instrument, self.nonlinearity,
self.Projection, sigma=self.sigma,
sample_list=self.sample_list,
samples=self.samples, inverter=self.inverter)
......@@ -80,36 +84,40 @@ class NonlinearPowerEnergy(Energy):
likelihood = 0.
for sample in self.sample_list:
likelihood += self._likelihood(sample)
return 0.5 * self.position.vdot(self.T(self.position)) + likelihood / float(len(self.sample_list))
likelihood /= float(len(self.sample_list))
return 0.5*self.position.vdot(self.T(self.position)) + likelihood
def _likelihood(self, m):
residual = self.d - \
self.Instrument(self.nonlinearity(
self.FFT.adjoint_times(self.power * m)))
energy = 0.5 * residual.vdot(self.N.inverse_times(residual))
return energy.real
self.FFT.adjoint_times(self.power*m)))
return 0.5 * residual.vdot(self.N.inverse_times(residual)).real
@property
@memo
def gradient(self):
likelihood_gradient = Field(self.position.domain, val=0.)
likelihood_gradient = None
for sample in self.sample_list:
likelihood_gradient += self._likelihood_gradient(sample)
return -likelihood_gradient / float(len(self.sample_list)) + self.T(self.position)
if likelihood_gradient is None:
likelihood_gradient = self._likelihood_gradient(sample)
else:
likelihood_gradient += self._likelihood_gradient(sample)
return -likelihood_gradient / float(len(self.sample_list)) + \
self.T(self.position)
def _likelihood_gradient(self, m):
residual = self.d - \
self.Instrument(self.nonlinearity(
self.FFT.adjoint_times(self.power * m)))
LinearizedResponse = LinearizedPowerResponse(self.Instrument, self.nonlinearity,
self.FFT, self.Projection, self.position, m)
gradient = LinearizedResponse.adjoint_times(
self.N.inverse_times(residual))
return gradient
self.FFT.adjoint_times(self.power*m)))
LinearizedResponse = LinearizedPowerResponse(
self.Instrument, self.nonlinearity, self.FFT, self.Projection,
self.position, m)
return LinearizedResponse.adjoint_times(self.N.inverse_times(residual))
@property
@memo
def curvature(self):
curvature = NonlinearPowerCurvature(self.position, self.FFT, self.Instrument, self.nonlinearity,
self.Projection, self.N, self.T, self.sample_list)
curvature = NonlinearPowerCurvature(
self.position, self.FFT, self.Instrument, self.nonlinearity,
self.Projection, self.N, self.T, self.sample_list)
return InversionEnabler(curvature, inverter=self.inverter)
from .wiener_filter_curvature import WienerFilterCurvature
from .. import Field, exp
from ..utilities import memo
from ..sugar import generate_posterior_sample
from ..minimization.energy import Energy
from ..operators.inversion_enabler import InversionEnabler
from .response_operators import LinearizedSignalResponse
......@@ -21,10 +19,7 @@ class NonlinearWienerFilterEnergy(Energy):
self.position)
position_map = FFT.adjoint_times(self.power * self.position)
# position_map = (Field(FFT.domain,val=position_map.val.real+0j))
self.residual = d - Instrument(nonlinearity(position_map))
# Field(d.domain,
# val=self.residual.val.get_full_data().view(np.complex128).conjugate().view(np.float64))
self.N = N
self.S = S
self.inverter = inverter
......@@ -33,14 +28,14 @@ class NonlinearWienerFilterEnergy(Energy):
def at(self, position):
return self.__class__(position, self.d, self.Instrument,
self.nonlinearity, self.FFT, self.power, self.N, self.S, inverter=self.inverter)
self.nonlinearity, self.FFT, self.power, self.N,
self.S, inverter=self.inverter)
@property
@memo
def value(self):
energy = 0.5 * self.position.vdot(self._t1)
energy += 0.5 * self.residual.vdot(self._t2)
return energy.real
energy = self.position.vdot(self._t1) + self.residual.vdot(self._t2)
return 0.5 * energy.real
@property
@memo
......
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