Commit cfeedd21 authored by Philipp Arras's avatar Philipp Arras

Some pep8 changes

parent 4df6cfa6
......@@ -22,4 +22,4 @@ from ..operators.diagonal_operator import DiagonalOperator
def CriticalPowerCurvature(theta, T, inverter):
theta = DiagonalOperator(theta)
return InversionEnabler(T+theta, inverter, theta.inverse_times)
return InversionEnabler(T + theta, inverter, theta.inverse_times)
......@@ -94,21 +94,21 @@ class CriticalPowerEnergy(Energy):
sample = self.D.generate_posterior_sample() + self.m
w += P(abs(sample)**2)
w *= 1./self.samples
w *= 1. / self.samples
else:
w = P(abs(self.m)**2)
self._w = w
self._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 = self._theta.sum()
energy += self.position.sum() * (self.alpha-0.5)
energy += 0.5*self.position.vdot(Tt)
energy += self.position.sum() * (self.alpha - 0.5)
energy += 0.5 * self.position.vdot(Tt)
self._value = energy.real
gradient = -self._theta
gradient += self.alpha-0.5
gradient += self.alpha - 0.5
gradient += Tt
self._gradient = gradient
......
......@@ -63,7 +63,7 @@ class LogNormalWienerFilterEnergy(Energy):
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))
self._value = 0.5 * (self.position.vdot(Sp) + Rexppd.vdot(NRexppd))
exppRNRexppd = self._ht.adjoint_times(
self._expp_sspace * self._ht(self.R.adjoint_times(NRexppd)))
self._gradient = Sp + exppRNRexppd
......
......@@ -50,7 +50,7 @@ class NoiseEnergy(Energy):
xi_sample_list = [xi]
else:
xi_sample_list = [D.generate_posterior_sample() + xi
for _ in range(samples)]
for _ in range(samples)]
self.xi_sample_list = xi_sample_list
self.inverter = inverter
......
......@@ -20,13 +20,23 @@ from ..operators.inversion_enabler import InversionEnabler
from .response_operators import LinearizedPowerResponse
def NonlinearPowerCurvature(tau, ht, Instrument, nonlinearity,
Projection, N, T, xi_sample_list, inverter, munit=1., sunit=1.):
def NonlinearPowerCurvature(
tau,
ht,
Instrument,
nonlinearity,
Projection,
N,
T,
xi_sample_list,
inverter,
munit=1.,
sunit=1.):
result = None
for xi_sample in xi_sample_list:
LinearizedResponse = LinearizedPowerResponse(
Instrument, nonlinearity, ht, Projection, tau, xi_sample, munit, sunit)
op = LinearizedResponse.adjoint*N.inverse*LinearizedResponse
op = LinearizedResponse.adjoint * N.inverse * LinearizedResponse
result = op if result is None else result + op
result = result * (1. / len(xi_sample_list)) + T
return InversionEnabler(result, inverter)
......@@ -73,7 +73,7 @@ class NonlinearPowerEnergy(Energy):
xi_sample_list = [xi]
else:
xi_sample_list = [D.generate_posterior_sample() + xi
for _ in range(samples)]
for _ in range(samples)]
self.xi_sample_list = xi_sample_list
self.inverter = inverter
......
......@@ -33,8 +33,8 @@ class NonlinearWienerFilterEnergy(Energy):
self.ht = ht
self.power = power
position_map = self.ht(self.power * self.position)
self.LinearizedResponse = LinearizedSignalResponse(Instrument, nonlinearity, ht, power,
position_map, sunit)
self.LinearizedResponse = LinearizedSignalResponse(
Instrument, nonlinearity, ht, power, position_map, sunit)
position_map = ht(self.power * self.position)
residual = d - Instrument(sunit * nonlinearity(position_map))
......
......@@ -23,8 +23,16 @@ def LinearizedSignalResponse(Instrument, nonlinearity, ht, power, m, sunit):
return sunit * (Instrument * nonlinearity.derivative(m) * ht * power)
def LinearizedPowerResponse(Instrument, nonlinearity, ht, Projection, tau, xi, munit, sunit):
power = exp(0.5*tau) * munit
def LinearizedPowerResponse(
Instrument,
nonlinearity,
ht,
Projection,
tau,
xi,
munit,
sunit):
power = exp(0.5 * tau) * munit
position = ht(Projection.adjoint_times(power) * xi)
linearization = nonlinearity.derivative(position)
return sunit * (0.5 * Instrument * linearization * ht * xi *
......
......@@ -52,7 +52,8 @@ class WienerFilterEnergy(Energy):
_j = self.R.adjoint_times(self.N.inverse_times(d))
self._j = _j
Dx = self._curvature(self.position)
self._value = 0.5*self.position.vdot(Dx) - self._j.vdot(self.position)
self._value = 0.5 * \
self.position.vdot(Dx) - self._j.vdot(self.position)
self._gradient = Dx - self._j
def at(self, position):
......
......@@ -75,7 +75,7 @@ class Energy_Tests(unittest.TestCase):
@expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[4, 78, 23]))
[4, 78, 23]))
def testLognormalMap(self, space, seed):
np.random.seed(seed)
dim = len(space.shape)
......@@ -121,7 +121,7 @@ class Energy_Tests(unittest.TestCase):
@expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[ift.library.Exponential, ift.library.Linear],
[4, 78, 23]))
[4, 78, 23]))
def testNonlinearMap(self, space, nonlinearity, seed):
np.random.seed(seed)
f = nonlinearity()
......@@ -194,7 +194,8 @@ class Curvature_Tests(unittest.TestCase):
inverter = ift.ConjugateGradient(IC)
S = ift.create_power_operator(hspace, power_spectrum=lambda k: 1.)
energy0 = ift.library.WienerFilterEnergy(position=s0, d=d, R=R, N=N, S=S, inverter=inverter)
energy0 = ift.library.WienerFilterEnergy(
position=s0, d=d, R=R, N=N, S=S, inverter=inverter)
gradient0 = energy0.gradient
gradient1 = energy0.at(s1).gradient
......@@ -205,7 +206,7 @@ class Curvature_Tests(unittest.TestCase):
@expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[4, 78, 23]))
[4, 78, 23]))
def testLognormalMapCurvature(self, space, seed):
np.random.seed(seed)
dim = len(space.shape)
......@@ -252,7 +253,7 @@ class Curvature_Tests(unittest.TestCase):
@expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[ift.library.Exponential, ift.library.Linear],
[4, 78, 23]))
[4, 78, 23]))
def testNonlinearMapCurvature(self, space, nonlinearity, seed):
np.random.seed(seed)
f = nonlinearity()
......@@ -285,7 +286,15 @@ class Curvature_Tests(unittest.TestCase):
tol_abs_gradnorm=1e-7)
inverter = ift.ConjugateGradient(IC)
energy0 = ift.library.NonlinearWienerFilterEnergy(
position=xi0, d=d, Instrument=R, nonlinearity=f, ht=ht, power=A, N=N, S=S, inverter=inverter)
position=xi0,
d=d,
Instrument=R,
nonlinearity=f,
ht=ht,
power=A,
N=N,
S=S,
inverter=inverter)
gradient0 = energy0.gradient
gradient1 = energy0.at(xi1).gradient
......
......@@ -47,7 +47,10 @@ class Noise_Energy_Tests(unittest.TestCase):
tau = ift.PS_field(pspace, pspec)
A = P.adjoint_times(ift.sqrt(tau))
var = 1.
n = ift.Field.from_random(domain=space, random_type='normal', std=np.sqrt(var))
n = ift.Field.from_random(
domain=space,
random_type='normal',
std=np.sqrt(var))
var = ift.Field(n.domain, val=var)
N = ift.DiagonalOperator(var)
eta0 = ift.log(var)
......@@ -80,7 +83,9 @@ class Noise_Energy_Tests(unittest.TestCase):
S=S,
inverter=inverter).curvature
Nsamples = 10
xi_sample_list = [D.generate_posterior_sample() + xi for i in range(Nsamples)]
xi_sample_list = [
D.generate_posterior_sample() +
xi for i in range(Nsamples)]
energy0 = ift.library.NoiseEnergy(
position=eta0, d=d, xi=xi, D=D, t=tau, Instrument=R,
......
......@@ -132,7 +132,9 @@ class Energy_Tests(unittest.TestCase):
ht=ht,
inverter=inverter).curvature
Nsamples = 10
xi_sample_list = [D.generate_posterior_sample() + xi for _ in range(Nsamples)]
xi_sample_list = [
D.generate_posterior_sample() +
xi for _ in range(Nsamples)]
energy0 = ift.library.NonlinearPowerEnergy(
position=tau0,
......@@ -271,7 +273,9 @@ class Curvature_Tests(unittest.TestCase):
ht=ht,
inverter=inverter).curvature
Nsamples = 10
xi_sample_list = [D.generate_posterior_sample() + xi for _ in range(Nsamples)]
xi_sample_list = [
D.generate_posterior_sample() +
xi for _ in range(Nsamples)]
energy0 = ift.library.NonlinearPowerEnergy(
position=tau0,
......
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