Commit 41d3585e authored by Martin Reinecke's avatar Martin Reinecke
Browse files

PEP8 work, revert some overly strict PEP8

parent b753ad27
Pipeline #24220 passed with stage
in 4 minutes and 58 seconds
...@@ -55,11 +55,11 @@ from .sugar import * ...@@ -55,11 +55,11 @@ from .sugar import *
from .plotting.plot import plot from .plotting.plot import plot
from . import library from . import library
__all__= ["DomainObject", "FieldArray", "Space", "RGSpace", "LMSpace", __all__ = ["DomainObject", "FieldArray", "Space", "RGSpace", "LMSpace",
"HPSpace", "GLSpace", "DOFSpace", "PowerSpace", "DomainTuple", "HPSpace", "GLSpace", "DOFSpace", "PowerSpace", "DomainTuple",
"LinearOperator", "EndomorphicOperator", "ScalingOperator", "LinearOperator", "EndomorphicOperator", "ScalingOperator",
"DiagonalOperator", "FFTOperator", "FFTSmoothingOperator", "DiagonalOperator", "FFTOperator", "FFTSmoothingOperator",
"DirectSmoothingOperator", "ResponseOperator", "LaplaceOperator", "DirectSmoothingOperator", "ResponseOperator", "LaplaceOperator",
"PowerProjectionOperator", "InversionEnabler", "PowerProjectionOperator", "InversionEnabler",
"Field", "sqrt", "exp", "log", "Field", "sqrt", "exp", "log",
"Prober", "DiagonalProberMixin", "TraceProberMixin"] "Prober", "DiagonalProberMixin", "TraceProberMixin"]
...@@ -22,4 +22,4 @@ from ..operators.diagonal_operator import DiagonalOperator ...@@ -22,4 +22,4 @@ from ..operators.diagonal_operator import DiagonalOperator
def CriticalPowerCurvature(theta, T, inverter): def CriticalPowerCurvature(theta, T, inverter):
theta = DiagonalOperator(theta) 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): ...@@ -94,21 +94,21 @@ class CriticalPowerEnergy(Energy):
sample = self.D.generate_posterior_sample() + self.m sample = self.D.generate_posterior_sample() + self.m
w += P(abs(sample)**2) w += P(abs(sample)**2)
w *= 1. / self.samples w *= 1./self.samples
else: else:
w = P(abs(self.m)**2) w = P(abs(self.m)**2)
self._w = w 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) Tt = self.T(self.position)
energy = self._theta.sum() energy = self._theta.sum()
energy += self.position.sum() * (self.alpha - 0.5) energy += self.position.sum() * (self.alpha-0.5)
energy += 0.5 * self.position.vdot(Tt) energy += 0.5*self.position.vdot(Tt)
self._value = energy.real self._value = energy.real
gradient = -self._theta gradient = -self._theta
gradient += self.alpha - 0.5 gradient += self.alpha-0.5
gradient += Tt gradient += Tt
self._gradient = gradient self._gradient = gradient
......
...@@ -20,7 +20,6 @@ from ..operators.inversion_enabler import InversionEnabler ...@@ -20,7 +20,6 @@ from ..operators.inversion_enabler import InversionEnabler
def LogNormalWienerFilterCurvature(R, N, S, ht, expp_sspace, inverter): def LogNormalWienerFilterCurvature(R, N, S, ht, expp_sspace, inverter):
part1 = S.inverse
LinResp = R * ht.adjoint * expp_sspace * ht LinResp = R * ht.adjoint * expp_sspace * ht
part3 = (LinResp.adjoint * N.inverse * LinResp) t_op = LinResp.adjoint * N.inverse * LinResp
return InversionEnabler(part1 + part3, inverter) return InversionEnabler(S.inverse + t_op, inverter)
...@@ -19,7 +19,7 @@ ...@@ -19,7 +19,7 @@
from ..minimization.energy import Energy from ..minimization.energy import Energy
from ..utilities import memo from ..utilities import memo
from .log_normal_wiener_filter_curvature import LogNormalWienerFilterCurvature 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 from ..field import exp
...@@ -52,8 +52,7 @@ class LogNormalWienerFilterEnergy(Energy): ...@@ -52,8 +52,7 @@ class LogNormalWienerFilterEnergy(Energy):
self._inverter = inverter self._inverter = inverter
if ht is None: if ht is None:
self._ht = create_composed_fft_operator(self.S.domain, self._ht = create_composed_ht_operator(self.S.domain)
all_to='position')
else: else:
self._ht = ht self._ht = ht
...@@ -63,7 +62,7 @@ class LogNormalWienerFilterEnergy(Energy): ...@@ -63,7 +62,7 @@ class LogNormalWienerFilterEnergy(Energy):
expp = self._ht.adjoint_times(self._expp_sspace) expp = self._ht.adjoint_times(self._expp_sspace)
Rexppd = self.R(expp) - self.d Rexppd = self.R(expp) - self.d
NRexppd = self.N.inverse_times(Rexppd) 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( exppRNRexppd = self._ht.adjoint_times(
self._expp_sspace * self._ht(self.R.adjoint_times(NRexppd))) self._expp_sspace * self._ht(self.R.adjoint_times(NRexppd)))
self._gradient = Sp + exppRNRexppd self._gradient = Sp + exppRNRexppd
......
...@@ -51,7 +51,7 @@ class NoiseEnergy(Energy): ...@@ -51,7 +51,7 @@ class NoiseEnergy(Energy):
self.xi_sample_list = xi_sample_list self.xi_sample_list = xi_sample_list
self.inverter = inverter self.inverter = inverter
A = Projection.adjoint_times(exp(.5 * self.t)) A = Projection.adjoint_times(exp(.5*self.t))
self._gradient = None self._gradient = None
for sample in self.xi_sample_list: for sample in self.xi_sample_list:
...@@ -60,7 +60,7 @@ class NoiseEnergy(Energy): ...@@ -60,7 +60,7 @@ class NoiseEnergy(Energy):
residual = self.d - \ residual = self.d - \
self.Instrument(self.nonlinearity(map_s)) self.Instrument(self.nonlinearity(map_s))
lh = .5 * residual.vdot(self.N.inverse_times(residual)) lh = .5 * residual.vdot(self.N.inverse_times(residual))
grad = -.5 * self.N.inverse_times(residual.conjugate() * residual) grad = -.5 * self.N.inverse_times(residual.conjugate()*residual)
if self._gradient is None: if self._gradient is None:
self._value = lh self._value = lh
...@@ -75,7 +75,7 @@ class NoiseEnergy(Energy): ...@@ -75,7 +75,7 @@ class NoiseEnergy(Energy):
self.q.vdot(exp(-self.position)) self.q.vdot(exp(-self.position))
self._gradient *= 1. / len(self.xi_sample_list) self._gradient *= 1. / len(self.xi_sample_list)
self._gradient += (self.alpha - 0.5) - self.q * (exp(-self.position)) self._gradient += (self.alpha-0.5) - self.q*(exp(-self.position))
def at(self, position): def at(self, position):
return self.__class__( return self.__class__(
......
...@@ -26,7 +26,7 @@ def NonlinearPowerCurvature(tau, ht, Instrument, nonlinearity, Projection, N, ...@@ -26,7 +26,7 @@ def NonlinearPowerCurvature(tau, ht, Instrument, nonlinearity, Projection, N,
for xi_sample in xi_sample_list: for xi_sample in xi_sample_list:
LinearizedResponse = LinearizedPowerResponse( LinearizedResponse = LinearizedPowerResponse(
Instrument, nonlinearity, ht, Projection, tau, xi_sample) Instrument, nonlinearity, ht, Projection, tau, xi_sample)
op = LinearizedResponse.adjoint * N.inverse * LinearizedResponse op = LinearizedResponse.adjoint*N.inverse*LinearizedResponse
result = op if result is None else result + op result = op if result is None else result + op
result = result * (1. / len(xi_sample_list)) + T result = result * (1./len(xi_sample_list)) + T
return InversionEnabler(result, inverter) return InversionEnabler(result, inverter)
...@@ -20,11 +20,11 @@ from ..field import exp ...@@ -20,11 +20,11 @@ from ..field import exp
def LinearizedSignalResponse(Instrument, nonlinearity, ht, power, m): def LinearizedSignalResponse(Instrument, nonlinearity, ht, power, m):
return (Instrument * nonlinearity.derivative(m) * ht * power) return Instrument * nonlinearity.derivative(m) * ht * power
def LinearizedPowerResponse(Instrument, nonlinearity, ht, Projection, tau, xi): def LinearizedPowerResponse(Instrument, nonlinearity, ht, Projection, tau, xi):
power = exp(0.5 * tau) power = exp(0.5*tau)
position = ht(Projection.adjoint_times(power) * xi) position = ht(Projection.adjoint_times(power)*xi)
linearization = nonlinearity.derivative(position) linearization = nonlinearity.derivative(position)
return (0.5 * Instrument * linearization * ht * xi * Projection.adjoint * power) return 0.5*Instrument*linearization*ht*xi*Projection.adjoint*power
...@@ -45,7 +45,7 @@ class WienerFilterCurvature(EndomorphicOperator): ...@@ -45,7 +45,7 @@ class WienerFilterCurvature(EndomorphicOperator):
self.R = R self.R = R
self.N = N self.N = N
self.S = S self.S = S
op = R.adjoint * N.inverse * R + S.inverse op = R.adjoint*N.inverse*R + S.inverse
self._op = InversionEnabler(op, inverter, S.times) self._op = InversionEnabler(op, inverter, S.times)
@property @property
...@@ -91,10 +91,9 @@ class WienerFilterCurvature(EndomorphicOperator): ...@@ -91,10 +91,9 @@ class WienerFilterCurvature(EndomorphicOperator):
def generate_posterior_sample2(self): def generate_posterior_sample2(self):
power = self.S.diagonal() power = self.S.diagonal()
mock_signal = Field.from_random( mock_signal = Field.from_random(random_type="normal",
random_type="normal", domain=self.S.domain,
domain=self.S.domain, dtype=power.dtype)
dtype=power.dtype)
mock_signal *= sqrt(power) mock_signal *= sqrt(power)
noise = self.N.diagonal() noise = self.N.diagonal()
......
...@@ -52,8 +52,7 @@ class WienerFilterEnergy(Energy): ...@@ -52,8 +52,7 @@ class WienerFilterEnergy(Energy):
_j = self.R.adjoint_times(self.N.inverse_times(d)) _j = self.R.adjoint_times(self.N.inverse_times(d))
self._j = _j self._j = _j
Dx = self._curvature(self.position) Dx = self._curvature(self.position)
self._value = 0.5 * \ self._value = 0.5*self.position.vdot(Dx) - self._j.vdot(self.position)
self.position.vdot(Dx) - self._j.vdot(self.position)
self._gradient = Dx - self._j self._gradient = Dx - self._j
def at(self, position): def at(self, position):
......
...@@ -69,7 +69,7 @@ class ScalingOperator(EndomorphicOperator): ...@@ -69,7 +69,7 @@ class ScalingOperator(EndomorphicOperator):
@property @property
def inverse(self): def inverse(self):
if self._factor!= 0.: if self._factor != 0.:
return ScalingOperator(1./self._factor, self._domain) return ScalingOperator(1./self._factor, self._domain)
from .inverse_operator import InverseOperator from .inverse_operator import InverseOperator
return InverseOperator(self) return InverseOperator(self)
...@@ -84,7 +84,7 @@ class ScalingOperator(EndomorphicOperator): ...@@ -84,7 +84,7 @@ class ScalingOperator(EndomorphicOperator):
@property @property
def capability(self): def capability(self):
if self._factor==0.: if self._factor == 0.:
return self.TIMES | self.ADJOINT_TIMES return self.TIMES | self.ADJOINT_TIMES
return (self.TIMES | self.ADJOINT_TIMES | return (self.TIMES | self.ADJOINT_TIMES |
self.INVERSE_TIMES | self.ADJOINT_INVERSE_TIMES) self.INVERSE_TIMES | self.ADJOINT_INVERSE_TIMES)
...@@ -189,13 +189,13 @@ def plot(f, **kwargs): ...@@ -189,13 +189,13 @@ def plot(f, **kwargs):
if fld.domain != dom: if fld.domain != dom:
raise ValueError("domain mismatch") raise ValueError("domain mismatch")
if not (isinstance(dom[0], PowerSpace) or 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") raise ValueError("PowerSpace or 1D RGSpace required")
label = _get_kw("label", None, **kwargs) label = _get_kw("label", None, **kwargs)
if label is None: if label is None:
label = [None] * len(f) label = [None] * len(f)
if not isinstance (label, list): if not isinstance(label, list):
label = [label] label = [label]
dom = dom[0] dom = dom[0]
...@@ -216,7 +216,7 @@ def plot(f, **kwargs): ...@@ -216,7 +216,7 @@ def plot(f, **kwargs):
xcoord = np.arange(npoints, dtype=np.float64)*dist xcoord = np.arange(npoints, dtype=np.float64)*dist
for i, fld in enumerate(f): for i, fld in enumerate(f):
ycoord = dobj.to_global_data(fld.val) ycoord = dobj.to_global_data(fld.val)
plt.plot(xcoord, ycoord,label=label[i]) plt.plot(xcoord, ycoord, label=label[i])
_limit_xy(**kwargs) _limit_xy(**kwargs)
if label != ([None]*len(f)): if label != ([None]*len(f)):
plt.legend() plt.legend()
......
...@@ -23,6 +23,7 @@ from .field import Field, sqrt ...@@ -23,6 +23,7 @@ from .field import Field, sqrt
from .operators.diagonal_operator import DiagonalOperator from .operators.diagonal_operator import DiagonalOperator
from .operators.power_projection_operator import PowerProjectionOperator from .operators.power_projection_operator import PowerProjectionOperator
from .operators.fft_operator import FFTOperator from .operators.fft_operator import FFTOperator
from .operators.harmonic_transform_operator import HarmonicTransformOperator
from .domain_tuple import DomainTuple from .domain_tuple import DomainTuple
from . import dobj, utilities from . import dobj, utilities
...@@ -32,6 +33,7 @@ __all__ = ['PS_field', ...@@ -32,6 +33,7 @@ __all__ = ['PS_field',
'power_synthesize_nonrandom', 'power_synthesize_nonrandom',
'create_power_field', 'create_power_field',
'create_power_operator', 'create_power_operator',
'create_composed_ht_operator',
'create_composed_fft_operator', 'create_composed_fft_operator',
'create_harmonic_smoothing_operator'] 'create_harmonic_smoothing_operator']
...@@ -241,6 +243,20 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None): ...@@ -241,6 +243,20 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None):
spaces=space) spaces=space)
def create_composed_ht_operator(domain, codomain=None):
if codomain is None:
codomain = [None]*len(domain)
res = None
for i, space in enumerate(domain):
if isinstance(space, Space) and space.harmonic:
tdom = domain if res is None else res.target
op = HarmonicTransformOperator(tdom, codomain[i], i)
res = op if res is None else op*res
if res is None:
raise ValueError("empty operator")
return res
def create_composed_fft_operator(domain, codomain=None, all_to='other'): def create_composed_fft_operator(domain, codomain=None, all_to='other'):
if codomain is None: if codomain is None:
codomain = [None]*len(domain) codomain = [None]*len(domain)
......
...@@ -147,7 +147,8 @@ class Energy_Tests(unittest.TestCase): ...@@ -147,7 +147,8 @@ class Energy_Tests(unittest.TestCase):
S = ift.create_power_operator(hspace, power_spectrum=_flat_PS) S = ift.create_power_operator(hspace, power_spectrum=_flat_PS)
energy0 = ift.library.NonlinearWienerFilterEnergy( energy0 = ift.library.NonlinearWienerFilterEnergy(
position=xi0, d=d, Instrument=R, nonlinearity=f, ht=ht, power=A, N=N, S=S) position=xi0, d=d, Instrument=R, nonlinearity=f, ht=ht, power=A,
N=N, S=S)
energy1 = energy0.at(xi1) energy1 = energy0.at(xi1)
a = (energy1.value - energy0.value) / eps a = (energy1.value - energy0.value) / eps
...@@ -204,7 +205,8 @@ class Curvature_Tests(unittest.TestCase): ...@@ -204,7 +205,8 @@ class Curvature_Tests(unittest.TestCase):
@expand(product([ift.RGSpace(64, distances=.789), @expand(product([ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)], ift.RGSpace([32, 32], distances=.789)],
# Only linear case due to approximation of Hessian in the case of nontrivial nonlinearities. # Only linear case due to approximation of Hessian in the
# case of nontrivial nonlinearities.
[ift.library.Linear], [ift.library.Linear],
[4, 78, 23])) [4, 78, 23]))
def testNonlinearMapCurvature(self, space, nonlinearity, seed): def testNonlinearMapCurvature(self, space, nonlinearity, seed):
......
...@@ -70,7 +70,8 @@ class Energy_Tests(unittest.TestCase): ...@@ -70,7 +70,8 @@ class Energy_Tests(unittest.TestCase):
inverter=inverter).curvature inverter=inverter).curvature
energy0 = ift.library.CriticalPowerEnergy( energy0 = ift.library.CriticalPowerEnergy(
position=tau0, m=s, inverter=inverter, D=D, samples=10, smoothness_prior=1.) position=tau0, m=s, inverter=inverter, D=D, samples=10,
smoothness_prior=1.)
energy1 = energy0.at(tau1) energy1 = energy0.at(tau1)
a = (energy1.value - energy0.value) / eps a = (energy1.value - energy0.value) / eps
......
...@@ -41,7 +41,8 @@ class Test_Minimizers(unittest.TestCase): ...@@ -41,7 +41,8 @@ class Test_Minimizers(unittest.TestCase):
covariance = ift.DiagonalOperator(covariance_diagonal) covariance = ift.DiagonalOperator(covariance_diagonal)
required_result = ift.Field.ones(space, dtype=np.float64) required_result = ift.Field.ones(space, dtype=np.float64)
IC = ift.GradientNormController(tol_abs_gradnorm=1e-5, iteration_limit=1000) IC = ift.GradientNormController(tol_abs_gradnorm=1e-5,
iteration_limit=1000)
try: try:
minimizer = minimizer_class(controller=IC) minimizer = minimizer_class(controller=IC)
energy = ift.QuadraticEnergy(A=covariance, b=required_result, energy = ift.QuadraticEnergy(A=covariance, b=required_result,
...@@ -57,4 +58,4 @@ class Test_Minimizers(unittest.TestCase): ...@@ -57,4 +58,4 @@ class Test_Minimizers(unittest.TestCase):
rtol=1e-3, atol=1e-3) rtol=1e-3, atol=1e-3)
#MR FIXME: add Rosenbrock test # MR FIXME: add Rosenbrock test
...@@ -39,6 +39,7 @@ def _check_adjointness(op, dtype=np.float64): ...@@ -39,6 +39,7 @@ def _check_adjointness(op, dtype=np.float64):
op.inverse_adjoint_times(f1).vdot(f2), op.inverse_adjoint_times(f1).vdot(f2),
rtol=1e-8) rtol=1e-8)
_h_RG_spaces = [ift.RGSpace(7, distances=0.2, harmonic=True), _h_RG_spaces = [ift.RGSpace(7, distances=0.2, harmonic=True),
ift.RGSpace((12, 46), distances=(.2, .3), harmonic=True)] ift.RGSpace((12, 46), distances=(.2, .3), harmonic=True)]
_h_spaces = _h_RG_spaces + [ift.LMSpace(17)] _h_spaces = _h_RG_spaces + [ift.LMSpace(17)]
......
...@@ -83,7 +83,6 @@ class FFTOperatorTests(unittest.TestCase): ...@@ -83,7 +83,6 @@ class FFTOperatorTests(unittest.TestCase):
assert_allclose(ift.dobj.to_global_data(inp.val), assert_allclose(ift.dobj.to_global_data(inp.val),
ift.dobj.to_global_data(out.val), rtol=tol, atol=tol) ift.dobj.to_global_data(out.val), rtol=tol, atol=tol)
@expand(product([0, 1, 2], @expand(product([0, 1, 2],
[np.float64, np.float32, np.complex64, np.complex128])) [np.float64, np.float32, np.complex64, np.complex128]))
def test_composed_fft(self, index, dtype): def test_composed_fft(self, index, dtype):
...@@ -98,9 +97,8 @@ class FFTOperatorTests(unittest.TestCase): ...@@ -98,9 +97,8 @@ class FFTOperatorTests(unittest.TestCase):
assert_allclose(ift.dobj.to_global_data(inp.val), assert_allclose(ift.dobj.to_global_data(inp.val),
ift.dobj.to_global_data(out.val), rtol=tol, atol=tol) ift.dobj.to_global_data(out.val), rtol=tol, atol=tol)
@expand(product([ift.RGSpace(128, distances=3.76, harmonic=True), @expand(product([ift.RGSpace(128, distances=3.76, harmonic=True),
ift.RGSpace((15,27), distances=(1.7,3.33), harmonic=True), ift.RGSpace((15, 27), distances=(.7, .33), harmonic=True),
ift.RGSpace(73, distances=0.5643)], ift.RGSpace(73, distances=0.5643)],
[np.float64, np.float32, np.complex64, np.complex128])) [np.float64, np.float32, np.complex64, np.complex128]))
def test_normalisation(self, space, tp): def test_normalisation(self, space, tp):
......
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