From a5b74cc3b315ab88b861a7c09f1f8ac8ed4a3f27 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Thu, 5 Oct 2017 13:54:41 +0200 Subject: [PATCH] WIP: improve repr by a lot - add helper `_separators` to determine string separators that keep a `repr` string within a given line width - split repr string rendering into two functions, one to determine the parts to be printed and one to compile the final string (with separators) - `attribute_repr_string` helper to class attributes - `element_repr_string` helper to print space elements - make repr printing of product space elements more consistent - improve default repr of `Operator` (now actually good by default, but unspecific) - use Numpy's print options as a single point to configure printing - add doctests to lots of `__repr__` methods Further miscellaneous changes: - fix bug of `uniform_discr_fromdiscr` ignoring dtype - sort imports using the `isort` tool (mode `-m 4`) - add `asweighted` to base_tensors to make weighted and unweighted spaces compatible in tensor operators - remove base classes for operator and adjoint in favor of inner adjoint classes (less code) - make most default operators take `domain` and `range` and implement the corresponding adjoints/inverses etc. - implement `noise_array` for product spaces that are not power spaces - remove the `almost_equal` test util - improve a number of doctests (use `odl.` calls, better names etc.) - make comparison of arrays in tests faster - implement broadcasting of shapes (M, N) and (N,) in product spaces - remove some remaining doc glitches, mostly `FnBase` TODO: - go through rest of the library - replace `almost_equal` by `pytest.approx` - tests for pspace broadcasting - make `uniform_discr_fromdiscr` respect weighting - fix failing unit tests --- .gitignore | 3 + odl/README.md | 2 +- odl/__init__.py | 3 + odl/contrib/fom/supervised.py | 36 +- odl/contrib/tensorflow/README.md | 2 +- odl/contrib/theano/layer.py | 4 +- odl/deform/linearized.py | 79 +- odl/discr/diff_ops.py | 223 +++--- odl/discr/discr_mappings.py | 138 +++- odl/discr/discr_ops.py | 210 +++-- odl/discr/discretization.py | 9 +- odl/discr/grid.py | 153 ++-- odl/discr/lp_discr.py | 101 ++- odl/discr/partition.py | 87 ++- odl/operator/default_ops.py | 723 ++++++++++++------ odl/operator/operator.py | 192 +++-- odl/operator/oputils.py | 7 +- odl/operator/pspace_ops.py | 535 +++++++------ odl/operator/tensor_ops.py | 511 +++++++------ odl/phantom/noise.py | 2 +- odl/set/domain.py | 111 +-- odl/set/sets.py | 95 ++- odl/solvers/functional/default_functionals.py | 2 +- odl/solvers/nonsmooth/proximal_operators.py | 2 +- odl/space/base_tensors.py | 27 +- odl/space/npy_tensors.py | 14 +- odl/space/pspace.py | 225 +++--- odl/test/operator/oputils_test.py | 12 +- odl/util/graphics.py | 7 +- odl/util/normalize.py | 4 +- odl/util/npy_compat.py | 4 +- odl/util/numerics.py | 4 +- odl/util/pytest_plugins.py | 6 +- odl/util/testutils.py | 131 ++-- odl/util/ufuncs.py | 7 +- odl/util/utility.py | 530 ++++++++++++- odl/util/vectorization.py | 5 +- 37 files changed, 2680 insertions(+), 1526 deletions(-) diff --git a/.gitignore b/.gitignore index 7ed03df106a..b6434b400f0 100644 --- a/.gitignore +++ b/.gitignore @@ -90,4 +90,7 @@ target/ *.ipynb .ipynb_checkpoints +# VS Code files +.vscode/ + # Documentation build files diff --git a/odl/README.md b/odl/README.md index 0ba62d860f5..7e05e316c45 100644 --- a/odl/README.md +++ b/odl/README.md @@ -21,7 +21,7 @@ This is a brief description of the content of each submodule, see the individual * [phantom](phantom) Standardized test images. Functions for generating standardized test examples such as `shepp_logan`. * [set](set) Sets of objects. Defines the abstract class `Set` and `LinearSpace` as well as some concrete implementations such as `RealNumbers`. * [solvers](solvers) Solution of equations and optimization. Contains both general solvers for problems of the form `A(x) = b` where `A` is an `Operator` as well as solvers of minimization problems. In addition, it defines the class `Functional` with several concrete implementations such as `L2Norm`. -* [space](space) Concrete vector spaces. Contains concrete implementations of `LinearSpace`, including `NumpyFnBase` and `ProductSpace`. +* [space](space) Concrete vector spaces. Contains concrete implementations of `LinearSpace`, including `NumpyTensorSpace` and `ProductSpace`. * [test](test) Tests for the ODL package. This contains automated tests for all other ODL functionality. In general, users should not be calling anything from this submoduel. * [tomo](tomo) Tomography. Defines the operator `RayTransform` as well as `Geometry` along with subclasses and utilities. Also defines problem dependent direct reconstruction such as `fbp_op`. * [trafos](trafos) Transformations between spaces. Defines `FourierTransform` and `WaveletTransform`. diff --git a/odl/__init__.py b/odl/__init__.py index 92b9436dea9..123379c9388 100644 --- a/odl/__init__.py +++ b/odl/__init__.py @@ -38,6 +38,9 @@ except TypeError: pass +# Set printing linewidth to 71 to allow method docstrings to not extend +# beyond 79 characters (2 times indent of 4) +np.set_printoptions(linewidth=71) # Propagate names defined in` __all__` of all "core" subpackages into # the top-level namespace diff --git a/odl/contrib/fom/supervised.py b/odl/contrib/fom/supervised.py index 5c2a899e4a4..0e3ca404553 100644 --- a/odl/contrib/fom/supervised.py +++ b/odl/contrib/fom/supervised.py @@ -27,9 +27,9 @@ def mean_squared_error(data, ground_truth, mask=None, normalized=False): Parameters ---------- - data : `FnBaseVector` + data : `Tensor` Input data to compare to the ground truth. - ground_truth : `FnBaseVector` + ground_truth : `Tensor` Reference to compare ``data`` to. mask : `array-like`, optional If given, ``data * mask`` is compared to ``ground_truth * mask``. @@ -85,9 +85,9 @@ def mean_absolute_error(data, ground_truth, mask=None, normalized=False): Parameters ---------- - data : `FnBaseVector` + data : `Tensor` Input data to compare to the ground truth. - ground_truth : `FnBaseVector` + ground_truth : `Tensor` Reference to compare ``data`` to. mask : `array-like`, optional If given, ``data * mask`` is compared to ``ground_truth * mask``. @@ -137,9 +137,9 @@ def mean_value_difference(data, ground_truth, mask=None, normalized=False): Parameters ---------- - data : `FnBaseVector` + data : `Tensor` Input data to compare to the ground truth. - ground_truth : `FnBaseVector` + ground_truth : `Tensor` Reference to compare ``data`` to. mask : `array-like`, optional If given, ``data * mask`` is compared to ``ground_truth * mask``. @@ -198,9 +198,9 @@ def standard_deviation_difference(data, ground_truth, mask=None, Parameters ---------- - data : `FnBaseVector` + data : `Tensor` Input data to compare to the ground truth. - ground_truth : `FnBaseVector` + ground_truth : `Tensor` Reference to compare ``data`` to. mask : `array-like`, optional If given, ``data * mask`` is compared to ``ground_truth * mask``. @@ -268,9 +268,9 @@ def range_difference(data, ground_truth, mask=None, normalized=False): Parameters ---------- - data : `FnBaseVector` + data : `Tensor` Input data to compare to the ground truth. - ground_truth : `FnBaseVector` + ground_truth : `Tensor` Reference to compare ``data`` to. mask : `array-like`, optional Binary mask or index array to define ROI in which FOM evaluation @@ -336,9 +336,9 @@ def blurring(data, ground_truth, mask=None, normalized=False, Parameters ---------- - data : `FnBaseVector` + data : `Tensor` Input data to compare to the ground truth. - ground_truth : `FnBaseVector` + ground_truth : `Tensor` Reference to compare ``data`` to. mask : `array-like`, optional Binary mask to define ROI in which FOM evaluation is performed. @@ -406,9 +406,9 @@ def false_structures(data, ground_truth, mask=None, normalized=False, Parameters ---------- - data : `FnBaseVector` + data : `Tensor` Input data to compare to the ground truth. - ground_truth : `FnBaseVector` + ground_truth : `Tensor` Reference to compare ``data`` to. mask : `array-like`, optional Binary mask to define ROI in which FOM evaluation is performed. @@ -471,9 +471,9 @@ def ssim(data, ground_truth, Parameters ---------- - data : `FnBaseVector` + data : `Tensor` Input data to compare to the ground truth. - ground_truth : `FnBaseVector` + ground_truth : `Tensor` Reference to compare ``data`` to. size : odd int Size in elements per axis of the Gaussian window that is used @@ -536,9 +536,9 @@ def psnr(data, ground_truth, normalized=False): Parameters ---------- - data : `FnBaseVector` + data : `Tensor` Input data to compare to the ground truth. - ground_truth : `FnBaseVector` + ground_truth : `Tensor` Reference to compare ``data`` to. normalized : bool If true, normalize ``data`` and ``ground_truth`` to have the diff --git a/odl/contrib/tensorflow/README.md b/odl/contrib/tensorflow/README.md index 43ef2864aa0..dc6bcfcad87 100644 --- a/odl/contrib/tensorflow/README.md +++ b/odl/contrib/tensorflow/README.md @@ -8,7 +8,7 @@ This package contains ODL functionality related to [TensorFlow](https://www.tens This allows using tensorflow neural networks as operators in ODL. * `as_tensorflow_layer` in [layer.py](layer.py) wraps an ODL operator into a tensorflow layer. This allows using arbitrary ODL operators inside tensorflow neural networks. -* `TensorflowSpace` in [space.py](space.py) is a `FnBase`-type space which uses tensorflow as a backend. +* `TensorflowSpace` in [space.py](space.py) is a `TensorSpace`-type space which uses tensorflow as a backend. ## Example usage diff --git a/odl/contrib/theano/layer.py b/odl/contrib/theano/layer.py index 52f0438e8fc..127ec06d909 100644 --- a/odl/contrib/theano/layer.py +++ b/odl/contrib/theano/layer.py @@ -37,8 +37,8 @@ def __init__(self, operator): Parameters ---------- operator : `Operator` - The operator that should be wrapped, must map `FnBase` spaces to - `FnBase` spaces. + The operator that should be wrapped, must map `TensorSpace`'s to + `TensorSpace`'s. Examples -------- diff --git a/odl/deform/linearized.py b/odl/deform/linearized.py index ef92091df87..62f9fb142f9 100644 --- a/odl/deform/linearized.py +++ b/odl/deform/linearized.py @@ -8,14 +8,14 @@ """Operators and functions for linearized deformation.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + import numpy as np -from odl.discr import DiscreteLp, Gradient, Divergence +from odl.discr import DiscreteLp, Divergence, Gradient from odl.operator import Operator, PointwiseInner from odl.space import ProductSpace -from odl.util import signature_string, indent - +from odl.util import repr_string, signature_string_parts __all__ = ('LinDeformFixedTempl', 'LinDeformFixedDisp', 'linear_deform') @@ -146,10 +146,10 @@ def __init__(self, template, domain=None): >>> space = odl.uniform_discr(0, 1, 5, interp='nearest') >>> template = space.element([0, 0, 1, 0, 0]) - >>> op = LinDeformFixedTempl(template) + >>> op = odl.deform.LinDeformFixedTempl(template) >>> disp_field = [[0, 0, 0, -0.2, 0]] - >>> print(op(disp_field)) - [ 0., 0., 1., 1., 0.] + >>> op(disp_field) + uniform_discr(0.0, 1.0, 5).element([ 0., 0., 1., 1., 0.]) The result depends on the chosen interpolation. With 'linear' interpolation and an offset equal to half the distance between two @@ -157,10 +157,12 @@ def __init__(self, template, domain=None): >>> space = odl.uniform_discr(0, 1, 5, interp='linear') >>> template = space.element([0, 0, 1, 0, 0]) - >>> op = LinDeformFixedTempl(template) + >>> op = odl.deform.LinDeformFixedTempl(template) >>> disp_field = [[0, 0, 0, -0.1, 0]] - >>> print(op(disp_field)) - [ 0. , 0. , 1. , 0.5, 0. ] + >>> op(disp_field) + uniform_discr(0.0, 1.0, 5, interp='linear').element( + [ 0. , 0. , 1. , 0.5, 0. ] + ) """ space = getattr(template, 'space', None) if not isinstance(space, DiscreteLp): @@ -231,11 +233,24 @@ def derivative(self, displacement): return PointwiseInner(self.domain, def_grad) def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> space = odl.uniform_discr(0, 1, 5) + >>> template = space.element([0, 0, 1, 0, 0]) + >>> op = odl.deform.LinDeformFixedTempl(template) + >>> op + LinDeformFixedTempl( + uniform_discr(0.0, 1.0, 5).element([ 0., 0., 1., 0., 0.]) + ) + """ posargs = [self.template] - optargs = [('domain', self.domain, self.template.space)] - inner_str = signature_string(posargs, optargs, mod='!r', sep=',\n') - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) + optargs = [('domain', self.domain, + self.template.space.real_space.tangent_bundle)] + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) class LinDeformFixedDisp(Operator): @@ -297,10 +312,10 @@ def __init__(self, displacement, templ_space=None): >>> space = odl.uniform_discr(0, 1, 5) >>> disp_field = space.tangent_bundle.element([[0, 0, 0, -0.2, 0]]) - >>> op = LinDeformFixedDisp(disp_field) + >>> op = odl.deform.LinDeformFixedDisp(disp_field) >>> template = [0, 0, 1, 0, 0] - >>> print(op([0, 0, 1, 0, 0])) - [ 0., 0., 1., 1., 0.] + >>> op(template) + uniform_discr(0.0, 1.0, 5).element([ 0., 0., 1., 1., 0.]) The result depends on the chosen interpolation. With 'linear' interpolation and an offset equal to half the distance between two @@ -308,10 +323,12 @@ def __init__(self, displacement, templ_space=None): >>> space = odl.uniform_discr(0, 1, 5, interp='linear') >>> disp_field = space.tangent_bundle.element([[0, 0, 0, -0.1, 0]]) - >>> op = LinDeformFixedDisp(disp_field) + >>> op = odl.deform.LinDeformFixedDisp(disp_field) >>> template = [0, 0, 1, 0, 0] - >>> print(op(template)) - [ 0. , 0. , 1. , 0.5, 0. ] + >>> op(template) + uniform_discr(0.0, 1.0, 5, interp='linear').element( + [ 0. , 0. , 1. , 0.5, 0. ] + ) """ space = getattr(displacement, 'space', None) if not isinstance(space, ProductSpace): @@ -374,11 +391,27 @@ def adjoint(self): return jacobian_det * self.inverse def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> space = odl.uniform_discr(0, 1, 5) + >>> disp_field = space.tangent_bundle.element([[0, 0, 0, -0.2, 0]]) + >>> op = odl.deform.LinDeformFixedDisp(disp_field) + >>> op + LinDeformFixedDisp( + ProductSpace(uniform_discr(0.0, 1.0, 5), 1).element( + [ + [ 0. , 0. , 0. , -0.2, 0. ] + ] + ) + ) + """ posargs = [self.displacement] optargs = [('templ_space', self.domain, self.displacement.space[0])] - inner_str = signature_string(posargs, optargs, mod='!r', sep=',\n') - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) if __name__ == '__main__': diff --git a/odl/discr/diff_ops.py b/odl/discr/diff_ops.py index fad6dbb939f..a821b8f928a 100644 --- a/odl/discr/diff_ops.py +++ b/odl/discr/diff_ops.py @@ -6,16 +6,18 @@ # v. 2.0. If a copy of the MPL was not distributed with this file, You can # obtain one at https://mozilla.org/MPL/2.0/. -"""Operators defined for tensor fields.""" +"""Differential operators.""" + +from __future__ import absolute_import, division, print_function -from __future__ import print_function, division, absolute_import import numpy as np from odl.discr.lp_discr import DiscreteLp from odl.operator.tensor_ops import PointwiseTensorFieldOperator from odl.space import ProductSpace -from odl.util import writable_array, signature_string, indent - +from odl.util import ( + REPR_PRECISION, npy_printoptions, repr_string, signature_string_parts, + writable_array) __all__ = ('PartialDerivative', 'Gradient', 'Divergence', 'Laplacian') @@ -45,11 +47,7 @@ class PartialDerivative(PointwiseTensorFieldOperator): - """Calculate the discrete partial derivative along a given axis. - - Calls helper function `finite_diff` to calculate finite difference. - Preserves the shape of the underlying grid. - """ + """Partial derivative operator along a given axis.""" def __init__(self, domain, axis, range=None, method='forward', pad_mode='constant', pad_const=0): @@ -98,9 +96,9 @@ def __init__(self, domain, axis, range=None, method='forward', -------- >>> f = np.array([[ 0., 1., 2., 3., 4.], ... [ 0., 2., 4., 6., 8.]]) - >>> discr = odl.uniform_discr([0, 0], [2, 1], f.shape) - >>> par_deriv = PartialDerivative(discr, axis=0, pad_mode='order1') - >>> par_deriv(f) + >>> space = odl.uniform_discr([0, 0], [2, 1], f.shape) + >>> pderiv = odl.PartialDerivative(space, axis=0, pad_mode='order1') + >>> pderiv(f) uniform_discr([ 0., 0.], [ 2., 1.], (2, 5)).element( [[ 0., 1., 2., 3., 4.], [ 0., 1., 2., 3., 4.]] @@ -177,40 +175,40 @@ def adjoint(self): self.pad_const) def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> space = odl.uniform_discr([0, 0], [2, 1], (2, 5)) + >>> pderiv = odl.PartialDerivative(space, axis=0, pad_mode='order1') + >>> pderiv + PartialDerivative( + uniform_discr([ 0., 0.], [ 2., 1.], (2, 5)), + axis=0, pad_mode='order1' + ) + """ posargs = [self.domain] optargs = [('axis', self.axis, None), ('range', self.range, self.domain), ('method', self.method, 'forward'), ('pad_mode', self.pad_mode, 'constant'), ('pad_const', self.pad_const, 0)] - inner_str = signature_string(posargs, optargs, - sep=',\n', mod=['!r', '']) - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) - - def __str__(self): - """Return ``str(self)``.""" - dom_ran_str = '\n-->\n'.join([repr(self.domain), repr(self.range)]) - return '{}:\n{}'.format(self.__class__.__name__, indent(dom_ran_str)) + optmod = ['', '!r', '', '', ''] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs, + mod=['!r', optmod]) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=True) class Gradient(PointwiseTensorFieldOperator): - """Spatial gradient operator for `DiscreteLp` spaces. - - Calls helper function `finite_diff` to calculate each component of the - resulting product space element. For the adjoint of the `Gradient` - operator, zero padding is assumed to match the negative `Divergence` - operator - """ + """Spatial gradient operator for `DiscreteLp` spaces.""" def __init__(self, domain=None, range=None, method='forward', pad_mode='constant', pad_const=0): """Initialize a new instance. - Zero padding is assumed for the adjoint of the `Gradient` - operator to match negative `Divergence` operator. - Parameters ---------- domain : `DiscreteLp`, optional @@ -253,25 +251,23 @@ def __init__(self, domain=None, range=None, method='forward', >>> dom = odl.uniform_discr([0, 0], [1, 1], (10, 20)) >>> ran = odl.ProductSpace(dom, dom.ndim) # 2-dimensional - >>> grad_op = Gradient(dom) - >>> grad_op.range == ran - True - >>> grad_op2 = Gradient(range=ran) - >>> grad_op2.domain == dom + >>> grad = odl.Gradient(dom) + >>> grad.range == ran True - >>> grad_op3 = Gradient(domain=dom, range=ran) - >>> grad_op3.domain == dom + >>> grad = odl.Gradient(range=ran) + >>> grad.domain == dom True - >>> grad_op3.range == ran + >>> grad = odl.Gradient(domain=dom, range=ran) + >>> grad.domain == dom and grad.range == ran True Calling the operator: >>> data = np.array([[ 0., 1., 2., 3., 4.], ... [ 0., 2., 4., 6., 8.]]) - >>> discr = odl.uniform_discr([0, 0], [2, 5], data.shape) - >>> f = discr.element(data) - >>> grad = Gradient(discr) + >>> space = odl.uniform_discr([0, 0], [2, 5], data.shape) + >>> f = space.element(data) + >>> grad = odl.Gradient(space) >>> grad_f = grad(f) >>> grad_f[0] uniform_discr([ 0., 0.], [ 2., 5.], (2, 5)).element( @@ -401,39 +397,39 @@ def adjoint(self): pad_const=self.pad_const) def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> space = odl.uniform_discr([0, 0], [2, 5], (2, 5)) + >>> grad = odl.Gradient(space, pad_mode='order1') + >>> grad + Gradient( + uniform_discr([ 0., 0.], [ 2., 5.], (2, 5)), + pad_mode='order1' + ) + """ posargs = [self.domain] optargs = [('range', self.range, self.domain ** self.domain.ndim), ('method', self.method, 'forward'), ('pad_mode', self.pad_mode, 'constant'), ('pad_const', self.pad_const, 0)] - inner_str = signature_string(posargs, optargs, - sep=[',\n', ', ', ',\n'], - mod=['!r', '']) - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) - - def __str__(self): - """Return ``str(self)``.""" - dom_ran_str = '\n-->\n'.join([repr(self.domain), repr(self.range)]) - return '{}:\n{}'.format(self.__class__.__name__, indent(dom_ran_str)) + optmod = ['!r', '', '', ''] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs, + mod=['!r', optmod]) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=True) class Divergence(PointwiseTensorFieldOperator): - """Divergence operator for `DiscreteLp` spaces. - - Calls helper function `finite_diff` for each component of the input - product space vector. For the adjoint of the `Divergence` operator to - match the negative `Gradient` operator implicit zero is assumed. - """ + """Divergence operator for `DiscreteLp` spaces.""" def __init__(self, domain=None, range=None, method='forward', pad_mode='constant', pad_const=0): """Initialize a new instance. - Zero padding is assumed for the adjoint of the `Divergence` - operator to match the negative `Gradient` operator. - Parameters ---------- domain : power space of `DiscreteLp`, optional @@ -475,16 +471,16 @@ def __init__(self, domain=None, range=None, method='forward', >>> ran = odl.uniform_discr([0, 0], [3, 5], (3, 5)) >>> dom = odl.ProductSpace(ran, ran.ndim) # 2-dimensional - >>> div = Divergence(dom) + >>> div = odl.Divergence(dom) >>> div.range == ran True - >>> div2 = Divergence(range=ran) - >>> div2.domain == dom + >>> div = odl.Divergence(range=ran) + >>> div.domain == dom True - >>> div3 = Divergence(domain=dom, range=ran) - >>> div3.domain == dom + >>> div = odl.Divergence(domain=dom, range=ran) + >>> div.domain == dom True - >>> div3.range == ran + >>> div.range == ran True Call the operator: @@ -493,11 +489,12 @@ def __init__(self, domain=None, range=None, method='forward', ... [1., 2., 3., 4., 5.], ... [2., 3., 4., 5., 6.]]) >>> f = div.domain.element([data, data]) - >>> div_f = div(f) - >>> print(div_f) - [[ 2., 2., 2., 2., -3.], - [ 2., 2., 2., 2., -4.], - [ -1., -2., -3., -4., -12.]] + >>> div(f) + uniform_discr([ 0., 0.], [ 3., 5.], (3, 5)).element( + [[ 2., 2., 2., 2., -3.], + [ 2., 2., 2., 2., -4.], + [ -1., -2., -3., -4., -12.]] + ) Verify adjoint: @@ -610,32 +607,34 @@ def adjoint(self): pad_mode=_ADJ_PADDING[self.pad_mode]) def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> ran = odl.uniform_discr([0, 0], [3, 5], (3, 5)) + >>> div = odl.Divergence(range=ran, pad_mode='order1') + >>> div + Divergence( + ProductSpace(uniform_discr([ 0., 0.], [ 3., 5.], (3, 5)), 2), + pad_mode='order1' + ) + """ posargs = [self.domain] optargs = [('range', self.range, self.domain[0]), ('method', self.method, 'forward'), ('pad_mode', self.pad_mode, 'constant'), ('pad_const', self.pad_const, 0)] - inner_str = signature_string(posargs, optargs, - sep=[',\n', ', ', ',\n'], - mod=['!r', '']) - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) - - def __str__(self): - """Return ``str(self)``.""" - dom_ran_str = '\n-->\n'.join([repr(self.domain), repr(self.range)]) - return '{}:\n{}'.format(self.__class__.__name__, indent(dom_ran_str)) + optmod = ['!r', '', '', ''] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs, + mod=['!r', optmod]) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=True) class Laplacian(PointwiseTensorFieldOperator): - """Spatial Laplacian operator for `DiscreteLp` spaces. - - Calls helper function `finite_diff` to calculate each component of the - resulting product space vector. - - Outside the domain zero padding is assumed. - """ + """Spatial Laplacian operator for `DiscreteLp` spaces.""" def __init__(self, domain, range=None, pad_mode='constant', pad_const=0): """Initialize a new instance. @@ -644,6 +643,9 @@ def __init__(self, domain, range=None, pad_mode='constant', pad_const=0): ---------- domain : `DiscreteLp` Space of elements which the operator is acting on. + range : `DiscreteLp`, optional + Space of elements to which the operator maps. + Default: ``domain`` pad_mode : string, optional The padding mode to use outside the domain. @@ -677,7 +679,7 @@ def __init__(self, domain, range=None, pad_mode='constant', pad_const=0): ... [ 0., 0., 0.]]) >>> space = odl.uniform_discr([0, 0], [3, 3], [3, 3]) >>> f = space.element(data) - >>> lap = Laplacian(space) + >>> lap = odl.Laplacian(space) >>> lap(f) uniform_discr([ 0., 0.], [ 3., 3.], (3, 3)).element( [[ 0., 1., 0.], @@ -692,9 +694,6 @@ def __init__(self, domain, range=None, pad_mode='constant', pad_const=0): if range is None: range = domain - super(Laplacian, self).__init__( - domain, range, base_space=domain, linear=True) - self.pad_mode, pad_mode_in = str(pad_mode).lower(), pad_mode if pad_mode not in _SUPPORTED_PAD_MODES: raise ValueError('`pad_mode` {} not understood' @@ -705,6 +704,10 @@ def __init__(self, domain, range=None, pad_mode='constant', pad_const=0): raise ValueError('`pad_mode` {} not implemented for Laplacian.' ''.format(pad_mode_in)) + linear = not (self.pad_mode == 'constant' and pad_const != 0) + super(Laplacian, self).__init__( + domain, range, base_space=domain, linear=linear) + self.pad_const = self.domain.field.element(pad_const) def _call(self, x, out=None): @@ -765,24 +768,36 @@ def adjoint(self): The laplacian is self-adjoint, so this returns ``self``. """ + if not self.is_linear: + raise ValueError('operator with nonzero pad_const ({}) is not' + ' linear and has no adjoint' + ''.format(self.pad_const)) return Laplacian(self.range, self.domain, pad_mode=self.pad_mode, pad_const=0) def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> space = odl.uniform_discr([0, 0], [3, 3], [3, 3]) + >>> lap = odl.Laplacian(space, pad_mode='symmetric') + >>> lap + Laplacian( + uniform_discr([ 0., 0.], [ 3., 3.], (3, 3)), + pad_mode='symmetric' + ) + """ posargs = [self.domain] - optargs = [('range', self.range, self.domain ** self.domain.ndim), + optargs = [('range', self.range, self.domain), ('pad_mode', self.pad_mode, 'constant'), ('pad_const', self.pad_const, 0)] - inner_str = signature_string(posargs, optargs, - sep=[',\n', ', ', ',\n'], - mod=['!r', '']) - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) - - def __str__(self): - """Return ``str(self)``.""" - dom_ran_str = '\n-->\n'.join([repr(self.domain), repr(self.range)]) - return '{}:\n{}'.format(self.__class__.__name__, indent(dom_ran_str)) + optmod = ['!r', '', ''] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs, + mod=['!r', optmod]) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=True) def finite_diff(f, axis, dx=1.0, method='forward', out=None, **kwargs): diff --git a/odl/discr/discr_mappings.py b/odl/discr/discr_mappings.py index 669f3b3c062..e182659d92b 100644 --- a/odl/discr/discr_mappings.py +++ b/odl/discr/discr_mappings.py @@ -12,20 +12,21 @@ operators. """ -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + from builtins import object from itertools import product + import numpy as np -from odl.operator import Operator from odl.discr.partition import RectPartition -from odl.space.base_tensors import TensorSpace +from odl.operator import Operator from odl.space import FunctionSpace +from odl.space.base_tensors import TensorSpace from odl.util import ( - is_valid_input_meshgrid, out_shape_from_array, out_shape_from_meshgrid, - is_string, is_numeric_dtype, signature_string, indent, dtype_repr, - writable_array) - + dtype_repr, is_numeric_dtype, is_string, is_valid_input_meshgrid, + out_shape_from_array, out_shape_from_meshgrid, repr_string, + signature_string_parts, writable_array) __all__ = ('FunctionSpaceMapping', 'PointCollocation', 'NearestInterpolation', 'LinearInterpolation', @@ -178,14 +179,13 @@ def __init__(self, fspace, partition, tspace): Partition the rectangle by a rectilinear grid: - >>> rect = odl.IntervalProd([1, 3], [2, 5]) >>> grid = odl.RectGrid([1, 2], [3, 4, 5]) >>> partition = odl.RectPartition(rect, grid) >>> tspace = odl.rn(grid.shape) Finally create the operator and test it on a function: - >>> coll_op = PointCollocation(fspace, partition, tspace) + >>> coll_op = odl.PointCollocation(fspace, partition, tspace) >>> >>> # Properly vectorized function >>> func_elem = fspace.element(lambda x: x[0] - x[1]) @@ -264,12 +264,26 @@ def _call(self, func, out=None, **kwargs): return out def __repr__(self): - """Return ``repr(self)``.""" - posargs = [self.range, self.grid, self.domain] - inner_str = signature_string(posargs, [], - sep=[',\n', ', ', ',\n'], - mod=['!r', '']) - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) + """Return ``repr(self)``. + + Examples + -------- + >>> rect = odl.IntervalProd([1, 3], [2, 5]) + >>> fspace = odl.FunctionSpace(rect) + >>> part = odl.uniform_partition_fromintv(rect, shape=(4, 2)) + >>> tspace = odl.rn(part.shape) + >>> coll_op = odl.PointCollocation(fspace, part, tspace) + >>> coll_op + PointCollocation( + FunctionSpace(IntervalProd([ 1., 3.], [ 2., 5.])), + uniform_partition([ 1., 3.], [ 2., 5.], (4, 2)), + rn((4, 2)) + ) + """ + posargs = [self.domain, self.partition, self.range] + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) class NearestInterpolation(FunctionSpaceMapping): @@ -335,7 +349,7 @@ def __init__(self, fspace, partition, tspace, variant='left'): Now we initialize the operator and test it with some points: >>> tspace = odl.tensor_space(part.shape, dtype='U1') - >>> interp_op = NearestInterpolation(fspace, part, tspace) + >>> interp_op = odl.NearestInterpolation(fspace, part, tspace) >>> values = np.array([['m', 'y'], ... ['s', 't'], ... ['r', 'i'], @@ -399,13 +413,27 @@ def nearest(arg, out=None): return self.range.element(nearest, vectorized=True) def __repr__(self): - """Return ``repr(self)``.""" - posargs = [self.range, self.grid, self.domain] + """Return ``repr(self)``. + + Examples + -------- + >>> rect = odl.IntervalProd([0, 0], [1, 1]) + >>> fspace = odl.FunctionSpace(rect) + >>> part = odl.uniform_partition_fromintv(rect, shape=(4, 2)) + >>> tspace = odl.rn(part.shape) + >>> interp_op = odl.NearestInterpolation(fspace, part, tspace) + >>> interp_op + NearestInterpolation( + FunctionSpace(IntervalProd([ 0., 0.], [ 1., 1.])), + uniform_partition([ 0., 0.], [ 1., 1.], (4, 2)), + rn((4, 2)) + ) + """ + posargs = [self.range, self.partition, self.domain] optargs = [('variant', self.variant, 'left')] - inner_str = signature_string(posargs, optargs, - sep=[',\n', ', ', ',\n'], - mod=['!r', '']) - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) + inner_parts = signature_string_parts(posargs, optargs, mod=['!r', '']) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) class LinearInterpolation(FunctionSpaceMapping): @@ -453,12 +481,26 @@ def linear(arg, out=None): return self.range.element(linear, vectorized=True) def __repr__(self): - """Return ``repr(self)``.""" - posargs = [self.range, self.grid, self.domain] - inner_str = signature_string(posargs, [], - sep=[',\n', ', ', ',\n'], - mod=['!r', '']) - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) + """Return ``repr(self)``. + + Examples + -------- + >>> rect = odl.IntervalProd([0, 0], [1, 1]) + >>> fspace = odl.FunctionSpace(rect) + >>> part = odl.uniform_partition_fromintv(rect, shape=(4, 2)) + >>> tspace = odl.rn(part.shape) + >>> interp_op = odl.LinearInterpolation(fspace, part, tspace) + >>> interp_op + LinearInterpolation( + FunctionSpace(IntervalProd([ 0., 0.], [ 1., 1.])), + uniform_partition([ 0., 0.], [ 1., 1.], (4, 2)), + rn((4, 2)) + ) + """ + posargs = [self.range, self.partition, self.domain] + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) class PerAxisInterpolation(FunctionSpaceMapping): @@ -543,8 +585,8 @@ def __init__(self, fspace, partition, tspace, schemes, nn_variants='left'): 'with `interp={!r}' ''.format(i, schemes_in[i])) - self.__schemes = schemes - self.__nn_variants = nn_variants + self.__schemes = tuple(schemes) + self.__nn_variants = tuple(nn_variants) @property def schemes(self): @@ -590,19 +632,34 @@ def per_axis_interp(arg, out=None): return self.range.element(per_axis_interp, vectorized=True) def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> rect = odl.IntervalProd([0, 0], [1, 1]) + >>> fspace = odl.FunctionSpace(rect) + >>> part = odl.uniform_partition_fromintv(rect, shape=(4, 2)) + >>> tspace = odl.rn(part.shape) + >>> interp_op = odl.PerAxisInterpolation(fspace, part, tspace, + ... schemes=['linear', 'nearest']) + >>> interp_op + PerAxisInterpolation( + FunctionSpace(IntervalProd([ 0., 0.], [ 1., 1.])), + uniform_grid([ 0.125, 0.25 ], [ 0.875, 0.75 ], (4, 2)), + rn((4, 2)), + schemes=('linear', 'nearest') + ) + """ if all(scm == self.schemes[0] for scm in self.schemes): schemes = self.schemes[0] else: schemes = self.schemes - posargs = [self.range, self.grid, self.domain, schemes] + posargs = [self.range, self.grid, self.domain] + optargs = [('schemes', schemes, None)] nn_relevant = [x for x in self.nn_variants if x is not None] - if not nn_relevant: - # No NN axes, ignore nn_variants - optargs = [] - else: + if nn_relevant: # Use single string if all are equal, one per axis otherwise first_relevant = nn_relevant[0] @@ -611,12 +668,11 @@ def __repr__(self): else: variants = self.nn_variants - optargs = [('nn_variants', variants, 'left')] + optargs.append(('nn_variants', variants, 'left')) - inner_str = signature_string(posargs, optargs, - sep=[',\n', ', ', ',\n'], - mod=['!r', '']) - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) + inner_parts = signature_string_parts(posargs, optargs, mod=['!r', '']) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) class _Interpolator(object): diff --git a/odl/discr/discr_ops.py b/odl/discr/discr_ops.py index cf2f2d1c858..6f954d27d81 100644 --- a/odl/discr/discr_ops.py +++ b/odl/discr/discr_ops.py @@ -8,7 +8,8 @@ """Operators defined on `DiscreteLp`.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + import numpy as np from odl.discr import DiscreteLp, uniform_partition @@ -16,10 +17,11 @@ from odl.set import IntervalProd from odl.space import FunctionSpace, tensor_space from odl.util import ( - normalized_scalar_param_list, safe_int_conv, writable_array, resize_array) + REPR_PRECISION, attribute_repr_string, normalized_scalar_param_list, + npy_printoptions, repr_string, resize_array, safe_int_conv, + signature_string_parts, writable_array) from odl.util.numerics import _SUPPORTED_RESIZE_PAD_MODES - __all__ = ('Resampling', 'ResizingOperator') @@ -61,16 +63,18 @@ def __init__(self, domain, range): Apply the corresponding resampling operator to an element: - >>> print(resampling([0, 1, 0])) - [ 0., 0., 1., 1., 0., 0.] + >>> resampling([0, 1, 0]) + uniform_discr(0.0, 1.0, 6).element([ 0., 0., 1., 1., 0., 0.]) The result depends on the interpolation chosen for the underlying spaces: >>> coarse_discr = odl.uniform_discr(0, 1, 3, interp='linear') >>> linear_resampling = odl.Resampling(coarse_discr, fine_discr) - >>> print(linear_resampling([0, 1, 0])) - [ 0. , 0.25, 0.75, 0.75, 0.25, 0. ] + >>> linear_resampling([0, 1, 0]) + uniform_discr(0.0, 1.0, 6).element( + [ 0. , 0.25, 0.75, 0.75, 0.25, 0. ] + ) """ if domain.fspace != range.fspace: raise ValueError('`domain.fspace` ({}) does not match ' @@ -98,6 +102,26 @@ def inverse(self): The returned operator is resampling defined in the opposite direction. + Examples + -------- + Create resampling operator and inverse: + + >>> coarse_discr = odl.uniform_discr(0, 1, 3) + >>> fine_discr = odl.uniform_discr(0, 1, 6) + >>> resampling = odl.Resampling(coarse_discr, fine_discr) + >>> resampling_inv = resampling.inverse + + The inverse is proper left inverse if the resampling goes from a + coarser to a finer sampling: + + >>> resampling_inv(resampling([0.0, 1.0, 0.0])) + uniform_discr(0.0, 1.0, 3).element([ 0., 1., 0.]) + + However, it can fail in the other direction: + + >>> resampling(resampling_inv([0.0, 0.0, 0.0, 1.0, 0.0, 0.0])) + uniform_discr(0.0, 1.0, 6).element([ 0., 0., 0., 0., 0., 0.]) + See Also -------- adjoint : resampling is unitary, so the adjoint is the inverse. @@ -115,38 +139,43 @@ def adjoint(self): ------- adjoint : Resampling Resampling operator defined in the opposite direction. + """ + return self.inverse + + def __repr__(self): + """Return ``repr(self)``. Examples -------- - Create resampling operator and inverse: - >>> coarse_discr = odl.uniform_discr(0, 1, 3) >>> fine_discr = odl.uniform_discr(0, 1, 6) >>> resampling = odl.Resampling(coarse_discr, fine_discr) - >>> resampling_inv = resampling.inverse - - The inverse is proper left inverse if the resampling goes from a - coarser to a finer sampling: - - >>> x = [0.0, 1.0, 0.0] - >>> print(resampling_inv(resampling(x))) - [ 0., 1., 0.] + >>> resampling + Resampling(uniform_discr(0.0, 1.0, 3), uniform_discr(0.0, 1.0, 6)) + """ + posargs = [self.domain, self.range] + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) - However, it can fail in the other direction: - >>> y = [0.0, 0.0, 0.0, 1.0, 0.0, 0.0] - >>> print(resampling(resampling_inv(y))) - [ 0., 0., 0., 0., 0., 0.] - """ - return self.inverse +class ResizingOperator(Operator): + """Operator mapping a discretized function to a new domain. -class ResizingOperatorBase(Operator): + This operator is a mapping between uniformly discretized + `DiscreteLp` spaces with the same `DiscreteLp.cell_sides`, + but different `DiscreteLp.shape`. The underlying operation is array + resizing, i.e. no resampling is performed. + In axes where the domain is enlarged, the new entries are filled + ("padded") according to a provided parameter ``pad_mode``. - """Base class for `ResizingOperator` and its adjoint. + All resizing operator variants are linear, except constant padding + with constant != 0. - This is an abstract class used to share code between the forward and - adjoint variants of the resizing operator. + See `the online documentation + `_ + on resizing operators for mathematical details. """ def __init__(self, domain, range=None, ran_shp=None, **kwargs): @@ -154,13 +183,13 @@ def __init__(self, domain, range=None, ran_shp=None, **kwargs): Parameters ---------- - domain : uniform `DiscreteLp` - Uniformly discretized space, the operator can be applied - to its elements. - range : uniform `DiscreteLp`, optional - Uniformly discretized space in which the result of the - application of this operator lies. - For the default ``None``, a space with the same attributes + domain : `DiscreteLp` + Space of discretized functions to which the operator can be + applied. It must be uniformly discretized in axes where + resizing is applied. + range : `DiscreteLp`, optional + Space in which the result of the application of this operator + lies. For the default ``None``, a space with the same attributes as ``domain`` is used, except for its shape, which is set to ``ran_shp``. ran_shp : sequence of ints, optional @@ -218,29 +247,35 @@ def __init__(self, domain, range=None, ran_shp=None, **kwargs): >>> x = [[1, 2, 3, 4], ... [5, 6, 7, 8]] >>> resize_op = odl.ResizingOperator(space, ran_shp=(4, 4)) - >>> print(resize_op(x)) - [[ 0., 0., 0., 0.], - [ 1., 2., 3., 4.], - [ 5., 6., 7., 8.], - [ 0., 0., 0., 0.]] + >>> resize_op(x) + uniform_discr([-0.5, 0. ], [ 1.5, 1. ], (4, 4)).element( + [[ 0., 0., 0., 0.], + [ 1., 2., 3., 4.], + [ 5., 6., 7., 8.], + [ 0., 0., 0., 0.]] + ) >>> >>> resize_op = odl.ResizingOperator(space, ran_shp=(4, 4), ... offset=(0, 0), ... pad_mode='periodic') - >>> print(resize_op(x)) - [[ 1., 2., 3., 4.], - [ 5., 6., 7., 8.], - [ 1., 2., 3., 4.], - [ 5., 6., 7., 8.]] + >>> resize_op(x) + uniform_discr([ 0., 0.], [ 2., 1.], (4, 4)).element( + [[ 1., 2., 3., 4.], + [ 5., 6., 7., 8.], + [ 1., 2., 3., 4.], + [ 5., 6., 7., 8.]] + ) >>> >>> resize_op = odl.ResizingOperator(space, ran_shp=(4, 4), ... offset=(0, 0), ... pad_mode='order0') - >>> print(resize_op(x)) - [[ 1., 2., 3., 4.], - [ 5., 6., 7., 8.], - [ 5., 6., 7., 8.], - [ 5., 6., 7., 8.]] + >>> resize_op(x) + uniform_discr([ 0., 0.], [ 2., 1.], (4, 4)).element( + [[ 1., 2., 3., 4.], + [ 5., 6., 7., 8.], + [ 5., 6., 7., 8.], + [ 5., 6., 7., 8.]] + ) Alternatively, the range of the operator can be provided directly. This requires that the partitions match, i.e. that the cell sizes @@ -250,11 +285,13 @@ def __init__(self, domain, range=None, ran_shp=None, **kwargs): >>> large_spc = odl.uniform_discr([-0.5, 0], [1.5, 1], (4, 4)) >>> resize_op = odl.ResizingOperator(space, large_spc, ... pad_mode='periodic') - >>> print(resize_op(x)) - [[ 5., 6., 7., 8.], - [ 1., 2., 3., 4.], - [ 5., 6., 7., 8.], - [ 1., 2., 3., 4.]] + >>> resize_op(x) + uniform_discr([-0.5, 0. ], [ 1.5, 1. ], (4, 4)).element( + [[ 5., 6., 7., 8.], + [ 1., 2., 3., 4.], + [ 5., 6., 7., 8.], + [ 1., 2., 3., 4.]] + ) """ # Swap names to be able to use the range iterator without worries import builtins @@ -313,8 +350,7 @@ def __init__(self, domain, range=None, ran_shp=None, **kwargs): # padding mode 'constant' with `pad_const != 0` is not linear linear = (self.pad_mode != 'constant' or self.pad_const == 0.0) - super(ResizingOperatorBase, self).__init__( - domain, ran, linear=linear) + super(ResizingOperator, self).__init__(domain, ran, linear=linear) @property def offset(self): @@ -337,26 +373,6 @@ def axes(self): return tuple(i for i in range(self.domain.ndim) if self.domain.shape[i] != self.range.shape[i]) - -class ResizingOperator(ResizingOperatorBase): - - """Operator mapping a discretized function to a new domain. - - This operator is a mapping between uniformly discretized - `DiscreteLp` spaces with the same `DiscreteLp.cell_sides`, - but different `DiscreteLp.shape`. The underlying operation is array - resizing, i.e. no resampling is performed. - In axes where the domain is enlarged, the new entries are filled - ("padded") according to a provided parameter ``pad_mode``. - - All resizing operator variants are linear, except constant padding - with constant != 0. - - See `the online documentation - `_ - on resizing operators for mathematical details. - """ - def _call(self, x, out): """Implement ``self(x, out)``.""" with writable_array(out) as out_arr: @@ -386,9 +402,9 @@ def adjoint(self): raise NotImplementedError('this operator is not linear and ' 'thus has no adjoint') - forward_op = self + op = self - class ResizingOperatorAdjoint(ResizingOperatorBase): + class ResizingOperatorAdjoint(Operator): """Adjoint of `ResizingOperator`. @@ -397,18 +413,23 @@ class ResizingOperatorAdjoint(ResizingOperatorBase): on resizing operators for mathematical details. """ + def __init__(self): + """Initialize a new instance.""" + super(ResizingOperatorAdjoint, self).__init__( + op.range, op.domain, linear=True) + def _call(self, x, out): """Implement ``self(x, out)``.""" with writable_array(out) as out_arr: - resize_array(x.asarray(), self.range.shape, - offset=self.offset, pad_mode=self.pad_mode, + resize_array(x.asarray(), op.domain.shape, + offset=op.offset, pad_mode=op.pad_mode, pad_const=0, direction='adjoint', out=out_arr) @property def adjoint(self): """Adjoint of the adjoint, i.e. the original operator.""" - return forward_op + return op @property def inverse(self): @@ -422,8 +443,11 @@ def inverse(self): domain=self.range, range=self.domain, pad_mode=self.pad_mode) - return ResizingOperatorAdjoint(domain=self.range, range=self.domain, - pad_mode=self.pad_mode) + def __repr__(self): + """Return ``repr(self)``.""" + return attribute_repr_string(repr(op), 'adjoint') + + return ResizingOperatorAdjoint() @property def inverse(self): @@ -437,6 +461,26 @@ def inverse(self): pad_mode=self.pad_mode, pad_const=self.pad_const) + def __repr__(self): + """Return ``repr(self)``. + + >>> space = odl.uniform_discr([0, 0], [1, 1], (2, 4)) + >>> resize_op = odl.ResizingOperator(space, ran_shp=(4, 4)) + >>> resize_op + ResizingOperator( + uniform_discr([ 0., 0.], [ 1., 1.], (2, 4)), + uniform_discr([-0.5, 0. ], [ 1.5, 1. ], (4, 4)) + ) + """ + posargs = [self.domain, self.range] + optargs = [('pad_mode', self.pad_mode, 'constant'), + ('pad_const', self.pad_const, 0)] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs, + mod=['!r', '']) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) + def _offset_from_spaces(dom, ran): """Return index offset corresponding to given spaces.""" diff --git a/odl/discr/discretization.py b/odl/discr/discretization.py index fb7915db7d8..731514b6bb6 100644 --- a/odl/discr/discretization.py +++ b/odl/discr/discretization.py @@ -8,16 +8,15 @@ """Base classes for discretization.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function from odl.operator import Operator +from odl.set import ComplexNumbers, RealNumbers from odl.set.sets import Set -from odl.space.base_tensors import TensorSpace, Tensor +from odl.space.base_tensors import Tensor, TensorSpace from odl.space.entry_points import tensor_space_impl -from odl.set import RealNumbers, ComplexNumbers from odl.util import ( - is_real_floating_dtype, is_complex_floating_dtype, is_numeric_dtype) - + is_complex_floating_dtype, is_numeric_dtype, is_real_floating_dtype) __all__ = ('DiscretizedSpace',) diff --git a/odl/discr/grid.py b/odl/discr/grid.py index d2aec484c19..bc735a609a6 100644 --- a/odl/discr/grid.py +++ b/odl/discr/grid.py @@ -12,14 +12,15 @@ space with a certain structure which is exploited to minimize storage. """ -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + import numpy as np -from odl.set import Set, IntervalProd +from odl.set import IntervalProd, Set from odl.util import ( - normalized_index_expression, normalized_scalar_param_list, safe_int_conv, - array_str, signature_string, indent, npy_printoptions) - + REPR_PRECISION, array_str, normalized_index_expression, + normalized_scalar_param_list, npy_printoptions, repr_string, safe_int_conv, + signature_string_parts) __all__ = ('RectGrid', 'uniform_grid', 'uniform_grid_fromintv') @@ -86,12 +87,9 @@ def __init__(self, *coord_vectors): Examples -------- - >>> g = RectGrid([1, 2, 5], [-2, 1.5, 2]) + >>> g = odl.RectGrid([1, 2, 5], [-2, 1.5, 2]) >>> g - RectGrid( - [ 1., 2., 5.], - [-2. , 1.5, 2. ] - ) + RectGrid([ 1., 2., 5.], [-2. , 1.5, 2. ]) >>> g.ndim # number of axes 2 >>> g.shape # points per axis @@ -102,26 +100,16 @@ def __init__(self, *coord_vectors): Grid points can be extracted with index notation (NOTE: This is slow, do not loop over the grid using indices!): - >>> g = RectGrid([-1, 0, 3], [2, 4, 5], [5], [2, 4, 7]) + >>> g = odl.RectGrid([-1, 0, 3], [2, 4, 5], [5], [2, 4, 7]) >>> g[0, 0, 0, 0] array([-1., 2., 5., 2.]) Slices and ellipsis are also supported: >>> g[:, 0, 0, 0] - RectGrid( - [-1., 0., 3.], - [ 2.], - [ 5.], - [ 2.] - ) + RectGrid([-1., 0., 3.], [ 2.], [ 5.], [ 2.]) >>> g[0, ..., 1:] - RectGrid( - [-1.], - [ 2., 4., 5.], - [ 5.], - [ 4., 7.] - ) + RectGrid([-1.], [ 2., 4., 5.], [ 5.], [ 4., 7.]) Notes ----- @@ -207,7 +195,7 @@ def coord_vectors(self): Examples -------- - >>> g = RectGrid([0, 1], [-1, 0, 2]) + >>> g = odl.RectGrid([0, 1], [-1, 0, 2]) >>> x, y = g.coord_vectors >>> x array([ 0., 1.]) @@ -253,7 +241,7 @@ def __len__(self): Examples -------- - >>> g = RectGrid([0, 1], [-1, 0, 2], [4, 5, 6]) + >>> g = odl.RectGrid([0, 1], [-1, 0, 2], [4, 5, 6]) >>> len(g) 2 @@ -269,7 +257,7 @@ def min_pt(self): Examples -------- - >>> g = RectGrid([1, 2, 5], [-2, 1.5, 2]) + >>> g = odl.RectGrid([1, 2, 5], [-2, 1.5, 2]) >>> g.min_pt array([ 1., -2.]) """ @@ -281,7 +269,7 @@ def max_pt(self): Examples -------- - >>> g = RectGrid([1, 2, 5], [-2, 1.5, 2]) + >>> g = odl.RectGrid([1, 2, 5], [-2, 1.5, 2]) >>> g.max_pt array([ 5., 2.]) """ @@ -325,7 +313,7 @@ def min(self, **kwargs): Examples -------- - >>> g = RectGrid([1, 2, 5], [-2, 1.5, 2]) + >>> g = odl.RectGrid([1, 2, 5], [-2, 1.5, 2]) >>> g.min() array([ 1., -2.]) @@ -356,7 +344,7 @@ def max(self, **kwargs): Examples -------- - >>> g = RectGrid([1, 2, 5], [-2, 1.5, 2]) + >>> g = odl.RectGrid([1, 2, 5], [-2, 1.5, 2]) >>> g.max() array([ 5., 2.]) @@ -406,13 +394,13 @@ def stride(self): NaN returned for non-uniform dimension: - >>> g = RectGrid([0, 1, 2], [0, 1, 4]) + >>> g = odl.RectGrid([0, 1, 2], [0, 1, 4]) >>> g.stride array([ 1., nan]) 0.0 returned for degenerate dimension: - >>> g = RectGrid([0, 1, 2], [0]) + >>> g = odl.RectGrid([0, 1, 2], [0]) >>> g.stride array([ 1., 0.]) """ @@ -436,7 +424,7 @@ def extent(self): Examples -------- - >>> g = RectGrid([1, 2, 5], [-2, 1.5, 2]) + >>> g = odl.RectGrid([1, 2, 5], [-2, 1.5, 2]) >>> g.extent array([ 4., 4.]) """ @@ -457,7 +445,7 @@ def convex_hull(self): Examples -------- - >>> g = RectGrid([-1, 0, 3], [2, 4], [5], [2, 4, 7]) + >>> g = odl.RectGrid([-1, 0, 3], [2, 4], [5], [2, 4, 7]) >>> g.convex_hull() IntervalProd([-1., 2., 5., 2.], [ 3., 4., 5., 7.]) """ @@ -487,8 +475,8 @@ def approx_equals(self, other, atol): Examples -------- - >>> g1 = RectGrid([0, 1], [-1, 0, 2]) - >>> g2 = RectGrid([-0.1, 1.1], [-1, 0.1, 2]) + >>> g1 = odl.RectGrid([0, 1], [-1, 0, 2]) + >>> g2 = odl.RectGrid([-0.1, 1.1], [-1, 0.1, 2]) >>> g1.approx_equals(g2, atol=0) False >>> g1.approx_equals(g2, atol=0.15) @@ -536,7 +524,7 @@ def approx_contains(self, other, atol): Examples -------- - >>> g = RectGrid([0, 1], [-1, 0, 2]) + >>> g = odl.RectGrid([0, 1], [-1, 0, 2]) >>> g.approx_contains([0, 0], atol=0.0) True >>> [0, 0] in g # equivalent @@ -662,15 +650,10 @@ def insert(self, index, *grids): Examples -------- - >>> g1 = RectGrid([0, 1], [-1, 0, 2]) - >>> g2 = RectGrid([1], [-6, 15]) + >>> g1 = odl.RectGrid([0, 1], [-1, 0, 2]) + >>> g2 = odl.RectGrid([1], [-6, 15]) >>> g1.insert(1, g2) - RectGrid( - [ 0., 1.], - [ 1.], - [ -6., 15.], - [-1., 0., 2.] - ) + RectGrid([ 0., 1.], [ 1.], [ -6., 15.], [-1., 0., 2.]) >>> g1.insert(1, g2, g2) RectGrid( [ 0., 1.], @@ -724,15 +707,10 @@ def append(self, *grids): Examples -------- - >>> g1 = RectGrid([0, 1], [-1, 0, 2]) - >>> g2 = RectGrid([1], [-6, 15]) + >>> g1 = odl.RectGrid([0, 1], [-1, 0, 2]) + >>> g2 = odl.RectGrid([1], [-6, 15]) >>> g1.append(g2) - RectGrid( - [ 0., 1.], - [-1., 0., 2.], - [ 1.], - [ -6., 15.] - ) + RectGrid([ 0., 1.], [-1., 0., 2.], [ 1.], [ -6., 15.]) >>> g1.append(g2, g2) RectGrid( [ 0., 1.], @@ -764,12 +742,9 @@ def squeeze(self, axis=None): Examples -------- - >>> g = RectGrid([0, 1], [-1], [-1, 0, 2]) + >>> g = odl.RectGrid([0, 1], [-1], [-1, 0, 2]) >>> g.squeeze() - RectGrid( - [ 0., 1.], - [-1., 0., 2.] - ) + RectGrid([ 0., 1.], [-1., 0., 2.]) """ if axis is None: rng = range(self.ndim) @@ -797,7 +772,7 @@ def points(self, order='C'): Examples -------- - >>> g = RectGrid([0, 1], [-1, 0, 2]) + >>> g = odl.RectGrid([0, 1], [-1, 0, 2]) >>> g.points() array([[ 0., -1.], [ 0., 0.], @@ -840,7 +815,7 @@ def corner_grid(self): Examples -------- - >>> g = RectGrid([0, 1], [-1, 0, 2]) + >>> g = odl.RectGrid([0, 1], [-1, 0, 2]) >>> g.corner_grid() uniform_grid([ 0., -1.], [ 1., 2.], (2, 2)) """ @@ -870,7 +845,7 @@ def corners(self, order='C'): Examples -------- - >>> g = RectGrid([0, 1], [-1, 0, 2]) + >>> g = odl.RectGrid([0, 1], [-1, 0, 2]) >>> g.corners() array([[ 0., -1.], [ 0., 2.], @@ -901,7 +876,7 @@ def meshgrid(self): Examples -------- - >>> g = RectGrid([0, 1], [-1, 0, 2]) + >>> g = odl.RectGrid([0, 1], [-1, 0, 2]) >>> x, y = g.meshgrid >>> x array([[ 0.], @@ -930,43 +905,23 @@ def __getitem__(self, indices): -------- Indexing with integers along all axes produces an array (a point): - >>> g = RectGrid([-1, 0, 3], [2, 4, 5], [5], [2, 4, 7]) + >>> g = odl.RectGrid([-1, 0, 3], [2, 4, 5], [5], [2, 4, 7]) >>> g[0, 0, 0, 0] array([-1., 2., 5., 2.]) Otherwise, a new RectGrid is returned: >>> g[:, 0, 0, 0] - RectGrid( - [-1., 0., 3.], - [ 2.], - [ 5.], - [ 2.] - ) + RectGrid([-1., 0., 3.], [ 2.], [ 5.], [ 2.]) >>> g[0, ..., 1:] - RectGrid( - [-1.], - [ 2., 4., 5.], - [ 5.], - [ 4., 7.] - ) + RectGrid([-1.], [ 2., 4., 5.], [ 5.], [ 4., 7.]) >>> g[::2, ..., ::2] - RectGrid( - [-1., 3.], - [ 2., 4., 5.], - [ 5.], - [ 2., 7.] - ) + RectGrid([-1., 3.], [ 2., 4., 5.], [ 5.], [ 2., 7.]) Too few indices are filled up with an ellipsis from the right: >>> g[0] - RectGrid( - [-1.], - [ 2., 4., 5.], - [ 5.], - [ 2., 4., 7.] - ) + RectGrid([-1.], [ 2., 4., 5.], [ 5.], [ 2., 4., 7.]) >>> g[0] == g[0, :, :, :] == g[0, ...] True """ @@ -1004,7 +959,7 @@ def __array__(self, dtype=None): Examples -------- - >>> g = RectGrid([0, 1], [-2, 0, 2]) + >>> g = odl.RectGrid([0, 1], [-2, 0, 2]) Convert to an array: @@ -1024,23 +979,31 @@ def __array__(self, dtype=None): return self.points().astype(dtype) def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> g = odl.RectGrid([-1, 0, 3], [2, 4, 5], [5], [2, 4, 7]) + >>> g + RectGrid([-1., 0., 3.], [ 2., 4., 5.], [ 5.], [ 2., 4., 7.]) + """ if self.is_uniform: ctor = 'uniform_grid' posargs = [self.min_pt, self.max_pt, self.shape] posmod = [array_str, array_str, ''] - with npy_printoptions(precision=4): - inner_str = signature_string(posargs, [], mod=[posmod, '']) - return '{}({})'.format(ctor, inner_str) else: ctor = self.__class__.__name__ posargs = self.coord_vectors posmod = array_str - inner_str = signature_string(posargs, [], sep=[',\n', ', ', ', '], - mod=[posmod, '']) - return '{}(\n{}\n)'.format(ctor, indent(inner_str)) - __str__ = __repr__ + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, [], + mod=[posmod, '']) + return repr_string(ctor, inner_parts) + + def __str__(self): + """Return ``str(self)``.""" + return repr(self) def uniform_grid_fromintv(intv_prod, shape, nodes_on_bdry=True): diff --git a/odl/discr/lp_discr.py b/odl/discr/lp_discr.py index f475fa88953..1382b62dd22 100644 --- a/odl/discr/lp_discr.py +++ b/odl/discr/lp_discr.py @@ -8,26 +8,29 @@ """Lebesgue L^p type discretizations of function spaces.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + from numbers import Integral + import numpy as np +from odl.discr.discr_mappings import ( + LinearInterpolation, NearestInterpolation, PerAxisInterpolation, + PointCollocation) from odl.discr.discretization import ( DiscretizedSpace, DiscretizedSpaceElement, tspace_type) -from odl.discr.discr_mappings import ( - PointCollocation, NearestInterpolation, LinearInterpolation, - PerAxisInterpolation) from odl.discr.partition import ( - RectPartition, uniform_partition_fromintv, uniform_partition) -from odl.set import RealNumbers, ComplexNumbers, IntervalProd + RectPartition, uniform_partition, uniform_partition_fromintv) +from odl.set import ComplexNumbers, IntervalProd, RealNumbers from odl.space import FunctionSpace, ProductSpace from odl.space.entry_points import tensor_space_impl from odl.space.weighting import ConstWeighting from odl.util import ( - apply_on_boundary, is_real_dtype, is_complex_floating_dtype, is_string, - is_floating_dtype, is_numeric_dtype, - dtype_str, array_str, signature_string, indent, npy_printoptions, - normalized_scalar_param_list, safe_int_conv, normalized_nodes_on_bdry) + REPR_PRECISION, apply_on_boundary, array_str, attribute_repr_string, + dtype_str, is_complex_floating_dtype, is_floating_dtype, is_numeric_dtype, + is_real_dtype, is_string, normalized_nodes_on_bdry, + normalized_scalar_param_list, npy_printoptions, repr_string, safe_int_conv, + signature_string_parts) __all__ = ('DiscreteLp', 'DiscreteLpElement', 'uniform_discr_frompartition', 'uniform_discr_fromspace', @@ -378,6 +381,13 @@ def _astype(self, dtype): fspace, self.partition, tspace, interp=self.interp, axis_labels=self.axis_labels) + def _asweighted(self, weighting): + """Internal helper for ``asweighted``.""" + tspace = self.tspace.asweighted(weighting) + return type(self)( + self.fspace, self.partition, tspace, interp=self.interp, + axis_labels=self.axis_labels) + # Overrides for space functions depending on partition # # The inherited methods by default use a weighting by a constant @@ -502,12 +512,46 @@ def __getitem__(self, indices): def __repr__(self): """Return ``repr(self)``.""" - return repr(space) + '.byaxis_in' + return attribute_repr_string(repr(space), 'byaxis_in') return DiscreteLpByaxisIn() def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + Uniform: + + >>> space = odl.uniform_discr(0, 1, 4, dtype='float32') + >>> space + uniform_discr(0.0, 1.0, 4, dtype='float32') + >>> space = odl.uniform_discr([0, 0, 0], [1, 1, 1], (2, 4, 8), + ... dtype=complex) + >>> space + uniform_discr( + [ 0., 0., 0.], [ 1., 1., 1.], (2, 4, 8), dtype=complex + ) + + Non-uniform: + + >>> rect = odl.IntervalProd([0, 0], [1, 1]) + >>> fspace = odl.FunctionSpace(rect) + >>> grid = odl.RectGrid([0.1, 0.25, 0.9], [0.25, 0.75]) + >>> part = odl.RectPartition(rect, grid) + >>> tspace = odl.rn(part.shape) + >>> space = odl.DiscreteLp(fspace, part, tspace, interp='linear') + >>> space + DiscreteLp( + FunctionSpace(IntervalProd([ 0., 0.], [ 1., 1.])), + nonuniform_partition( + [ 0.1 , 0.25, 0.9 ], [ 0.25, 0.75], + min_pt=[ 0., 0.], max_pt=[ 1., 1.] + ), + rn((3, 2)), + interp='linear' + ) + """ # Clunky check if the factory repr can be used if (uniform_partition_fromintv( self.fspace.domain, self.shape, @@ -526,7 +570,7 @@ def __repr__(self): ctor = 'uniform_discr' if self.ndim == 1: posargs = [self.min_pt[0], self.max_pt[0], self.shape[0]] - posmod = [':.4', ':.4', ''] + posmod = '' else: posargs = [self.min_pt, self.max_pt, self.shape] posmod = [array_str, array_str, ''] @@ -558,22 +602,19 @@ def __repr__(self): if self.dtype in (float, complex, int, bool): optmod[3] = '!s' - with npy_printoptions(precision=4): - inner_str = signature_string(posargs, optargs, - mod=[posmod, optmod]) - return '{}({})'.format(ctor, inner_str) + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs, + mod=[posmod, optmod]) else: ctor = self.__class__.__name__ posargs = [self.fspace, self.partition, self.tspace] optargs = [('interp', self.interp, 'nearest')] - with npy_printoptions(precision=4): - inner_str = signature_string(posargs, optargs, - sep=[',\n', ', ', ',\n'], - mod=['!r', '!s']) + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs) - return '{}({})'.format(ctor, indent(inner_str)) + return repr_string(ctor, inner_parts, allow_mixed_seps=True) def __str__(self): """Return ``str(self)``.""" @@ -641,7 +682,7 @@ def conj(self, out=None): Examples -------- - >>> discr = uniform_discr(0, 1, 4, dtype=complex) + >>> discr = odl.uniform_discr(0, 1, 4, dtype=complex) >>> x = discr.element([5+1j, 3, 2-2j, 1j]) >>> y = x.conj() >>> print(y) @@ -1303,7 +1344,7 @@ def uniform_discr_frompartition(partition, dtype=None, impl='numpy', **kwargs): Examples -------- >>> part = odl.uniform_partition(0, 1, 10) - >>> uniform_discr_frompartition(part) + >>> odl.uniform_discr_frompartition(part) uniform_discr(0.0, 1.0, 10) See Also @@ -1369,7 +1410,7 @@ def uniform_discr_fromspace(fspace, shape, dtype=None, impl='numpy', **kwargs): -------- >>> intv = odl.IntervalProd(0, 1) >>> space = odl.FunctionSpace(intv) - >>> uniform_discr_fromspace(space, 10) + >>> odl.uniform_discr_fromspace(space, 10) uniform_discr(0.0, 1.0, 10) See Also @@ -1444,7 +1485,7 @@ def uniform_discr_fromintv(intv_prod, shape, dtype=None, impl='numpy', Examples -------- >>> intv = IntervalProd(0, 1) - >>> uniform_discr_fromintv(intv, 10) + >>> odl.uniform_discr_fromintv(intv, 10) uniform_discr(0.0, 1.0, 10) See Also @@ -1518,7 +1559,7 @@ def uniform_discr(min_pt, max_pt, shape, dtype=None, impl='numpy', **kwargs): -------- Create real space: - >>> space = uniform_discr([0, 0], [1, 1], (10, 10)) + >>> space = odl.uniform_discr([0, 0], [1, 1], (10, 10)) >>> space uniform_discr([ 0., 0.], [ 1., 1.], (10, 10)) >>> space.cell_sides @@ -1530,7 +1571,7 @@ def uniform_discr(min_pt, max_pt, shape, dtype=None, impl='numpy', **kwargs): Create complex space by giving a dtype: - >>> space = uniform_discr([0, 0], [1, 1], (10, 10), dtype=complex) + >>> space = odl.uniform_discr([0, 0], [1, 1], (10, 10), dtype=complex) >>> space uniform_discr([ 0., 0.], [ 1., 1.], (10, 10), dtype=complex) >>> space.is_complex @@ -1789,9 +1830,11 @@ def uniform_discr_fromdiscr(discr, min_pt=None, max_pt=None, cell_sides=new_csides, nodes_on_bdry=nodes_on_bdry) + dtype = kwargs.pop('dtype', discr.dtype) + # TODO: weighting return uniform_discr_frompartition(new_part, exponent=discr.exponent, interp=discr.interp, impl=discr.impl, - **kwargs) + dtype=dtype, **kwargs) def _scaling_func_list(bdry_fracs, exponent): diff --git a/odl/discr/partition.py b/odl/discr/partition.py index 15c6070ed4d..be3f46c8444 100644 --- a/odl/discr/partition.py +++ b/odl/discr/partition.py @@ -14,17 +14,19 @@ of partitions of intervals. """ -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + from builtins import object + import numpy as np from odl.discr.grid import RectGrid, uniform_grid_fromintv from odl.set import IntervalProd from odl.util import ( + REPR_PRECISION, array_str, attribute_repr_string, normalized_index_expression, normalized_nodes_on_bdry, - normalized_scalar_param_list, safe_int_conv, - signature_string, indent, array_str, npy_printoptions) - + normalized_scalar_param_list, npy_printoptions, repr_string, safe_int_conv, + signature_string_parts) __all__ = ('RectPartition', 'uniform_partition_fromintv', 'uniform_partition_fromgrid', 'uniform_partition', @@ -150,11 +152,11 @@ def nodes_on_bdry(self): nodes_on_bdry = [] for on_bdry in self.nodes_on_bdry_byaxis: - l, r = on_bdry - if l == r: - nodes_on_bdry.append(l) + left, right = on_bdry + if left == right: + nodes_on_bdry.append(left) else: - nodes_on_bdry.append((l, r)) + nodes_on_bdry.append((left, right)) if all(on_bdry == nodes_on_bdry[0] for on_bdry in nodes_on_bdry[1:]): return nodes_on_bdry[0] else: @@ -499,11 +501,8 @@ def __getitem__(self, indices): Take every second grid point. Note that is is in general non-uniform: >>> partition = odl.uniform_partition(0, 10, 10) - >>> partition[::2] - nonuniform_partition( - [ 0.5, 2.5, 4.5, 6.5, 8.5], - min_pt=0.0, max_pt=10.0 - ) + >>> partition[2::2] + nonuniform_partition([ 2.5, 4.5, 6.5, 8.5], min_pt=2.0, max_pt=10.0) A more advanced example is: @@ -847,12 +846,27 @@ def __repr__(self): >>> p.byaxis uniform_partition(0, 1, 5).byaxis """ - return '{!r}.byaxis'.format(partition) + return attribute_repr_string(repr(partition), 'byaxis') return RectPartitionByAxis() def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> p = odl.uniform_partition([0, 1, 2], [1, 3, 5], (3, 5, 6)) + >>> p + uniform_partition([ 0., 1., 2.], [ 1., 3., 5.], (3, 5, 6)) + + >>> p = odl.nonuniform_partition([0, 1, 2], [1, 1.5, 2], + ... min_pt=[0, 0], max_pt=[2.5, 3]) + >>> p + nonuniform_partition( + [ 0., 1., 2.], [ 1. , 1.5, 2. ], + min_pt=[ 0., 0.], max_pt=[ 2.5, 3. ] + ) + """ if self.ndim == 0: return 'uniform_partition([], [], ())' @@ -876,16 +890,16 @@ def __repr__(self): ctor = 'uniform_partition' if self.ndim == 1: posargs = [self.min_pt[0], self.max_pt[0], self.shape[0]] - posmod = [':.4', ':.4', ''] + posmod = '' else: posargs = [self.min_pt, self.max_pt, self.shape] posmod = [array_str, array_str, ''] optargs = [('nodes_on_bdry', self.nodes_on_bdry, False)] - with npy_printoptions(precision=4): - sig_str = signature_string(posargs, optargs, mod=[posmod, '']) - return '{}({})'.format(ctor, sig_str) + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs, + mod=[posmod, '']) else: ctor = 'nonuniform_partition' posargs = self.coord_vectors @@ -913,9 +927,9 @@ def __repr__(self): if not np.allclose(self.min_pt, def_min_pt): if self.ndim == 1: optargs.append(('min_pt', self.min_pt[0], None)) - optmod.append(':.4') + optmod.append('') else: - with npy_printoptions(precision=4): + with npy_printoptions(precision=REPR_PRECISION): optargs.append( ('min_pt', array_str(self.min_pt), '')) optmod.append('!s') @@ -923,16 +937,18 @@ def __repr__(self): if not np.allclose(self.max_pt, def_max_pt): if self.ndim == 1: optargs.append(('max_pt', self.max_pt[0], None)) - optmod.append(':.4') + optmod.append('') else: - with npy_printoptions(precision=4): + with npy_printoptions(precision=REPR_PRECISION): optargs.append( ('max_pt', array_str(self.max_pt), '')) optmod.append('!s') - sig_str = signature_string(posargs, optargs, mod=[posmod, optmod], - sep=[',\n', ', ', ',\n']) - return '{}(\n{}\n)'.format(ctor, indent(sig_str)) + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs, + mod=[posmod, optmod]) + + return repr_string(ctor, inner_parts, allow_mixed_seps=True) def __str__(self): """Return ``str(self)``.""" @@ -1338,18 +1354,13 @@ def nonuniform_partition(*coord_vecs, **kwargs): that the points are in the middle of the sub-intervals: >>> odl.nonuniform_partition([0, 1, 3]) - nonuniform_partition( - [ 0., 1., 3.] - ) + nonuniform_partition([ 0., 1., 3.]) Higher dimensional partitions are created by specifying the gridpoints along each dimension: >>> odl.nonuniform_partition([0, 1, 3], [1, 2]) - nonuniform_partition( - [ 0., 1., 3.], - [ 1., 2.] - ) + nonuniform_partition([ 0., 1., 3.], [ 1., 2.]) Partitions with a single element are by default degenerate @@ -1360,19 +1371,13 @@ def nonuniform_partition(*coord_vecs, **kwargs): can be used: >>> odl.nonuniform_partition([0, 1, 3], nodes_on_bdry=True) - nonuniform_partition( - [ 0., 1., 3.], - nodes_on_bdry=True - ) + nonuniform_partition([ 0., 1., 3.], nodes_on_bdry=True) Users can also manually specify the containing intervals dimensions by using the ``min_pt`` and ``max_pt`` arguments: >>> odl.nonuniform_partition([0, 1, 3], min_pt=-2, max_pt=3) - nonuniform_partition( - [ 0., 1., 3.], - min_pt=-2.0, max_pt=3.0 - ) + nonuniform_partition([ 0., 1., 3.], min_pt=-2.0, max_pt=3.0) """ # Get parameters from kwargs min_pt = kwargs.pop('min_pt', None) diff --git a/odl/operator/default_ops.py b/odl/operator/default_ops.py index dd6107e620e..920bee6217d 100644 --- a/odl/operator/default_ops.py +++ b/odl/operator/default_ops.py @@ -8,15 +8,17 @@ """Default operators defined on any (reasonable) space.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + from copy import copy + import numpy as np from odl.operator.operator import Operator -from odl.set import LinearSpace, Field, RealNumbers -from odl.set.space import LinearSpaceElement +from odl.set import Field, LinearSpace, RealNumbers from odl.space import ProductSpace - +from odl.util import ( + REPR_PRECISION, npy_printoptions, repr_string, signature_string_parts) __all__ = ('ScalingOperator', 'ZeroOperator', 'IdentityOperator', 'LinCombOperator', 'MultiplyOperator', 'PowerOperator', @@ -34,7 +36,7 @@ class ScalingOperator(Operator): ScalingOperator(s)(x) == s * x """ - def __init__(self, domain, scalar): + def __init__(self, domain, scalar, range=None): """Initialize a new instance. Parameters @@ -43,13 +45,16 @@ def __init__(self, domain, scalar): Set of elements on which this operator acts. scalar : ``domain.field`` element Fixed scaling factor of this operator. + range : `LinearSpace` or `Field`, optional + Set to which this operator maps. + Default: ``domain``. Examples -------- >>> r3 = odl.rn(3) >>> vec = r3.element([1, 2, 3]) >>> out = r3.element() - >>> op = ScalingOperator(r3, 2.0) + >>> op = odl.ScalingOperator(r3, 2.0) >>> op(vec, out) # In-place, Returns out rn(3).element([ 2., 4., 6.]) >>> out @@ -58,10 +63,17 @@ def __init__(self, domain, scalar): rn(3).element([ 2., 4., 6.]) """ if not isinstance(domain, (LinearSpace, Field)): - raise TypeError('`space` {!r} not a `LinearSpace` or `Field` ' + raise TypeError('`domain` {!r} not a `LinearSpace` or `Field` ' 'instance'.format(domain)) - super(ScalingOperator, self).__init__(domain, domain, linear=True) + if range is None: + range = domain + else: + if not isinstance(range, (LinearSpace, Field)): + raise TypeError('`range` {!r} not a `LinearSpace` or `Field` ' + 'instance'.format(range)) + + super(ScalingOperator, self).__init__(domain, range, linear=True) self.__scalar = domain.field.element(scalar) @property @@ -84,18 +96,18 @@ def inverse(self): Examples -------- >>> r3 = odl.rn(3) - >>> vec = r3.element([1, 2, 3]) - >>> op = ScalingOperator(r3, 2.0) + >>> x = r3.element([1, 2, 3]) + >>> op = odl.ScalingOperator(r3, 2.0) >>> inv = op.inverse - >>> inv(op(vec)) == vec + >>> inv(op(x)) == x True - >>> op(inv(vec)) == vec + >>> op(inv(x)) == x True """ if self.scalar == 0.0: raise ZeroDivisionError('scaling operator not invertible for ' 'scalar==0') - return ScalingOperator(self.domain, 1.0 / self.scalar) + return ScalingOperator(self.range, 1.0 / self.scalar, self.domain) @property def adjoint(self): @@ -107,9 +119,10 @@ def adjoint(self): ``self`` if `scalar` is real, else `scalar` is conjugated. """ if complex(self.scalar).imag == 0.0: - return self + return ScalingOperator(self.range, self.scalar, self.domain) else: - return ScalingOperator(self.domain, self.scalar.conjugate()) + return ScalingOperator(self.range, self.scalar.conjugate(), + self.domain) def norm(self, estimate=False, **kwargs): """Return the operator norm of this operator. @@ -128,19 +141,27 @@ def norm(self, estimate=False, **kwargs): -------- >>> spc = odl.rn(3) >>> scaling = odl.ScalingOperator(spc, 3.0) - >>> scaling.norm(True) + >>> scaling.norm(estimate=True) 3.0 """ return np.abs(self.scalar) def __repr__(self): - """Return ``repr(self)``.""" - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.domain, self.scalar) + """Return ``repr(self)``. - def __str__(self): - """Return ``str(self)``.""" - return '{} * I'.format(self.scalar) + Examples + -------- + >>> r3 = odl.rn(3) + >>> op = odl.ScalingOperator(r3, 2.0) + >>> op + ScalingOperator(rn(3), scalar=2.0) + """ + posargs = [self.domain] + optargs = [('scalar', self.scalar, None), + ('range', self.range, self.domain)] + with npy_printoptions(precision=4): + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class IdentityOperator(ScalingOperator): @@ -152,23 +173,38 @@ class IdentityOperator(ScalingOperator): IdentityOperator()(x) == x """ - def __init__(self, space): + def __init__(self, domain, range=None): """Initialize a new instance. Parameters ---------- - space : `LinearSpace` - Space of elements which the operator is acting on. + domain : `LinearSpace` or `Field` + Set of elements on which this operator acts. + range : `LinearSpace` or `Field`, optional + Set to which this operator maps. + Default: ``domain``. """ - super(IdentityOperator, self).__init__(space, 1) + super(IdentityOperator, self).__init__(domain, 1, range) + + @property + def adjoint(self): + """Adjoint of the identity operator.""" + return IdentityOperator(self.range, self.domain) def __repr__(self): - """Return ``repr(self)``.""" - return '{}({!r})'.format(self.__class__.__name__, self.domain) + """Return ``repr(self)``. - def __str__(self): - """Return ``str(self)``.""" - return "I" + Examples + -------- + >>> r3 = odl.rn(3) + >>> op = odl.IdentityOperator(r3) + >>> op + IdentityOperator(rn(3)) + """ + posargs = [self.domain] + optargs = [('range', self.range, self.domain)] + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class LinCombOperator(Operator): @@ -196,7 +232,7 @@ def __init__(self, space, a, b): >>> r3xr3 = odl.ProductSpace(r3, r3) >>> xy = r3xr3.element([[1, 2, 3], [1, 2, 3]]) >>> z = r3.element() - >>> op = LinCombOperator(r3, 1.0, 1.0) + >>> op = odl.LinCombOperator(r3, 1.0, 1.0) >>> op(xy, out=z) # Returns z rn(3).element([ 2., 4., 6.]) >>> z @@ -244,48 +280,53 @@ def __init__(self, multiplicand, domain=None, range=None): Parameters ---------- - multiplicand : `LinearSpaceElement` or scalar - Value to multiply by. + multiplicand : `LinearSpaceElement` or ``domain`` `element-like` + Vector or scalar with which should be multiplied. If ``domain`` + is provided, this parameter can be an `element-like` object + for ``domain``. Otherwise it must be a `LinearSpaceElement`. domain : `LinearSpace` or `Field`, optional - Set to which the operator can be applied. + Set to which the operator can be applied. Mandatory if + ``multiplicand`` is not a `LinearSpaceElement`. Default: ``multiplicand.space``. range : `LinearSpace` or `Field`, optional - Set to which the operator maps. Default: ``multiplicand.space``. + Set to which the operator maps. + Default: ``domain`` if given, otherwise ``multiplicand.space``. Examples -------- + If a `LinearSpaceElement` is provided, domain and range are + inferred: + >>> r3 = odl.rn(3) >>> x = r3.element([1, 2, 3]) - - Multiply by vector: - - >>> op = MultiplyOperator(x) - >>> op(x) - rn(3).element([ 1., 4., 9.]) + >>> op = odl.MultiplyOperator(x) + >>> op([2, 4, 6]) + rn(3).element([ 2., 8., 18.]) >>> out = r3.element() >>> op(x, out) rn(3).element([ 1., 4., 9.]) - Multiply by scalar: + For a scalar or `element-like` multiplicand, ``domain`` (and + ``range``) should be given: - >>> op2 = MultiplyOperator(x, domain=r3.field) - >>> op2(3) + >>> op = odl.MultiplyOperator(x, domain=r3.field, range=r3) + >>> op(3) rn(3).element([ 3., 6., 9.]) >>> out = r3.element() - >>> op2(3, out) + >>> op(3, out) rn(3).element([ 3., 6., 9.]) """ if domain is None: domain = multiplicand.space if range is None: - range = multiplicand.space + range = domain super(MultiplyOperator, self).__init__(domain, range, linear=True) self.__multiplicand = multiplicand - self.__domain_is_field = isinstance(domain, Field) - self.__range_is_field = isinstance(range, Field) + self.__domain_is_field = isinstance(self.domain, Field) + self.__range_is_field = isinstance(self.range, Field) @property def multiplicand(self): @@ -317,42 +358,51 @@ def adjoint(self): Examples -------- - >>> r3 = odl.rn(3) - >>> x = r3.element([1, 2, 3]) - Multiply by a space element: - >>> op = MultiplyOperator(x) - >>> out = r3.element() + >>> r3 = odl.rn(3) + >>> x = r3.element([1, 2, 3]) + >>> op = odl.MultiplyOperator(x) >>> op.adjoint(x) rn(3).element([ 1., 4., 9.]) Multiply scalars with a fixed vector: - >>> op2 = MultiplyOperator(x, domain=r3.field) - >>> op2.adjoint(x) + >>> op = odl.MultiplyOperator(x, domain=r3.field, range=r3) + >>> op.adjoint(x) 14.0 Multiply vectors with a fixed scalar: - >>> op2 = MultiplyOperator(3.0, domain=r3, range=r3) - >>> op2.adjoint(x) + >>> op = odl.MultiplyOperator(3.0, domain=r3, range=r3) + >>> op.adjoint(x) rn(3).element([ 3., 6., 9.]) """ if self.__domain_is_field: - return InnerProductOperator(self.multiplicand) + return InnerProductOperator(self.multiplicand, self.range, + self.domain) else: # TODO: complex case - return MultiplyOperator(self.multiplicand, - domain=self.range, range=self.domain) + return MultiplyOperator(self.multiplicand, self.range, self.domain) def __repr__(self): - """Return ``repr(self)``.""" - return '{}({!r})'.format(self.__class__.__name__, self.multiplicand) + """Return ``repr(self)``. - def __str__(self): - """Return ``str(self)``.""" - return "x * {}".format(self.y) + Examples + -------- + >>> r3 = odl.rn(3) + >>> x = r3.element([1, 2, 3]) + >>> op = odl.MultiplyOperator(x) + >>> op + MultiplyOperator(rn(3).element([ 1., 2., 3.])) + """ + posargs = [self.multiplicand] + optargs = [('domain', self.domain, + getattr(self.multiplicand, 'space', None)), + ('range', self.range, self.domain)] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class PowerOperator(Operator): @@ -368,7 +418,7 @@ class PowerOperator(Operator): `LinearSpace` or on a `Field`. """ - def __init__(self, domain, exponent): + def __init__(self, domain, exponent, range=None): """Initialize a new instance. Parameters @@ -377,25 +427,30 @@ def __init__(self, domain, exponent): Set of elements on which the operator can be applied. exponent : float Exponent parameter of the power function applied to an element. + range : `LinearSpace` or `Field`, optional + Set to which this operator maps. + Default: ``domain``. Examples -------- - Use with vectors + Usage on a space of vectors: - >>> op = PowerOperator(odl.rn(3), exponent=2) + >>> op = odl.PowerOperator(odl.rn(3), exponent=2) >>> op([1, 2, 3]) rn(3).element([ 1., 4., 9.]) - or scalars + For a scalar space: - >>> op = PowerOperator(odl.RealNumbers(), exponent=2) + >>> op = odl.PowerOperator(odl.RealNumbers(), exponent=2) >>> op(2.0) 4.0 """ + if range is None: + range = domain super(PowerOperator, self).__init__( - domain, domain, linear=(exponent == 1)) + domain, range, linear=(exponent == 1)) self.__exponent = float(exponent) - self.__domain_is_field = isinstance(domain, Field) + self.__range_is_field = isinstance(range, Field) @property def exponent(self): @@ -406,7 +461,7 @@ def _call(self, x, out=None): """Take the power of ``x`` and write to ``out`` if given.""" if out is None: return x ** self.exponent - elif self.__domain_is_field: + elif self.__range_is_field: raise ValueError('cannot use `out` with field') else: out.assign(x) @@ -431,30 +486,38 @@ def derivative(self, point): -------- Use on vector spaces: - >>> op = PowerOperator(odl.rn(3), exponent=2) + >>> op = odl.PowerOperator(odl.rn(3), exponent=2) >>> dop = op.derivative(op.domain.element([1, 2, 3])) >>> dop([1, 1, 1]) rn(3).element([ 2., 4., 6.]) Use with scalars: - >>> op = PowerOperator(odl.RealNumbers(), exponent=2) + >>> op = odl.PowerOperator(odl.RealNumbers(), exponent=2) >>> dop = op.derivative(2.0) >>> dop(2.0) 8.0 """ - return self.exponent * MultiplyOperator(point ** (self.exponent - 1), - domain=self.domain, - range=self.range) + return (self.exponent * + MultiplyOperator(point ** (self.exponent - 1), + domain=self.domain, range=self.range) + ) def __repr__(self): - """Return ``repr(self)``.""" - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.domain, self.exponent) + """Return ``repr(self)``. - def __str__(self): - """Return ``str(self)``.""" - return "x ** {}".format(self.exponent) + Examples + -------- + >>> op = odl.PowerOperator(odl.rn(3), exponent=2) + >>> op + PowerOperator(rn(3), exponent=2.0) + """ + posargs = [self.domain] + optargs = [('exponent', self.exponent, None), + ('range', self.range, self.domain)] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class InnerProductOperator(Operator): @@ -472,25 +535,62 @@ class InnerProductOperator(Operator): NormOperator : Vector space norm as operator. """ - def __init__(self, vector): + def __init__(self, vector, domain=None, range=None): """Initialize a new instance. Parameters ---------- - vector : `LinearSpaceElement` - The element to take the inner product with. + vector : `LinearSpaceElement` or ``domain`` `element-like` + The element with which the inner product is taken. If ``domain`` + is given, this can be an `element-like` object for ``domain``, + otherwise it must be a `LinearSpaceElement`. + domain : `LinearSpace` or `Field`, optional + Set of elements on which the operator can be applied. Optional + if ``vector`` is a `LinearSpaceElement`, in which case + ``vector.space`` is taken as default. Otherwise this parameter + is mandatory. + range : `LinearSpace` or `Field`, optional + Set to which this operator maps. + Default: ``domain.field``. Examples -------- + By default, ``domain`` and ``range`` are inferred from the + given ``vector``: + >>> r3 = odl.rn(3) >>> x = r3.element([1, 2, 3]) - >>> op = InnerProductOperator(x) - >>> op(r3.element([1, 2, 3])) + >>> op = odl.InnerProductOperator(x) + >>> op.range + RealNumbers() + >>> op([1, 2, 3]) 14.0 + + With an explicit domain, we do not need a `LinearSpaceElement` + as ``vector``: + + >>> op = odl.InnerProductOperator([1, 2, 3], domain=r3) + >>> op([1, 2, 3]) + 14.0 + + We can also specify an explicit range, which should be able to hold + a single scalar: + + >>> r1 = odl.rn(1) + >>> op = odl.InnerProductOperator([1, 2, 3], domain=r3, range=r1) + >>> op([1, 2, 3]) + rn(1).element([ 14.]) + >>> r = odl.rn(()) + >>> op = odl.InnerProductOperator([1, 2, 3], domain=r3, range=r) + >>> op([1, 2, 3]) + rn(()).element(14.0) """ - super(InnerProductOperator, self).__init__( - vector.space, vector.space.field, linear=True) - self.__vector = vector + if domain is None: + domain = vector.space + if range is None: + range = domain.field + super(InnerProductOperator, self).__init__(domain, range, linear=True) + self.__vector = self.domain.element(vector) @property def vector(self): @@ -514,11 +614,12 @@ def adjoint(self): -------- >>> r3 = odl.rn(3) >>> x = r3.element([1, 2, 3]) - >>> op = InnerProductOperator(x) + >>> op = odl.InnerProductOperator(x) >>> op.adjoint(2.0) rn(3).element([ 2., 4., 6.]) """ - return MultiplyOperator(self.vector, self.vector.space.field) + return MultiplyOperator(self.vector, domain=self.range, + range=self.domain) @property def T(self): @@ -541,12 +642,23 @@ def T(self): return self.vector def __repr__(self): - """Return ``repr(self)``.""" - return '{}({!r})'.format(self.__class__.__name__, self.vector) + """Return ``repr(self)``. - def __str__(self): - """Return ``str(self)``.""" - return '{}.T'.format(self.vector) + Examples + -------- + >>> r3 = odl.rn(3) + >>> x = r3.element([1, 2, 3]) + >>> op = odl.InnerProductOperator(x) + >>> op + InnerProductOperator(rn(3).element([ 1., 2., 3.])) + """ + posargs = [self.vector] + optargs = [('domain', self.domain, + getattr(self.vector, 'space', None)), + ('range', self.range, self.domain.field)] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class NormOperator(Operator): @@ -566,22 +678,28 @@ class NormOperator(Operator): DistOperator : Distance to a fixed space element. """ - def __init__(self, space): + def __init__(self, domain, range=None): """Initialize a new instance. Parameters ---------- - space : `LinearSpace` - Space to take the norm in. + domain : `LinearSpace` + Set of elements on which the operator can be applied. Needs + to implement ``space.norm``. + range : `LinearSpace` or `Field`, optional + Set to which this operator maps. + Default: ``RealNumbers``. Examples -------- >>> r2 = odl.rn(2) - >>> op = NormOperator(r2) + >>> op = odl.NormOperator(r2) >>> op([3, 4]) 5.0 """ - super(NormOperator, self).__init__(space, RealNumbers(), linear=False) + if range is None: + range = RealNumbers() + super(NormOperator, self).__init__(domain, range, linear=False) def _call(self, x): """Return the norm of ``x``.""" @@ -620,7 +738,7 @@ def derivative(self, point): Examples -------- >>> r3 = odl.rn(3) - >>> op = NormOperator(r3) + >>> op = odl.NormOperator(r3) >>> derivative = op.derivative([1, 0, 0]) >>> derivative([1, 0, 0]) 1.0 @@ -630,15 +748,22 @@ def derivative(self, point): if norm == 0: raise ValueError('not differentiable in 0') - return InnerProductOperator(point / norm) + return InnerProductOperator(point / norm, self.domain, self.range) def __repr__(self): - """Return ``repr(self)``.""" - return '{}({!r})'.format(self.__class__.__name__, self.domain) + """Return ``repr(self)``. - def __str__(self): - """Return ``str(self)``.""" - return '{}({})'.format(self.__class__.__name__, self.domain) + Examples + -------- + >>> r2 = odl.rn(2) + >>> op = odl.NormOperator(r2) + >>> op + NormOperator(rn(2)) + """ + posargs = [self.domain] + optargs = [('range', self.range, RealNumbers())] + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class DistOperator(Operator): @@ -658,25 +783,38 @@ class DistOperator(Operator): NormOperator : Vector space norm as an operator. """ - def __init__(self, vector): + def __init__(self, vector, domain=None, range=None): """Initialize a new instance. Parameters ---------- - vector : `LinearSpaceElement` - Point to calculate the distance to. + vector : `LinearSpaceElement` or ``domain`` `element-like` + Point to which to calculate the distance. If ``domain`` is + given, this can be `element-like` for ``domain``, otherwise + it must be a `LinearSpaceElement`. + domain : `LinearSpace`, optional + Set of elements on which the operator can be applied. Needs + to implement ``space.dist``. Optional if ``vector`` is a + `LinearSpaceElement`, in which case ``vector.space`` is taken + as default. Otherwise this parameter is mandatory. + range : `LinearSpace` or `Field`, optional + Set to which this operator maps. + Default: ``RealNumbers``. Examples -------- >>> r2 = odl.rn(2) >>> x = r2.element([1, 1]) - >>> op = DistOperator(x) + >>> op = odl.DistOperator(x) >>> op([4, 5]) 5.0 """ - super(DistOperator, self).__init__( - vector.space, RealNumbers(), linear=False) - self.__vector = vector + if domain is None: + domain = vector.space + if range is None: + range = RealNumbers() + super(DistOperator, self).__init__(domain, range, linear=False) + self.__vector = self.domain.element(vector) @property def vector(self): @@ -722,7 +860,7 @@ def derivative(self, point): -------- >>> r2 = odl.rn(2) >>> x = r2.element([1, 1]) - >>> op = DistOperator(x) + >>> op = odl.DistOperator(x) >>> derivative = op.derivative([2, 1]) >>> derivative([1, 0]) 1.0 @@ -734,15 +872,26 @@ def derivative(self, point): raise ValueError('not differentiable at the reference vector {!r}' ''.format(self.vector)) - return InnerProductOperator(diff / dist) + return InnerProductOperator(diff / dist, self.domain, self.range) def __repr__(self): - """Return ``repr(self)``.""" - return '{}({!r})'.format(self.__class__.__name__, self.vector) + """Return ``repr(self)``. - def __str__(self): - """Return ``str(self)``.""" - return '{}({})'.format(self.__class__.__name__, self.vector) + Examples + -------- + >>> r2 = odl.rn(2) + >>> x = r2.element([1, 1]) + >>> op = odl.DistOperator(x) + >>> op + DistOperator(rn(2).element([ 1., 1.])) + """ + posargs = [self.vector] + optargs = [('domain', self.domain, + getattr(self.vector, 'space', None)), + ('range', self.range, RealNumbers())] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class ConstantOperator(Operator): @@ -760,37 +909,32 @@ def __init__(self, constant, domain=None, range=None): Parameters ---------- constant : `LinearSpaceElement` or ``range`` `element-like` - The constant space element to be returned. If ``range`` is not - provided, ``constant`` must be a `LinearSpaceElement` since the - operator range is then inferred from it. + Constant vector that should be returned. If ``domain`` is + given, this can be `element-like` for ``domain``, otherwise + it must be a `LinearSpaceElement`. domain : `LinearSpace`, optional - Domain of the operator. Default: ``vector.space`` - range : `LinearSpace`, optional - Range of the operator. Default: ``vector.space`` + Set of elements on which the operator can be applied. + Default: ``range`` if provided, otherwise ``constant.space``. + range : `LinearSpace` or `Field`, optional + Set to which this operator maps. Optional if ``constant`` is a + `LinearSpaceElement`, in which case ``constant.space`` is taken + as default. Otherwise this parameter is mandatory. Examples -------- >>> r3 = odl.rn(3) >>> x = r3.element([1, 2, 3]) - >>> op = ConstantOperator(x) + >>> op = odl.ConstantOperator(x) >>> op(x, out=r3.element()) rn(3).element([ 1., 2., 3.]) """ - - if ((domain is None or range is None) and - not isinstance(constant, LinearSpaceElement)): - raise TypeError('If either domain or range is unspecified ' - '`constant` must be LinearSpaceVector, got ' - '{!r}.'.format(constant)) - - if domain is None: - domain = constant.space if range is None: range = constant.space + if domain is None: + domain = range - self.__constant = range.element(constant) - linear = self.constant.norm() == 0 - super(ConstantOperator, self).__init__(domain, range, linear=linear) + super(ConstantOperator, self).__init__(domain, range, linear=False) + self.__constant = self.range.element(constant) @property def constant(self): @@ -804,13 +948,6 @@ def _call(self, x, out=None): else: out.assign(self.constant) - @property - def adjoint(self): - """Adjoint of the operator. - - Only defined if the operator is the constant operator. - """ - def derivative(self, point): """Derivative of this operator, always zero. @@ -822,20 +959,31 @@ def derivative(self, point): -------- >>> r3 = odl.rn(3) >>> x = r3.element([1, 2, 3]) - >>> op = ConstantOperator(x) + >>> op = odl.ConstantOperator(x) >>> deriv = op.derivative([1, 1, 1]) >>> deriv([2, 2, 2]) rn(3).element([ 0., 0., 0.]) """ - return ZeroOperator(domain=self.domain, range=self.range) + return ZeroOperator(self.domain, self.range) def __repr__(self): - """Return ``repr(self)``.""" - return '{}({!r})'.format(self.__class__.__name__, self.constant) + """Return ``repr(self)``. - def __str__(self): - """Return ``str(self)``.""" - return "{}".format(self.constant) + Examples + -------- + >>> r3 = odl.rn(3) + >>> x = r3.element([1, 2, 3]) + >>> op = odl.ConstantOperator(x) + >>> op + ConstantOperator(rn(3).element([ 1., 2., 3.])) + """ + posargs = [self.constant] + optargs = [('domain', self.domain, self.range), + ('range', self.range, + getattr(self.constant, 'space', None))] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class ZeroOperator(Operator): @@ -852,10 +1000,11 @@ def __init__(self, domain, range=None): Parameters ---------- - domain : `LinearSpace` - Domain of the operator. - range : `LinearSpace`, optional - Range of the operator. Default: ``domain`` + domain : `LinearSpace`, optional + Set of elements on which the operator can be applied. + range : `LinearSpace` or `Field`, optional + Set to which this operator maps. + Default: ``domain`` Examples -------- @@ -896,15 +1045,21 @@ def adjoint(self): If ``self.domain == self.range`` the zero operator is self-adjoint, otherwise it is the `ZeroOperator` from `range` to `domain`. """ - return ZeroOperator(domain=self.range, range=self.domain) + return ZeroOperator(self.range, self.domain) def __repr__(self): - """Return ``repr(self)``.""" - return '{}({!r})'.format(self.__class__.__name__, self.domain) + """Return ``repr(self)``. - def __str__(self): - """Return ``str(self)``.""" - return '0' + Examples + -------- + >>> op = odl.ZeroOperator(odl.rn(3)) + >>> op + ZeroOperator(rn(3)) + """ + posargs = [self.domain] + optargs = [('range', self.range, self.domain)] + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class RealPart(Operator): @@ -916,28 +1071,31 @@ class RealPart(Operator): RealPart(x) == x.real """ - def __init__(self, space): + def __init__(self, domain, range=None): """Initialize a new instance. Parameters ---------- - space : `TensorSpace` - Space in which the real part should be taken, needs to implement - ``space.real_space``. + domain : `TensorSpace` or `Field` + Space in which the real part should be taken. Needs to + implement ``domain.real_space`` and ``domain.real``. + range : `TensorSpace` or `Field`, optional + Set to which this operator maps. + Default: ``domain.real_space`` Examples -------- Take the real part of complex vector: >>> c3 = odl.cn(3) - >>> op = RealPart(c3) + >>> op = odl.RealPart(c3) >>> op([1 + 2j, 2, 3 - 1j]) rn(3).element([ 1., 2., 3.]) The operator is the identity on real spaces: >>> r3 = odl.rn(3) - >>> op = RealPart(r3) + >>> op = odl.RealPart(r3) >>> op([1, 2, 3]) rn(3).element([ 1., 2., 3.]) @@ -945,13 +1103,14 @@ def __init__(self, space): `DiscreteLp` spaces: >>> r3 = odl.uniform_discr(0, 1, 3, dtype=complex) - >>> op = RealPart(r3) + >>> op = odl.RealPart(r3) >>> op([1, 2, 3]) uniform_discr(0.0, 1.0, 3).element([ 1., 2., 3.]) """ - real_space = space.real_space - linear = (space == real_space) - super(RealPart, self).__init__(space, real_space, linear=linear) + if range is None: + range = domain.real_space + linear = (domain == range) + super(RealPart, self).__init__(domain, range, linear=linear) def _call(self, x): """Return ``self(x)``.""" @@ -966,7 +1125,7 @@ def inverse(self): The inverse is its own inverse if its domain is real: >>> r3 = odl.rn(3) - >>> op = RealPart(r3) + >>> op = odl.RealPart(r3) >>> op.inverse(op([1, 2, 3])) rn(3).element([ 1., 2., 3.]) @@ -974,14 +1133,14 @@ def inverse(self): will by necessity be lost. >>> c3 = odl.cn(3) - >>> op = RealPart(c3) + >>> op = odl.RealPart(c3) >>> op.inverse(op([1 + 2j, 2, 3 - 1j])) cn(3).element([ 1.+0.j, 2.+0.j, 3.+0.j]) """ if self.is_linear: return self else: - return ComplexEmbedding(self.domain, scalar=1) + return ComplexEmbedding(self.range, self.domain, scalar=1) @property def adjoint(self): @@ -1005,7 +1164,7 @@ def adjoint(self): The adjoint satisfies the adjoint equation for real spaces: >>> r3 = odl.rn(3) - >>> op = RealPart(r3) + >>> op = odl.RealPart(r3) >>> x = op.domain.element([1, 2, 3]) >>> y = op.range.element([3, 2, 1]) >>> x.inner(op.adjoint(y)) == op(x).inner(y) @@ -1014,7 +1173,7 @@ def adjoint(self): If the domain is complex, it only satisfies the weaker definition: >>> c3 = odl.cn(3) - >>> op = RealPart(c3) + >>> op = odl.RealPart(c3) >>> x = op.range.element([1, 2, 3]) >>> y = op.range.element([3, 2, 1]) >>> AtAxy = op(op.adjoint(x)).inner(y) @@ -1023,9 +1182,24 @@ def adjoint(self): True """ if self.is_linear: - return self + return RealPart(self.range, self.domain) else: - return ComplexEmbedding(self.domain, scalar=1) + return ComplexEmbedding(self.range, self.domain, scalar=1) + + def __repr__(self): + """Return ``repr(self)``. + + Examples + -------- + >>> c3 = odl.cn(3) + >>> op = odl.RealPart(c3) + >>> op + RealPart(cn(3)) + """ + posargs = [self.domain] + optargs = [('range', self.range, self.domain.real_space)] + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class ImagPart(Operator): @@ -1037,34 +1211,38 @@ class ImagPart(Operator): ImagPart(x) == x.imag """ - def __init__(self, space): + def __init__(self, domain, range=None): """Initialize a new instance. Parameters ---------- - space : `TensorSpace` - Space in which the imaginary part should be taken, needs to - implement ``space.real_space``. + domain : `TensorSpace` or `Field` + Space in which the imaginary part should be taken. Needs to + implement ``domain.real_space`` and ``domain.imag``. + range : `TensorSpace` or `Field`, optional + Set to which this operator maps. + Default: ``domain.real_space`` Examples -------- Take the imaginary part of complex vector: >>> c3 = odl.cn(3) - >>> op = ImagPart(c3) + >>> op = odl.ImagPart(c3) >>> op([1 + 2j, 2, 3 - 1j]) rn(3).element([ 2., 0., -1.]) The operator is the zero operator on real spaces: >>> r3 = odl.rn(3) - >>> op = ImagPart(r3) + >>> op = odl.ImagPart(r3) >>> op([1, 2, 3]) rn(3).element([ 0., 0., 0.]) """ - real_space = space.real_space - linear = (space == real_space) - super(ImagPart, self).__init__(space, real_space, linear=linear) + if range is None: + range = domain.real_space + linear = (domain == range) + super(ImagPart, self).__init__(domain, range, linear=linear) def _call(self, x): """Return ``self(x)``.""" @@ -1079,7 +1257,7 @@ def inverse(self): The inverse is the zero operator if the domain is real: >>> r3 = odl.rn(3) - >>> op = ImagPart(r3) + >>> op = odl.ImagPart(r3) >>> op.inverse(op([1, 2, 3])) rn(3).element([ 0., 0., 0.]) @@ -1087,14 +1265,14 @@ def inverse(self): will by necessity be lost. >>> c3 = odl.cn(3) - >>> op = ImagPart(c3) + >>> op = odl.ImagPart(c3) >>> op.inverse(op([1 + 2j, 2, 3 - 1j])) cn(3).element([ 0.+2.j, 0.+0.j, -0.-1.j]) """ if self.is_linear: - return ZeroOperator(self.domain) + return ZeroOperator(self.range, self.domain) else: - return ComplexEmbedding(self.domain, scalar=1j) + return ComplexEmbedding(self.range, self.domain, scalar=1j) @property def adjoint(self): @@ -1118,7 +1296,7 @@ def adjoint(self): The adjoint satisfies the adjoint equation for real spaces: >>> r3 = odl.rn(3) - >>> op = ImagPart(r3) + >>> op = odl.ImagPart(r3) >>> x = op.domain.element([1, 2, 3]) >>> y = op.range.element([3, 2, 1]) >>> x.inner(op.adjoint(y)) == op(x).inner(y) @@ -1127,7 +1305,7 @@ def adjoint(self): If the domain is complex, it only satisfies the weaker definition: >>> c3 = odl.cn(3) - >>> op = ImagPart(c3) + >>> op = odl.ImagPart(c3) >>> x = op.range.element([1, 2, 3]) >>> y = op.range.element([3, 2, 1]) >>> AtAxy = op(op.adjoint(x)).inner(y) @@ -1136,9 +1314,24 @@ def adjoint(self): True """ if self.is_linear: - return ZeroOperator(self.domain) + return ZeroOperator(self.range, self.domain) else: - return ComplexEmbedding(self.domain, scalar=1j) + return ComplexEmbedding(self.range, self.domain, scalar=1j) + + def __repr__(self): + """Return ``repr(self)``. + + Examples + -------- + >>> c3 = odl.cn(3) + >>> op = odl.ImagPart(c3) + >>> op + ImagPart(cn(3)) + """ + posargs = [self.domain] + optargs = [('range', self.range, self.domain.real_space)] + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class ComplexEmbedding(Operator): @@ -1150,14 +1343,17 @@ class ComplexEmbedding(Operator): ComplexEmbedding(space)(x) <==> space.complex_space.element(x) """ - def __init__(self, space, scalar=1.0): + def __init__(self, domain, range=None, scalar=1.0): """Initialize a new instance. Parameters ---------- - space : `TensorSpace` + domain : `TensorSpace` or `Field` Space that should be embedded into its complex counterpart. - It must implement `TensorSpace.complex_space`. + Needs to implement ``domain.complex_space``. + range : `TensorSpace` or `Field`, optional + Set to which this operator maps. + Default: ``domain.complex_space`` scalar : ``space.complex_space.field`` element, optional Scalar to be multiplied with incoming vectors in order to get the complex vector. @@ -1167,13 +1363,13 @@ def __init__(self, space, scalar=1.0): Embed real vector into complex space: >>> r3 = odl.rn(3) - >>> op = ComplexEmbedding(r3) + >>> op = odl.ComplexEmbedding(r3) >>> op([1, 2, 3]) cn(3).element([ 1.+0.j, 2.+0.j, 3.+0.j]) Embed real vector as imaginary part into complex space: - >>> op = ComplexEmbedding(r3, scalar=1j) + >>> op = odl.ComplexEmbedding(r3, scalar=1j) >>> op([1, 2, 3]) cn(3).element([ 0.+1.j, 0.+2.j, 0.+3.j]) @@ -1181,14 +1377,14 @@ def __init__(self, space, scalar=1.0): scalar: >>> c3 = odl.cn(3) - >>> op = ComplexEmbedding(c3, scalar=1 + 2j) + >>> op = odl.ComplexEmbedding(c3, scalar=1 + 2j) >>> op([1 + 1j, 2 + 2j, 3 + 3j]) cn(3).element([-1.+3.j, -2.+6.j, -3.+9.j]) """ - complex_space = space.complex_space - self.scalar = complex_space.field.element(scalar) - super(ComplexEmbedding, self).__init__( - space, complex_space, linear=True) + if range is None: + range = domain.complex_space + super(ComplexEmbedding, self).__init__(domain, range, linear=True) + self.scalar = self.range.field.element(scalar) def _call(self, x, out): """Return ``self(x)``.""" @@ -1210,7 +1406,7 @@ def inverse(self): Examples -------- >>> r3 = odl.rn(3) - >>> op = ComplexEmbedding(r3, scalar=1) + >>> op = odl.ComplexEmbedding(r3, scalar=1) >>> op.inverse(op([1, 2, 4])) rn(3).element([ 1., 2., 4.]) """ @@ -1218,17 +1414,20 @@ def inverse(self): # Real domain # Optimizations for simple cases. if self.scalar.real == self.scalar: - return (1 / self.scalar.real) * RealPart(self.range) + return (1 / self.scalar.real) * RealPart(self.range, + self.domain) elif 1j * self.scalar.imag == self.scalar: - return (1 / self.scalar.imag) * ImagPart(self.range) + return (1 / self.scalar.imag) * ImagPart(self.range, + self.domain) else: # General case inv_scalar = (1 / self.scalar).conjugate() - return ((inv_scalar.real) * RealPart(self.range) + - (inv_scalar.imag) * ImagPart(self.range)) + return ((inv_scalar.real) * RealPart(self.range, self.domain) + + (inv_scalar.imag) * ImagPart(self.range, self.domain)) else: # Complex domain - return ComplexEmbedding(self.range, self.scalar.conjugate()) + return ComplexEmbedding(self.range, self.domain, + self.scalar.conjugate()) @property def adjoint(self): @@ -1275,30 +1474,55 @@ def adjoint(self): # Real domain # Optimizations for simple cases. if self.scalar.real == self.scalar: - return self.scalar.real * RealPart(self.range) + return self.scalar.real * ComplexEmbedding(self.range, + self.domain) elif 1j * self.scalar.imag == self.scalar: - return self.scalar.imag * ImagPart(self.range) + return self.scalar.imag * ImagPart(self.range, self.domain) else: # General case - return (self.scalar.real * RealPart(self.range) + - self.scalar.imag * ImagPart(self.range)) + return (self.scalar.real * RealPart(self.range, self.domain) + + self.scalar.imag * ImagPart(self.range, self.domain)) else: # Complex domain - return ComplexEmbedding(self.range, self.scalar.conjugate()) + return ComplexEmbedding(self.range, self.domain, + self.scalar.conjugate()) + + def __repr__(self): + """Return ``repr(self)``. + + Examples + -------- + >>> r3 = odl.rn(3) + >>> op = odl.ComplexEmbedding(r3) + >>> op + ComplexEmbedding(rn(3)) + >>> op = odl.ComplexEmbedding(r3, scalar=1j) + >>> op + ComplexEmbedding(rn(3), scalar=1j) + """ + posargs = [self.domain] + optargs = [('range', self.range, self.domain.complex_space), + ('scalar', self.scalar, 1.0)] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class ComplexModulus(Operator): - """Operator that computes the modolus (absolute value) of a vector.""" + """Operator that computes the modulus (absolute value) of a vector.""" - def __init__(self, space): + def __init__(self, domain, range=None): """Initialize a new instance. Parameters ---------- - space : `FnBase` - Space whose real part should be taken, needs to implement - ``space.real_space``. + domain : `TensorSpace` or `Field` + Space in which the complex modulus should be taken. Needs to + implement ``domain.real_space``. + range : `TensorSpace` or `Field`, optional + Set to which this operator maps. + Default: ``domain.real_space`` Examples -------- @@ -1316,21 +1540,26 @@ def __init__(self, space): >>> op([1, -2]) rn(2).element([ 1., 2.]) - The operator also works on other `FnBase` spaces such as - `DiscreteLp` spaces: + The operator also works on other `TensorSpace`'s such as + `DiscreteLp`: >>> r2 = odl.uniform_discr(0, 1, 2, dtype=complex) >>> op = odl.ComplexModulus(r2) >>> op([3 + 4j, 2]) uniform_discr(0.0, 1.0, 2).element([ 5., 2.]) """ - real_space = space.real_space - linear = (space == real_space) - super(ComplexModulus, self).__init__(space, real_space, linear=linear) + if range is None: + range = domain.real_space + linear = (domain == range) + super(ComplexModulus, self).__init__(domain, range, linear=linear) def _call(self, x): """Return ``self(x)``.""" - return (x.real ** 2 + x.imag ** 2).ufuncs.sqrt() + squared_mod = x.real ** 2 + x.imag ** 2 + if hasattr(squared_mod, 'ufuncs'): + return squared_mod.ufuncs.sqrt() + else: + return np.sqrt(squared_mod) @property def inverse(self): @@ -1341,7 +1570,7 @@ def inverse(self): The (pseudo-)inverse in the real case is the identity: >>> r2 = odl.rn(2) - >>> op = ComplexModulus(r2) + >>> op = odl.ComplexModulus(r2) >>> op.inverse(op([1, -2])) rn(2).element([ 1., 2.]) @@ -1349,14 +1578,30 @@ def inverse(self): positive weights to the real and complex parts: >>> c2 = odl.cn(2) - >>> op = ComplexModulus(c2) + >>> op = odl.ComplexModulus(c2) >>> op.inverse(op([np.sqrt(2), 2 + 2j])) cn(2).element([ 1.+1.j, 2.+2.j]) """ if self.is_linear: - return IdentityOperator(self.domain) + return IdentityOperator(self.range, self.domain) else: - return ComplexEmbedding(self.range, scalar=(1 + 1j) / np.sqrt(2)) + return ComplexEmbedding(self.range, self.domain, + scalar=(1 + 1j) / np.sqrt(2)) + + def __repr__(self): + """Return ``repr(self)``. + + Examples + -------- + >>> c2 = odl.cn(2) + >>> op = odl.ComplexModulus(c2) + >>> op + ComplexModulus(cn(2)) + """ + posargs = [self.domain] + optargs = [('range', self.range, self.domain.real_space)] + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) if __name__ == '__main__': diff --git a/odl/operator/operator.py b/odl/operator/operator.py index 1a5e812f27c..86e42008be8 100644 --- a/odl/operator/operator.py +++ b/odl/operator/operator.py @@ -8,16 +8,20 @@ """Abstract mathematical operators.""" -from __future__ import print_function, division, absolute_import -from builtins import object +from __future__ import absolute_import, division, print_function + import inspect -from numbers import Number, Integral import sys +from builtins import object +from numbers import Integral, Number -from odl.set import LinearSpace, Set, Field -from odl.set.space import LinearSpaceElement -from odl.util import cache_arguments +import numpy as np +from odl.set import Field, LinearSpace, Set +from odl.set.space import LinearSpaceElement +from odl.util import ( + REPR_PRECISION, cache_arguments, npy_printoptions, repr_string, + signature_string_parts) __all__ = ('Operator', 'OperatorComp', 'OperatorSum', 'OperatorVectorSum', 'OperatorLeftScalarMult', 'OperatorRightScalarMult', @@ -576,9 +580,8 @@ def adjoint(self): OpNotImplementedError Since the adjoint cannot be default implemented. """ - raise OpNotImplementedError('adjoint not implemented ' - 'for operator {!r}' - ''.format(self)) + raise OpNotImplementedError( + 'adjoint not implemented for operator {!r}'.format(self)) def derivative(self, point): """Return the operator derivative at ``point``. @@ -592,9 +595,8 @@ def derivative(self, point): if self.is_linear: return self else: - raise OpNotImplementedError('derivative not implemented ' - 'for operator {!r}' - ''.format(self)) + raise OpNotImplementedError( + 'derivative not implemented for operator {!r}'.format(self)) @property def inverse(self): @@ -605,8 +607,8 @@ def inverse(self): OpNotImplementedError Since the inverse cannot be default implemented. """ - raise OpNotImplementedError('inverse not implemented for operator {!r}' - ''.format(self)) + raise OpNotImplementedError( + 'inverse not implemented for operator {!r}'.format(self)) def __call__(self, x, out=None, **kwargs): """Return ``self(x[, out, **kwargs])``. @@ -723,8 +725,8 @@ def norm(self, estimate=False, **kwargs): Some operators know their own operator norm and do not need an estimate >>> spc = odl.rn(3) - >>> id = odl.IdentityOperator(spc) - >>> id.norm(True) + >>> op = odl.IdentityOperator(spc) + >>> op.norm(estimate=True) 1.0 For others, there is no closed form expression and an estimate is @@ -735,9 +737,10 @@ def norm(self, estimate=False, **kwargs): >>> opnorm = grad.norm(estimate=True) """ if not estimate: - raise NotImplementedError('`Operator.norm()` not implemented, use ' - '`Operator.norm(estimate=True)` to ' - 'obtain an estimate.') + raise OpNotImplementedError( + '`norm` not implemented for operator {!r}, use ' + '`Operator.norm(estimate=True)` to obtain an estimate.' + ''.format(self)) else: from odl.operator.oputils import power_method_opnorm return power_method_opnorm(self, **kwargs) @@ -841,8 +844,8 @@ def __mul__(self, other): >>> x = rn.element([1, 2, 3]) >>> op(x) rn(3).element([ 1., 2., 3.]) - >>> Scaled = op * 3 - >>> Scaled(x) + >>> scaled = op * 3 + >>> scaled(x) rn(3).element([ 3., 6., 9.]) """ if isinstance(other, Operator): @@ -923,8 +926,8 @@ def __rmul__(self, other): >>> x = rn.element([1, 2, 3]) >>> op(x) rn(3).element([ 1., 2., 3.]) - >>> Scaled = 3 * op - >>> Scaled(x) + >>> scaled = 3 * op + >>> scaled(x) rn(3).element([ 3., 6., 9.]) """ if isinstance(other, Operator): @@ -1016,8 +1019,8 @@ def __truediv__(self, other): >>> x = rn.element([3, 6, 9]) >>> op(x) rn(3).element([ 3., 6., 9.]) - >>> Scaled = op / 3.0 - >>> Scaled(x) + >>> scaled = op / 3.0 + >>> scaled(x) rn(3).element([ 1., 2., 3.]) """ if isinstance(other, Number): @@ -1044,8 +1047,15 @@ def __repr__(self): The default `repr` implementation. Should be overridden by subclasses. """ - return '{}: {!r} -> {!r}'.format(self.__class__.__name__, self.domain, - self.range) + linewidth = np.get_printoptions()['linewidth'] + oneline_str = '{}: {!r} -> {!r}'.format(self.__class__.__name__, + self.domain, self.range) + if len(oneline_str) <= linewidth: + return oneline_str + else: + return '{}:\n {!r} ->\n {!r}'.format(self.__class__.__name__, + self.domain, + self.range) def __str__(self): """Return ``str(self)``. @@ -1055,7 +1065,7 @@ def __str__(self): """ return self.__class__.__name__ - # Give a `Operator` a higher priority than any NumPy array type. This + # Give an `Operator` a higher priority than any NumPy array type. This # forces the usage of `__op__` of `Operator` if the other operand # is a NumPy object (applies also to scalars!). # Set higher than LinearSpaceElement.__array_priority__ to handle @@ -1098,11 +1108,11 @@ def __init__(self, left, right, tmp_ran=None, tmp_dom=None): >>> op = odl.IdentityOperator(r3) >>> x = r3.element([1, 2, 3]) >>> out = r3.element() - >>> OperatorSum(op, op)(x, out) # In-place, returns out + >>> odl.OperatorSum(op, op)(x, out) # In-place, returns out rn(3).element([ 2., 4., 6.]) >>> out rn(3).element([ 2., 4., 6.]) - >>> OperatorSum(op, op)(x) + >>> odl.OperatorSum(op, op)(x) rn(3).element([ 2., 4., 6.]) """ if left.range != right.range: @@ -1199,8 +1209,10 @@ def adjoint(self): def __repr__(self): """Return ``repr(self)``.""" - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.left, self.right) + posargs = [self.left, self.right] + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): """Return ``str(self)``.""" @@ -1229,8 +1241,8 @@ def __init__(self, operator, vector): -------- >>> r3 = odl.rn(3) >>> y = r3.element([1, 2, 3]) - >>> ident_op = odl.IdentityOperator(r3) - >>> sum_op = odl.OperatorVectorSum(ident_op, y) + >>> I = odl.IdentityOperator(r3) + >>> sum_op = odl.OperatorVectorSum(I, y) >>> x = r3.element([4, 5, 6]) >>> sum_op(x) rn(3).element([ 5., 7., 9.]) @@ -1293,12 +1305,15 @@ def derivative(self, point): def __repr__(self): """Return ``repr(self)``.""" - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.operator, self.vector) + posargs = [self.operator, self.vector] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): """Return ``str(self)``.""" - return '({} + {})'.format(self.left, self.right) + return '({} + {})'.format(self.operator, self.vector) class OperatorComp(Operator): @@ -1316,7 +1331,7 @@ def __init__(self, left, right, tmp=None): Parameters ---------- left : `Operator` - The left ("outer") operator + The left ("outer") operator. right : `Operator` The right ("inner") operator. Its range must coincide with the domain of ``left``. @@ -1431,8 +1446,10 @@ def adjoint(self): def __repr__(self): """Return ``repr(self)``.""" - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.left, self.right) + posargs = [self.left, self.right] + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): """Return ``str(self)``.""" @@ -1451,11 +1468,9 @@ def __init__(self, left, right): Parameters ---------- - left : `Operator` - The first factor - right : `Operator` - The second factor. Must have the same domain and range as - ``left``. + left, right : `Operator` + Factors in the pointwise product. They must have the same domain + and range. """ if left.range != right.range: raise OpTypeError('operator ranges {!r} and {!r} do not match' @@ -1505,8 +1520,10 @@ def derivative(self, x): def __repr__(self): """Return ``repr(self)``.""" - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.left, self.right) + posargs = [self.left, self.right] + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): """Return ``str(self)``.""" @@ -1538,8 +1555,8 @@ def __init__(self, operator, scalar): Examples -------- >>> space = odl.rn(3) - >>> operator = odl.IdentityOperator(space) - >>> left_mul_op = OperatorLeftScalarMult(operator, 3) + >>> I = odl.IdentityOperator(space) + >>> left_mul_op = odl.OperatorLeftScalarMult(I, 3) >>> left_mul_op([1, 2, 3]) rn(3).element([ 3., 6., 9.]) """ @@ -1596,8 +1613,8 @@ def inverse(self): Examples -------- >>> space = odl.rn(3) - >>> operator = odl.IdentityOperator(space) - >>> left_mul_op = OperatorLeftScalarMult(operator, 3) + >>> I = odl.IdentityOperator(space) + >>> left_mul_op = odl.OperatorLeftScalarMult(I, 3) >>> left_mul_op.inverse([3, 3, 3]) rn(3).element([ 1., 1., 1.]) """ @@ -1626,8 +1643,8 @@ def derivative(self, x): Examples -------- >>> space = odl.rn(3) - >>> operator = odl.IdentityOperator(space) - space.element([1, 1, 1]) - >>> left_mul_op = OperatorLeftScalarMult(operator, 3) + >>> op = odl.IdentityOperator(space) - space.element([1, 1, 1]) + >>> left_mul_op = odl.OperatorLeftScalarMult(op, 3) >>> derivative = left_mul_op.derivative([0, 0, 0]) >>> derivative([1, 1, 1]) rn(3).element([ 3., 3., 3.]) @@ -1655,8 +1672,8 @@ def adjoint(self): Examples -------- >>> space = odl.rn(3) - >>> operator = odl.IdentityOperator(space) - >>> left_mul_op = OperatorLeftScalarMult(operator, 3) + >>> I = odl.IdentityOperator(space) + >>> left_mul_op = odl.OperatorLeftScalarMult(I, 3) >>> left_mul_op.adjoint([1, 2, 3]) rn(3).element([ 3., 6., 9.]) """ @@ -1668,8 +1685,11 @@ def adjoint(self): def __repr__(self): """Return ``repr(self)``.""" - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.operator, self.scalar) + posargs = [self.operator, self.scalar] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): """Return ``str(self)``.""" @@ -1704,9 +1724,9 @@ def __init__(self, operator, scalar, tmp=None): Examples -------- >>> space = odl.rn(3) - >>> operator = odl.IdentityOperator(space) - >>> left_mul_op = OperatorRightScalarMult(operator, 3) - >>> left_mul_op([1, 2, 3]) + >>> I = odl.IdentityOperator(space) + >>> right_mul_op = odl.OperatorRightScalarMult(I, 3) + >>> right_mul_op([1, 2, 3]) rn(3).element([ 3., 6., 9.]) """ if not isinstance(operator.domain, (LinearSpace, Field)): @@ -1783,9 +1803,9 @@ def inverse(self): Examples -------- >>> space = odl.rn(3) - >>> operator = odl.IdentityOperator(space) - >>> left_mul_op = OperatorRightScalarMult(operator, 3) - >>> left_mul_op.inverse([3, 3, 3]) + >>> I = odl.IdentityOperator(space) + >>> right_mul_op = odl.OperatorRightScalarMult(I, 3) + >>> right_mul_op.inverse([3, 3, 3]) rn(3).element([ 1., 1., 1.]) """ if self.scalar == 0.0: @@ -1810,9 +1830,9 @@ def derivative(self, x): Examples -------- >>> space = odl.rn(3) - >>> operator = odl.IdentityOperator(space) - space.element([1, 1, 1]) - >>> left_mul_op = OperatorRightScalarMult(operator, 3) - >>> derivative = left_mul_op.derivative([0, 0, 0]) + >>> op = odl.IdentityOperator(space) - space.element([1, 1, 1]) + >>> right_mul_op = odl.OperatorRightScalarMult(op, 3) + >>> derivative = right_mul_op.derivative([0, 0, 0]) >>> derivative([1, 1, 1]) rn(3).element([ 3., 3., 3.]) """ @@ -1836,9 +1856,9 @@ def adjoint(self): Examples -------- >>> space = odl.rn(3) - >>> operator = odl.IdentityOperator(space) - >>> left_mul_op = OperatorRightScalarMult(operator, 3) - >>> left_mul_op.adjoint([1, 2, 3]) + >>> I = odl.IdentityOperator(space) + >>> right_mul_op = odl.OperatorRightScalarMult(I, 3) + >>> right_mul_op.adjoint([1, 2, 3]) rn(3).element([ 3., 6., 9.]) """ @@ -1849,8 +1869,11 @@ def adjoint(self): def __repr__(self): """Return ``repr(self)``.""" - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.operator, self.scalar) + posargs = [self.operator, self.scalar] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): """Return ``str(self)``.""" @@ -1888,7 +1911,7 @@ def __init__(self, functional, vector): >>> space = odl.rn(3) >>> y = space.element([1, 2, 3]) >>> functional = odl.InnerProductOperator(y) - >>> left_mul_op = FunctionalLeftVectorMult(functional, y) + >>> left_mul_op = odl.FunctionalLeftVectorMult(functional, y) >>> left_mul_op([1, 2, 3]) rn(3).element([ 14., 28., 42.]) """ @@ -1965,8 +1988,11 @@ def adjoint(self): def __repr__(self): """Return ``repr(self)``.""" - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.functional, self.vector) + posargs = [self.functional, self.vector] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): """Return ``str(self)``.""" @@ -1986,9 +2012,9 @@ def __init__(self, operator, vector): Parameters ---------- operator : `Operator` - The range of ``op`` must be a `LinearSpace`. + The operator in the product. Its range must be a `LinearSpace`. vector : `LinearSpaceElement` in ``op.range`` - The vector to multiply by + The vector to multiply with. """ if vector not in operator.range: raise OpRangeError('`vector` {!r} not in operator.range {!r}' @@ -2078,8 +2104,11 @@ def adjoint(self): def __repr__(self): """Return ``repr(self)``.""" - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.operator, self.vector) + posargs = [self.operator, self.vector] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): """Return ``str(self)``.""" @@ -2198,8 +2227,11 @@ def adjoint(self): def __repr__(self): """Return ``repr(self)``.""" - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.operator, self.vector) + posargs = [self.operator, self.vector] + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): """Return ``str(self)``.""" diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index 4a584f23a6f..8a08a5fb1d9 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -8,12 +8,15 @@ """Convenience functions for operators.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function from future.utils import native + import numpy as np from odl.space.base_tensors import TensorSpace -from odl.space import ProductSpace +from odl.space.pspace import ProductSpace, ProductSpaceElement +from odl.space.weighting import ArrayWeighting, ConstWeighting +from odl.util import OptionalArgDecorator __all__ = ('matrix_representation', 'power_method_opnorm', 'as_scipy_operator', 'as_scipy_functional', 'as_proximal_lang_operator') diff --git a/odl/operator/pspace_ops.py b/odl/operator/pspace_ops.py index 1cbbda36ec5..1627bcdd041 100644 --- a/odl/operator/pspace_ops.py +++ b/odl/operator/pspace_ops.py @@ -8,17 +8,19 @@ """Default operators defined on any `ProductSpace`.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + from numbers import Integral + import numpy as np -from odl.operator.operator import Operator from odl.operator.default_ops import ZeroOperator +from odl.operator.operator import Operator from odl.space import ProductSpace +from odl.util import ( + array_str, attribute_repr_string, repr_string, signature_string_parts) - -__all__ = ('ProductSpaceOperator', - 'ComponentProjection', 'ComponentProjectionAdjoint', +__all__ = ('ProductSpaceOperator', 'ComponentProjection', 'BroadcastOperator', 'ReductionOperator', 'DiagonalOperator') @@ -114,30 +116,38 @@ def __init__(self, operators, domain=None, range=None): Sum of elements: - >>> prod_op = ProductSpaceOperator([I, I]) + >>> prod_op = odl.ProductSpaceOperator([I, I]) >>> prod_op(x) - ProductSpace(rn(3), 1).element([ - [ 5., 7., 9.] - ]) + ProductSpace(rn(3), 1).element( + [ + [ 5., 7., 9.] + ] + ) Diagonal operator -- 0 or ``None`` means ignore, or the implicit zero operator: - >>> prod_op = ProductSpaceOperator([[I, 0], [0, I]]) + >>> prod_op = odl.ProductSpaceOperator([[I, 0], + ... [0, I]]) >>> prod_op(x) - ProductSpace(rn(3), 2).element([ - [ 1., 2., 3.], - [ 4., 5., 6.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 1., 2., 3.], + [ 4., 5., 6.] + ] + ) Complicated combinations: - >>> prod_op = ProductSpaceOperator([[I, I], [I, 0]]) + >>> prod_op = odl.ProductSpaceOperator([[I, I], + ... [I, 0]]) >>> prod_op(x) - ProductSpace(rn(3), 2).element([ - [ 5., 7., 9.], - [ 1., 2., 3.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 5., 7., 9.], + [ 1., 2., 3.] + ] + ) """ # Lazy import to improve `import odl` time import scipy.sparse @@ -245,7 +255,7 @@ def derivative(self, x): Parameters ---------- x : `domain` element - The point to take the derivative in + Point in which to take the derivative. Returns ------- @@ -257,44 +267,55 @@ def derivative(self, x): >>> r3 = odl.rn(3) >>> pspace = odl.ProductSpace(r3, r3) >>> I = odl.IdentityOperator(r3) - >>> x = pspace.element([[1, 2, 3], [4, 5, 6]]) + >>> x = pspace.element([[1, 2, 3], + ... [4, 5, 6]]) - Example with linear operator (derivative is itself) + Example with linear operator (derivative is itself): - >>> prod_op = ProductSpaceOperator([[0, I], [0, 0]], - ... domain=pspace, range=pspace) + >>> prod_op = odl.ProductSpaceOperator([[0, I], + ... [0, 0]], + ... domain=pspace, range=pspace) >>> prod_op(x) - ProductSpace(rn(3), 2).element([ - [ 4., 5., 6.], - [ 0., 0., 0.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 4., 5., 6.], + [ 0., 0., 0.] + ] + ) >>> prod_op.derivative(x)(x) - ProductSpace(rn(3), 2).element([ - [ 4., 5., 6.], - [ 0., 0., 0.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 4., 5., 6.], + [ 0., 0., 0.] + ] + ) Example with affine operator >>> residual_op = I - r3.element([1, 1, 1]) - >>> op = ProductSpaceOperator([[0, residual_op], [0, 0]], - ... domain=pspace, range=pspace) + >>> op = odl.ProductSpaceOperator([[0, residual_op], + ... [0, 0]], + ... domain=pspace, range=pspace) - Calling operator gives offset by [1, 1, 1] + Calling operator gives offset by [1, 1, 1]: >>> op(x) - ProductSpace(rn(3), 2).element([ - [ 3., 4., 5.], - [ 0., 0., 0.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 3., 4., 5.], + [ 0., 0., 0.] + ] + ) Derivative of affine operator does not have this offset >>> op.derivative(x)(x) - ProductSpace(rn(3), 2).element([ - [ 4., 5., 6.], - [ 0., 0., 0.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 4., 5., 6.], + [ 0., 0., 0.] + ] + ) """ # Lazy import to improve `import odl` time import scipy.sparse @@ -323,7 +344,7 @@ def adjoint(self): Returns ------- adjoint : `ProductSpaceOperator` - The adjoint + The adjoint operator. Examples -------- @@ -335,18 +356,23 @@ def adjoint(self): Matrix is transposed: - >>> prod_op = ProductSpaceOperator([[0, I], [0, 0]], - ... domain=pspace, range=pspace) + >>> prod_op = odl.ProductSpaceOperator([[0, I], + ... [0, 0]], + ... domain=pspace, range=pspace) >>> prod_op(x) - ProductSpace(rn(3), 2).element([ - [ 4., 5., 6.], - [ 0., 0., 0.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 4., 5., 6.], + [ 0., 0., 0.] + ] + ) >>> prod_op.adjoint(x) - ProductSpace(rn(3), 2).element([ - [ 0., 0., 0.], - [ 1., 2., 3.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 0., 0., 0.], + [ 1., 2., 3.] + ] + ) """ # Lazy import to improve `import odl` time import scipy.sparse @@ -379,9 +405,9 @@ def __getitem__(self, index): >>> r3 = odl.rn(3) >>> pspace = odl.ProductSpace(r3, r3) >>> I = odl.IdentityOperator(r3) - >>> prod_op = ProductSpaceOperator([[0, I], - ... [0, 0]], - ... domain=pspace, range=pspace) + >>> prod_op = odl.ProductSpaceOperator([[0, I], + ... [0, 0]], + ... domain=pspace, range=pspace) >>> prod_op[0, 0] 0 >>> prod_op[0, 1] @@ -435,11 +461,65 @@ def size(self): return np.prod(self.shape) def __repr__(self): - """Return ``repr(self)``.""" - aslist = [[0] * len(self.domain) for _ in range(len(self.range))] + """Return ``repr(self)``. + + Examples + -------- + All rows and columns with an operator, i.e., domain and range + are unambiguous: + + >>> r3 = odl.rn(3) + >>> pspace = odl.ProductSpace(r3, r3) + >>> I = odl.IdentityOperator(r3) + >>> prod_op = odl.ProductSpaceOperator([[I, 0], + ... [0, I]]) + >>> prod_op + ProductSpaceOperator( + [[IdentityOperator(rn(3)), 0], + [0, IdentityOperator(rn(3))]] + ) + + Domain not completely determined (column with all zeros): + + >>> prod_op = odl.ProductSpaceOperator([[I, 0], + ... [I, 0]], domain=pspace) + >>> prod_op + ProductSpaceOperator( + [[IdentityOperator(rn(3)), 0], + [IdentityOperator(rn(3)), 0]], + domain=ProductSpace(rn(3), 2) + ) + + Range not completely determined (row with all zeros): + + >>> prod_op = odl.ProductSpaceOperator([[I, I], + ... [0, 0]], range=pspace) + >>> prod_op + ProductSpaceOperator( + [[IdentityOperator(rn(3)), IdentityOperator(rn(3))], + [0, 0]], + range=ProductSpace(rn(3), 2) + ) + """ + op_arr = np.zeros(self.ops.shape, dtype=object) for i, j, op in zip(self.ops.row, self.ops.col, self.ops.data): - aslist[i][j] = op - return '{}({!r})'.format(self.__class__.__name__, aslist) + op_arr[i, j] = op + + posargs = [op_arr] + posmod = array_str + optargs = [] + + if any(j not in self.ops.col for j in range(self.ops.shape[1])): + # Not all columns have operators, need explicit domain + optargs.append(('domain', self.domain, None)) + if any(i not in self.ops.row for i in range(self.ops.shape[0])): + # Not all rows have operators, need explicit range + optargs.append(('range', self.range, None)) + + inner_parts = signature_string_parts(posargs, optargs, + mod=[posmod, '!r']) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) class ComponentProjection(Operator): @@ -495,18 +575,25 @@ def __init__(self, space, index): >>> proj = odl.ComponentProjection(pspace, [0, 2]) >>> proj(x) - ProductSpace(rn(1), rn(3)).element([ - [ 1.], - [ 4., 5., 6.] - ]) + ProductSpace(rn(1), rn(3)).element( + [ + [ 1.], + [ 4., 5., 6.] + ] + ) """ - self.__index = index + if np.isscalar(index): + self.__index = int(index) + elif isinstance(index, slice): + self.__index = index + else: + self.__index = list(index) super(ComponentProjection, self).__init__( space, space[index], linear=True) @property def index(self): - """Index of the subspace.""" + """Index, indices or slice defining the subspace.""" return self.__index def _call(self, x, out=None): @@ -519,113 +606,76 @@ def _call(self, x, out=None): @property def adjoint(self): - """The adjoint operator. - - The adjoint is given by extending along `ComponentProjection.index`, - and setting zero along the others. - - See Also - -------- - ComponentProjectionAdjoint - """ - return ComponentProjectionAdjoint(self.domain, self.index) - - def __repr__(self): - """Return ``repr(self)``. - - Examples - -------- - >>> pspace = odl.ProductSpace(odl.rn(1), odl.rn(2)) - >>> odl.ComponentProjection(pspace, 0) - ComponentProjection(ProductSpace(rn(1), rn(2)), 0) - """ - return '{}({!r}, {})'.format(self.__class__.__name__, - self.domain, self.index) - - -class ComponentProjectionAdjoint(Operator): - - """Adjoint operator to `ComponentProjection`. - - As a special case of the adjoint of a `ProductSpaceOperator`, - this operator is given as a column vector of identity operators - and zero operators, with the identities placed in the positions - defined by `ComponentProjectionAdjoint.index`. - - In weighted product spaces, the adjoint needs to take the - weightings into account. This is currently not supported. - """ - - def __init__(self, space, index): - """Initialize a new instance - - Parameters - ---------- - space : `ProductSpace` - Space to project to. - index : int, slice, or list - Indexes to project from. - - Examples - -------- - >>> r1 = odl.rn(1) - >>> r2 = odl.rn(2) - >>> r3 = odl.rn(3) - >>> pspace = odl.ProductSpace(r1, r2, r3) - >>> x = pspace.element([[1], - ... [2, 3], - ... [4, 5, 6]]) - - Projection on the 0-th component: - - >>> proj_adj = odl.ComponentProjectionAdjoint(pspace, 0) - >>> proj_adj(x[0]) - ProductSpace(rn(1), rn(2), rn(3)).element([ - [ 1.], - [ 0., 0.], - [ 0., 0., 0.] - ]) - - Projection on a sub-space corresponding to indices 0 and 2: - - >>> proj_adj = odl.ComponentProjectionAdjoint(pspace, [0, 2]) - >>> proj_adj(x[[0, 2]]) - ProductSpace(rn(1), rn(2), rn(3)).element([ - [ 1.], - [ 0., 0.], - [ 4., 5., 6.] - ]) - """ - self.__index = index - super(ComponentProjectionAdjoint, self).__init__( - space[index], space, linear=True) - - @property - def index(self): - """Index of the subspace.""" - return self.__index + """Adjoint operator extending by zeros from subspace to full space.""" + op = self + + class ComponentProjectionAdjoint(Operator): + + """Adjoint operator to `ComponentProjection`.""" + + def __init__(self): + """Initialize a new instance. + + Examples + -------- + >>> r1 = odl.rn(1) + >>> r2 = odl.rn(2) + >>> r3 = odl.rn(3) + >>> pspace = odl.ProductSpace(r1, r2, r3) + >>> x = pspace.element([[1], + ... [2, 3], + ... [4, 5, 6]]) + + Projection on the 0-th component: + + >>> proj_adj = odl.ComponentProjection(pspace, 0).adjoint + >>> proj_adj(x[0]) + ProductSpace(rn(1), rn(2), rn(3)).element([ + [ 1.], + [ 0., 0.], + [ 0., 0., 0.] + ]) + + Projection on a subspace corresponding to indices 0 and 2: + + >>> proj_adj = odl.ComponentProjection(pspace, [0, 2]).adjoint + >>> proj_adj(x[[0, 2]]) + ProductSpace(rn(1), rn(2), rn(3)).element([ + [ 1.], + [ 0., 0.], + [ 4., 5., 6.] + ]) + """ + super(ComponentProjectionAdjoint, self).__init__( + op.range, op.domain, linear=True) + + @property + def index(self): + """Index, indices or slice defining the subspace.""" + return op.index + + def _call(self, x, out=None): + """Extend ``x`` to the full space.""" + if out is None: + out = self.range.zero() + else: + out.set_zero() - def _call(self, x, out=None): - """Extend ``x`` from the subspace.""" - if out is None: - out = self.range.zero() - else: - out.set_zero() + out[self.index] = x + return out - out[self.index] = x - return out + def __repr__(self): + """Return ``repr(self)``. - @property - def adjoint(self): - """Adjoint of this operator. + Examples + -------- + >>> pspace = odl.ProductSpace(odl.rn(1), odl.rn(2)) + >>> odl.ComponentProjection(pspace, 0).adjoint + ComponentProjection(ProductSpace(rn(1), rn(2)), 0).adjoint + """ + return attribute_repr_string(repr(op), 'adjoint') - Returns - ------- - adjoint : `ComponentProjection` - The adjoint is given by the `ComponentProjection` related to this - operator's `index`. - """ - return ComponentProjection(self.range, self.index) + return ComponentProjectionAdjoint() def __repr__(self): """Return ``repr(self)``. @@ -633,11 +683,12 @@ def __repr__(self): Examples -------- >>> pspace = odl.ProductSpace(odl.rn(1), odl.rn(2)) - >>> odl.ComponentProjectionAdjoint(pspace, 0) - ComponentProjectionAdjoint(ProductSpace(rn(1), rn(2)), 0) + >>> odl.ComponentProjection(pspace, 0) + ComponentProjection(ProductSpace(rn(1), rn(2)), 0) """ - return '{}({!r}, {})'.format(self.__class__.__name__, - self.range, self.index) + posargs = [self.domain, self.index] + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts) class BroadcastOperator(Operator): @@ -679,10 +730,12 @@ def __init__(self, *operators): >>> x = [1, 2, 3] >>> op(x) - ProductSpace(rn(3), 2).element([ - [ 1., 2., 3.], - [ 2., 4., 6.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 1., 2., 3.], + [ 2., 4., 6.] + ] + ) Can also initialize by calling an operator repeatedly: @@ -736,12 +789,12 @@ def derivative(self, x): Parameters ---------- x : `domain` element - The point to take the derivative in + The point in which to take the derivative. Returns ------- - adjoint : linear `BroadcastOperator` - The derivative + derivative : linear `BroadcastOperator` + The derivative operator. Examples -------- @@ -751,22 +804,26 @@ def derivative(self, x): >>> residual_op = I - I.domain.element([1, 1, 1]) >>> op = BroadcastOperator(residual_op, 2 * residual_op) - Calling operator offsets by ``[1, 1, 1]``: + Calling the affine operator: >>> x = [1, 2, 3] >>> op(x) - ProductSpace(rn(3), 2).element([ - [ 0., 1., 2.], - [ 0., 2., 4.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 0., 1., 2.], + [ 0., 2., 4.] + ] + ) The derivative of this affine operator does not have an offset: >>> op.derivative(x)(x) - ProductSpace(rn(3), 2).element([ - [ 1., 2., 3.], - [ 2., 4., 6.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 1., 2., 3.], + [ 2., 4., 6.] + ] + ) """ return BroadcastOperator(*[op.derivative(x) for op in self.operators]) @@ -794,19 +851,24 @@ def __repr__(self): Examples -------- >>> spc = odl.rn(3) - >>> id = odl.IdentityOperator(spc) - >>> odl.BroadcastOperator(id, 3) + >>> ident = odl.IdentityOperator(spc) + >>> odl.BroadcastOperator(ident, 3) BroadcastOperator(IdentityOperator(rn(3)), 3) - >>> scale = odl.ScalingOperator(spc, 3) - >>> odl.BroadcastOperator(id, scale) - BroadcastOperator(IdentityOperator(rn(3)), ScalingOperator(rn(3), 3.0)) + >>> scaling = odl.ScalingOperator(spc, 3) + >>> odl.BroadcastOperator(ident, scaling) + BroadcastOperator( + IdentityOperator(rn(3)), + ScalingOperator(rn(3), scalar=3.0) + ) """ - if all(op == self[0] for op in self): - return '{}({!r}, {})'.format(self.__class__.__name__, - self[0], len(self)) + if all(op is self[0] for op in self): + posargs = [self[0], len(self)] else: - op_repr = ', '.join(repr(op) for op in self) - return '{}({})'.format(self.__class__.__name__, op_repr) + posargs = self.operators + + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) class ReductionOperator(Operator): @@ -968,10 +1030,12 @@ def adjoint(self): >>> I = odl.IdentityOperator(odl.rn(3)) >>> op = ReductionOperator(I, 2 * I) >>> op.adjoint([1, 2, 3]) - ProductSpace(rn(3), 2).element([ - [ 1., 2., 3.], - [ 2., 4., 6.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 1., 2., 3.], + [ 2., 4., 6.] + ] + ) """ return BroadcastOperator(*[op.adjoint for op in self.operators]) @@ -981,19 +1045,24 @@ def __repr__(self): Examples -------- >>> spc = odl.rn(3) - >>> id = odl.IdentityOperator(spc) - >>> odl.ReductionOperator(id, 3) + >>> I = odl.IdentityOperator(spc) + >>> odl.ReductionOperator(I, 3) ReductionOperator(IdentityOperator(rn(3)), 3) - >>> scale = odl.ScalingOperator(spc, 3) - >>> odl.ReductionOperator(id, scale) - ReductionOperator(IdentityOperator(rn(3)), ScalingOperator(rn(3), 3.0)) + >>> S = odl.ScalingOperator(spc, 3) + >>> odl.ReductionOperator(I, S) + ReductionOperator( + IdentityOperator(rn(3)), + ScalingOperator(rn(3), scalar=3.0) + ) """ - if all(op == self[0] for op in self): - return '{}({!r}, {})'.format(self.__class__.__name__, - self[0], len(self)) + if all(op is self[0] for op in self): + posargs = [self[0], len(self)] else: - op_repr = ', '.join(repr(op) for op in self) - return '{}({})'.format(self.__class__.__name__, op_repr) + posargs = self.operators + + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) class DiagonalOperator(ProductSpaceOperator): @@ -1043,12 +1112,15 @@ def __init__(self, *operators, **kwargs): >>> op([[1, 2, 3], ... [4, 5, 6]]) - ProductSpace(rn(3), 2).element([ - [ 1., 2., 3.], - [ 8., 10., 12.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 1., 2., 3.], + [ 8., 10., 12.] + ] + ) - Can also be created using a multiple of a single operator + The diagonal operator can also be created using a multiple of a + single operator: >>> op = DiagonalOperator(I, 2) >>> op.operators @@ -1186,19 +1258,24 @@ def __repr__(self): Examples -------- >>> spc = odl.rn(3) - >>> id = odl.IdentityOperator(spc) - >>> odl.DiagonalOperator(id, 3) + >>> I = odl.IdentityOperator(spc) + >>> odl.DiagonalOperator(I, 3) DiagonalOperator(IdentityOperator(rn(3)), 3) - >>> scale = odl.ScalingOperator(spc, 3) - >>> odl.DiagonalOperator(id, scale) - DiagonalOperator(IdentityOperator(rn(3)), ScalingOperator(rn(3), 3.0)) + >>> S = odl.ScalingOperator(spc, 3) + >>> odl.DiagonalOperator(I, S) + DiagonalOperator( + IdentityOperator(rn(3)), + ScalingOperator(rn(3), scalar=3.0) + ) """ - if all(op == self[0] for op in self): - return '{}({!r}, {})'.format(self.__class__.__name__, - self[0], len(self)) + if all(op is self[0] for op in self): + posargs = [self[0], len(self)] else: - op_repr = ', '.join(repr(op) for op in self) - return '{}({})'.format(self.__class__.__name__, op_repr) + posargs = self.operators + + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) if __name__ == '__main__': diff --git a/odl/operator/tensor_ops.py b/odl/operator/tensor_ops.py index 0e909dcb2b0..fc1122bf37b 100644 --- a/odl/operator/tensor_ops.py +++ b/odl/operator/tensor_ops.py @@ -8,18 +8,20 @@ """Operators defined for tensor fields.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + from numbers import Integral + import numpy as np from odl.operator.operator import Operator -from odl.set import RealNumbers, ComplexNumbers +from odl.set import ComplexNumbers, RealNumbers from odl.space import ProductSpace, tensor_space from odl.space.base_tensors import TensorSpace from odl.space.weighting import ArrayWeighting from odl.util import ( - signature_string, indent, dtype_repr, moveaxis, writable_array) - + REPR_PRECISION, array_str, attribute_repr_string, dtype_repr, moveaxis, + npy_printoptions, repr_string, signature_string_parts, writable_array) __all__ = ('PointwiseNorm', 'PointwiseInner', 'PointwiseSum', 'MatrixOperator', 'SamplingOperator', 'WeightedSumSamplingOperator', @@ -109,7 +111,7 @@ class PointwiseNorm(PointwiseTensorFieldOperator): ``d``. """ - def __init__(self, vfspace, exponent=None, weighting=None): + def __init__(self, vfspace, range=None, exponent=None, weighting=None): """Initialize a new instance. Parameters @@ -118,6 +120,10 @@ def __init__(self, vfspace, exponent=None, weighting=None): Space of vector fields on which the operator acts. It has to be a product space of identical spaces, i.e. a power space. + range : `TensorSpace`, optional + Scalar space to which the operator maps, must be compatible + with ``vfspace``'s shape. + Default: ``vfspace[0]`` exponent : non-zero float, optional Exponent of the norm in each point. Values between 0 and 1 are currently not supported due to numerical @@ -138,7 +144,7 @@ def __init__(self, vfspace, exponent=None, weighting=None): maps a vector field to a scalar function: >>> spc = odl.uniform_discr([-1, -1], [1, 1], (1, 2)) - >>> vfspace = odl.ProductSpace(spc, 2) + >>> vfspace = spc ** 2 >>> pw_norm = odl.PointwiseNorm(vfspace) >>> pw_norm.range == spc True @@ -147,27 +153,29 @@ def __init__(self, vfspace, exponent=None, weighting=None): >>> x = vfspace.element([[[1, -4]], ... [[0, 3]]]) - >>> print(pw_norm(x)) - [[ 1., 5.]] + >>> pw_norm(x) + uniform_discr([-1., -1.], [ 1., 1.], (1, 2)).element([[ 1., 5.]]) We can change the exponent either in the vector field space or in the operator directly: >>> vfspace = odl.ProductSpace(spc, 2, exponent=1) - >>> pw_norm = PointwiseNorm(vfspace) - >>> print(pw_norm(x)) - [[ 1., 7.]] + >>> pw_norm = odl.PointwiseNorm(vfspace) + >>> pw_norm(x) + uniform_discr([-1., -1.], [ 1., 1.], (1, 2)).element([[ 1., 7.]]) >>> vfspace = odl.ProductSpace(spc, 2) >>> pw_norm = PointwiseNorm(vfspace, exponent=1) - >>> print(pw_norm(x)) - [[ 1., 7.]] + >>> pw_norm(x) + uniform_discr([-1., -1.], [ 1., 1.], (1, 2)).element([[ 1., 7.]]) """ + # TODO: allow `Tensorspace` and `DiscreteLp` with shaped dtype if not isinstance(vfspace, ProductSpace): raise TypeError('`vfspace` {!r} is not a ProductSpace ' 'instance'.format(vfspace)) + if range is None: + range = vfspace[0] super(PointwiseNorm, self).__init__( - domain=vfspace, range=vfspace[0], base_space=vfspace[0], - linear=False) + domain=vfspace, range=range, base_space=vfspace[0], linear=False) # Need to check for product space shape once higher order tensors # are implemented @@ -176,13 +184,13 @@ def __init__(self, vfspace, exponent=None, weighting=None): if self.domain.exponent is None: raise ValueError('cannot determine `exponent` from {}' ''.format(self.domain)) - self._exponent = self.domain.exponent + self.__exponent = self.domain.exponent elif exponent < 1: raise ValueError('`exponent` smaller than 1 not allowed') else: - self._exponent = float(exponent) + self.__exponent = float(exponent) - # Handle weighting, including sanity checks + # Weighting checks if weighting is None: # TODO: find a more robust way of getting the weights as an array if hasattr(self.domain.weighting, 'array'): @@ -210,7 +218,7 @@ def __init__(self, vfspace, exponent=None, weighting=None): @property def exponent(self): """Exponent ``p`` of this norm.""" - return self._exponent + return self.__exponent @property def weights(self): @@ -353,25 +361,55 @@ def derivative(self, vf): return PointwiseInner(self.domain, inner_vf, weighting=self.weights) + def __repr__(self): + """Return ``repr(self)``. + + Examples + -------- + >>> spc = odl.uniform_discr([-1, -1], [1, 1], (1, 2)) + >>> vfspace = spc ** 2 + >>> pw_norm = odl.PointwiseNorm(vfspace) + >>> pw_norm + PointwiseNorm( + ProductSpace(uniform_discr([-1., -1.], [ 1., 1.], (1, 2)), 2) + ) + """ + posargs = [self.domain] + optargs = [('range', self.range, self.domain[0]), + ('exponent', self.exponent, self.domain.exponent)] + optmod = ['!r', '!r'] + if self.is_weighted: + optargs.append(('weighting', array_str(self.weights), '')) + optmod.append('!s') + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs, + mod=['!r', optmod]) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) + + +class PointwiseInner(PointwiseTensorFieldOperator): -class PointwiseInnerBase(PointwiseTensorFieldOperator): - """Base class for `PointwiseInner` and `PointwiseInnerAdjoint`. + """Take the point-wise inner product with a given vector field. + + This operator takes the (weighted) inner product :: + + P[F](x) = = sum_j ( w_j * F_j(x) * conj(G_j(x)) ) + + for a given vector field ``G``, where ``F`` is the vector field + acting as a variable to this operator. - Implemented to allow code reuse between the classes. + This implies that the `Operator.domain` is a power space of a + discretized function space. For example, if ``X`` is a `DiscreteLp` + space, then ``ProductSpace(X, d)`` is a valid domain for any + positive integer ``d``. """ - def __init__(self, adjoint, vfspace, vecfield, weighting=None): + def __init__(self, vfspace, vecfield, weighting=None): """Initialize a new instance. - All parameters are given according to the specifics of the "usual" - operator. The ``adjoint`` parameter is used to control conversions - for the inverse transform. - Parameters ---------- - adjoint : bool - ``True`` if the operator should be the adjoint, ``False`` - otherwise. vfspace : `ProductSpace` Space of vector fields on which the operator acts. It has to be a product space of identical spaces, i.e. a @@ -381,22 +419,39 @@ def __init__(self, adjoint, vfspace, vecfield, weighting=None): product of an input vector field weighting : `array-like` or float, optional Weighting array or constant for the norm. If an array is - given, its length must be equal to ``len(domain)``. + given, its length must be equal to ``len(domain)``, and + all entries must be positive. By default, the weights are is taken from ``domain.weighting``. Note that this excludes unusual weightings with custom inner product, norm or dist. + + Examples + -------- + We make a tiny vector field space in 2D and create the + point-wise inner product operator with a fixed vector field. + The operator maps a vector field to a scalar function: + + >>> spc = odl.uniform_discr([-1, -1], [1, 1], (1, 2)) + >>> vfspace = spc ** 2 + >>> fixed_vf = np.array([[[0, 1]], + ... [[1, -1]]]) + >>> pw_inner = odl.PointwiseInner(vfspace, fixed_vf) + >>> pw_inner.range == spc + True + + Now we can calculate the inner product in each point: + + >>> x = vfspace.element([[[1, -4]], + ... [[0, 3]]]) + >>> pw_inner(x) + uniform_discr([-1., -1.], [ 1., 1.], (1, 2)).element([[ 0., -7.]]) """ if not isinstance(vfspace, ProductSpace): raise TypeError('`vfsoace` {!r} is not a ProductSpace ' 'instance'.format(vfspace)) - if adjoint: - super(PointwiseInnerBase, self).__init__( - domain=vfspace[0], range=vfspace, base_space=vfspace[0], - linear=True) - else: - super(PointwiseInnerBase, self).__init__( - domain=vfspace, range=vfspace[0], base_space=vfspace[0], - linear=True) + # TODO: custom range + super(PointwiseInner, self).__init__( + vfspace, vfspace[0], base_space=vfspace[0], linear=True) # Bail out if the space is complex but we cannot take the complex # conjugate. @@ -407,7 +462,7 @@ def __init__(self, adjoint, vfspace, vecfield, weighting=None): 'method required for complex inner products' ''.format(self.base_space.element_type)) - self._vecfield = vfspace.element(vecfield) + self.__vecfield = vfspace.element(vecfield) # Handle weighting, including sanity checks if weighting is None: @@ -429,7 +484,7 @@ def __init__(self, adjoint, vfspace, vecfield, weighting=None): @property def vecfield(self): """Fixed vector field ``G`` of this inner product.""" - return self._vecfield + return self.__vecfield @property def weights(self): @@ -441,85 +496,12 @@ def is_weighted(self): """``True`` if weighting is not 1 or all ones.""" return self.__is_weighted - @property - def adjoint(self): - """Adjoint operator.""" - raise NotImplementedError('abstract method') - - -class PointwiseInner(PointwiseInnerBase): - - """Take the point-wise inner product with a given vector field. - - This operator takes the (weighted) inner product :: - - = sum_j ( w_j * F_j(x) * conj(G_j(x)) ) - - for a given vector field ``G``, where ``F`` is the vector field - acting as a variable to this operator. - - This implies that the `Operator.domain` is a power space of a - discretized function space. For example, if ``X`` is a `DiscreteLp` - space, then ``ProductSpace(X, d)`` is a valid domain for any - positive integer ``d``. - """ - - def __init__(self, vfspace, vecfield, weighting=None): - """Initialize a new instance. - - Parameters - ---------- - vfspace : `ProductSpace` - Space of vector fields on which the operator acts. - It has to be a product space of identical spaces, i.e. a - power space. - vecfield : ``vfspace`` `element-like` - Vector field with which to calculate the point-wise inner - product of an input vector field - weighting : `array-like` or float, optional - Weighting array or constant for the norm. If an array is - given, its length must be equal to ``len(domain)``, and - all entries must be positive. - By default, the weights are is taken from - ``domain.weighting``. Note that this excludes unusual - weightings with custom inner product, norm or dist. - - Examples - -------- - We make a tiny vector field space in 2D and create the - point-wise inner product operator with a fixed vector field. - The operator maps a vector field to a scalar function: - - >>> spc = odl.uniform_discr([-1, -1], [1, 1], (1, 2)) - >>> vfspace = odl.ProductSpace(spc, 2) - >>> fixed_vf = np.array([[[0, 1]], - ... [[1, -1]]]) - >>> pw_inner = PointwiseInner(vfspace, fixed_vf) - >>> pw_inner.range == spc - True - - Now we can calculate the inner product in each point: - - >>> x = vfspace.element([[[1, -4]], - ... [[0, 3]]]) - >>> print(pw_inner(x)) - [[ 0., -7.]] - """ - super(PointwiseInner, self).__init__( - adjoint=False, vfspace=vfspace, vecfield=vecfield, - weighting=weighting) - - @property - def vecfield(self): - """Fixed vector field ``G`` of this inner product.""" - return self._vecfield - def _call(self, vf, out): """Implement ``self(vf, out)``.""" if self.domain.field == ComplexNumbers(): - vf[0].multiply(self._vecfield[0].conj(), out=out) + vf[0].multiply(self.vecfield[0].conj(), out=out) else: - vf[0].multiply(self._vecfield[0], out=out) + vf[0].multiply(self.vecfield[0], out=out) if self.is_weighted: out *= self.weights[0] @@ -544,101 +526,90 @@ def _call(self, vf, out): def adjoint(self): """Adjoint of this operator. - Returns - ------- - adjoint : `PointwiseInnerAdjoint` - """ - return PointwiseInnerAdjoint( - sspace=self.base_space, vecfield=self.vecfield, - vfspace=self.domain, weighting=self.weights) - - -class PointwiseInnerAdjoint(PointwiseInnerBase): + The adjoint operator is given by mapping a scalar-valued function + ``h`` and multiplying it with each component of the vector field + ``G``, scaled inversely by the weights:: - """Adjoint of the point-wise inner product operator. + P^*[h]_j(x) = h(x) * G_j(x) / w_j - The adjoint of the inner product operator is a mapping :: - - A^* : X --> X^d - - If the vector field space ``X^d`` is weighted by a vector ``v``, - the adjoint, applied to a function ``h`` from ``X`` is the vector - field :: + Examples + -------- + >>> spc = odl.uniform_discr([-1, -1], [1, 1], (1, 2)) + >>> vfspace = spc ** 2 + >>> fixed_vf = np.array([[[0, 1]], + ... [[1, -1]]]) + >>> pw_inner = odl.PointwiseInner(vfspace, fixed_vf) + >>> y = spc.element([[1, 2]]) + >>> pw_inner.adjoint(y) + ProductSpace(uniform_discr([-1., -1.], [ 1., 1.], (1, 2)), 2).element( + [ + [[ 0., 2.]], + [[ 1., -2.]] + ] + ) + """ + op = self - x --> h(x) * (w / v) * G(x), + class PointwiseInnerAdjoint(PointwiseTensorFieldOperator): - where ``G`` and ``w`` are the vector field and weighting from the - inner product operator, resp., and all multiplications are understood - component-wise. - """ + """Adjoint of the point-wise inner product operator.""" - def __init__(self, sspace, vecfield, vfspace=None, weighting=None): - """Initialize a new instance. + def __init__(self): + """Initialize a new instance.""" + super(PointwiseInnerAdjoint, self).__init__( + domain=op.range, range=op.domain, + base_space=op.base_space, linear=True) + + def _call(self, f, out): + """Implement ``self(vf, out)``.""" + for vfi, oi, wi in zip(op.vecfield, out, op.weights): + vfi.multiply(f, out=oi) + if not np.isclose(wi, 1.0): + oi *= wi + + @property + def adjoint(self): + """Adjoint of this operator. - Parameters - ---------- - sspace : `LinearSpace` - "Scalar" space on which the operator acts - vecfield : range `element-like` - Vector field of the point-wise inner product operator - vfspace : `ProductSpace`, optional - Space of vector fields to which the operator maps. It must - be a power space with ``sspace`` as base space. - This option is intended to enforce an operator range - with a certain weighting. - Default: ``ProductSpace(space, len(vecfield), - weighting=weighting)`` - weighting : `array-like` or float, optional - Weighting array or constant of the inner product operator. - If an array is given, its length must be equal to - ``len(vecfield)``. - By default, the weights are is taken from - ``range.weighting`` if applicable. Note that this excludes - unusual weightings with custom inner product, norm or dist. - """ - if vfspace is None: - vfspace = ProductSpace(sspace, len(vecfield), weighting=weighting) - else: - if not isinstance(vfspace, ProductSpace): - raise TypeError('`vfspace` {!r} is not a ' - 'ProductSpace instance'.format(vfspace)) - if vfspace[0] != sspace: - raise ValueError('base space of the range is different from ' - 'the given scalar space ({!r} != {!r})' - ''.format(vfspace[0], sspace)) - super(PointwiseInnerAdjoint, self).__init__( - adjoint=True, vfspace=vfspace, vecfield=vecfield, - weighting=weighting) - - # Get weighting from range - if hasattr(self.range.weighting, 'array'): - self.__ran_weights = self.range.weighting.array - elif hasattr(self.range.weighting, 'const'): - self.__ran_weights = (self.range.weighting.const * - np.ones(len(self.range))) - else: - raise ValueError('weighting scheme {!r} of the range does ' - 'not define a weighting array or constant' - ''.format(self.range.weighting)) + Returns + ------- + adjoint : `PointwiseInner` + """ + return op - def _call(self, f, out): - """Implement ``self(vf, out)``.""" - for vfi, oi, ran_wi, dom_wi in zip(self.vecfield, out, - self.__ran_weights, self.weights): - vfi.multiply(f, out=oi) - if not np.isclose(ran_wi, dom_wi): - oi *= dom_wi / ran_wi + return PointwiseInnerAdjoint() - @property - def adjoint(self): - """Adjoint of this operator. + def __repr__(self): + """Return ``repr(self)``. - Returns - ------- - adjoint : `PointwiseInner` + Examples + -------- + >>> spc = odl.rn((1, 2)) + >>> vfspace = spc ** 2 + >>> fixed_vf = np.array([[[0, 1]], + ... [[1, -1]]]) + >>> pw_inner = odl.PointwiseInner(vfspace, fixed_vf) + >>> pw_inner + PointwiseInner( + ProductSpace(rn((1, 2)), 2), + ProductSpace(rn((1, 2)), 2).element( + [ + [[ 0., 1.]], + [[ 1., -1.]] + ] + ) + ) """ - return PointwiseInner(vfspace=self.range, vecfield=self.vecfield, - weighting=self.weights) + posargs = [self.domain, self.vecfield] + optargs = [] + optmod = [] + if self.is_weighted: + optargs.append(('weighting', array_str(self.weights), '')) + optmod.append('!s') + inner_parts = signature_string_parts(posargs, optargs, + mod=['!r', optmod]) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) # TODO: Make this an optimized operator on its own. @@ -681,7 +652,7 @@ def __init__(self, vfspace, weighting=None): >>> spc = odl.uniform_discr([-1, -1], [1, 1], (1, 2)) >>> vfspace = odl.ProductSpace(spc, 2) - >>> pw_sum = PointwiseSum(vfspace) + >>> pw_sum = odl.PointwiseSum(vfspace) >>> pw_sum.range == spc True @@ -696,9 +667,32 @@ def __init__(self, vfspace, weighting=None): raise TypeError('`vfspace` {!r} is not a ProductSpace ' 'instance'.format(vfspace)) - ones = vfspace.one() super(PointwiseSum, self).__init__( - vfspace, vecfield=ones, weighting=weighting) + vfspace, vecfield=vfspace.one(), weighting=weighting) + + def __repr__(self): + """Return ``repr(self)``. + + Examples + -------- + >>> spc = odl.uniform_discr([-1, -1], [1, 1], (1, 2)) + >>> vfspace = odl.ProductSpace(spc, 2) + >>> pw_sum = odl.PointwiseSum(vfspace) + >>> pw_sum + PointwiseSum( + ProductSpace(uniform_discr([-1., -1.], [ 1., 1.], (1, 2)), 2) + ) + """ + posargs = [self.domain] + optargs = [] + optmod = [] + if self.is_weighted: + optargs.append(('weighting', array_str(self.weights), '')) + optmod.append('!s') + inner_parts = signature_string_parts(posargs, optargs, + mod=['!r', optmod]) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) class MatrixOperator(Operator): @@ -741,7 +735,7 @@ def __init__(self, matrix, domain=None, range=None, axis=0): By default, ``domain`` and ``range`` are spaces of with one axis: >>> m = np.ones((3, 4)) - >>> op = MatrixOperator(m) + >>> op = odl.MatrixOperator(m) >>> op.domain rn(4) >>> op.range @@ -755,7 +749,7 @@ def __init__(self, matrix, domain=None, range=None, axis=0): domain shape entry in the given axis: >>> dom = odl.rn((5, 4, 4)) # can use axis=1 or axis=2 - >>> op = MatrixOperator(m, domain=dom, axis=1) + >>> op = odl.MatrixOperator(m, domain=dom, axis=1) >>> op(dom.one()).shape (5, 3, 4) >>> op = MatrixOperator(m, domain=dom, axis=2) @@ -941,7 +935,22 @@ def _call(self, x, out=None): return out def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> m = np.ones((3, 4)) + >>> dom = odl.rn((5, 4, 4)) + >>> op = odl.MatrixOperator(m, domain=dom, axis=1) + >>> op + MatrixOperator( + [[ 1., 1., 1., 1.], + [ 1., 1., 1., 1.], + [ 1., 1., 1., 1.]], + domain=rn((5, 4, 4)), + axis=1 + ) + """ # Lazy import to improve `import odl` time import scipy.sparse @@ -964,9 +973,10 @@ def __repr__(self): ('axis', self.axis, 0) ] - inner_str = signature_string(posargs, optargs, sep=[', ', ', ', ',\n'], - mod=[['!s'], ['!r', '!r', '']]) - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) + inner_parts = signature_string_parts(posargs, optargs, + mod=[['!s'], ['!r', '!r', '']]) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): """Return ``str(self)``.""" @@ -1193,12 +1203,25 @@ def __repr__(self): """Return ``repr(self)``.""" posargs = [self.domain, self.sampling_points] optargs = [('variant', self.variant, 'point_eval')] - sig_str = signature_string(posargs, optargs, mod=['!r', ''], - sep=[',\n', '', ',\n']) - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(sig_str)) + inner_parts = signature_string_parts(posargs, optargs, mod=['!r', '']) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): - """Return ``str(self)``.""" + """Return ``str(self)``. + + Examples + -------- + >>> space = odl.uniform_discr(0, 1, 4) + >>> op = odl.SamplingOperator(space, sampling_points=[1, 2, 1], + ... variant='integrate') + >>> op + SamplingOperator( + uniform_discr(0.0, 1.0, 4), + [array([1, 2, 1])], + variant='integrate' + ) + """ return repr(self) @@ -1388,12 +1411,26 @@ def adjoint(self): return SamplingOperator(self.range, self.sampling_points, variant) def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> space = odl.uniform_discr(0, 1, 4) + >>> op = odl.WeightedSumSamplingOperator(space, sampling_points=1) + >>> op = odl.WeightedSumSamplingOperator( + ... space, sampling_points=[1, 2, 1], variant='dirac') + >>> op + WeightedSumSamplingOperator( + uniform_discr(0.0, 1.0, 4), + [array([1, 2, 1])], + variant='dirac' + ) + """ posargs = [self.range, self.sampling_points] optargs = [('variant', self.variant, 'char_fun')] - sig_str = signature_string(posargs, optargs, mod=['!r', ''], - sep=[',\n', '', ',\n']) - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(sig_str)) + inner_parts = signature_string_parts(posargs, optargs, mod=['!r', '']) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): """Return ``str(self)``.""" @@ -1479,8 +1516,7 @@ def adjoint(self): >>> abs(op.adjoint(op(x)).inner(x) - op(x).inner(op(x))) < 1e-10 True """ - scaling = getattr(self.domain, 'cell_volume', 1.0) - return 1 / scaling * self.inverse + return self.inverse @property def inverse(self): @@ -1506,7 +1542,6 @@ def inverse(self): True """ op = self - scaling = getattr(self.domain, 'cell_volume', 1.0) class FlatteningOperatorInverse(Operator): @@ -1528,16 +1563,16 @@ def _call(self, x): order=op.order) def adjoint(self): - """Adjoint of this operator, a scaled `FlatteningOperator`.""" - return scaling * op + """Adjoint of this operator, the original operator.""" + return op def inverse(self): - """Inverse of this operator.""" + """Inverse of this operator, the original operator.""" return op def __repr__(self): """Return ``repr(self)``.""" - return '{!r}.inverse'.format(op) + return attribute_repr_string(repr(op), 'inverse') def __str__(self): """Return ``str(self)``.""" @@ -1546,12 +1581,20 @@ def __str__(self): return FlatteningOperatorInverse() def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> space = odl.uniform_discr([-1, -1], [1, 1], shape=(2, 3)) + >>> op = odl.FlatteningOperator(space) + >>> op + FlatteningOperator(uniform_discr([-1., -1.], [ 1., 1.], (2, 3))) + """ posargs = [self.domain] optargs = [('order', self.order, 'C')] - sig_str = signature_string(posargs, optargs, mod=['!r', ''], - sep=['', '', ',\n']) - return '{}(\n{}\n)'.format(self.__class__.__name__, indent(sig_str)) + inner_parts = signature_string_parts(posargs, optargs, mod=['!r', '']) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def __str__(self): """Return ``str(self)``.""" @@ -1617,12 +1660,18 @@ def is_compatible_space(space, base_space): return is_compatible_space(space[0], base_space) else: if hasattr(space, 'astype') and hasattr(base_space, 'dtype'): - # TODO: maybe only the shape should play a role? comp_space = space.astype(base_space.dtype) + if (hasattr(comp_space, 'asweighted') and + hasattr(base_space, 'asweighted')): + comp_space = comp_space.asweighted(None) + comp_base = base_space.asweighted(None) + else: + comp_base = base_space else: comp_space = space + comp_base = base_space - return comp_space == base_space + return comp_space == comp_base if __name__ == '__main__': diff --git a/odl/phantom/noise.py b/odl/phantom/noise.py index 6a886cfd7ae..bf19d766235 100644 --- a/odl/phantom/noise.py +++ b/odl/phantom/noise.py @@ -72,7 +72,7 @@ def uniform_noise(space, low=0, high=1, seed=None): Parameters ---------- - space : `FnBase` or `ProductSpace` + space : `TensorSpace` or `ProductSpace` The space in which the noise is created. low : ``space.field`` element or ``space`` `element-like`, optional The lower bound of the uniform noise. If a scalar, it is interpreted as diff --git a/odl/set/domain.py b/odl/set/domain.py index 64ab5a17542..7e04e6f67d7 100644 --- a/odl/set/domain.py +++ b/odl/set/domain.py @@ -13,7 +13,8 @@ from odl.set.sets import Set from odl.util import ( - array_str, is_valid_input_array, is_valid_input_meshgrid, safe_int_conv) + array_str, is_valid_input_array, is_valid_input_meshgrid, safe_int_conv, + signature_string_parts, repr_string, npy_printoptions, REPR_PRECISION) __all__ = ('IntervalProd',) @@ -42,8 +43,8 @@ def __init__(self, min_pt, max_pt): Examples -------- >>> min_pt, max_pt = [-1, 2.5, 70], [-0.5, 10, 75] - >>> rbox = odl.IntervalProd(min_pt, max_pt) - >>> rbox + >>> intv = odl.IntervalProd(min_pt, max_pt) + >>> intv IntervalProd([ -1. , 2.5, 70. ], [ -0.5, 10. , 75. ]) """ super(IntervalProd, self).__init__() @@ -199,11 +200,11 @@ def approx_equals(self, other, atol): Examples -------- - >>> rbox1 = IntervalProd(0, 0.5) - >>> rbox2 = IntervalProd(0, np.sqrt(0.5)**2) - >>> rbox1.approx_equals(rbox2, atol=0) # Numerical error + >>> intv1 = IntervalProd(0, 0.5) + >>> intv2 = IntervalProd(0, np.sqrt(0.5)**2) + >>> intv1.approx_equals(intv2, atol=0) # Numerical error False - >>> rbox1.approx_equals(rbox2, atol=1e-15) + >>> intv1.approx_equals(intv2, atol=1e-15) True """ if other is self: @@ -244,11 +245,11 @@ def approx_contains(self, point, atol): Examples -------- >>> min_pt, max_pt = [-1, 0, 2], [-0.5, 0, 3] - >>> rbox = IntervalProd(min_pt, max_pt) + >>> intv = IntervalProd(min_pt, max_pt) >>> # Numerical error - >>> rbox.approx_contains([-1 + np.sqrt(0.5)**2, 0., 2.9], atol=0) + >>> intv.approx_contains([-1 + np.sqrt(0.5)**2, 0., 2.9], atol=0) False - >>> rbox.approx_contains([-1 + np.sqrt(0.5)**2, 0., 2.9], atol=1e-9) + >>> intv.approx_contains([-1 + np.sqrt(0.5)**2, 0., 2.9], atol=1e-9) True """ try: @@ -306,12 +307,12 @@ def contains_set(self, other, atol=0.0): Examples -------- >>> min_pt1, max_pt1 = [-1, 0, 2], [-0.5, 0, 3] - >>> rbox1 = IntervalProd(min_pt1, max_pt1) + >>> intv1 = IntervalProd(min_pt1, max_pt1) >>> min_pt2, max_pt2 = [-0.6, 0, 2.1], [-0.5, 0, 2.5] - >>> rbox2 = IntervalProd(min_pt2, max_pt2) - >>> rbox1.contains_set(rbox2) + >>> intv2 = IntervalProd(min_pt2, max_pt2) + >>> intv1.contains_set(intv2) True - >>> rbox2.contains_set(rbox1) + >>> intv2.contains_set(intv1) False """ if self is other: @@ -345,13 +346,13 @@ def contains_all(self, other, atol=0.0): Examples -------- >>> min_pt, max_pt = [-1, 0, 2], [-0.5, 0, 3] - >>> rbox = IntervalProd(min_pt, max_pt) + >>> intv = IntervalProd(min_pt, max_pt) Arrays are expected in ``(ndim, npoints)`` shape: >>> arr = np.array([[-1, 0, 2], # defining one point at a time ... [-0.5, 0, 2]]) - >>> rbox.contains_all(arr.T) + >>> intv.contains_all(arr.T) True Implicit meshgrids defined by coordinate vectors: @@ -361,20 +362,20 @@ def contains_all(self, other, atol=0.0): >>> vec2 = (0, 0, 0) >>> vec3 = (2.5, 2.75, 3) >>> mg = sparse_meshgrid(vec1, vec2, vec3) - >>> rbox.contains_all(mg) + >>> intv.contains_all(mg) True Works also with an arbitrary iterable: - >>> rbox.contains_all([[-1, -0.5], # define points by axis + >>> intv.contains_all([[-1, -0.5], # define points by axis ... [0, 0], ... [2, 2]]) True Grids are also accepted as input: - >>> agrid = odl.uniform_grid(rbox.min_pt, rbox.max_pt, [3, 1, 3]) - >>> rbox.contains_all(agrid) + >>> agrid = odl.uniform_grid(intv.min_pt, intv.max_pt, [3, 1, 3]) + >>> intv.contains_all(agrid) True """ atol = float(atol) @@ -417,16 +418,16 @@ def measure(self, ndim=None): Examples -------- >>> min_pt, max_pt = [-1, 2.5, 0], [-0.5, 10, 0] - >>> rbox = IntervalProd(min_pt, max_pt) - >>> rbox.measure() + >>> intv = IntervalProd(min_pt, max_pt) + >>> intv.measure() 3.75 - >>> rbox.measure(ndim=3) + >>> intv.measure(ndim=3) 0.0 - >>> rbox.measure(ndim=3) == rbox.volume + >>> intv.measure(ndim=3) == intv.volume True - >>> rbox.measure(ndim=1) + >>> intv.measure(ndim=1) inf - >>> rbox.measure() == rbox.squeeze().volume + >>> intv.measure() == intv.squeeze().volume True """ if self.true_ndim == 0: @@ -466,10 +467,10 @@ def dist(self, point, exponent=2.0): Examples -------- >>> min_pt, max_pt = [-1, 0, 2], [-0.5, 0, 3] - >>> rbox = IntervalProd(min_pt, max_pt) - >>> rbox.dist([-5, 3, 2]) + >>> intv = IntervalProd(min_pt, max_pt) + >>> intv.dist([-5, 3, 2]) 5.0 - >>> rbox.dist([-5, 3, 2], exponent=float('inf')) + >>> intv.dist([-5, 3, 2], exponent=float('inf')) 4.0 """ point = np.atleast_1d(point) @@ -514,10 +515,10 @@ def collapse(self, indices, values): Examples -------- >>> min_pt, max_pt = [-1, 0, 2], [-0.5, 1, 3] - >>> rbox = IntervalProd(min_pt, max_pt) - >>> rbox.collapse(1, 0) + >>> intv = IntervalProd(min_pt, max_pt) + >>> intv.collapse(1, 0) IntervalProd([-1., 0., 2.], [-0.5, 0. , 3. ]) - >>> rbox.collapse([1, 2], [0, 2.5]) + >>> intv.collapse([1, 2], [0, 2.5]) IntervalProd([-1. , 0. , 2.5], [-0.5, 0. , 2.5]) """ indices = np.atleast_1d(indices).astype('int64', casting='safe') @@ -563,12 +564,12 @@ def squeeze(self): Examples -------- >>> min_pt, max_pt = [-1, 0, 2], [-0.5, 1, 3] - >>> rbox = IntervalProd(min_pt, max_pt) - >>> rbox.collapse(1, 0).squeeze() + >>> intv = IntervalProd(min_pt, max_pt) + >>> intv.collapse(1, 0).squeeze() IntervalProd([-1., 2.], [-0.5, 3. ]) - >>> rbox.collapse([1, 2], [0, 2.5]).squeeze() + >>> intv.collapse([1, 2], [0, 2.5]).squeeze() IntervalProd(-1.0, -0.5) - >>> rbox.collapse([0, 1, 2], [-1, 0, 2.5]).squeeze() + >>> intv.collapse([0, 1, 2], [-1, 0, 2.5]).squeeze() IntervalProd([], []) """ b_new = self.min_pt[self.nondegen_byaxis] @@ -742,23 +743,23 @@ def __getitem__(self, indices): Examples -------- - >>> rbox = IntervalProd([-1, 2, 0], [-0.5, 3, 0.5]) + >>> intv = IntervalProd([-1, 2, 0], [-0.5, 3, 0.5]) Indexing by integer selects single axes: - >>> rbox[0] + >>> intv[0] IntervalProd(-1.0, -0.5) With slices, multiple axes can be selected: - >>> rbox[:] + >>> intv[:] IntervalProd([-1., 2., 0.], [-0.5, 3. , 0.5]) - >>> rbox[::2] + >>> intv[::2] IntervalProd([-1., 0.], [-0.5, 0.5]) A list of integers can be used for free combinations of axes: - >>> rbox[[0, 1, 0]] + >>> intv[[0, 1, 0]] IntervalProd([-1., 2., -1.], [-0.5, 3. , -0.5]) """ return IntervalProd(self.min_pt[indices], self.max_pt[indices]) @@ -836,14 +837,30 @@ def __rdiv__(self, other): __rtruediv__ = __rdiv__ def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> intv = odl.IntervalProd(0, 1) + >>> intv + IntervalProd(0.0, 1.0) + >>> intv = odl.IntervalProd([0] * 10, [1] * 10) + >>> intv + IntervalProd( + [ 0., 0., 0., ..., 0., 0., 0.], + [ 1., 1., 1., ..., 1., 1., 1.] + ) + """ if self.ndim == 1: - return '{}({:.4}, {:.4})'.format(self.__class__.__name__, - self.min_pt[0], self.max_pt[0]) + posargs = [self.min_pt[0], self.max_pt[0]] + posmod = '' else: - return '{}({}, {})'.format(self.__class__.__name__, - array_str(self.min_pt), - array_str(self.max_pt)) + posargs = [self.min_pt, self.max_pt] + posmod = array_str + + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, [], mod=[posmod, '']) + return repr_string(self.__class__.__name__, inner_parts) def __str__(self): """Return ``str(self)``.""" diff --git a/odl/set/sets.py b/odl/set/sets.py index ad9287afdeb..03e46b9182f 100644 --- a/odl/set/sets.py +++ b/odl/set/sets.py @@ -14,7 +14,9 @@ from past.builtins import basestring import numpy as np -from odl.util import is_int_dtype, is_real_dtype, is_numeric_dtype, unique +from odl.util import ( + is_int_dtype, is_real_dtype, is_numeric_dtype, unique, + signature_string_parts, repr_string) __all__ = ('Set', 'EmptySet', 'UniversalSet', 'Field', 'Integers', @@ -283,8 +285,17 @@ def examples(self): return [('hello', hello_str), ('world', world_str)] def __repr__(self): - """Return ``repr(self)``.""" - return 'Strings({})'.format(self.length) + """Return ``repr(self)``. + + Examples + -------- + >>> s8 = odl.Strings(8) + >>> s8 + Strings(8) + """ + posargs = [self.length] + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts) class Field(Set): @@ -348,6 +359,16 @@ def contains_all(self, other): dtype = np.result_type(*other) return is_numeric_dtype(dtype) + @property + def real_space(self): + """``RealNumbers``, for compatibility with tensor spaces.""" + return RealNumbers() + + @property + def complex_space(self): + """``ComplexNumbers``, for compatibility with tensor spaces.""" + return ComplexNumbers() + def __eq__(self, other): """Return ``self == other``.""" if other is self: @@ -415,6 +436,16 @@ def contains_all(self, array): dtype = np.result_type(*array) return is_real_dtype(dtype) + @property + def real_space(self): + """``RealNumbers``, for compatibility with tensor spaces.""" + return RealNumbers() + + @property + def complex_space(self): + """``ComplexNumbers``, for compatibility with tensor spaces.""" + return ComplexNumbers() + def __eq__(self, other): """Return ``self == other``.""" if other is self: @@ -613,9 +644,20 @@ def __str__(self): return ' x '.join(str(set_) for set_ in self.sets) def __repr__(self): - """Return ``repr(self)``.""" - sets_str = ', '.join(repr(set_) for set_ in self.sets) - return '{}({})'.format(self.__class__.__name__, sets_str) + """Return ``repr(self)``. + + Examples + -------- + >>> rs, cs = odl.RealNumbers(), odl.ComplexNumbers() + >>> prod = odl.CartesianProduct(rs, cs, cs, rs) + >>> prod + CartesianProduct( + RealNumbers(), ComplexNumbers(), ComplexNumbers(), RealNumbers() + ) + """ + posargs = self.sets + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts) class SetUnion(Set): @@ -733,12 +775,14 @@ def __repr__(self): Examples -------- - >>> reals, complexnrs = odl.RealNumbers(), odl.ComplexNumbers() - >>> odl.SetUnion(reals, complexnrs) + >>> rs, cs = odl.RealNumbers(), odl.ComplexNumbers() + >>> union = odl.SetUnion(rs, cs) + >>> union SetUnion(RealNumbers(), ComplexNumbers()) """ - sets_str = ', '.join(repr(set_) for set_ in self.sets) - return '{}({})'.format(self.__class__.__name__, sets_str) + posargs = self.sets + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts) class SetIntersection(Set): @@ -763,7 +807,7 @@ def __init__(self, *sets): Examples -------- >>> reals, complexnrs = odl.RealNumbers(), odl.ComplexNumbers() - >>> union = odl.SetIntersection(reals, complexnrs) + >>> intersec = odl.SetIntersection(reals, complexnrs) """ for set_ in sets: if not isinstance(set_, Set): @@ -844,8 +888,9 @@ def __repr__(self): >>> odl.SetIntersection(reals, complexnrs) SetIntersection(RealNumbers(), ComplexNumbers()) """ - sets_str = ', '.join(repr(set_) for set_ in self.sets) - return '{}({})'.format(self.__class__.__name__, sets_str) + posargs = self.sets + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts) class FiniteSet(Set): @@ -862,7 +907,7 @@ def __init__(self, *elements): Examples -------- - >>> set = odl.FiniteSet(1, 'string') + >>> finite_set = odl.FiniteSet(1, 'string') """ self.__elements = tuple(unique(elements)) @@ -882,12 +927,12 @@ def __contains__(self, other): Examples -------- - >>> set = odl.FiniteSet(1, 'string') - >>> 1 in set + >>> finite_set = odl.FiniteSet(1, 'string') + >>> 1 in finite_set True - >>> 2 in set + >>> 2 in finite_set False - >>> 'string' in set + >>> 'string' in finite_set True """ return other in self.elements @@ -929,10 +974,10 @@ def __getitem__(self, indcs): Examples -------- - >>> set = odl.FiniteSet(1, 2, 3, 'string') - >>> set[:3] + >>> finite_set = odl.FiniteSet(1, 2, 3, 'string') + >>> finite_set[:3] FiniteSet(1, 2, 3) - >>> set[3] + >>> finite_set[3] 'string' """ if isinstance(indcs, slice): @@ -945,11 +990,13 @@ def __repr__(self): Examples -------- - >>> odl.FiniteSet(1, 'string') + >>> finite_set = odl.FiniteSet(1, 'string') + >>> finite_set FiniteSet(1, 'string') """ - elements_str = ', '.join(repr(el) for el in self.elements) - return '{}({})'.format(self.__class__.__name__, elements_str) + posargs = self.elements + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts) if __name__ == '__main__': diff --git a/odl/solvers/functional/default_functionals.py b/odl/solvers/functional/default_functionals.py index b54d32a6bbc..4c0144656d8 100644 --- a/odl/solvers/functional/default_functionals.py +++ b/odl/solvers/functional/default_functionals.py @@ -2311,7 +2311,7 @@ def __init__(self, space, gamma): Parameters ---------- - space : `FnBase` + space : `TensorSpace` Domain of the functional. gamma : float Smoothing parameter of the Huber functional. If ``gamma = 0``, diff --git a/odl/solvers/nonsmooth/proximal_operators.py b/odl/solvers/nonsmooth/proximal_operators.py index d8955b2a7dc..d8ee620a999 100644 --- a/odl/solvers/nonsmooth/proximal_operators.py +++ b/odl/solvers/nonsmooth/proximal_operators.py @@ -1697,7 +1697,7 @@ def proximal_huber(space, gamma): Parameters ---------- - space : `FnBase` + space : `TensorSpace` The domain of the functional gamma : float The smoothing parameter of the Huber norm functional. diff --git a/odl/space/base_tensors.py b/odl/space/base_tensors.py index 67092cf1a17..bc9afca2673 100644 --- a/odl/space/base_tensors.py +++ b/odl/space/base_tensors.py @@ -18,7 +18,8 @@ from odl.util import ( is_numeric_dtype, is_real_dtype, is_floating_dtype, is_real_floating_dtype, is_complex_floating_dtype, safe_int_conv, - array_str, dtype_str, signature_string, indent, writable_array) + array_str, dtype_str, signature_string_parts, repr_string, + element_repr_string, writable_array) from odl.util.ufuncs import TensorSpaceUfuncs from odl.util.utility import TYPE_MAP_R2C, TYPE_MAP_C2R @@ -258,6 +259,21 @@ def astype(self, dtype): else: return self._astype(dtype) + def asweighted(self, weighting): + """Return a copy of this space with new ``weighting``.""" + if weighting == self.weighting: # let `AttributeError` bubble up + return self + else: + return self._asweighted(weighting) + + def _asweighted(self, weighting): + """Internal helper for `asweighted`. + + Subclasses with differing init parameters should overload this + method. + """ + return type(self)(self.shape, dtype=self.dtype, weighting=weighting) + @property def default_order(self): """Default storage order for new elements in this space. @@ -398,8 +414,8 @@ def __hash__(self): def __repr__(self): """Return ``repr(self)``.""" posargs = [self.shape, dtype_str(self.dtype)] - return "{}({})".format(self.__class__.__name__, - signature_string(posargs, [])) + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts) def __str__(self): """Return ``str(self)``.""" @@ -629,10 +645,7 @@ def __repr__(self): """Return ``repr(self)``.""" maxsize_full_print = 2 * np.get_printoptions()['edgeitems'] self_str = array_str(self, nprint=maxsize_full_print) - if self.ndim == 1 and self.size <= maxsize_full_print: - return '{!r}.element({})'.format(self.space, self_str) - else: - return '{!r}.element(\n{}\n)'.format(self.space, indent(self_str)) + return element_repr_string(repr(self.space), self_str) def __str__(self): """Return ``str(self)``.""" diff --git a/odl/space/npy_tensors.py b/odl/space/npy_tensors.py index e7e45a0b17d..24da6c13405 100644 --- a/odl/space/npy_tensors.py +++ b/odl/space/npy_tensors.py @@ -22,8 +22,8 @@ Weighting, ArrayWeighting, ConstWeighting, CustomInner, CustomNorm, CustomDist) from odl.util import ( - dtype_str, signature_string, is_real_dtype, is_numeric_dtype, - writable_array, is_floating_dtype) + dtype_str, signature_string_parts, repr_string, attribute_repr_string, + is_real_dtype, is_numeric_dtype, writable_array, is_floating_dtype) __all__ = ('NumpyTensorSpace',) @@ -807,7 +807,7 @@ def __getitem__(self, indices): def __repr__(self): """Return ``repr(self)``.""" - return repr(space) + '.byaxis' + return attribute_repr_string(repr(space), 'byaxis') return NpyTensorSpacebyaxis() @@ -837,12 +837,14 @@ def __repr__(self): optargs = [] optmod = '' - inner_str = signature_string(posargs, optargs, mod=['', optmod]) + inner_parts = signature_string_parts(posargs, optargs, + mod=['', optmod]) + inner_parts = [list(p) for p in inner_parts] weight_str = self.weighting.repr_part if weight_str: - inner_str += ', ' + weight_str + inner_parts[1].append(weight_str) - return '{}({})'.format(ctor_name, inner_str) + return repr_string(ctor_name, inner_parts) @property def element_type(self): diff --git a/odl/space/pspace.py b/odl/space/pspace.py index f57f19f4926..7c0d20c1d86 100644 --- a/odl/space/pspace.py +++ b/odl/space/pspace.py @@ -8,20 +8,23 @@ """Cartesian products of `LinearSpace` instances.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + from itertools import product from numbers import Integral + import numpy as np from odl.set import LinearSpace from odl.set.space import LinearSpaceElement from odl.space.weighting import ( - Weighting, ArrayWeighting, ConstWeighting, - CustomInner, CustomNorm, CustomDist) -from odl.util import is_real_dtype, signature_string, indent + ArrayWeighting, ConstWeighting, CustomDist, CustomInner, CustomNorm, + Weighting) +from odl.util import ( + dedent, element_repr_string, indent, is_real_dtype, repr_string, + signature_string_parts) from odl.util.ufuncs import ProductSpaceUfuncs - __all__ = ('ProductSpace',) @@ -127,11 +130,11 @@ class instance can be retrieved from the space by the -------- Product of two rn spaces - >>> r2x3 = ProductSpace(odl.rn(2), odl.rn(3)) + >>> r2x3 = odl.ProductSpace(odl.rn(2), odl.rn(3)) Powerspace of rn space - >>> r2x2x2 = ProductSpace(odl.rn(2), 3) + >>> r2x2x2 = odl.ProductSpace(odl.rn(2), 3) Notes ----- @@ -407,7 +410,7 @@ def element(self, inp=None, cast=True): -------- >>> r2, r3 = odl.rn(2), odl.rn(3) >>> vec_2, vec_3 = r2.element(), r3.element() - >>> r2x3 = ProductSpace(r2, r3) + >>> r2x3 = odl.ProductSpace(r2, r3) >>> vec_2x3 = r2x3.element() >>> vec_2.space == vec_2x3[0].space True @@ -417,15 +420,17 @@ def element(self, inp=None, cast=True): Create an element of the product space >>> r2, r3 = odl.rn(2), odl.rn(3) - >>> prod = ProductSpace(r2, r3) + >>> prod = odl.ProductSpace(r2, r3) >>> x2 = r2.element([1, 2]) >>> x3 = r3.element([1, 2, 3]) >>> x = prod.element([x2, x3]) >>> x - ProductSpace(rn(2), rn(3)).element([ - [ 1., 2.], - [ 1., 2., 3.] - ]) + ProductSpace(rn(2), rn(3)).element( + [ + [ 1., 2.], + [ 1., 2., 3.] + ] + ) """ # If data is given as keyword arg, prefer it over arg list if inp is None: @@ -478,7 +483,7 @@ def zero(self): -------- >>> r2, r3 = odl.rn(2), odl.rn(3) >>> zero_2, zero_3 = r2.zero(), r3.zero() - >>> r2x3 = ProductSpace(r2, r3) + >>> r2x3 = odl.ProductSpace(r2, r3) >>> zero_2x3 = r2x3.zero() >>> zero_2 == zero_2x3[0] True @@ -506,7 +511,7 @@ def one(self): -------- >>> r2, r3 = odl.rn(2), odl.rn(3) >>> one_2, one_3 = r2.one(), r3.one() - >>> r2x3 = ProductSpace(r2, r3) + >>> r2x3 = odl.ProductSpace(r2, r3) >>> one_2x3 = r2x3.one() >>> one_2 == one_2x3[0] True @@ -558,13 +563,13 @@ def __eq__(self, other): -------- >>> r2, r3 = odl.rn(2), odl.rn(3) >>> rn, rm = odl.rn(2), odl.rn(3) - >>> r2x3, rnxm = ProductSpace(r2, r3), ProductSpace(rn, rm) + >>> r2x3, rnxm = odl.ProductSpace(r2, r3), odl.ProductSpace(rn, rm) >>> r2x3 == rnxm True - >>> r3x2 = ProductSpace(r3, r2) + >>> r3x2 = odl.ProductSpace(r3, r2) >>> r2x3 == r3x2 False - >>> r5 = ProductSpace(*[odl.rn(1)]*5) + >>> r5 = odl.ProductSpace(*[odl.rn(1)]*5) >>> r2x3 == r5 False >>> r5 = odl.rn(5) @@ -684,7 +689,7 @@ def __str__(self): elif self.is_power_space: return '({}) ** {}'.format(self.spaces[0], len(self)) else: - return ' x '.join(str(space) for space in self.spaces) + return ' * '.join(str(space) for space in self.spaces) def __repr__(self): """Return ``repr(self)``.""" @@ -694,40 +699,29 @@ def __repr__(self): posargs = [] posmod = '' optargs = [('field', self.field, None)] - oneline = True + optmod = [''] elif self.is_power_space: posargs = [self.spaces[0], len(self)] posmod = '!r' - optargs = [] - oneline = True + optargs = optmod = [] elif self.size <= 2 * edgeitems: posargs = self.spaces posmod = '!r' - optargs = [] - argstr = ', '.join(repr(s) for s in self.spaces) - oneline = (len(argstr + weight_str) <= 40 and - '\n' not in argstr + weight_str) + optargs = optmod = [] else: posargs = (self.spaces[:edgeitems] + ('...',) + self.spaces[-edgeitems:]) posmod = ['!r'] * edgeitems + ['!s'] + ['!r'] * edgeitems - optargs = [] - oneline = False - - if oneline: - inner_str = signature_string(posargs, optargs, sep=', ', - mod=[posmod, '!r']) - if weight_str: - inner_str = ', '.join([inner_str, weight_str]) - return '{}({})'.format(self.__class__.__name__, inner_str) - else: - inner_str = signature_string(posargs, optargs, sep=',\n', - mod=[posmod, '!r']) - if weight_str: - inner_str = ',\n'.join([inner_str, weight_str]) - return '{}(\n{}\n)'.format(self.__class__.__name__, - indent(inner_str)) + optargs = optmod = [] + + inner_parts = signature_string_parts(posargs, optargs, + mod=[posmod, optmod]) + inner_parts = [list(p) for p in inner_parts] + if weight_str: + inner_parts[1].append(weight_str) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) @property def element_type(self): @@ -1018,24 +1012,30 @@ def ufuncs(self): >>> r22 = odl.ProductSpace(odl.rn(2), 2) >>> x = r22.element([[1, -2], [-3, 4]]) >>> x.ufuncs.absolute() - ProductSpace(rn(2), 2).element([ - [ 1., 2.], - [ 3., 4.] - ]) + ProductSpace(rn(2), 2).element( + [ + [ 1., 2.], + [ 3., 4.] + ] + ) These functions can also be used with non-vector arguments and support broadcasting, per component and even recursively: >>> x.ufuncs.add([1, 2]) - ProductSpace(rn(2), 2).element([ - [ 2., 0.], - [-2., 6.] - ]) + ProductSpace(rn(2), 2).element( + [ + [ 2., 0.], + [-2., 6.] + ] + ) >>> x.ufuncs.subtract(1) - ProductSpace(rn(2), 2).element([ - [ 0., -3.], - [-4., 3.] - ]) + ProductSpace(rn(2), 2).element( + [ + [ 0., -3.], + [-4., 3.] + ] + ) There is also support for various reductions (sum, prod, min, max): @@ -1047,10 +1047,12 @@ def ufuncs(self): >>> y = r22.element() >>> result = x.ufuncs.absolute(out=y) >>> result - ProductSpace(rn(2), 2).element([ - [ 1., 2.], - [ 3., 4.] - ]) + ProductSpace(rn(2), 2).element( + [ + [ 1., 2.], + [ 3., 4.] + ] + ) >>> result is y True @@ -1073,53 +1075,59 @@ def __repr__(self): Examples -------- - >>> from odl import rn # need to import rn into namespace >>> r2, r3 = odl.rn(2), odl.rn(3) - >>> r2x3 = ProductSpace(r2, r3) + >>> r2x3 = odl.ProductSpace(r2, r3) >>> x = r2x3.element([[1, 2], [3, 4, 5]]) - >>> eval(repr(x)) == x - True The result is readable: >>> x - ProductSpace(rn(2), rn(3)).element([ - [ 1., 2.], - [ 3., 4., 5.] - ]) + ProductSpace(rn(2), rn(3)).element( + [ + [ 1., 2.], + [ 3., 4., 5.] + ] + ) - Nestled spaces work as well: + Nested spaces work as well: - >>> X = ProductSpace(r2x3, r2x3) - >>> x = X.element([[[1, 2], [3, 4, 5]],[[1, 2], [3, 4, 5]]]) - >>> eval(repr(x)) == x - True + >>> pspace = odl.ProductSpace(r2x3, r2x3) + >>> x = pspace.element([[[1, 2], + ... [3, 4, 5]], + ... [[1, 2], + ... [3, 4, 5]]]) >>> x - ProductSpace(ProductSpace(rn(2), rn(3)), 2).element([ - [ - [ 1., 2.], - [ 3., 4., 5.] - ], + ProductSpace(ProductSpace(rn(2), rn(3)), 2).element( [ - [ 1., 2.], - [ 3., 4., 5.] + [ + [ 1., 2.], + [ 3., 4., 5.] + ], + [ + [ 1., 2.], + [ 3., 4., 5.] + ] ] - ]) + ) """ - inner_str = '[\n' - if len(self) < 5: - inner_str += ',\n'.join('{}'.format( - _indent(_strip_space(part))) for part in self.parts) + data_str = '[\n' + edgeitems = np.get_printoptions()['edgeitems'] + if len(self) <= 2 * edgeitems: + data_str += ',\n'.join( + '{}'.format(indent(_strip_space(part), ' ')) + for part in self.parts) else: - inner_str += ',\n'.join('{}'.format( - _indent(_strip_space(part))) for part in self.parts[:3]) - inner_str += ',\n ...\n' - inner_str += ',\n'.join('{}'.format( - _indent(_strip_space(part))) for part in self.parts[-1:]) + data_str += ',\n'.join( + '{}'.format(indent(_strip_space(part), ' ')) + for part in self.parts[:edgeitems]) + data_str += ',\n ...\n' + data_str += ',\n'.join( + '{}'.format(indent(_strip_space(part), ' ')) + for part in self.parts[-edgeitems:]) - inner_str += '\n]' + data_str += '\n]' - return '{!r}.element({})'.format(self.space, inner_str) + return element_repr_string(repr(self.space), data_str) def show(self, title=None, indices=None, **kwargs): """Display the parts of this product space element graphically. @@ -1232,7 +1240,9 @@ def show(self, title=None, indices=None, **kwargs): return tuple(figs) -# --- Add arithmetic operators that broadcast --- +# --- Add arithmetic operators that broadcast --- # + + def _broadcast_arithmetic(op): """Return ``op(self, other)`` with broadcasting. @@ -1259,25 +1269,33 @@ def _broadcast_arithmetic(op): layer" broadcasting, i.e., we do not support broadcasting over several product spaces at once. """ - def _broadcast_arithmetic_impl(self, other): - if (self.space.is_power_space and other in self.space[0]): + def broadcast_arithmetic_wrapper(self, other): + """Wrapper function for the arithmetic operation.""" + if other in self.space: + return getattr(LinearSpaceElement, op)(self, other) + + elif self.space.is_power_space and other in self.space[0]: + # Implement broadcasting along implicit axes "to the left" -- + # corresponding to Numpy broadcasting of the shapes + # (M, N) and (N,) results = [] - for xi in self: - res = getattr(xi, op)(other) + for self_i in self: + res = getattr(self_i, op)(other) if res is NotImplemented: return NotImplemented else: results.append(res) return self.space.element(results) + else: return getattr(LinearSpaceElement, op)(self, other) # Set docstring - docstring = """Broadcasted {op}.""".format(op=op) - _broadcast_arithmetic_impl.__doc__ = docstring + docstring = """Broadcasting ``{}``.""".format(op) + broadcast_arithmetic_wrapper.__doc__ = docstring - return _broadcast_arithmetic_impl + return broadcast_arithmetic_wrapper for op in ['add', 'sub', 'mul', 'div', 'truediv']: @@ -1610,17 +1628,12 @@ def _strip_space(x): space_repr = '{!r}.element('.format(x.space) if r.startswith(space_repr) and r.endswith(')'): r = r[len(space_repr):-1] + if r.startswith('\n'): + r = r.strip('\n') + r = dedent(r, max_levels=1) return r -def _indent(x): - """Indent a string by 4 characters.""" - lines = x.splitlines() - for i, line in enumerate(lines): - lines[i] = ' ' + line - return '\n'.join(lines) - - if __name__ == '__main__': from odl.util.testutils import run_doctests run_doctests() diff --git a/odl/test/operator/oputils_test.py b/odl/test/operator/oputils_test.py index 79e941f2b2b..49f98df7d1d 100644 --- a/odl/test/operator/oputils_test.py +++ b/odl/test/operator/oputils_test.py @@ -11,10 +11,16 @@ import pytest import odl -from odl.operator.oputils import matrix_representation, power_method_opnorm -from odl.space.pspace import ProductSpace +from odl.operator.oputils import ( + matrix_representation, power_method_opnorm) from odl.operator.pspace_ops import ProductSpaceOperator -from odl.util.testutils import almost_equal +from odl.space.pspace import ProductSpace +from odl.util.testutils import almost_equal, simple_fixture + + +optimize_weighting = simple_fixture('optimize', [True, False]) +call_variant = simple_fixture('call_variant', ['oop', 'ip', 'dual']) +weighting = simple_fixture('weighting', [1.0, 2.0, [1.0, 2.0]]) def test_matrix_representation(): diff --git a/odl/util/graphics.py b/odl/util/graphics.py index 240937bf54b..5868686b0f1 100644 --- a/odl/util/graphics.py +++ b/odl/util/graphics.py @@ -8,14 +8,15 @@ """Functions for graphical output.""" -from __future__ import print_function, division, absolute_import -import numpy as np +from __future__ import absolute_import, division, print_function + import warnings +import numpy as np + from odl.util.testutils import run_doctests from odl.util.utility import is_real_dtype - __all__ = ('show_discrete_data',) diff --git a/odl/util/normalize.py b/odl/util/normalize.py index daf09e48130..f436c7e58dd 100644 --- a/odl/util/normalize.py +++ b/odl/util/normalize.py @@ -8,9 +8,9 @@ """Utilities for normalization of user input.""" -from __future__ import print_function, division, absolute_import -import numpy as np +from __future__ import absolute_import, division, print_function +import numpy as np __all__ = ('normalized_scalar_param_list', 'normalized_index_expression', 'normalized_nodes_on_bdry', 'normalized_axes_tuple', diff --git a/odl/util/npy_compat.py b/odl/util/npy_compat.py index 8b40310b833..608fd20239a 100644 --- a/odl/util/npy_compat.py +++ b/odl/util/npy_compat.py @@ -8,9 +8,9 @@ """Numpy functions not available in the minimal required version.""" -from __future__ import print_function, division, absolute_import -import numpy as np +from __future__ import absolute_import, division, print_function +import numpy as np __all__ = ('broadcast_to', 'moveaxis') diff --git a/odl/util/numerics.py b/odl/util/numerics.py index 9366b02625f..e348c2494db 100644 --- a/odl/util/numerics.py +++ b/odl/util/numerics.py @@ -8,12 +8,12 @@ """Numerical helper functions for convenience or speed.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + import numpy as np from odl.util.normalize import normalized_scalar_param_list, safe_int_conv - __all__ = ('apply_on_boundary', 'fast_1d_tensor_mult', 'resize_array', 'zscore') diff --git a/odl/util/pytest_plugins.py b/odl/util/pytest_plugins.py index 331fa70e051..8fdef174a32 100644 --- a/odl/util/pytest_plugins.py +++ b/odl/util/pytest_plugins.py @@ -8,11 +8,13 @@ """Test configuration file.""" -from __future__ import print_function, division, absolute_import -import numpy as np +from __future__ import absolute_import, division, print_function + import operator import os +import numpy as np + import odl from odl.space.entry_points import tensor_space_impl_names from odl.trafos.backends import PYFFTW_AVAILABLE, PYWT_AVAILABLE diff --git a/odl/util/testutils.py b/odl/util/testutils.py index 78f0f12a992..a4e4b549b77 100644 --- a/odl/util/testutils.py +++ b/odl/util/testutils.py @@ -20,23 +20,23 @@ from odl.util.utility import run_from_ipython, is_string -__all__ = ('almost_equal', 'all_equal', 'all_almost_equal', 'never_skip', - 'skip_if_no_stir', 'skip_if_no_pywavelets', +__all__ = ('all_equal', 'all_almost_equal', 'is_subdict', + 'never_skip', 'skip_if_no_stir', 'skip_if_no_pywavelets', 'skip_if_no_pyfftw', 'skip_if_no_largescale', - 'noise_array', 'noise_element', 'noise_elements', - 'Timer', 'timeit', 'ProgressBar', 'ProgressRange', - 'test', 'run_doctests', 'test_file') + 'skip_if_no_benchmark', 'simple_fixture', 'noise_array', + 'noise_element', 'noise_elements', 'FailCounter', 'Timer', + 'timeit', 'ProgressBar', 'ProgressRange', 'test', 'run_doctests', + 'test_file') def _places(a, b, default=None): """Return number of expected correct digits between ``a`` and ``b``. Returned numbers if one of ``a.dtype`` and ``b.dtype`` is as below: - 1 -- for ``np.float16`` - 3 -- for ``np.float32`` or ``np.complex64`` - - 5 -- for all other cases + - 1 -- for ``np.float16`` + - 3 -- for ``np.float32`` or ``np.complex64`` + - 5 -- for all other cases """ dtype1 = getattr(a, 'dtype', object) dtype2 = getattr(b, 'dtype', object) @@ -47,11 +47,10 @@ def dtype_places(dtype, default=None): """Return number of correct digits expected for given dtype. Returned numbers: - 1 -- for ``np.float16`` - - 3 -- for ``np.float32`` or ``np.complex64`` - 5 -- for all other cases + - 1 -- for ``np.float16`` + - 3 -- for ``np.float32`` or ``np.complex64`` + - 5 -- for all other cases """ small_dtypes = [np.float32, np.complex64] tiny_dtypes = [np.float16] @@ -64,34 +63,6 @@ def dtype_places(dtype, default=None): return default if default is not None else 5 -def almost_equal(a, b, places=None): - """Return ``True`` if the scalars ``a`` and ``b`` are almost equal.""" - if a is None and b is None: - return True - - if places is None: - places = _places(a, b) - - eps = 10 ** -places - - try: - complex(a) - complex(b) - except TypeError: - return False - - if np.isnan(a) and np.isnan(b): - return True - - if np.isinf(a) and np.isinf(b): - return a == b - - if abs(complex(b)) < eps: - return abs(complex(a) - complex(b)) < eps - else: - return abs(a / b - 1) < eps - - def all_equal(iter1, iter2): """Return ``True`` if all elements in ``a`` and ``b`` are equal.""" # Direct comparison for scalars, tuples or lists @@ -105,6 +76,10 @@ def all_equal(iter1, iter2): if iter1 is None and iter2 is None: return True + # Do it faster for arrays + if hasattr(iter1, '__array__') and hasattr(iter2, '__array__'): + return np.array_equal(iter1, iter2) + # If one nested iterator is exhausted, go to direct comparison try: it1 = iter(iter1) @@ -132,18 +107,12 @@ def all_equal(iter1, iter2): return True -def all_almost_equal_array(v1, v2, places): - return np.allclose(v1, v2, - rtol=10 ** (-places), atol=10 ** (-places), - equal_nan=True) - - def all_almost_equal(iter1, iter2, places=None): """Return ``True`` if all elements in ``a`` and ``b`` are almost equal.""" try: if iter1 is iter2 or iter1 == iter2: return True - except ValueError: + except (ValueError, TypeError): pass if iter1 is None and iter2 is None: @@ -154,13 +123,23 @@ def all_almost_equal(iter1, iter2, places=None): # otherwise for recursive calls. if places is None: places = _places(iter1, iter2, None) - return all_almost_equal_array(iter1, iter2, places) + return np.allclose(iter1, iter2, + rtol=10 ** (-places), atol=10 ** (-places), + equal_nan=True) try: it1 = iter(iter1) it2 = iter(iter2) except TypeError: - return almost_equal(iter1, iter2, places) + if places is None: + places = _places(iter1, iter2, None) + isclose = np.isclose(iter1, iter2, + atol=10 ** (-places), + equal_nan=True) + try: + return bool(isclose) + except ValueError: + return False diff_length_sentinel = object() for [ip1, ip2] in zip_longest(it1, it2, @@ -304,7 +283,7 @@ def noise_array(space): Returns ------- - noise_array : `numpy.ndarray` element + noise_array : `numpy.ndarray` Array with white noise such that ``space.element``'s can be created from it. @@ -322,9 +301,21 @@ def noise_array(space): odl.set.space.LinearSpace.examples : Examples of elements typical to the space. """ - from odl.space import ProductSpace + from odl.space.pspace import ProductSpace + if isinstance(space, ProductSpace): - return np.array([noise_array(si) for si in space]) + arr_list = [noise_array(spc_i) for spc_i in space] + if space.is_power_space: + arr = np.empty((len(arr_list),) + arr_list[0].shape, + dtype=space[0].dtype) + else: + arr = np.empty((len(arr_list),) + arr_list[0].shape, + dtype=object) + for i in range(len(arr)): + arr[i] = arr_list[i] + + return arr + else: if space.dtype == bool: arr = np.random.randint(0, 2, size=space.shape, dtype=bool) @@ -433,13 +424,21 @@ def noise_elements(space, n=1): noise_array noise_element """ - arrs = tuple(noise_array(space) for _ in range(n)) + from odl.space.pspace import ProductSpace + + if isinstance(space, ProductSpace) and not space.is_power_space: + raise ValueError('`space` cannot be a non-power product space') + + if isinstance(space, ProductSpace): + impl = space[0].impl + else: + impl = space.impl - # Make space elements from arrays - elems = tuple(space.element(arr.copy()) for arr in arrs) + arrs = tuple(noise_array(space, impl=impl) for _ in range(n)) + elems = tuple(space.element(arr) for arr in arrs) if n == 1: - return tuple(arrs + elems) + return arrs + elems else: return arrs, elems @@ -555,23 +554,23 @@ class ProgressBar(object): Usage: >>> progress = ProgressBar('Reading data', 10) - \rReading data: [ ] Starting + Reading data: [ ] Starting >>> progress.update(4) #halfway, zero indexing - \rReading data: [############### ] 50.0% + Reading data: [############### ] 50.0% Multi-indices, from slowest to fastest: >>> progress = ProgressBar('Reading data', 10, 10) - \rReading data: [ ] Starting + Reading data: [ ] Starting >>> progress.update(9, 8) - \rReading data: [############################# ] 99.0% + Reading data: [############################# ] 99.0% Supports simply calling update, which moves the counter forward: >>> progress = ProgressBar('Reading data', 10, 10) - \rReading data: [ ] Starting + Reading data: [ ] Starting >>> progress.update() - \rReading data: [ ] 1.0% + Reading data: [ ] 1.0% """ def __init__(self, text='progress', *njobs): @@ -685,12 +684,14 @@ def run_doctests(skip_if=False, **kwargs): Extra keyword arguments passed on to the ``doctest.testmod`` function. """ - from doctest import testmod, NORMALIZE_WHITESPACE, SKIP + from doctest import ( + testmod, NORMALIZE_WHITESPACE, SKIP, IGNORE_EXCEPTION_DETAIL) from pkg_resources import parse_version import odl import numpy as np - optionflags = kwargs.pop('optionflags', NORMALIZE_WHITESPACE) + optionflags = kwargs.pop('optionflags', + NORMALIZE_WHITESPACE | IGNORE_EXCEPTION_DETAIL) if skip_if: optionflags |= SKIP diff --git a/odl/util/ufuncs.py b/odl/util/ufuncs.py index caa6f1fbb0d..975b8e6e9b8 100644 --- a/odl/util/ufuncs.py +++ b/odl/util/ufuncs.py @@ -23,11 +23,12 @@ arrays.classes.rst#special-attributes-and-methods>`_. """ -from __future__ import print_function, division, absolute_import -from builtins import object -import numpy as np +from __future__ import absolute_import, division, print_function + import re +from builtins import object +import numpy as np __all__ = ('TensorSpaceUfuncs', 'ProductSpaceUfuncs') diff --git a/odl/util/utility.py b/odl/util/utility.py index 25fc3b268d7..660b347af7d 100644 --- a/odl/util/utility.py +++ b/odl/util/utility.py @@ -8,22 +8,28 @@ """Utilities mainly for internal use.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + +import inspect +import sys from builtins import object from collections import OrderedDict from functools import wraps -import inspect + import numpy as np -import sys +__all__ = ( + 'array_str', 'dtype_str', 'dtype_repr', 'npy_printoptions', + 'signature_string', 'signature_string_parts', 'repr_string', + 'indent', 'dedent', 'attribute_repr_string', 'element_repr_string', + 'is_numeric_dtype', 'is_int_dtype', 'is_floating_dtype', 'is_real_dtype', + 'is_real_floating_dtype', 'is_complex_floating_dtype', + 'real_dtype', 'complex_dtype', 'is_string', 'conj_exponent', + 'writable_array', 'run_from_ipython', 'NumpyRandomSeed', + 'cache_arguments', 'unique', + 'REPR_PRECISION') -__all__ = ('array_str', 'dtype_str', 'dtype_repr', 'npy_printoptions', - 'signature_string', 'indent', - 'is_numeric_dtype', 'is_int_dtype', 'is_floating_dtype', - 'is_real_dtype', 'is_real_floating_dtype', - 'is_complex_floating_dtype', 'real_dtype', 'complex_dtype', - 'is_string', 'conj_exponent', 'writable_array', - 'run_from_ipython', 'NumpyRandomSeed', 'cache_arguments', 'unique') +REPR_PRECISION = 4 # For printing scalars and array entries TYPE_MAP_R2C = {np.dtype(dtype): np.result_type(dtype, 1j) for dtype in np.sctypes['float']} @@ -70,12 +76,80 @@ def indent(string, indent_str=' '): Indenting by random stuff: - >>> print(indent(text, indent_str='|----|')) - |----|This is line 1. - |----|Next line. - |----|And another one. + >>> print(indent(text, indent_str='<->')) + <->This is line 1. + <->Next line. + <->And another one. + """ + return '\n'.join(indent_str + row for row in string.splitlines()) + + +def dedent(string, indent_str=' ', max_levels=None): + """Revert the effect of indentation. + + Examples + -------- + Remove a simple one-level indentation: + + >>> text = '''<->This is line 1. + ... <->Next line. + ... <->And another one.''' + >>> print(text) + <->This is line 1. + <->Next line. + <->And another one. + >>> print(dedent(text, '<->')) + This is line 1. + Next line. + And another one. + + Multiple levels of indentation: + + >>> text = '''<->Level 1. + ... <-><->Level 2. + ... <-><-><->Level 3.''' + >>> print(text) + <->Level 1. + <-><->Level 2. + <-><-><->Level 3. + >>> print(dedent(text, '<->')) + Level 1. + <->Level 2. + <-><->Level 3. + >>> text = '''<-><->Level 2. + ... <-><-><->Level 3.''' + >>> print(text) + <-><->Level 2. + <-><-><->Level 3. + >>> print(dedent(text, '<->')) + Level 2. + <->Level 3. + >>> print(dedent(text, '<->', max_levels=1)) + <->Level 2. + <-><->Level 3. """ - return '\n'.join(indent_str + row for row in string.split('\n')) + lines = string.splitlines() + + # Determine common (minumum) number of indentation levels, capped at + # `max_levels` if given + def num_indents(line): + nindents = 0 + while True: + if line.startswith(indent_str): + nindents += 1 + line = line[len(indent_str):] + else: + break + + return nindents + + num_levels = num_indents(min(lines, key=num_indents)) + if max_levels is not None: + num_levels = min(num_levels, max_levels) + + # Dedent + dedent_len = num_levels * len(indent_str) + return '\n'.join(line[dedent_len:] for line in lines) class npy_printoptions(object): @@ -623,7 +697,7 @@ def __exit__(self, type, value, traceback): self.arr = None -def signature_string(posargs, optargs, sep=', ', mod=''): +def signature_string(posargs, optargs, sep=', ', mod='!r'): """Return a stringified signature from given arguments. Parameters @@ -673,6 +747,11 @@ def signature_string(posargs, optargs, sep=', ', mod=''): Finally, if ``mod`` is a string or callable, it is applied to all arguments. + The default behavior is to apply the "{!r}" (``repr``) conversion. + For floating point scalars, the number of digits printed is + determined by the ``precision`` value in NumPy's printing options, + which can be temporarily modified with `npy_printoptions`. + Returns ------- signature : string @@ -738,6 +817,18 @@ def signature_string(posargs, optargs, sep=', ', mod=''): >>> optargs = [] >>> signature_string(posargs, optargs, mod=[['', array_str], []]) "'arg1', [ 1., 1., 1.]" + + The number of printed digits in floating point numbers can be changed + with `npy_printoptions`: + + >>> posargs = ['hello', 0.123456789012345] + >>> optargs = [('extent', 1.234567890123456, 1.0)] + >>> signature_string(posargs, optargs) # default is 8 digits + "'hello', 0.12345679, extent=1.2345679" + >>> with npy_printoptions(precision=2): + ... sig_str = signature_string(posargs, optargs) + >>> sig_str + "'hello', 0.12, extent=1.2" """ # Define the separators for the two possible cases if is_string(sep): @@ -745,6 +836,72 @@ def signature_string(posargs, optargs, sep=', ', mod=''): else: pos_sep, opt_sep, part_sep = sep + # Get the stringified parts + posargs_conv, optargs_conv = signature_string_parts(posargs, optargs, mod) + + # Join the arguments using the separators + parts = [] + if posargs_conv: + parts.append(pos_sep.join(argstr for argstr in posargs_conv)) + if optargs_conv: + parts.append(opt_sep.join(optargs_conv)) + + return part_sep.join(parts) + + +def signature_string_parts(posargs, optargs, mod='!r'): + """Return stringified arguments as tuples. + + Parameters + ---------- + posargs : sequence + Positional argument values, always included in the returned string + tuple. + optargs : sequence of 3-tuples + Optional arguments with names and defaults, given in the form:: + + [(name1, value1, default1), (name2, value2, default2), ...] + + Only those parameters that are different from the given default + are included as ``name=value`` keyword pairs. + + **Note:** The comparison is done by using ``if value == default:``, + which is not valid for, e.g., NumPy arrays. + + mod : string or callable or sequence, optional + Format modifier(s) for the argument strings. + In its most general form, ``mod`` is a sequence of 2 sequences + ``pos_mod, opt_mod`` with ``len(pos_mod) == len(posargs)`` and + ``len(opt_mod) == len(optargs)``. Each entry ``m`` in those sequences + can be a string, resulting in the following stringification + of ``arg``:: + + arg_fmt = {{{}}}.format(m) + arg_str = arg_fmt.format(arg) + + For a callable ``to_str``, the stringification is simply + ``arg_str = to_str(arg)``. + + The entries ``pos_mod, opt_mod`` of ``mod`` can also be strings + or callables instead of sequences, in which case the modifier + applies to all corresponding arguments. + + Finally, if ``mod`` is a string or callable, it is applied to + all arguments. + + The default behavior is to apply the "{!r}" (``repr``) conversion. + For floating point scalars, the number of digits printed is + determined by the ``precision`` value in NumPy's printing options, + which can be temporarily modified with `npy_printoptions`. + + Returns + ------- + pos_strings : tuple of str + The stringified positional arguments. + opt_strings : tuple of str + The stringified optional arguments, not including the ones + equal to their respective defaults. + """ # Convert modifiers to 2-sequence of sequence of strings if is_string(mod) or callable(mod): pos_mod = opt_mod = mod @@ -765,9 +922,7 @@ def signature_string(posargs, optargs, sep=', ', mod=''): 'len({}) != len({})'.format(m, args)) pos_mod, opt_mod = mods - - # Convert the arguments to strings - parts = [] + precision = np.get_printoptions()['precision'] # Stringify values, treating strings specially posargs_conv = [] @@ -781,14 +936,17 @@ def signature_string(posargs, optargs, sep=', ', mod=''): else: fmt = "'{}'" posargs_conv.append(fmt.format(arg)) + elif (np.isscalar(arg) and + np.array(arg).real.astype('int64') != arg and + modifier in ('', '!s', '!r')): + # Floating point value, use Numpy print option 'precision' + fmt = '{{:.{}}}'.format(precision) + posargs_conv.append(fmt.format(arg)) else: # All non-string types are passed through a format conversion fmt = '{{{}}}'.format(modifier) posargs_conv.append(fmt.format(arg)) - if posargs_conv: - parts.append(pos_sep.join(argstr for argstr in posargs_conv)) - # Build 'key=value' strings for values that are not equal to default optargs_conv = [] for (name, value, default), modifier in zip(optargs, opt_mod): @@ -806,15 +964,337 @@ def signature_string(posargs, optargs, sep=', ', mod=''): fmt = "'{}'" value_str = fmt.format(value) optargs_conv.append('{}={}'.format(name, value_str)) + elif (np.isscalar(value) and + np.array(value).real.astype('int64') != value and + modifier in ('', '!s', '!r')): + fmt = '{{:.{}}}'.format(precision) + value_str = fmt.format(value) + optargs_conv.append('{}={}'.format(name, value_str)) else: fmt = '{{{}}}'.format(modifier) value_str = fmt.format(value) optargs_conv.append('{}={}'.format(name, value_str)) - if optargs_conv: - parts.append(opt_sep.join(optargs_conv)) + return tuple(posargs_conv), tuple(optargs_conv) - return part_sep.join(parts) + +def _separators(strings, linewidth): + """Return separators that keep joined strings within the line width.""" + if len(strings) <= 1: + return () + + indent_len = 4 + separators = [] + cur_line_len = indent_len + len(strings[0]) + 1 + if cur_line_len + 2 <= linewidth and '\n' not in strings[0]: + # Next string might fit on same line + separators.append(', ') + cur_line_len += 1 # for the extra space + else: + # Use linebreak if string contains newline or doesn't fit + separators.append(',\n') + cur_line_len = indent_len + + for i, s in enumerate(strings[1:-1]): + + cur_line_len += len(s) + 1 + + if '\n' in s: + # Use linebreak before and after if string contains newline + separators[i] = ',\n' + cur_line_len = indent_len + separators.append(',\n') + + elif cur_line_len + 2 <= linewidth: + # This string fits, next one might also fit on same line + separators.append(', ') + cur_line_len += 1 # for the extra space + + elif cur_line_len <= linewidth: + # This string fits, but next one won't + separators.append(',\n') + cur_line_len = indent_len + + else: + # This string doesn't fit but has no newlines in it + separators[i] = ',\n' + cur_line_len = indent_len + len(s) + 1 + + # Need to determine again what should come next + if cur_line_len + 2 <= linewidth: + # Next string might fit on same line + separators.append(', ') + else: + separators.append(',\n') + + cur_line_len += len(strings[-1]) + if cur_line_len + 1 > linewidth or '\n' in strings[-1]: + # This string and a comma don't fit on this line + separators[-1] = ',\n' + + return tuple(separators) + + +def repr_string(outer_string, inner_strings, allow_mixed_seps=True): + """Return a pretty string for ``repr``. + + The returned string is formatted such that it does not extend + beyond the line boundary if avoidable. The line width is taken from + NumPy's printing options that can be retrieved with + ``np.get_printoptions()``. They can be temporarily overridden + using the `npy_printoptions` context manager. See Examples for details. + + Parameters + ---------- + outer_str : str + Name of the class or function that should be printed outside + the parentheses. + inner_strings : sequence of sequence of str + Stringifications of the positional and optional arguments. + This is usually the return value of `signature_string_parts`. + allow_mixed_seps : bool, optional + If ``False`` and the string does not fit on one line, use + ``',\\n'`` to separate all strings. + By default, a mixture of ``', '`` and ``,\\n`` is used to fit + as much on one line as possible. + + Returns + ------- + repr_string : str + Full string that can be returned by a class' ``repr`` method. + + Examples + -------- + Things that fit into one line are printed on one line: + + >>> outer_string = 'MyClass' + >>> inner_strings = [('1', "'hello'", 'None'), + ... ("dtype='float32'",)] + >>> print(repr_string(outer_string, inner_strings)) + MyClass(1, 'hello', None, dtype='float32') + + Otherwise, if a part of ``inner_strings`` fits on a line of its own, + it is printed on one line, but separated from the other part with + a line break: + + >>> outer_string = 'MyClass' + >>> inner_strings = [('2.0', "'this_is_a_very_long_argument_string'"), + ... ("long_opt_arg='another_quite_long_string'",)] + >>> print(repr_string(outer_string, inner_strings)) + MyClass( + 2.0, 'this_is_a_very_long_argument_string', + long_opt_arg='another_quite_long_string' + ) + + If those parts are themselves too long, they are broken down into + several lines: + + >>> outer_string = 'MyClass' + >>> inner_strings = [("'this_is_a_very_long_argument_string'", + ... "'another_very_long_argument_string'"), + ... ("long_opt_arg='another_quite_long_string'", + ... "long_opt2_arg='this_wont_fit_on_one_line_either'")] + >>> print(repr_string(outer_string, inner_strings)) + MyClass( + 'this_is_a_very_long_argument_string', + 'another_very_long_argument_string', + long_opt_arg='another_quite_long_string', + long_opt2_arg='this_wont_fit_on_one_line_either' + ) + + The usage of mixed separators to optimally use horizontal space can + be disabled by setting ``allow_mixed_seps=False``: + + >>> outer_string = 'MyClass' + >>> inner_strings = [('2.0', "'this_is_a_very_long_argument_string'"), + ... ("long_opt_arg='another_quite_long_string'",)] + >>> print(repr_string(outer_string, inner_strings, allow_mixed_seps=False)) + MyClass( + 2.0, + 'this_is_a_very_long_argument_string', + long_opt_arg='another_quite_long_string' + ) + + With the ``npy_printoptions`` context manager, the available line + width can be changed: + + >>> outer_string = 'MyClass' + >>> inner_strings = [('1', "'hello'", 'None'), + ... ("dtype='float32'",)] + >>> with npy_printoptions(linewidth=20): + ... print(repr_string(outer_string, inner_strings)) + MyClass( + 1, 'hello', + None, + dtype='float32' + ) + """ + linewidth = np.get_printoptions()['linewidth'] + pos_strings, opt_strings = inner_strings + # Length of the positional and optional argument parts of the signature, + # including separators `', '` + pos_sig_len = (sum(len(pstr) for pstr in pos_strings) + + 2 * max((len(pos_strings) - 1), 0)) + opt_sig_len = (sum(len(pstr) for pstr in opt_strings) + + 2 * max((len(opt_strings) - 1), 0)) + + # Length of the one-line string, including 2 for the parentheses and + # 2 for the joining ', ' + repr_len = len(outer_string) + 2 + pos_sig_len + 2 + opt_sig_len + + if repr_len <= linewidth and not any('\n' in s + for s in pos_strings + opt_strings): + # Everything fits on one line + fmt = '{}({})' + pos_str = ', '.join(pos_strings) + opt_str = ', '.join(opt_strings) + parts_sep = ', ' + else: + # Need to split lines in some way + fmt = '{}(\n{}\n)' + + if not allow_mixed_seps: + pos_separators = [',\n'] * (len(pos_strings) - 1) + else: + pos_separators = _separators(pos_strings, linewidth) + if len(pos_strings) == 0: + pos_str = '' + else: + pos_str = pos_strings[0] + for s, sep in zip(pos_strings[1:], pos_separators): + pos_str = sep.join([pos_str, s]) + + if not allow_mixed_seps: + opt_separators = [',\n'] * (len(opt_strings) - 1) + else: + opt_separators = _separators(opt_strings, linewidth) + if len(opt_strings) == 0: + opt_str = '' + else: + opt_str = opt_strings[0] + for s, sep in zip(opt_strings[1:], opt_separators): + opt_str = sep.join([opt_str, s]) + + # Check if we can put both parts on one line. This requires their + # concatenation including 4 for indentation and 2 for ', ' to + # be less than the line width. And they should contain no newline. + if pos_str and opt_str: + inner_len = 4 + len(pos_str) + 2 + len(opt_str) + elif (pos_str and not opt_str) or (opt_str and not pos_str): + inner_len = 4 + len(pos_str) + len(opt_str) + else: + inner_len = 0 + + if (not allow_mixed_seps or + any('\n' in s for s in [pos_str, opt_str]) or + inner_len > linewidth): + parts_sep = ',\n' + pos_str = indent(pos_str) + opt_str = indent(opt_str) + else: + parts_sep = ', ' + pos_str = indent(pos_str) + # Don't indent `opt_str` + + parts = [s for s in [pos_str, opt_str] if s.strip()] # ignore empty + inner_string = parts_sep.join(parts) + return fmt.format(outer_string, inner_string) + + +def attribute_repr_string(instance_str, attr_str): + """Return a repr string for an attribute that respects line width. + + Parameters + ---------- + instance_str : str + Stringification of a class instance. + attr_str : str + Name of the attribute (not including the ``'.'``). + + Returns + ------- + attr_repr_str : str + Concatenation of the two strings in a way that the line width + is respected. + + Examples + -------- + >>> inst_str = 'rn((2, 3))' + >>> attr_str = 'byaxis' + >>> print(attribute_repr_string(inst_str, attr_str)) + rn((2, 3)).byaxis + >>> long_inst_str = ( + ... "MyClass('long string that will definitely trigger a line break')" + ... ) + >>> long_attr_str = 'long_attribute_name' + >>> print(attribute_repr_string(long_inst_str, long_attr_str)) + MyClass( + 'long string that will definitely trigger a line break' + ).long_attribute_name + """ + linewidth = np.get_printoptions()['linewidth'] + if len(instance_str) + len(attr_str) + 1 <= linewidth: + parts = [instance_str, attr_str] + else: + left, rest = instance_str.split('(', maxsplit=1) + right, middle = rest[::-1].split(')', maxsplit=1) + middle, right = middle[::-1], right[::-1] + new_instance_str = '(\n'.join([left, indent(middle)]) + '\n)' + right + parts = [new_instance_str, attr_str] + + return '.'.join(parts) + + +def element_repr_string(space_str, data_str): + """Return a repr string for a space element that respects line width. + + Parameters + ---------- + space_str : str + Stringification of the space in which the element lives. + data_str : str + Stringification of the data that the element contains. + + Returns + ------- + elem_repr_str : str + Stringification of the element in a way that the line width + is respected. + + Examples + -------- + >>> space_str = 'rn(3)' + >>> data_str = array_str([1.0, 0.0, 1.0]) + >>> print(element_repr_string(space_str, data_str)) + rn(3).element([ 1., 0., 1.]) + >>> space_str = ( + ... 'uniform_discr([ 0., 0.], [ 1., 1.], (2, 2), dtype=complex)' + ... ) + >>> long_data_str = array_str([[1.0, 2.0], + ... [3.0, 4.0]]) + >>> print(element_repr_string(space_str, long_data_str)) + uniform_discr([ 0., 0.], [ 1., 1.], (2, 2), dtype=complex).element( + [[ 1., 2.], + [ 3., 4.]] + ) + """ + linewidth = np.get_printoptions()['linewidth'] + + spc_len = len(space_str) + len('.element(') + if spc_len <= linewidth: + # Can fit the 'space.element(' part on one line + spc_part = space_str + '.element(' + else: + left, rest = space_str.split('(', maxsplit=1) + right, middle = rest[::-1].split(')', maxsplit=1) + middle, right = middle[::-1], right[::-1] + spc_part = ('(\n'.join([left, indent(middle)]) + '\n)' + right + + '.element(') + + if '\n' in data_str or len(spc_part) + len(data_str) + 1 > linewidth: + return spc_part + '\n' + indent(data_str) + '\n)' + else: + return spc_part + data_str + ')' def run_from_ipython(): diff --git a/odl/util/vectorization.py b/odl/util/vectorization.py index ffe1648e1d3..89b6385397e 100644 --- a/odl/util/vectorization.py +++ b/odl/util/vectorization.py @@ -8,11 +8,12 @@ """Utilities for internal functionality connected to vectorization.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + from builtins import object from functools import wraps -import numpy as np +import numpy as np __all__ = ('is_valid_input_array', 'is_valid_input_meshgrid', 'out_shape_from_meshgrid', 'out_shape_from_array',