Skip to content

GitLab

Commit 177c1976 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'energy_tests' into 'NIFTy_4' 

Energy tests

See merge request ift/NIFTy!215
parents e7d4b853 8a09eadd
Pipeline #24225 passed with stage
in 5 minutes and 20 seconds
Hide whitespace changes
Inline Side-by-side
Showing with 721 additions and 161 deletions
nifty4/__init__.py
View file @ 177c1976
......@@ -55,11 +55,11 @@ from .sugar import *
from .plotting.plot import plot
from . import library
__all__= ["DomainObject", "FieldArray", "Space", "RGSpace", "LMSpace",
"HPSpace", "GLSpace", "DOFSpace", "PowerSpace", "DomainTuple",
"LinearOperator", "EndomorphicOperator", "ScalingOperator",
"DiagonalOperator", "FFTOperator", "FFTSmoothingOperator",
"DirectSmoothingOperator", "ResponseOperator", "LaplaceOperator",
"PowerProjectionOperator", "InversionEnabler",
"Field", "sqrt", "exp", "log",
"Prober", "DiagonalProberMixin", "TraceProberMixin"]
__all__ = ["DomainObject", "FieldArray", "Space", "RGSpace", "LMSpace",
"HPSpace", "GLSpace", "DOFSpace", "PowerSpace", "DomainTuple",
"LinearOperator", "EndomorphicOperator", "ScalingOperator",
"DiagonalOperator", "FFTOperator", "FFTSmoothingOperator",
"DirectSmoothingOperator", "ResponseOperator", "LaplaceOperator",
"PowerProjectionOperator", "InversionEnabler",
"Field", "sqrt", "exp", "log",
"Prober", "DiagonalProberMixin", "TraceProberMixin"]
nifty4/library/critical_power_energy.py
View file @ 177c1976
......@@ -35,9 +35,10 @@ class CriticalPowerEnergy(Energy):
Parameters
----------
position : Field,
The current position of this energy.
The current position of this energy. (Logarithm of power spectrum)
m : Field,
The map whose power spectrum has to be inferred
The map whose power spectrum is inferred. Needs to live in harmonic
signal space.
D : EndomorphicOperator,
The curvature of the Gaussian encoding the posterior covariance.
If not specified, the map is assumed to be no reconstruction.
......@@ -101,8 +102,8 @@ class CriticalPowerEnergy(Energy):
self._theta = exp(-self.position) * (self.q + self._w*0.5)
Tt = self.T(self.position)
energy = self._theta.integrate()
energy += self.position.integrate()*(self.alpha-0.5)
energy = self._theta.sum()
energy += self.position.sum() * (self.alpha-0.5)
energy += 0.5*self.position.vdot(Tt)
self._value = energy.real
......
nifty4/library/log_normal_wiener_filter_curvature.py
View file @ 177c1976
......@@ -19,9 +19,7 @@
from ..operators.inversion_enabler import InversionEnabler
def LogNormalWienerFilterCurvature(R, N, S, fft, expp_sspace, inverter):
part1 = S.inverse
part3 = (fft.adjoint * expp_sspace * fft *
R. adjoint * N.inverse * R *
fft.adjoint * expp_sspace * fft)
return InversionEnabler(part1 + part3, inverter)
def LogNormalWienerFilterCurvature(R, N, S, ht, expp_sspace, inverter):
LinResp = R * ht.adjoint * expp_sspace * ht
t_op = LinResp.adjoint * N.inverse * LinResp
return InversionEnabler(S.inverse + t_op, inverter)
nifty4/library/log_normal_wiener_filter_energy.py
View file @ 177c1976
......@@ -19,7 +19,7 @@
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 ..sugar import create_composed_ht_operator
from ..field import exp
......@@ -43,7 +43,7 @@ class LogNormalWienerFilterEnergy(Energy):
The prior signal covariance in harmonic space.
"""
def __init__(self, position, d, R, N, S, inverter, fft=None):
def __init__(self, position, d, R, N, S, inverter, ht=None):
super(LogNormalWienerFilterEnergy, self).__init__(position=position)
self.d = d
self.R = R
......@@ -51,26 +51,25 @@ class LogNormalWienerFilterEnergy(Energy):
self.S = S
self._inverter = inverter
if fft is None:
self._fft = create_composed_fft_operator(self.S.domain,
all_to='position')
if ht is None:
self._ht = create_composed_ht_operator(self.S.domain)
else:
self._fft = fft
self._ht = ht
self._expp_sspace = exp(self._fft(self.position))
self._expp_sspace = exp(self._ht(self.position))
Sp = self.S.inverse_times(self.position)
expp = self._fft.adjoint_times(self._expp_sspace)
expp = self._ht.adjoint_times(self._expp_sspace)
Rexppd = self.R(expp) - self.d
NRexppd = self.N.inverse_times(Rexppd)
self._value = 0.5*(self.position.vdot(Sp) + Rexppd.vdot(NRexppd))
exppRNRexppd = self._fft.adjoint_times(
self._expp_sspace * self._fft(self.R.adjoint_times(NRexppd)))
exppRNRexppd = self._ht.adjoint_times(
self._expp_sspace * self._ht(self.R.adjoint_times(NRexppd)))
self._gradient = Sp + exppRNRexppd
def at(self, position):
return self.__class__(position, self.d, self.R, self.N, self.S,
self._inverter, self._fft)
self._inverter, self._ht)
@property
def value(self):
......@@ -84,5 +83,5 @@ class LogNormalWienerFilterEnergy(Energy):
@memo
def curvature(self):
return LogNormalWienerFilterCurvature(
R=self.R, N=self.N, S=self.S, fft=self._fft,
R=self.R, N=self.N, S=self.S, ht=self._ht,
expp_sspace=self._expp_sspace, inverter=self._inverter)
nifty4/library/noise_energy.py
View file @ 177c1976
......@@ -20,19 +20,21 @@ from .. import Field, exp
from ..operators.diagonal_operator import DiagonalOperator
from ..minimization.energy import Energy
# TODO Take only residual_sample_list as argument
class NoiseEnergy(Energy):
def __init__(self, position, d, m, D, t, FFT, Instrument, nonlinearity,
alpha, q, Projection, samples=3, sample_list=None,
inverter=None):
def __init__(self, position, d, xi, D, t, ht, Instrument,
nonlinearity, alpha, q, Projection, samples=3,
xi_sample_list=None, inverter=None):
super(NoiseEnergy, self).__init__(position=position)
self.m = m
self.xi = xi
self.D = D
self.d = d
self.N = DiagonalOperator(diagonal=exp(self.position))
self.t = t
self.samples = samples
self.FFT = FFT
self.ht = ht
self.Instrument = Instrument
self.nonlinearity = nonlinearity
......@@ -40,48 +42,47 @@ class NoiseEnergy(Energy):
self.q = q
self.Projection = Projection
self.power = self.Projection.adjoint_times(exp(0.5 * self.t))
self.one = Field(self.position.domain, val=1.)
if sample_list is None:
if xi_sample_list is None:
if samples is None or samples == 0:
sample_list = [m]
xi_sample_list = [xi]
else:
sample_list = [D.generate_posterior_sample() + m
for _ in range(samples)]
self.sample_list = sample_list
xi_sample_list = [D.generate_posterior_sample() + xi
for _ in range(samples)]
self.xi_sample_list = xi_sample_list
self.inverter = inverter
self._value, self._gradient = self._value_and_grad()
def at(self, position):
return self.__class__(
position, self.d, self.m, self.D, self.t, self.FFT,
self.Instrument, self.nonlinearity, self.alpha, self.q,
self.Projection, sample_list=self.sample_list,
samples=self.samples, inverter=self.inverter)
A = Projection.adjoint_times(exp(.5*self.t))
self._gradient = None
for sample in self.xi_sample_list:
map_s = self.ht(A * sample)
def _value_and_grad(self):
likelihood_gradient = None
for sample in self.sample_list:
residual = self.d - \
self.Instrument(self.nonlinearity(
self.FFT.adjoint_times(self.power*sample)))
lh = 0.5 * residual.vdot(self.N.inverse_times(residual))
grad = -0.5 * self.N.inverse_times(residual.conjugate() * residual)
if likelihood_gradient is None:
likelihood = lh
likelihood_gradient = grad.copy()
self.Instrument(self.nonlinearity(map_s))
lh = .5 * residual.vdot(self.N.inverse_times(residual))
grad = -.5 * self.N.inverse_times(residual.conjugate()*residual)
if self._gradient is None:
self._value = lh
self._gradient = grad.copy()
else:
likelihood += lh
likelihood_gradient += grad
self._value += lh
self._gradient += grad
likelihood = ((likelihood / float(len(self.sample_list))) +
0.5 * self.position.integrate() +
(self.alpha - 1.).vdot(self.position) +
self.q.vdot(exp(-self.position)))
likelihood_gradient = (
likelihood_gradient / float(len(self.sample_list)) +
(self.alpha-0.5) -
self.q * (exp(-self.position)))
return likelihood, likelihood_gradient
self._value *= 1. / len(self.xi_sample_list)
self._value += .5 * self.position.sum()
self._value += (self.alpha - 1.).vdot(self.position) + \
self.q.vdot(exp(-self.position))
self._gradient *= 1. / len(self.xi_sample_list)
self._gradient += (self.alpha-0.5) - self.q*(exp(-self.position))
def at(self, position):
return self.__class__(
position, self.d, self.xi, self.D, self.t, self.ht,
self.Instrument, self.nonlinearity, self.alpha, self.q,
self.Projection, xi_sample_list=self.xi_sample_list,
samples=self.samples, inverter=self.inverter)
@property
def value(self):
......
nifty4/library/nonlinear_power_curvature.py
View file @ 177c1976
......@@ -20,13 +20,13 @@ from ..operators.inversion_enabler import InversionEnabler
from .response_operators import LinearizedPowerResponse
def NonlinearPowerCurvature(position, FFT, Instrument, nonlinearity,
Projection, N, T, sample_list, inverter):
def NonlinearPowerCurvature(tau, ht, Instrument, nonlinearity, Projection, N,
T, xi_sample_list, inverter):
result = None
for sample in sample_list:
for xi_sample in xi_sample_list:
LinearizedResponse = LinearizedPowerResponse(
Instrument, nonlinearity, FFT, Projection, position, sample)
Instrument, nonlinearity, ht, Projection, tau, xi_sample)
op = LinearizedResponse.adjoint*N.inverse*LinearizedResponse
result = op if result is None else result + op
result = result*(1./len(sample_list)) + T
result = result*(1./len(xi_sample_list)) + T
return InversionEnabler(result, inverter)
nifty4/library/nonlinear_power_energy.py
View file @ 177c1976
......@@ -17,11 +17,11 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from .. import exp
from ..utilities import memo
from ..minimization.energy import Energy
from ..operators.smoothness_operator import SmoothnessOperator
from ..utilities import memo
from .nonlinear_power_curvature import NonlinearPowerCurvature
from .response_operators import LinearizedPowerResponse
from ..minimization.energy import Energy
class NonlinearPowerEnergy(Energy):
......@@ -51,64 +51,65 @@ class NonlinearPowerEnergy(Energy):
default : 3
"""
def __init__(self, position, d, N, m, D, FFT, Instrument, nonlinearity,
Projection, sigma=0., samples=3, sample_list=None,
def __init__(self, position, d, N, xi, D, ht, Instrument, nonlinearity,
Projection, sigma=0., samples=3, xi_sample_list=None,
inverter=None):
super(NonlinearPowerEnergy, self).__init__(position)
self.m = m
self.xi = xi
self.D = D
self.d = d
self.N = N
self.T = SmoothnessOperator(domain=self.position.domain[0],
strength=sigma, logarithmic=True)
self.FFT = FFT
self.ht = ht
self.Instrument = Instrument
self.nonlinearity = nonlinearity
self.Projection = Projection
self._sigma = sigma
self.power = self.Projection.adjoint_times(exp(0.5*self.position))
if sample_list is None:
self.sigma = sigma
if xi_sample_list is None:
if samples is None or samples == 0:
sample_list = [m]
xi_sample_list = [xi]
else:
sample_list = [D.generate_posterior_sample() + m
for _ in range(samples)]
self.sample_list = sample_list
xi_sample_list = [D.generate_posterior_sample() + xi
for _ in range(samples)]
self.xi_sample_list = xi_sample_list
self.inverter = inverter
self._value, self._gradient = self._value_and_grad()
def at(self, position):
return self.__class__(position, self.d, self.N, self.m, self.D,
self.FFT, self.Instrument, self.nonlinearity,
self.Projection, sigma=self._sigma,
samples=len(self.sample_list),
sample_list=self.sample_list,
inverter=self.inverter)
A = Projection.adjoint_times(exp(.5 * position))
map_s = self.ht(A * xi)
Tpos = self.T(position)
self._gradient = None
for xi_sample in self.xi_sample_list:
map_s = self.ht(A * xi_sample)
LinR = LinearizedPowerResponse(
self.Instrument, self.nonlinearity, self.ht, self.Projection,
self.position, xi_sample)
def _value_and_grad(self):
likelihood_gradient = None
for sample in self.sample_list:
residual = self.d - \
self.Instrument(self.nonlinearity(
self.FFT.adjoint_times(self.power*sample)))
self.Instrument(self.nonlinearity(map_s))
lh = 0.5 * residual.vdot(self.N.inverse_times(residual))
LinR = LinearizedPowerResponse(
self.Instrument, self.nonlinearity, self.FFT, self.Projection,
self.position, sample)
grad = LinR.adjoint_times(self.N.inverse_times(residual))
if likelihood_gradient is None:
likelihood = lh
likelihood_gradient = grad.copy()
if self._gradient is None:
self._value = lh
self._gradient = grad.copy()
else:
likelihood += lh
likelihood_gradient += grad
Tpos = self.T(self.position)
likelihood *= 1./len(self.sample_list)
likelihood += 0.5*self.position.vdot(Tpos)
likelihood_gradient *= -1./len(self.sample_list)
likelihood_gradient += Tpos
return likelihood, likelihood_gradient
self._value += lh
self._gradient += grad
self._value *= 1. / len(self.xi_sample_list)
self._value += 0.5 * self.position.vdot(Tpos)
self._gradient *= -1. / len(self.xi_sample_list)
self._gradient += Tpos
def at(self, position):
return self.__class__(position, self.d, self.N, self.xi, self.D,
self.ht, self.Instrument, self.nonlinearity,
self.Projection, sigma=self.sigma,
samples=len(self.xi_sample_list),
xi_sample_list=self.xi_sample_list,
inverter=self.inverter)
@property
def value(self):
......@@ -122,6 +123,6 @@ class NonlinearPowerEnergy(Energy):
@memo
def curvature(self):
return NonlinearPowerCurvature(
self.position, self.FFT, self.Instrument, self.nonlinearity,
self.Projection, self.N, self.T, self.sample_list,
inverter=self.inverter)
self.position, self.ht, self.Instrument, self.nonlinearity,
self.Projection, self.N, self.T, self.xi_sample_list,
self.inverter)
nifty4/library/nonlinear_wiener_filter_energy.py
View file @ 177c1976
......@@ -23,20 +23,19 @@ from .response_operators import LinearizedSignalResponse
class NonlinearWienerFilterEnergy(Energy):
def __init__(self, position, d, Instrument, nonlinearity, FFT, power, N, S,
def __init__(self, position, d, Instrument, nonlinearity, ht, power, N, S,
inverter=None):
super(NonlinearWienerFilterEnergy, self).__init__(position=position)
self.d = d
self.Instrument = Instrument
self.nonlinearity = nonlinearity
self.FFT = FFT
self.ht = ht
self.power = power
self.LinearizedResponse = \
LinearizedSignalResponse(Instrument, nonlinearity, FFT, power,
self.position)
m = self.ht(self.power * self.position)
self.LinearizedResponse = LinearizedSignalResponse(
Instrument, nonlinearity, ht, power, m)
position_map = FFT.adjoint_times(self.power * self.position)
residual = d - Instrument(nonlinearity(position_map))
residual = d - Instrument(nonlinearity(m))
self.N = N
self.S = S
self.inverter = inverter
......@@ -48,8 +47,8 @@ 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.ht, self.power, self.N,
self.S, self.inverter)
@property
def value(self):
......
nifty4/library/response_operators.py
View file @ 177c1976
......@@ -19,14 +19,12 @@
from ..field import exp
def LinearizedSignalResponse(Instrument, nonlinearity, FFT, power, m):
position = FFT.adjoint_times(power*m)
return (Instrument * nonlinearity.derivative(position) *
FFT.adjoint * power)
def LinearizedSignalResponse(Instrument, nonlinearity, ht, power, m):
return Instrument * nonlinearity.derivative(m) * ht * power
def LinearizedPowerResponse(Instrument, nonlinearity, FFT, Projection, t, m):
power = exp(0.5*t)
position = FFT.adjoint_times(Projection.adjoint_times(power) * m)
def LinearizedPowerResponse(Instrument, nonlinearity, ht, Projection, tau, xi):
power = exp(0.5*tau)
position = ht(Projection.adjoint_times(power)*xi)
linearization = nonlinearity.derivative(position)
return (0.5 * Instrument * linearization * FFT.adjoint * m *
Projection.adjoint * power)
return 0.5*Instrument*linearization*ht*xi*Projection.adjoint*power
nifty4/library/wiener_filter_curvature.py
View file @ 177c1976
......@@ -92,7 +92,8 @@ class WienerFilterCurvature(EndomorphicOperator):
def generate_posterior_sample2(self):
power = self.S.diagonal()
mock_signal = Field.from_random(random_type="normal",
domain=self.S.domain, dtype=power.dtype)
domain=self.S.domain,
dtype=power.dtype)
mock_signal *= sqrt(power)
noise = self.N.diagonal()
......
nifty4/library/wiener_filter_energy.py
View file @ 177c1976
......@@ -29,11 +29,12 @@ class WienerFilterEnergy(Energy):
Parameters
----------
position: Field,
The current position.
The current map in harmonic space.
d: Field,
the data
R: LinearOperator,
The response operator, description of the measurement process.
The response operator, description of the measurement process. It needs
to map from harmonic signal space to data space.
N: EndomorphicOperator,
The noise covariance in data space.
S: EndomorphicOperator,
......
nifty4/minimization/gradient_norm_controller.py
View file @ 177c1976
......@@ -57,8 +57,8 @@ class GradientNormController(IterationController):
if self._name is not None:
msg = self._name+":"
msg += " Iteration #" + str(self._itcount)
msg += " energy=" + str(energy.value)
msg += " gradnorm=" + str(energy.gradient_norm)
msg += " energy={:.6E}".format(energy.value)
msg += " gradnorm={:.2E}".format(energy.gradient_norm)
msg += " clvl=" + str(self._ccount)
dobj.mprint(msg)
# self.logger.info(msg)
......
nifty4/operators/scaling_operator.py
View file @ 177c1976
......@@ -69,7 +69,7 @@ class ScalingOperator(EndomorphicOperator):
@property
def inverse(self):
if self._factor!= 0.:
if self._factor != 0.:
return ScalingOperator(1./self._factor, self._domain)
from .inverse_operator import InverseOperator
return InverseOperator(self)
......@@ -84,7 +84,7 @@ class ScalingOperator(EndomorphicOperator):
@property
def capability(self):
if self._factor==0.:
if self._factor == 0.:
return self.TIMES | self.ADJOINT_TIMES
return (self.TIMES | self.ADJOINT_TIMES |
self.INVERSE_TIMES | self.ADJOINT_INVERSE_TIMES)
nifty4/plotting/plot.py
View file @ 177c1976
......@@ -189,13 +189,13 @@ def plot(f, **kwargs):
if fld.domain != dom:
raise ValueError("domain mismatch")
if not (isinstance(dom[0], PowerSpace) or
(isinstance(dom[0], RGSpace) and len(dom[0].shape)==1)):
(isinstance(dom[0], RGSpace) and len(dom[0].shape) == 1)):
raise ValueError("PowerSpace or 1D RGSpace required")
label = _get_kw("label", None, **kwargs)
if label is None:
label = [None] * len(f)
if not isinstance (label, list):
if not isinstance(label, list):
label = [label]
dom = dom[0]
......@@ -216,7 +216,7 @@ def plot(f, **kwargs):
xcoord = np.arange(npoints, dtype=np.float64)*dist
for i, fld in enumerate(f):
ycoord = dobj.to_global_data(fld.val)
plt.plot(xcoord, ycoord,label=label[i])
plt.plot(xcoord, ycoord, label=label[i])
_limit_xy(**kwargs)
if label != ([None]*len(f)):
plt.legend()
......
nifty4/probing/diagonal_prober_mixin.py
View file @ 177c1976
......@@ -18,7 +18,7 @@
from __future__ import division
from builtins import object
from ..sugar import create_composed_fft_operator
from ..sugar import create_composed_ht_operator
class DiagonalProberMixin(object):
......@@ -36,8 +36,8 @@ class DiagonalProberMixin(object):
def finish_probe(self, probe, pre_result):
if self.__evaluate_probe_in_signal_space:
fft = create_composed_fft_operator(self._domain, all_to='position')
result = fft(probe[1]).conjugate()*fft(pre_result)
ht = create_composed_ht_operator(self._domain)
result = ht(probe[1]).conjugate()*ht(pre_result)
else:
result = probe[1].conjugate()*pre_result
self.__sum_of_probings += result
......
nifty4/probing/trace_prober_mixin.py
View file @ 177c1976
......@@ -18,7 +18,7 @@
from __future__ import division
from builtins import object
from ..sugar import create_composed_fft_operator
from ..sugar import create_composed_ht_operator
class TraceProberMixin(object):
......@@ -36,8 +36,8 @@ class TraceProberMixin(object):
def finish_probe(self, probe, pre_result):
if self.__evaluate_probe_in_signal_space:
fft = create_composed_fft_operator(self._domain, all_to='position')
result = fft(probe[1]).vdot(fft(pre_result))
ht = create_composed_ht_operator(self._domain)
result = ht(probe[1]).vdot(ht(pre_result))
else:
result = probe[1].vdot(pre_result)
......
nifty4/sugar.py
View file @ 177c1976
......@@ -22,7 +22,7 @@ from .spaces.power_space import PowerSpace
from .field import Field, sqrt
from .operators.diagonal_operator import DiagonalOperator
from .operators.power_projection_operator import PowerProjectionOperator
from .operators.fft_operator import FFTOperator
from .operators.harmonic_transform_operator import HarmonicTransformOperator
from .domain_tuple import DomainTuple
from . import dobj, utilities
......@@ -32,7 +32,7 @@ __all__ = ['PS_field',
'power_synthesize_nonrandom',
'create_power_field',
'create_power_operator',