From 002fc6fe823d07299af2769e08904144c5c886ac Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Sun, 14 Apr 2019 23:54:33 +0200 Subject: [PATCH] Misc fixes --- examples/solvers/douglas_rachford_pd_heron.py | 13 +++++---- examples/solvers/douglas_rachford_pd_mri.py | 10 +++---- examples/space/simple_rn.py | 24 ---------------- odl/solvers/functional/default_functionals.py | 5 ++-- odl/solvers/nonsmooth/proximal_operators.py | 28 +++++++++---------- odl/test/operator/operator_test.py | 15 ++++++---- odl/util/testutils.py | 2 ++ 7 files changed, 40 insertions(+), 57 deletions(-) diff --git a/examples/solvers/douglas_rachford_pd_heron.py b/examples/solvers/douglas_rachford_pd_heron.py index 3a799b9eaa9..78ccc1f8a27 100644 --- a/examples/solvers/douglas_rachford_pd_heron.py +++ b/examples/solvers/douglas_rachford_pd_heron.py @@ -63,14 +63,15 @@ def print_objective(x): tau=tau, sigma=sigma, niter=20, lam=lam, callback=print_objective, l=l) -# plot the result +# Plot the result +fig, ax = plt.subplots() for minp, maxp in rectangles: xp = [minp[0], maxp[0], maxp[0], minp[0], minp[0]] yp = [minp[1], minp[1], maxp[1], maxp[1], minp[1]] - plt.plot(xp, yp) + ax.plot(xp, yp) -plt.scatter(x[0], x[1]) +ax.scatter(x[0], x[1]) -plt.xlim(-1, 4) -plt.ylim(-1, 4) -plt.show() +ax.set_xlim(-1, 4) +ax.set_ylim(-1, 4) +fig.show() diff --git a/examples/solvers/douglas_rachford_pd_mri.py b/examples/solvers/douglas_rachford_pd_mri.py index adf3f5448c0..460bde16e38 100644 --- a/examples/solvers/douglas_rachford_pd_mri.py +++ b/examples/solvers/douglas_rachford_pd_mri.py @@ -28,8 +28,8 @@ # Create noisy MRI data phantom = odl.phantom.shepp_logan(space, modified=True) noisy_data = mri_op(phantom) + odl.phantom.white_noise(mri_op.range) * 0.1 -phantom.show('Phantom') -noisy_data.show('Noisy MRI Data') +space.show(phantom, 'Phantom') +ft.range.show(noisy_data, 'Noisy MRI Data') # Gradient for TV regularization gradient = odl.Gradient(space) @@ -44,11 +44,11 @@ # Solve x = mri_op.domain.zero() -callback = (odl.solvers.CallbackShow(step=5, clim=[0, 1]) & +callback = (odl.solvers.CallbackShow(space, step=5, clim=[0, 1]) & odl.solvers.CallbackPrintIteration()) odl.solvers.douglas_rachford_pd(x, f, g, lin_ops, tau=2.0, sigma=[1.0, 0.1], niter=500, callback=callback) -x.show('Douglas-Rachford Result') -ft.inverse(noisy_data).show('Fourier Inversion Result', force_show=True) +space.show(x, 'TV-regularized Result (Douglas-Rachford)') +space.show(ft.inverse(noisy_data), 'Fourier Inversion Result', force_show=True) diff --git a/examples/space/simple_rn.py b/examples/space/simple_rn.py index a089c5e1a6d..77beeb347b9 100644 --- a/examples/space/simple_rn.py +++ b/examples/space/simple_rn.py @@ -61,7 +61,6 @@ def element(self, *args, **kwargs): # Do some tests to compare n = 10 ** 7 iterations = 10 -cuda_supported = 'cuda' in odl.space.entry_points.tensor_space_impl_names() # Perform some benchmarks with rn opt_space = odl.rn(n) @@ -71,10 +70,6 @@ def element(self, *args, **kwargs): ox, oy, oz = (opt_space.copy(a) for a in (x, y, z)) sx, sy, sz = (simple_space.copy(a) for a in (x, y, z)) -if cuda_supported: - cu_space = odl.rn(n, impl='cuda') - cx, cy, cz = (cu_space.element(a.copy()) for a in (x, y, z)) - print(" lincomb:") with timer("SimpleRn"): for _ in range(iterations): @@ -86,13 +81,6 @@ def element(self, *args, **kwargs): opt_space.lincomb(2.13, ox, 3.14, oy, out=oz) print("result: {}".format(oz[1:5])) -if cuda_supported: - with timer("odl cuda"): - for _ in range(iterations): - cu_space.lincomb(2.13, cx, 3.14, cy, out=cz) - print("result: {}".format(cz[1:5])) - - print("\n Norm:") with timer("SimpleRn"): for _ in range(iterations): @@ -104,12 +92,6 @@ def element(self, *args, **kwargs): result = opt_space.norm(oz) print("result: {}".format(result)) -if cuda_supported: - with timer("odl cuda"): - for _ in range(iterations): - result = cu_space.norm(cz) - print("result: {}".format(result)) - print("\n Inner:") with timer("SimpleRn"): @@ -121,9 +103,3 @@ def element(self, *args, **kwargs): for _ in range(iterations): result = opt_space.inner(ox, oz) print("result: {}".format(result)) - -if cuda_supported: - with timer("odl cuda"): - for _ in range(iterations): - result = cu_space.inner(cx, cz) - print("result: {}".format(result)) diff --git a/odl/solvers/functional/default_functionals.py b/odl/solvers/functional/default_functionals.py index f0cd122ecad..1cea6dfeb68 100644 --- a/odl/solvers/functional/default_functionals.py +++ b/odl/solvers/functional/default_functionals.py @@ -555,8 +555,9 @@ def proximal(self): elif self.exponent == 1: return proximal_convex_conj_linfty(space=self.domain) else: - raise NotImplementedError('`proximal` only implemented for p=1, ' - 'p=2 or p=inf') + raise NotImplementedError( + '`proximal` only implemented for p=2 and p=inf' + ) def __repr__(self): """Return ``repr(self)``.""" diff --git a/odl/solvers/nonsmooth/proximal_operators.py b/odl/solvers/nonsmooth/proximal_operators.py index 974a6801291..95672289033 100644 --- a/odl/solvers/nonsmooth/proximal_operators.py +++ b/odl/solvers/nonsmooth/proximal_operators.py @@ -169,10 +169,10 @@ def convex_conj_prox_factory(sigma): # prox_factory accepts stepsize objects of the type given by sigma. space = prox_factory(sigma).domain - mult_inner = MultiplyOperator(space, 1.0 / sigma) + mult_inner = MultiplyOperator(space, 1 / sigma) mult_outer = MultiplyOperator(space, sigma) result = (IdentityOperator(space) - - mult_outer * prox_factory(1.0 / sigma) * mult_inner) + mult_outer * prox_factory(1 / sigma) * mult_inner) return result return convex_conj_prox_factory @@ -292,7 +292,7 @@ def proximal_arg_scaling(prox_factory, scaling): # unconditionally, but only if the scaling factor is a scalar: if np.isscalar(scaling): if scaling == 0: - return proximal_const_func(prox_factory(1.0).domain) + return proximal_const_func(prox_factory(1).domain) elif scaling.imag != 0: raise ValueError("Complex scaling not supported.") else: @@ -327,7 +327,7 @@ def arg_scaling_prox_factory(sigma): def proximal_quadratic_perturbation(prox_factory, a, u=None): - r"""Calculate the proximal of function F(x) + a * \|x\|^2 + . + r"""Calculate the proximal of function F(x) + a * ||x||^2 + . Parameters ---------- @@ -377,8 +377,9 @@ def proximal_quadratic_perturbation(prox_factory, a, u=None): """ a = float(a) if a < 0: - raise ValueError('scaling parameter muts be non-negative, got {}' - ''.format(a)) + raise ValueError( + 'scaling parameter must be non-negative, got {}'.format(a) + ) def quadratic_perturbation_prox_factory(sigma): r"""Create proximal for the quadratic perturbation with a given sigma. @@ -399,7 +400,7 @@ def quadratic_perturbation_prox_factory(sigma): else: sigma = np.asarray(sigma) - const = 1.0 / np.sqrt(sigma * 2.0 * a + 1) + const = 1 / np.sqrt(2 * sigma * a + 1) prox = proximal_arg_scaling(prox_factory, const)(sigma) space = prox.domain if u is not None: @@ -486,8 +487,7 @@ def proximal_composition_factory(sigma): Id = IdentityOperator(operator.domain) Ir = IdentityOperator(operator.range) prox_muf = proximal(mu * sigma) - return (Id + - (1.0 / mu) * operator.adjoint * ((prox_muf - Ir) * operator)) + return Id + (1 / mu) * operator.adjoint * ((prox_muf - Ir) * operator) return proximal_composition_factory @@ -788,7 +788,7 @@ def _call(self, x, out): else: step = np.infty - if step < 1.0: + if step < 1: self.range.lincomb(1 - step, x, out=out) else: self.range.lincomb(0, out, out=out) @@ -800,7 +800,7 @@ def _call(self, x, out): else: step = np.infty - if step < 1.0: + if step < 1: self.range.lincomb(1 - step, x, step, g, out=out) else: self.range.assign(out, g) @@ -1012,7 +1012,7 @@ def _call(self, x, out): space.lincomb(1, x, 1, tmp, out=out) else: F.multiply(sig, 2 * lam * g, out=out) - space.lincomb.lincomb(1, x, 1, out, out=out) + space.lincomb(1, x, 1, out, out=out) F.divide(out, 1 + 2 * sig * lam, out=out) return ProximalL2Squared @@ -1363,7 +1363,7 @@ def _call(self, x, out): F.divide(diff, denom, out=out) # out = x - ... - space.lincomb(1, x, -1, out, out=out) + space.lincomb(1, x_old, -1, out, out=out) return ProximalL1 @@ -1810,7 +1810,7 @@ def _call(self, x, out): # out = ... + 4*lam*sigma*g # If g is None, it is taken as the one element if g is None: - out += 4.0 * lam * self.sigma + out += 4 * lam * self.sigma else: space.lincomb(1, out, 4 * lam * self.sigma, g, out=out) diff --git a/odl/test/operator/operator_test.py b/odl/test/operator/operator_test.py index 40c93a1494b..2ff4d080d91 100644 --- a/odl/test/operator/operator_test.py +++ b/odl/test/operator/operator_test.py @@ -779,7 +779,10 @@ def test_nonlinear_functional_operators(): assert C(x) == pytest.approx(mat(x / 2.0)) -# test functions to dispatch +# Test functions to dispatch +# First doc line is the true signature +# Second doc line contains `has_out` and `out_optional` booleans +# Third doc line indicates whether the signature is OK for Operator._call def f1(x): """f1(x) False, False @@ -898,7 +901,6 @@ def func(request): def test_function_signature(func): - true_sig = func.__doc__.splitlines()[0].strip() sig = _function_signature(func) assert true_sig == sig @@ -906,17 +908,17 @@ def test_function_signature(func): def test_dispatch_call_args(func): # Unbound functions - true_has, true_opt = eval(func.__doc__.splitlines()[1].strip()) + true_has_out, true_out_opt = eval(func.__doc__.splitlines()[1].strip()) good = func.__doc__.splitlines()[2].strip() == 'good' if good: truespec = getargspec(func) truespec.args.insert(0, 'self') - has, opt, spec = _dispatch_call_args(unbound_call=func) + has_out, out_opt, spec = _dispatch_call_args(unbound_call=func) - assert has == true_has - assert opt == true_opt + assert has_out == true_has_out + assert out_opt == true_out_opt assert spec == truespec else: with pytest.raises(ValueError): @@ -926,6 +928,7 @@ def test_dispatch_call_args(func): def test_dispatch_call_args_class(): # Two sneaky classes whose _call method would pass the signature check + # because it looks okay from the second argument on class WithStaticMethod(object): @staticmethod def _call(x, y, out): diff --git a/odl/util/testutils.py b/odl/util/testutils.py index d38d3d95ded..489863a1705 100644 --- a/odl/util/testutils.py +++ b/odl/util/testutils.py @@ -308,6 +308,8 @@ def simple_fixture(name, params, fmt=None): # Helpers to generate data + +# TODO(kohr-h): rename to noise_np_array def noise_array(space): """Generate a white noise array for ``space``.