From 84cec7e46bb275b94f4870ac28a92959cf11fe46 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Tue, 19 Jun 2018 23:28:11 +0200 Subject: [PATCH 1/9] ENH: update repr and docs of operator and space --- odl/discr/diff_ops.py | 2 +- odl/operator/__init__.py | 2 +- odl/operator/default_ops.py | 848 +++++++++++++++++++++++++----------- odl/operator/operator.py | 204 +++++---- odl/operator/oputils.py | 9 +- odl/operator/pspace_ops.py | 476 +++++++++++--------- odl/operator/tensor_ops.py | 546 ++++++++++++----------- odl/space/__init__.py | 2 +- odl/space/base_tensors.py | 58 ++- odl/space/entry_points.py | 2 +- odl/space/fspace.py | 29 +- odl/space/npy_tensors.py | 37 +- odl/space/pspace.py | 231 +++++----- odl/space/space_utils.py | 8 +- odl/space/weighting.py | 32 +- 15 files changed, 1486 insertions(+), 1000 deletions(-) diff --git a/odl/discr/diff_ops.py b/odl/discr/diff_ops.py index fad6dbb939f..49f5cee287b 100644 --- a/odl/discr/diff_ops.py +++ b/odl/discr/diff_ops.py @@ -503,7 +503,7 @@ def __init__(self, domain=None, range=None, method='forward', >>> g = div.range.element(data ** 2) >>> adj_div_g = div.adjoint(g) - >>> g.inner(div_f) / f.inner(adj_div_g) + >>> div(f).inner(g) / f.inner(div.adjoint(g)) 1.0 """ if domain is None and range is None: diff --git a/odl/operator/__init__.py b/odl/operator/__init__.py index b23d113b281..cc9ea55f788 100644 --- a/odl/operator/__init__.py +++ b/odl/operator/__init__.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # diff --git a/odl/operator/default_ops.py b/odl/operator/default_ops.py index 68fd8f9c069..7cb960b7813 100644 --- a/odl/operator/default_ops.py +++ b/odl/operator/default_ops.py @@ -1,6 +1,6 @@ # coding=utf-8 -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -10,15 +10,18 @@ """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, ComplexNumbers -from odl.set.space import LinearSpaceElement +from odl.set import ComplexNumbers, Field, LinearSpace, RealNumbers from odl.space import ProductSpace - +from odl.util import ( + REPR_PRECISION, attribute_repr_string, method_repr_string, + npy_printoptions, repr_string, signature_string_parts) __all__ = ('ScalingOperator', 'ZeroOperator', 'IdentityOperator', 'LinCombOperator', 'MultiplyOperator', 'PowerOperator', @@ -27,6 +30,45 @@ 'ComplexModulus', 'ComplexModulusSquared') +def _scale_op(operator, scalar): + """Scale an operator, optimizing for ``scalar=0`` and ``scalar=1``.""" + if scalar == 0: + return ZeroOperator(operator.domain, operator.range) + elif scalar == 1: + return operator + else: + return scalar * operator + + +def _lico_ops(a, op1, b, op2): + """Linear combination of operators, optimizing trivial cases.""" + if op1.domain != op2.domain or op1.range != op2.range: + raise ValueError('domain/range mismatch between {!r} and {!r}' + .format(op1, op2)) + dom, ran = op1.domain, op1.range + if a == 0: + if b == 0: + return ZeroOperator(dom, ran) + elif b == 1: + return op2 + else: + return b * op2 + elif a == 1: + if b == 0: + return op1 + elif b == 1: + return op1 + op2 + else: + return op1 + b * op2 + else: + if b == 0: + return a * op1 + elif b == 1: + return a * op1 + op2 + else: + return a * op1 + b * op2 + + class ScalingOperator(Operator): """Operator of multiplication with a scalar. @@ -36,7 +78,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 @@ -45,13 +87,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 @@ -63,7 +108,14 @@ def __init__(self, domain, scalar): 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 @@ -86,18 +138,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): @@ -132,9 +184,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. @@ -153,19 +206,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=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs) + return repr_string(self.__class__.__name__, inner_parts) class IdentityOperator(ScalingOperator): @@ -177,23 +238,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): @@ -221,7 +297,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 @@ -253,7 +329,7 @@ class MultiplyOperator(Operator): """Operator multiplying by a fixed space or field element. - Implements:: + Implements :: MultiplyOperator(y)(x) == x * y @@ -269,48 +345,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): @@ -342,26 +423,24 @@ 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.]) Multiplication operator with complex space: @@ -388,12 +467,23 @@ def adjoint(self): domain=self.range, 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): @@ -409,7 +499,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 @@ -418,25 +508,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): @@ -447,7 +542,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) @@ -472,30 +567,43 @@ 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) + if self.exponent == 1: + # Trivial case + return IdentityOperator(self.range, self.domain) + else: + 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): @@ -513,25 +621,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): @@ -555,11 +700,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): @@ -582,12 +728,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): @@ -607,29 +764,35 @@ 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``.""" return x.norm() def derivative(self, point): - """Derivative of this operator in ``point``. + r"""Derivative of this operator in ``point``. ``NormOperator().derivative(y)(x) == (y / y.norm()).inner(x)`` @@ -656,12 +819,12 @@ def derivative(self, point): spaces, in which case it is given by .. math:: - (D \|\cdot\|)(y)(x) = \langle y / \|y\|, x \\rangle + (D \|\cdot\|)(y)(x) = \langle y / \|y\|, x \rangle Examples -------- >>> r3 = odl.rn(3) - >>> op = NormOperator(r3) + >>> op = odl.NormOperator(r3) >>> derivative = op.derivative([1, 0, 0]) >>> derivative([1, 0, 0]) 1.0 @@ -671,15 +834,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): @@ -699,25 +869,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): @@ -729,7 +912,7 @@ def _call(self, x): return self.vector.dist(x) def derivative(self, point): - """The derivative operator. + r"""The derivative operator. ``DistOperator(y).derivative(z)(x) == ((y - z) / y.dist(z)).inner(x)`` @@ -757,13 +940,13 @@ def derivative(self, point): spaces, in which case it is given by .. math:: - (D d(\cdot, y))(z)(x) = \\langle (y-z) / d(y, z), x \\rangle + (D d(\cdot, y))(z)(x) = \langle (y-z) / d(y, z), x \rangle Examples -------- >>> 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 @@ -775,15 +958,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): @@ -801,37 +995,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): @@ -845,13 +1034,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. @@ -863,20 +1045,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): @@ -893,10 +1086,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 -------- @@ -937,15 +1131,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): @@ -957,28 +1157,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.]) @@ -986,14 +1189,18 @@ 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 - self.space_is_real = (space == real_space) - linear = True - super(RealPart, self).__init__(space, real_space, linear=linear) + if range is None: + range = domain.real_space + + # `linear=True` is a compromise in favor of efficiency, since + # operator compositions can optimize `derivative` and similar . + # See https://odlgroup.github.io/odl/math/derivatives_guide.html, + # subsection "Complex <-> Real mappings" for details. + super(RealPart, self).__init__(domain, range, linear=True) def _call(self, x): """Return ``self(x)``.""" @@ -1019,10 +1226,10 @@ def inverse(self): Examples -------- - The inverse is its own inverse if its domain is real: + The operator 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.]) @@ -1030,18 +1237,22 @@ 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.space_is_real: - return self + if self.domain.is_real and self.range.is_real: + # Identity case + return IdentityOperator(self.range, self.domain) + elif self.range.is_complex: + # Self case or odd case corresponding to ComplexEmbedding + return RealPart(self.range, self.domain) else: - return ComplexEmbedding(self.domain, scalar=1) + return ComplexEmbedding(self.range, self.domain, scalar=1) @property def adjoint(self): - """Return the (left) adjoint. + r"""Return the (left) adjoint. Notes ----- @@ -1049,19 +1260,19 @@ def adjoint(self): space, this does not satisfy the usual adjoint equation: .. math:: - \langle Ax, y \\rangle = \langle x, A^*y \\rangle + \langle Ax, y \rangle = \langle x, A^*y \rangle Instead it is an adjoint in a weaker sense as follows: .. math:: - \langle AA^*x, y \\rangle = \langle A^*x, A^*y \\rangle + \langle AA^*x, y \rangle = \langle A^*x, A^*y \rangle Examples -------- 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) @@ -1070,7 +1281,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) @@ -1078,10 +1289,22 @@ def adjoint(self): >>> AtAxy == AtxAty True """ - if self.space_is_real: - return self - else: - return ComplexEmbedding(self.domain, scalar=1) + return self.inverse + + 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): @@ -1093,35 +1316,42 @@ 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 - self.space_is_real = (space == real_space) - linear = True - super(ImagPart, self).__init__(space, real_space, linear=linear) + if range is None: + range = domain.real_space + + # `linear=True` is a compromise in favor of efficiency, since + # operator compositions can optimize `derivative` and similar . + # See https://odlgroup.github.io/odl/math/derivatives_guide.html, + # subsection "Complex <-> Real mappings" for details. + super(ImagPart, self).__init__(domain, range, linear=True) def _call(self, x): """Return ``self(x)``.""" @@ -1150,7 +1380,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.]) @@ -1158,18 +1388,22 @@ 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.space_is_real: - return ZeroOperator(self.domain) + if self.domain.is_real: + # Zero case + return ZeroOperator(self.range, self.domain) + elif self.domain.is_complex and self.range.is_complex: + # Self case + return ImagPart(self.range, self.domain) else: - return ComplexEmbedding(self.domain, scalar=1j) + return ComplexEmbedding(self.range, self.domain, scalar=1) @property def adjoint(self): - """Return the (left) adjoint. + r"""Return the (left) adjoint. Notes ----- @@ -1177,19 +1411,19 @@ def adjoint(self): space, this does not satisfy the usual adjoint equation: .. math:: - \langle Ax, y \\rangle = \langle x, A^*y \\rangle + \langle Ax, y \rangle = \langle x, A^*y \rangle Instead it is an adjoint in a weaker sense as follows: .. math:: - \langle AA^*x, y \\rangle = \langle A^*x, A^*y \\rangle + \langle AA^*x, y \rangle = \langle A^*x, A^*y \rangle Examples -------- 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) @@ -1198,7 +1432,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) @@ -1206,10 +1440,22 @@ def adjoint(self): >>> AtAxy == AtxAty True """ - if self.space_is_real: - return ZeroOperator(self.domain) - else: - return ComplexEmbedding(self.domain, scalar=1j) + return self.inverse + + 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): @@ -1221,14 +1467,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. @@ -1238,13 +1487,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]) @@ -1252,14 +1501,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)``.""" @@ -1281,29 +1530,36 @@ 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.]) """ + if self.scalar == 0: + return ZeroOperator(self.range, self.domain) + if self.domain.is_real: # Real domain # Optimizations for simple cases. if self.scalar.real == self.scalar: - return (1 / self.scalar.real) * RealPart(self.range) + return _scale_op(RealPart(self.range, self.domain), + 1 / self.scalar.real) elif 1j * self.scalar.imag == self.scalar: - return (1 / self.scalar.imag) * ImagPart(self.range) + return _scale_op(ImagPart(self.range, self.domain), + 1 / self.scalar.imag) else: # General case inv_scalar = (1 / self.scalar).conjugate() - return ((inv_scalar.real) * RealPart(self.range) + - (inv_scalar.imag) * ImagPart(self.range)) + return _lico_ops( + 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): - """Return the (right) adjoint. + r"""Return the (right) adjoint. Notes ----- @@ -1311,12 +1567,12 @@ def adjoint(self): space, this does not satisfy the usual adjoint equation: .. math:: - \langle Ax, y \\rangle = \langle x, A^*y \\rangle + \langle Ax, y \rangle = \langle x, A^*y \rangle Instead it is an adjoint in a weaker sense as follows: .. math:: - \langle A^*Ax, y \\rangle = \langle Ax, Ay \\rangle + \langle A^*Ax, y \rangle = \langle Ax, Ay \rangle Examples -------- @@ -1342,34 +1598,64 @@ def adjoint(self): >>> AtAxy == AxAy True """ + if self.scalar == 0: + return ZeroOperator(self.range, self.domain) + if self.domain.is_real: # Real domain # Optimizations for simple cases. if self.scalar.real == self.scalar: - return self.scalar.real * RealPart(self.range) + return _scale_op(self.scalar.real, + ComplexEmbedding(self.range, self.domain)) elif 1j * self.scalar.imag == self.scalar: - return self.scalar.imag * ImagPart(self.range) + return _scale_op(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 _lico_ops( + 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 modulus (absolute value) of a vector.""" - def __init__(self, space): + def __init__(self, domain, range=None): """Initialize a new instance. Parameters ---------- - space : `TensorSpace` - Space in which the modulus 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 -------- @@ -1395,12 +1681,17 @@ def __init__(self, space): >>> op([3 + 4j, 2]) uniform_discr(0.0, 1.0, 2).element([ 5., 2.]) """ - real_space = space.real_space - super(ComplexModulus, self).__init__(space, real_space, linear=False) + if range is None: + range = domain.real_space + super(ComplexModulus, self).__init__(domain, range, linear=False) 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) def derivative(self, x): r"""Return the derivative operator in the "C = R^2" sense. @@ -1551,18 +1842,41 @@ def adjoint(self): """Adjoint in the "C = R^2" sense.""" return deriv + def __repr__(self): + """Return ``repr(self)``.""" + return attribute_repr_string(repr(deriv), 'adjoint') + return ComplexModulusDerivativeAdjoint( deriv.range, deriv.domain, linear=deriv.domain.is_real) + def __repr__(self): + """Return ``repr(self)``.""" + return method_repr_string( + repr(op), 'derivative', arg_strs=[repr(x)]) + return ComplexModulusDerivative(op.domain, op.range, linear=op.domain.is_real) + 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) class ComplexModulusSquared(Operator): """Operator that computes the squared complex modulus (absolute value).""" - def __init__(self, space): + def __init__(self, domain, range=None): """Initialize a new instance. Parameters @@ -1595,9 +1909,10 @@ def __init__(self, space): >>> op([3 + 4j, 2]) uniform_discr(0.0, 1.0, 2).element([ 25., 4.]) """ - real_space = space.real_space + if range is None: + range = domain.real_space super(ComplexModulusSquared, self).__init__( - space, real_space, linear=False) + domain, range, linear=False) def _call(self, x): """Return ``self(x)``.""" @@ -1745,12 +2060,35 @@ def adjoint(self): """Adjoint in the "C = R^2" sense.""" return deriv + def __repr__(self): + """Return ``repr(self)``.""" + return attribute_repr_string(repr(deriv), 'adjoint') + return ComplexModulusSquaredDerivAdj( deriv.range, deriv.domain, linear=deriv.domain.is_real) + def __repr__(self): + """Return ``repr(self)``.""" + return method_repr_string( + repr(op), 'derivative', arg_strs=[repr(x)]) + return ComplexModulusSquaredDerivative(op.domain, op.range, linear=op.domain.is_real) + 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__': from odl.util.testutils import run_doctests diff --git a/odl/operator/operator.py b/odl/operator/operator.py index 33d9b2acb52..3f8c4f0c567 100644 --- a/odl/operator/operator.py +++ b/odl/operator/operator.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -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: norm = getattr(self, '__norm', None) if norm is not None: @@ -846,8 +849,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): @@ -928,8 +931,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): @@ -1021,8 +1024,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): @@ -1049,18 +1052,25 @@ 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)``. - The default string implementation. Should be overridden by - subclasses. + By default, ``str(self)`` is the same as ``repr(self)``. Subclasses + that need a different behavior should override this method. """ - return self.__class__.__name__ + return repr(self) - # 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 @@ -1103,11 +1113,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: @@ -1204,8 +1214,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)``.""" @@ -1234,8 +1246,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.]) @@ -1298,12 +1310,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): @@ -1321,7 +1336,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``. @@ -1436,8 +1451,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)``.""" @@ -1456,11 +1473,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' @@ -1510,12 +1525,10 @@ def derivative(self, x): def __repr__(self): """Return ``repr(self)``.""" - return '{}({!r}, {!r})'.format(self.__class__.__name__, - self.left, self.right) - - def __str__(self): - """Return ``str(self)``.""" - return repr(self) + posargs = [self.left, self.right] + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) class OperatorLeftScalarMult(Operator): @@ -1543,8 +1556,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.]) """ @@ -1601,8 +1614,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.]) """ @@ -1631,8 +1644,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.]) @@ -1660,8 +1673,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.]) """ @@ -1673,8 +1686,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)``.""" @@ -1709,9 +1725,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)): @@ -1788,9 +1804,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: @@ -1815,9 +1831,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.]) """ @@ -1841,9 +1857,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.]) """ @@ -1854,8 +1870,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)``.""" @@ -1893,7 +1912,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.]) """ @@ -1970,8 +1989,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)``.""" @@ -1991,9 +2013,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}' @@ -2083,8 +2105,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)``.""" @@ -2203,8 +2228,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 1922af41861..39869c804b1 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -8,12 +8,13 @@ """Convenience functions for operators.""" -from __future__ import print_function, division, absolute_import -from future.utils import native +from __future__ import absolute_import, division, print_function + import numpy as np +from future.utils import native from odl.space.base_tensors import TensorSpace -from odl.space import ProductSpace +from odl.space.pspace import ProductSpace from odl.util import nd_iterator from odl.util.testutils import noise_element diff --git a/odl/operator/pspace_ops.py b/odl/operator/pspace_ops.py index 149b260c556..f04930e1067 100644 --- a/odl/operator/pspace_ops.py +++ b/odl/operator/pspace_ops.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -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') @@ -116,9 +118,11 @@ def __init__(self, operators, domain=None, range=None): >>> 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: @@ -126,10 +130,12 @@ def __init__(self, operators, domain=None, range=None): >>> 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.] + ] + ) If a column is empty, the operator domain must be specified. The same holds for an empty row and the range of the operator: @@ -137,17 +143,21 @@ def __init__(self, operators, domain=None, range=None): >>> prod_op = odl.ProductSpaceOperator([[I, 0], ... [I, 0]], domain=r3 ** 2) >>> prod_op(x) - ProductSpace(rn(3), 2).element([ - [ 1., 2., 3.], - [ 1., 2., 3.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 1., 2., 3.], + [ 1., 2., 3.] + ] + ) >>> prod_op = odl.ProductSpaceOperator([[I, I], ... [0, 0]], range=r3 ** 2) >>> prod_op(x) - ProductSpace(rn(3), 2).element([ - [ 5., 7., 9.], - [ 0., 0., 0.] - ]) + ProductSpace(rn(3), 2).element( + [ + [ 5., 7., 9.], + [ 0., 0., 0.] + ] + ) """ # Lazy import to improve `import odl` time import scipy.sparse @@ -318,7 +328,7 @@ def derivative(self, x): Parameters ---------- x : `domain` element - The point to take the derivative in + Point in which to take the derivative. Returns ------- @@ -330,44 +340,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 @@ -398,7 +419,7 @@ def adjoint(self): Returns ------- adjoint : `ProductSpaceOperator` - The adjoint + The adjoint operator. Examples -------- @@ -410,18 +431,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 @@ -456,9 +482,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] @@ -512,11 +538,65 @@ def size(self): return np.prod(self.shape, dtype='int64') 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): @@ -572,18 +652,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): @@ -596,52 +683,7 @@ 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. + """Adjoint operator extending by zeros from subspace to full space. Examples -------- @@ -655,7 +697,7 @@ def __init__(self, space, index): Projection on the 0-th component: - >>> proj_adj = odl.ComponentProjectionAdjoint(pspace, 0) + >>> proj_adj = odl.ComponentProjection(pspace, 0).adjoint >>> proj_adj(x[0]) ProductSpace(rn(1), rn(2), rn(3)).element([ [ 1.], @@ -663,9 +705,9 @@ def __init__(self, space, index): [ 0., 0., 0.] ]) - Projection on a sub-space corresponding to indices 0 and 2: + Projection on a subspace corresponding to indices 0 and 2: - >>> proj_adj = odl.ComponentProjectionAdjoint(pspace, [0, 2]) + >>> proj_adj = odl.ComponentProjection(pspace, [0, 2]).adjoint >>> proj_adj(x[[0, 2]]) ProductSpace(rn(1), rn(2), rn(3)).element([ [ 1.], @@ -673,36 +715,39 @@ def __init__(self, space, index): [ 4., 5., 6.] ]) """ - self.__index = index - super(ComponentProjectionAdjoint, self).__init__( - space[index], space, linear=True) + op = self - @property - def index(self): - """Index of the subspace.""" - return self.__index + class ComponentProjectionAdjoint(Operator): - def _call(self, x, out=None): - """Extend ``x`` from the subspace.""" - if out is None: - out = self.range.zero() - else: - out.set_zero() + """Adjoint operator to `ComponentProjection`.""" - out[self.index] = x - return out + def _call(self, x, out=None): + """Extend ``x`` to the full space.""" + if out is None: + out = self.range.zero() + else: + out.set_zero() - @property - def adjoint(self): - """Adjoint of this operator. + out[op.index] = x + return out - Returns - ------- - adjoint : `ComponentProjection` - The adjoint is given by the `ComponentProjection` related to this - operator's `index`. - """ - return ComponentProjection(self.range, self.index) + @property + def adjoint(self): + """Adjoint of the adjoint, the original operator.""" + return op + + def __repr__(self): + """Return ``repr(self)``. + + 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') + + return ComponentProjectionAdjoint(self.range, self.domain, linear=True) def __repr__(self): """Return ``repr(self)``. @@ -710,11 +755,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): @@ -756,10 +802,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: @@ -813,12 +861,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 -------- @@ -828,22 +876,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]) @@ -871,19 +923,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): @@ -1045,10 +1102,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]) @@ -1058,19 +1117,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): @@ -1120,12 +1184,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 @@ -1266,19 +1333,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 2ba85bce13a..f5a9b1f9f38 100644 --- a/odl/operator/tensor_ops.py +++ b/odl/operator/tensor_ops.py @@ -8,19 +8,21 @@ """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 packaging.version import parse as parse_version 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', @@ -110,7 +112,7 @@ class PointwiseNorm(PointwiseTensorFieldOperator): ``d``. """ - def __init__(self, vfspace, exponent=None, weighting=None): + def __init__(self, vfspace, exponent=None, range=None, weighting=None): """Initialize a new instance. Parameters @@ -124,6 +126,10 @@ def __init__(self, vfspace, exponent=None, weighting=None): 0 and 1 are currently not supported due to numerical instability. Default: ``vfspace.exponent`` + range : `TensorSpace`, optional + Scalar space to which the operator maps, must be compatible + with ``vfspace``'s shape. + Default: ``vfspace[0]`` weighting : `array-like` or positive float, optional Weighting array or constant for the norm. If an array is given, its length must be equal to ``len(domain)``, and @@ -139,7 +145,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 @@ -148,27 +154,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(kohr-h): 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 @@ -177,15 +185,16 @@ 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 + # TODO(kohr-h): find a more robust way of getting the weights as + # an array if hasattr(self.domain.weighting, 'array'): self.__weights = self.domain.weighting.array elif hasattr(self.domain.weighting, 'const'): @@ -211,7 +220,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): @@ -357,25 +366,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 @@ -385,22 +424,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. @@ -411,7 +467,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: @@ -433,7 +489,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): @@ -445,85 +501,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] @@ -548,101 +531,85 @@ 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 + 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 - 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 :: + class PointwiseInnerAdjoint(PointwiseTensorFieldOperator): - x --> h(x) * (w / v) * G(x), + """Adjoint of the point-wise inner product operator.""" - where ``G`` and ``w`` are the vector field and weighting from the - inner product operator, resp., and all multiplications are understood - component-wise. - """ + 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 - def __init__(self, sspace, vecfield, vfspace=None, weighting=None): - """Initialize a new instance. + @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(self.range, self.domain, + base_space=self.base_space, linear=True) - @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. @@ -685,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 @@ -700,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): @@ -745,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 @@ -759,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) @@ -936,7 +926,6 @@ def _call(self, x, out=None): self.matrix.dot(x, out=out_arr) else: # Could use einsum to have out, but it's damn slow - # TODO: investigate speed issue dot = np.tensordot(self.matrix, x, axes=(1, self.axis)) # New axis ends up as first, need to move it to its place out[:] = moveaxis(dot, 0, self.axis) @@ -944,7 +933,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 @@ -967,13 +971,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)) - - def __str__(self): - """Return ``str(self)``.""" - return repr(self) + inner_parts = signature_string_parts(posargs, optargs, + mod=[['!s'], ['!r', '!r', '']]) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def _normalize_sampling_points(sampling_points, ndim): @@ -1199,16 +1200,25 @@ def adjoint(self): variant) def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(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' + ) + """ 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)) - - def __str__(self): - """Return ``str(self)``.""" - return repr(self) + inner_parts = signature_string_parts(posargs, optargs, mod=['!r', '']) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) class WeightedSumSamplingOperator(Operator): @@ -1397,16 +1407,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)) - - def __str__(self): - """Return ``str(self)``.""" - return repr(self) + inner_parts = signature_string_parts(posargs, optargs, mod=['!r', '']) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) class FlatteningOperator(Operator): @@ -1488,8 +1508,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): @@ -1515,7 +1534,6 @@ def inverse(self): True """ op = self - scaling = getattr(self.domain, 'cell_volume', 1.0) class FlatteningOperatorInverse(Operator): @@ -1526,45 +1544,51 @@ class FlatteningOperatorInverse(Operator): FlatteningOperatorInverse(x) == reshape(x, orig_shape) """ - def __init__(self): - """Initialize a new instance.""" - super(FlatteningOperatorInverse, self).__init__( - op.range, op.domain, linear=True) - def _call(self, x): """Reshape ``x`` back to n-dim. shape.""" return np.reshape(x.asarray(), self.range.shape, 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 ``repr(self)``. - def __str__(self): - """Return ``str(self)``.""" - return repr(self) + Examples + -------- + >>> space = odl.uniform_discr([-1, -1], [1, 1], shape=(2, 3)) + >>> op = odl.FlatteningOperator(space) + >>> op.inverse + FlatteningOperator( + uniform_discr([-1., -1.], [ 1., 1.], (2, 3)) + ).inverse - return FlatteningOperatorInverse() + """ + return attribute_repr_string(repr(op), 'inverse') + + return FlatteningOperatorInverse(self.range, self.domain, linear=True) 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)) - - def __str__(self): - """Return ``str(self)``.""" - return repr(self) + inner_parts = signature_string_parts(posargs, optargs, mod=['!r', '']) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) def is_compatible_space(space, base_space): @@ -1626,12 +1650,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/space/__init__.py b/odl/space/__init__.py index 36e99548a05..e27b9986f5c 100644 --- a/odl/space/__init__.py +++ b/odl/space/__init__.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # diff --git a/odl/space/base_tensors.py b/odl/space/base_tensors.py index 5e620baaafb..aa48e834010 100644 --- a/odl/space/base_tensors.py +++ b/odl/space/base_tensors.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -8,20 +8,21 @@ """Base classes for implementations of tensor spaces.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + from builtins import object from numbers import Integral + import numpy as np -from odl.set.sets import RealNumbers, ComplexNumbers +from odl.set.sets import ComplexNumbers, RealNumbers from odl.set.space import LinearSpace, LinearSpaceElement 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, method_repr_string, is_complex_floating_dtype, + is_floating_dtype, is_numeric_dtype, is_real_dtype, is_real_floating_dtype, + repr_string, safe_int_conv, signature_string_parts, writable_array) from odl.util.ufuncs import TensorSpaceUfuncs -from odl.util.utility import TYPE_MAP_R2C, TYPE_MAP_C2R - +from odl.util.utility import TYPE_MAP_C2R, TYPE_MAP_R2C __all__ = ('TensorSpace',) @@ -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. @@ -399,12 +415,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, [])) - - def __str__(self): - """Return ``str(self)``.""" - return repr(self) + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts) @property def examples(self): @@ -427,7 +439,8 @@ def examples(self): self.element(np.random.uniform(size=self.shape) + np.random.uniform(size=self.shape) * 1j)) else: - # TODO: return something that always works, like zeros or ones? + # TODO(kohr-h): return something that always works, like zeros + # or ones? raise NotImplementedError('no examples available for non-numeric' 'data type') @@ -627,13 +640,18 @@ def astype(self, dtype): raise NotImplementedError('abstract method') def __repr__(self): - """Return ``repr(self)``.""" + """Return ``repr(self)``. + + Examples + -------- + >>> rn = odl.rn(3) + >>> x = rn.element([1, 2, 3]) + >>> x + rn(3).element([ 1., 2., 3.]) + """ 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 method_repr_string(repr(self.space), 'element', [self_str]) def __str__(self): """Return ``str(self)``.""" diff --git a/odl/space/entry_points.py b/odl/space/entry_points.py index fe1fc7644f8..62498ce6721 100644 --- a/odl/space/entry_points.py +++ b/odl/space/entry_points.py @@ -20,7 +20,7 @@ NumpyTensorSpace : Numpy-based implementation of `TensorSpace` """ -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function from odl.space.npy_tensors import NumpyTensorSpace diff --git a/odl/space/fspace.py b/odl/space/fspace.py index dfbe9042141..05733c318c5 100644 --- a/odl/space/fspace.py +++ b/odl/space/fspace.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -8,22 +8,23 @@ """Spaces of scalar-, vector- and tensor-valued functions on a given domain.""" -from __future__ import print_function, division, absolute_import -from builtins import object +from __future__ import absolute_import, division, print_function + import inspect -import numpy as np import sys import warnings +from builtins import object + +import numpy as np -from odl.set import RealNumbers, ComplexNumbers, Set, LinearSpace +from odl.set import ComplexNumbers, LinearSpace, RealNumbers, Set from odl.set.space import LinearSpaceElement from odl.util import ( - is_real_dtype, is_complex_floating_dtype, dtype_repr, dtype_str, - complex_dtype, real_dtype, signature_string, is_real_floating_dtype, - is_valid_input_array, is_valid_input_meshgrid, - out_shape_from_array, out_shape_from_meshgrid, vectorize, writable_array) -from odl.util.utility import preload_first_arg, getargspec - + complex_dtype, dtype_repr, dtype_str, is_complex_floating_dtype, + is_real_dtype, is_real_floating_dtype, is_valid_input_array, + is_valid_input_meshgrid, out_shape_from_array, out_shape_from_meshgrid, + real_dtype, signature_string, vectorize, writable_array) +from odl.util.utility import getargspec, preload_first_arg __all__ = ('FunctionSpace',) @@ -1018,7 +1019,7 @@ def examples(self): Gaussian Linear gradients """ - # TODO: adapt for tensor-valued functions + # TODO(kohr-h): adapt for tensor-valued functions # Get the points and calculate some statistics on them mins = self.domain.min() @@ -1096,10 +1097,6 @@ def __repr__(self): inner_str = signature_string(posargs, optargs, mod=['!r', optmod]) return '{}({})'.format(self.__class__.__name__, inner_str) - def __str__(self): - """Return ``str(self)``.""" - return repr(self) - class FunctionSpaceElement(LinearSpaceElement): diff --git a/odl/space/npy_tensors.py b/odl/space/npy_tensors.py index 5b122a720fd..f5a423c522e 100644 --- a/odl/space/npy_tensors.py +++ b/odl/space/npy_tensors.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -8,23 +8,24 @@ """NumPy implementation of tensor spaces.""" -from __future__ import print_function, division, absolute_import -from future.utils import native -from builtins import object +from __future__ import absolute_import, division, print_function + import ctypes +from builtins import object from functools import partial + import numpy as np +from future.utils import native -from odl.set.sets import RealNumbers, ComplexNumbers +from odl.set.sets import ComplexNumbers, RealNumbers from odl.set.space import LinearSpaceTypeError -from odl.space.base_tensors import TensorSpace, Tensor +from odl.space.base_tensors import Tensor, TensorSpace from odl.space.weighting import ( - Weighting, ArrayWeighting, ConstWeighting, - CustomInner, CustomNorm, CustomDist) + ArrayWeighting, ConstWeighting, CustomDist, CustomInner, CustomNorm, + Weighting) from odl.util import ( - dtype_str, signature_string, is_real_dtype, is_numeric_dtype, - writable_array, is_floating_dtype) - + attribute_repr_string, dtype_str, is_floating_dtype, is_numeric_dtype, + is_real_dtype, repr_string, signature_string_parts, writable_array) __all__ = ('NumpyTensorSpace',) @@ -807,7 +808,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 +838,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): @@ -1794,7 +1797,7 @@ def _blas_is_applicable(*args): return False elif any(x.size > np.iinfo('int32').max for x in args): # Temporary fix for 32 bit int overflow in BLAS - # TODO: use chunking instead + # TODO(#1200): use chunking instead return False else: return True @@ -2164,7 +2167,7 @@ def norm(self, x): The norm of the provided tensor. """ if self.exponent == 2.0: - norm_squared = self.inner(x, x).real # TODO: optimize?! + norm_squared = self.inner(x, x).real if norm_squared < 0: norm_squared = 0.0 # Compensate for numerical error return float(np.sqrt(norm_squared)) diff --git a/odl/space/pspace.py b/odl/space/pspace.py index 663b41c2aab..7630592f053 100644 --- a/odl/space/pspace.py +++ b/odl/space/pspace.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -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, method_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 ----- @@ -462,7 +465,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 @@ -472,15 +475,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: @@ -533,7 +538,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 @@ -561,7 +566,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 @@ -613,13 +618,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) @@ -739,7 +744,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)``.""" @@ -749,40 +754,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): @@ -1094,24 +1088,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): @@ -1123,10 +1123,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 @@ -1311,62 +1313,64 @@ def conj(self): complex_conj = [part.conj() for part in self.parts] return self.space.element(complex_conj) - def __str__(self): - """Return ``str(self)``.""" - return repr(self) - def __repr__(self): """Return ``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 method_repr_string(repr(self.space), 'element', [data_str]) def show(self, title=None, indices=None, **kwargs): """Display the parts of this product space element graphically. @@ -1508,25 +1512,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']: @@ -1635,7 +1647,7 @@ def norm(self, x): The norm of the provided element. """ if self.exponent == 2.0: - norm_squared = self.inner(x, x).real # TODO: optimize?! + norm_squared = self.inner(x, x).real return np.sqrt(norm_squared) else: norms = np.fromiter( @@ -1740,7 +1752,7 @@ def norm(self, x): The norm of the element. """ if self.exponent == 2.0: - norm_squared = self.inner(x, x).real # TODO: optimize?! + norm_squared = self.inner(x, x).real return np.sqrt(norm_squared) else: norms = np.fromiter( @@ -1859,17 +1871,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/space/space_utils.py b/odl/space/space_utils.py index 7beb9dff672..be1efe79d13 100644 --- a/odl/space/space_utils.py +++ b/odl/space/space_utils.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -8,13 +8,13 @@ """Utility functions for space implementations.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + import numpy as np -from odl.set import RealNumbers, ComplexNumbers +from odl.set import ComplexNumbers, RealNumbers from odl.space.entry_points import tensor_space_impl - __all__ = ('vector', 'tensor_space', 'cn', 'rn') diff --git a/odl/space/weighting.py b/odl/space/weighting.py index 8fcabbf50d8..898d7250a38 100644 --- a/odl/space/weighting.py +++ b/odl/space/weighting.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -8,13 +8,14 @@ """Weightings for finite-dimensional spaces.""" -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.space.base_tensors import TensorSpace -from odl.util import array_str, signature_string, indent - +from odl.util import array_str, indent, signature_string __all__ = ('MatrixWeighting', 'ArrayWeighting', 'ConstWeighting', 'CustomInner', 'CustomNorm', 'CustomDist') @@ -149,6 +150,10 @@ def dist(self, x1, x2): """ return self.norm(x1 - x2) + def __str__(self): + """Return ``str(self)``.""" + return repr(self) + class MatrixWeighting(Weighting): @@ -212,7 +217,7 @@ def __init__(self, matrix, impl, exponent=2.0, **kwargs): # Lazy import to improve `import odl` time import scipy.sparse - # TODO: fix dead link `scipy.sparse.spmatrix` + # TODO(kohr-h): fix dead link `scipy.sparse.spmatrix` precomp_mat_pow = kwargs.pop('precomp_mat_pow', False) self._cache_mat_pow = bool(kwargs.pop('cache_mat_pow', True)) self._cache_mat_decomp = bool(kwargs.pop('cache_mat_decomp', False)) @@ -317,7 +322,7 @@ def matrix_decomp(self, cache=None): import scipy.linalg import scipy.sparse - # TODO: fix dead link `scipy.linalg.decomp.eigh` + # TODO(kohr-h): fix dead link `scipy.linalg.decomp.eigh` if scipy.sparse.isspmatrix(self.matrix): raise NotImplementedError('sparse matrix not supported') @@ -450,10 +455,6 @@ def __repr__(self): mod=['!s', '']) return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) - def __str__(self): - """Return ``str(self)``.""" - return repr(self) - class ArrayWeighting(Weighting): @@ -567,10 +568,6 @@ def __repr__(self): mod=['!s', ':.4']) return '{}(\n{}\n)'.format(self.__class__.__name__, indent(inner_str)) - def __str__(self): - """Return ``str(self)``.""" - return repr(self) - class ConstWeighting(Weighting): @@ -654,10 +651,6 @@ def __repr__(self): return '{}({})'.format(self.__class__.__name__, signature_string(posargs, optargs)) - def __str__(self): - """Return ``str(self)``.""" - return repr(self) - class CustomInner(Weighting): @@ -719,8 +712,7 @@ def repr_part(self): def __repr__(self): """Return ``repr(self)``.""" posargs = [self.inner] - optargs = [] - inner_str = signature_string(posargs, optargs, mod='!r') + inner_str = signature_string(posargs, [], mod='!r') return '{}({})'.format(self.__class__.__name__, inner_str) From 8bf750ad44f081c241609124660962136dcadf80 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Wed, 20 Jun 2018 00:34:37 +0200 Subject: [PATCH 2/9] ENH: fix discr repr and operator tests --- odl/discr/__init__.py | 2 +- odl/discr/diff_ops.py | 226 ++++++++++++++------------ odl/discr/discr_mappings.py | 146 ++++++++++++----- odl/discr/discr_ops.py | 207 +++++++++++++---------- odl/discr/discretization.py | 11 +- odl/discr/grid.py | 157 +++++++----------- odl/discr/lp_discr.py | 113 ++++++++----- odl/discr/partition.py | 81 ++++----- odl/operator/default_ops.py | 6 +- odl/operator/pspace_ops.py | 217 +++++++++++-------------- odl/operator/tensor_ops.py | 39 +++-- odl/space/pspace.py | 85 ++++------ odl/test/discr/diff_ops_test.py | 10 +- odl/test/discr/discr_mappings_test.py | 9 +- odl/test/discr/discr_ops_test.py | 8 +- odl/test/discr/grid_test.py | 8 +- odl/test/discr/lp_discr_test.py | 4 +- odl/test/discr/partition_test.py | 8 +- odl/test/operator/operator_test.py | 21 +-- odl/test/operator/oputils_test.py | 5 +- odl/test/operator/pspace_ops_test.py | 4 +- odl/test/operator/tensor_ops_test.py | 10 +- odl/util/utility.py | 16 +- 23 files changed, 741 insertions(+), 652 deletions(-) diff --git a/odl/discr/__init__.py b/odl/discr/__init__.py index 380483fc761..c60c2b853f0 100644 --- a/odl/discr/__init__.py +++ b/odl/discr/__init__.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # diff --git a/odl/discr/diff_ops.py b/odl/discr/diff_ops.py index 49f5cee287b..17aab220562 100644 --- a/odl/discr/diff_ops.py +++ b/odl/discr/diff_ops.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -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.]] @@ -137,7 +135,6 @@ def _call(self, x, out=None): if out is None: out = self.range.element() - # TODO: this pipes CUDA arrays through NumPy. Write native operator. with writable_array(out) as out_arr: finite_diff(x.asarray(), axis=self.axis, dx=self.dx, method=self.method, pad_mode=self.pad_mode, @@ -177,40 +174,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 +250,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 +396,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 +470,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 +488,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 +606,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 +642,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 +678,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 +693,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 +703,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 +767,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 99a01258ac6..69a55541288 100644 --- a/odl/discr/discr_mappings.py +++ b/odl/discr/discr_mappings.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -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'], @@ -382,7 +396,7 @@ def variant(self): def _call(self, x, out=None): """Return ``self(x[, out])``.""" - # TODO: pass reasonable options on to the interpolator + # TODO(kohr-h): pass reasonable options on to the interpolator def nearest(arg, out=None): """Interpolating function with vectorization.""" if is_valid_input_meshgrid(arg, self.grid.ndim): @@ -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): @@ -437,7 +465,7 @@ def __init__(self, fspace, partition, tspace): def _call(self, x, out=None): """Return ``self(x[, out])``.""" - # TODO: pass reasonable options on to the interpolator + # TODO(kohr-h): pass reasonable options on to the interpolator def linear(arg, out=None): """Interpolating function with vectorization.""" if is_valid_input_meshgrid(arg, self.grid.ndim): @@ -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): @@ -929,7 +985,7 @@ def _evaluate(self, indices, norm_distances, out=None): for lo_hi, edge in zip(product(*([['l', 'h']] * len(indices))), product(*edge_indices)): weight = 1.0 - # TODO: determine best summation order from array strides + # TODO(kohr-h): determine best summation order from array strides for lh, w_lo, w_hi in zip(lo_hi, low_weights, high_weights): # We don't multiply in-place to exploit the cheap operations diff --git a/odl/discr/discr_ops.py b/odl/discr/discr_ops.py index cf2f2d1c858..6ec5a280e78 100644 --- a/odl/discr/discr_ops.py +++ b/odl/discr/discr_ops.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -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`. @@ -400,15 +416,15 @@ class ResizingOperatorAdjoint(ResizingOperatorBase): 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 +438,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(self.range, self.domain, linear=True) @property def inverse(self): @@ -437,6 +456,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 14c5ec8515d..51ea75b15da 100644 --- a/odl/discr/discretization.py +++ b/odl/discr/discretization.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -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 fea4aa3a6a7..0795e54eb10 100644 --- a/odl/discr/grid.py +++ b/odl/discr/grid.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -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.]) @@ -254,7 +242,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 @@ -270,7 +258,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.]) """ @@ -282,7 +270,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.]) """ @@ -326,7 +314,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.]) @@ -357,7 +345,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.]) @@ -407,13 +395,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.]) """ @@ -437,7 +425,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.]) """ @@ -458,7 +446,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.]) """ @@ -488,8 +476,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) @@ -520,7 +508,7 @@ def __eq__(self, other): def __hash__(self): """Return ``hash(self)``.""" - # TODO: update with #841 + # TODO(doable, #841): update (what exactly) coord_vec_str = tuple(cv.tobytes() for cv in self.coord_vectors) return hash((type(self), coord_vec_str)) @@ -537,7 +525,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 @@ -663,15 +651,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.], @@ -725,15 +708,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.], @@ -765,12 +743,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) @@ -798,7 +773,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.], @@ -841,7 +816,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)) """ @@ -871,7 +846,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.], @@ -902,7 +877,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.], @@ -931,43 +906,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 """ @@ -1005,7 +960,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: @@ -1025,23 +980,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 54a849b2ff4..4b13c061763 100644 --- a/odl/discr/lp_discr.py +++ b/odl/discr/lp_discr.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -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 @@ -388,7 +398,7 @@ def _astype(self, dtype): def _inner(self, x, y): """Return ``self.inner(x, y)``.""" if self.is_uniform and not self.is_uniformly_weighted: - # TODO: implement without copying x + # TODO(kohr-h): implement without copying x bdry_fracs = self.partition.boundary_cell_fractions func_list = _scaling_func_list(bdry_fracs, exponent=1.0) x_arr = apply_on_boundary(x, func=func_list, only_once=False) @@ -399,7 +409,7 @@ def _inner(self, x, y): def _norm(self, x): """Return ``self.norm(x)``.""" if self.is_uniform and not self.is_uniformly_weighted: - # TODO: implement without copying x + # TODO(kohr-h): implement without copying x bdry_fracs = self.partition.boundary_cell_fractions func_list = _scaling_func_list(bdry_fracs, exponent=self.exponent) x_arr = apply_on_boundary(x, func=func_list, only_once=False) @@ -420,7 +430,7 @@ def _dist(self, x, y): else: return super(DiscreteLp, self)._dist(x, y) - # TODO: add byaxis_out when discretized tensor-valued functions are + # TODO(kohr-h): add byaxis_out when discretized tensor-valued functions are # available @property @@ -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,26 +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']) - - return '{}({})'.format(ctor, indent(inner_str)) + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, optargs) - def __str__(self): - """Return ``str(self)``.""" - return repr(self) + return repr_string(ctor, inner_parts, allow_mixed_seps=True) @property def element_type(self): @@ -726,7 +763,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) @@ -1392,7 +1429,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 @@ -1458,7 +1495,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 @@ -1533,7 +1570,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 @@ -1607,7 +1644,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 @@ -1619,7 +1656,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 @@ -1878,9 +1915,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(kohr-h): 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 21cc6727fc1..831cc5b7c8e 100644 --- a/odl/discr/partition.py +++ b/odl/discr/partition.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -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', @@ -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: @@ -848,12 +847,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([], [], ())' @@ -877,16 +891,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 @@ -914,9 +928,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') @@ -924,16 +938,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)``.""" @@ -1339,18 +1355,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 @@ -1361,19 +1372,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 7cb960b7813..2be2611b88e 100644 --- a/odl/operator/default_ops.py +++ b/odl/operator/default_ops.py @@ -1396,10 +1396,10 @@ def inverse(self): # Zero case return ZeroOperator(self.range, self.domain) elif self.domain.is_complex and self.range.is_complex: - # Self case - return ImagPart(self.range, self.domain) + # "Self" case + return 1j * ImagPart(self.range, self.domain) else: - return ComplexEmbedding(self.range, self.domain, scalar=1) + return ComplexEmbedding(self.range, self.domain, scalar=1j) @property def adjoint(self): diff --git a/odl/operator/pspace_ops.py b/odl/operator/pspace_ops.py index f04930e1067..6a3eff749aa 100644 --- a/odl/operator/pspace_ops.py +++ b/odl/operator/pspace_ops.py @@ -26,7 +26,7 @@ class ProductSpaceOperator(Operator): - """A "matrix of operators" on product spaces. + r"""A "matrix of operators" on product spaces. For example a matrix of operators can act on a vector by @@ -39,49 +39,46 @@ class ProductSpaceOperator(Operator): as a linear combination of "sub-operators", e.g. .. math:: - \\left( - \\begin{array}{ccc} - A & B & 0 \\\\ - 0 & C & 0 \\\\ + \left( + \begin{array}{ccc} + A & B & 0 \\ + 0 & C & 0 \\ 0 & 0 & D - \end{array}\\right) - \\left( - \\begin{array}{c} - x \\\\ - y \\\\ + \end{array}\right) + \left( + \begin{array}{c} + x \\ + y \\ z - \end{array}\\right) + \end{array}\right) = - \\left( - \\begin{array}{c} - A(x) + B(y) \\\\ - C(y) \\\\ + \left( + \begin{array}{c} + A(x) + B(y) \\ + C(y) \\ D(z) - \end{array}\\right) + \end{array}\right) Mathematically, a `ProductSpaceOperator` is an operator .. math:: - \mathcal{A}: \mathcal{X} \\to \mathcal{Y} + A: X \to Y between product spaces - :math:`\mathcal{X}=\mathcal{X}_1 \\times\dots\\times \mathcal{X}_m` - and - :math:`\mathcal{Y}=\mathcal{Y}_1 \\times\dots\\times \mathcal{Y}_n` - which can be written in the form + :math:`X = X_1 \times\dots\times X_m` and + :math:`Y = Y_1 \times\dots\times Y_n`, which can be written in the form .. math:: - \mathcal{A} = (\mathcal{A}_{ij})_{i,j}, \quad - i = 1, \dots, n, \\ j = 1, \dots, m + A = (A_{ij})_{i,j}, \quad i = 1, \dots, n, \ j = 1, \dots, m with *component operators* - :math:`\mathcal{A}_{ij}: \mathcal{X}_j \\to \mathcal{Y}_i`. + :math:`A_{ij}: X_j \to Y_i`. Its action on a vector :math:`x = (x_1, \dots, x_m)` is defined as the matrix multiplication .. math:: - [\mathcal{A}(x)]_i = \sum_{j=1}^m \mathcal{A}_{ij}(x_j). + [A(x)]_i = \sum_{j=1}^m A_{ij}(x_j). See Also -------- @@ -118,11 +115,9 @@ def __init__(self, operators, domain=None, range=None): >>> 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: @@ -130,12 +125,10 @@ def __init__(self, operators, domain=None, range=None): >>> 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.] + ]) If a column is empty, the operator domain must be specified. The same holds for an empty row and the range of the operator: @@ -143,21 +136,17 @@ def __init__(self, operators, domain=None, range=None): >>> prod_op = odl.ProductSpaceOperator([[I, 0], ... [I, 0]], domain=r3 ** 2) >>> prod_op(x) - ProductSpace(rn(3), 2).element( - [ - [ 1., 2., 3.], - [ 1., 2., 3.] - ] - ) + ProductSpace(rn(3), 2).element([ + [ 1., 2., 3.], + [ 1., 2., 3.] + ]) >>> prod_op = odl.ProductSpaceOperator([[I, I], ... [0, 0]], range=r3 ** 2) >>> prod_op(x) - ProductSpace(rn(3), 2).element( - [ - [ 5., 7., 9.], - [ 0., 0., 0.] - ] - ) + ProductSpace(rn(3), 2).element([ + [ 5., 7., 9.], + [ 0., 0., 0.] + ]) """ # Lazy import to improve `import odl` time import scipy.sparse @@ -349,19 +338,15 @@ def derivative(self, x): ... [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 @@ -373,22 +358,18 @@ def derivative(self, x): 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 @@ -435,19 +416,15 @@ def adjoint(self): ... [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 @@ -601,20 +578,20 @@ def __repr__(self): class ComponentProjection(Operator): - """Projection onto the subspace identified by an index. + r"""Projection onto the subspace identified by an index. - For a product space :math:`\mathcal{X} = \mathcal{X}_1 \\times \dots - \\times \mathcal{X}_n`, the component projection + For a product space :math:`X = X_1 \times \dots + \times X_n`, the component projection .. math:: - \mathcal{P}_i: \mathcal{X} \\to \mathcal{X}_i + \mathcal{P}_i: X \to X_i is given by :math:`\mathcal{P}_i(x) = x_i` for an element - :math:`x = (x_1, \dots, x_n) \\in \mathcal{X}`. + :math:`x = (x_1, \dots, x_n) \in X`. More generally, for an index set :math:`I \subset \{1, \dots, n\}`, the projection operator :math:`\mathcal{P}_I` is defined by - :math:`\mathcal{P}_I(x) = (x_i)_{i \\in I}`. + :math:`\mathcal{P}_I(x) = (x_i)_{i \in I}`. Note that this is a special case of a product space operator where the "operator matrix" has only one row and contains only @@ -652,12 +629,10 @@ 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.] + ]) """ if np.isscalar(index): self.__index = int(index) @@ -802,12 +777,10 @@ 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: @@ -880,22 +853,18 @@ def derivative(self, x): >>> 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]) @@ -1102,12 +1071,10 @@ 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]) @@ -1184,12 +1151,10 @@ 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.] + ]) The diagonal operator can also be created using a multiple of a single operator: diff --git a/odl/operator/tensor_ops.py b/odl/operator/tensor_ops.py index f5a9b1f9f38..5879b0187e4 100644 --- a/odl/operator/tensor_ops.py +++ b/odl/operator/tensor_ops.py @@ -546,25 +546,34 @@ def adjoint(self): >>> 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.]] - ] - ) + ProductSpace(uniform_discr([-1., -1.], [ 1., 1.], (1, 2)), 2).element([ + [[ 0., 2.]], + [[ 1., -2.]] + ]) """ op = self + # Get weighting as array + if hasattr(op.domain.weighting, 'array'): + adj_ran_weights = op.domain.weighting.array + elif hasattr(op.domain.weighting, 'const'): + adj_ran_weights = [op.domain.weighting.const] * len(op.domain) + else: + raise ValueError('domain weighting scheme {!r} does not define a ' + 'weighting array or constant' + ''.format(op.domain.weighting)) + class PointwiseInnerAdjoint(PointwiseTensorFieldOperator): """Adjoint of the point-wise inner product operator.""" def _call(self, f, out): """Implement ``self(vf, out)``.""" - for vfi, oi, wi in zip(op.vecfield, out, op.weights): + for vfi, oi, ran_wi, dom_wi in zip( + op.vecfield, out, adj_ran_weights, op.weights): vfi.multiply(f, out=oi) - if not np.isclose(wi, 1.0): - oi *= wi + if not np.isclose(ran_wi, dom_wi): + oi *= dom_wi / ran_wi @property def adjoint(self): @@ -592,12 +601,10 @@ def __repr__(self): >>> pw_inner PointwiseInner( ProductSpace(rn((1, 2)), 2), - ProductSpace(rn((1, 2)), 2).element( - [ - [[ 0., 1.]], - [[ 1., -1.]] - ] - ) + ProductSpace(rn((1, 2)), 2).element([ + [[ 0., 1.]], + [[ 1., -1.]] + ]) ) """ posargs = [self.domain, self.vecfield] @@ -1508,7 +1515,7 @@ def adjoint(self): >>> abs(op.adjoint(op(x)).inner(x) - op(x).inner(op(x))) < 1e-10 True """ - return self.inverse + return (1 / self.domain.cell_volume) * self.inverse @property def inverse(self): diff --git a/odl/space/pspace.py b/odl/space/pspace.py index 7630592f053..b7eb03e5e39 100644 --- a/odl/space/pspace.py +++ b/odl/space/pspace.py @@ -480,12 +480,10 @@ def element(self, inp=None, cast=True): >>> 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: @@ -1088,30 +1086,24 @@ 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): @@ -1123,12 +1115,10 @@ 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 @@ -1325,12 +1315,10 @@ def __repr__(self): 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.] + ]) Nested spaces work as well: @@ -1340,18 +1328,16 @@ def __repr__(self): ... [[1, 2], ... [3, 4, 5]]]) >>> x - ProductSpace(ProductSpace(rn(2), rn(3)), 2).element( - [ - [ - [ 1., 2.], - [ 3., 4., 5.] - ], - [ - [ 1., 2.], - [ 3., 4., 5.] - ] - ] - ) + ProductSpace(ProductSpace(rn(2), rn(3)), 2).element([ + [ + [ 1., 2.], + [ 3., 4., 5.] + ], + [ + [ 1., 2.], + [ 3., 4., 5.] + ] + ]) """ data_str = '[\n' edgeitems = np.get_printoptions()['edgeitems'] @@ -1370,7 +1356,8 @@ def __repr__(self): data_str += '\n]' - return method_repr_string(repr(self.space), 'element', [data_str]) + return method_repr_string(repr(self.space), 'element', [data_str], + force_manual_arg_fmt=True) def show(self, title=None, indices=None, **kwargs): """Display the parts of this product space element graphically. diff --git a/odl/test/discr/diff_ops_test.py b/odl/test/discr/diff_ops_test.py index 5b37889d499..802e41cf57f 100644 --- a/odl/test/discr/diff_ops_test.py +++ b/odl/test/discr/diff_ops_test.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -9,15 +9,15 @@ """Unit tests for `diff_ops`.""" from __future__ import division -import pytest + import numpy as np +import pytest import odl from odl.discr.diff_ops import ( - finite_diff, PartialDerivative, Gradient, Divergence, Laplacian) + Divergence, Gradient, Laplacian, PartialDerivative, finite_diff) from odl.util.testutils import ( - all_equal, all_almost_equal, dtype_tol, noise_element, simple_fixture) - + all_almost_equal, all_equal, dtype_tol, noise_element, simple_fixture) # --- pytest fixtures --- # diff --git a/odl/test/discr/discr_mappings_test.py b/odl/test/discr/discr_mappings_test.py index 3cab6e85d80..0fe99f57244 100644 --- a/odl/test/discr/discr_mappings_test.py +++ b/odl/test/discr/discr_mappings_test.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -9,14 +9,15 @@ """Unit tests for `discr_mappings`.""" from __future__ import division + import numpy as np import pytest import odl -from odl.discr.grid import sparse_meshgrid from odl.discr.discr_mappings import ( - PointCollocation, NearestInterpolation, LinearInterpolation, - PerAxisInterpolation) + LinearInterpolation, NearestInterpolation, PerAxisInterpolation, + PointCollocation) +from odl.discr.grid import sparse_meshgrid from odl.util.testutils import all_almost_equal, all_equal diff --git a/odl/test/discr/discr_ops_test.py b/odl/test/discr/discr_ops_test.py index 061cccbdfdb..1c854e56eac 100644 --- a/odl/test/discr/discr_ops_test.py +++ b/odl/test/discr/discr_ops_test.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -9,15 +9,15 @@ """Unit tests for `discr_ops`.""" from __future__ import division -import pytest + import numpy as np +import pytest import odl from odl.discr.discr_ops import _SUPPORTED_RESIZE_PAD_MODES from odl.space.entry_points import tensor_space_impl from odl.util import is_numeric_dtype, is_real_floating_dtype -from odl.util.testutils import noise_element, dtype_tol - +from odl.util.testutils import dtype_tol, noise_element # --- pytest fixtures --- # diff --git a/odl/test/discr/grid_test.py b/odl/test/discr/grid_test.py index e91827038ba..84173af2b83 100644 --- a/odl/test/discr/grid_test.py +++ b/odl/test/discr/grid_test.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -7,14 +7,14 @@ # obtain one at https://mozilla.org/MPL/2.0/. from __future__ import division -import pytest + import numpy as np +import pytest import odl -from odl.discr.grid import RectGrid, uniform_grid, sparse_meshgrid +from odl.discr.grid import RectGrid, sparse_meshgrid, uniform_grid from odl.util.testutils import all_equal - # ---- RectGrid ---- # diff --git a/odl/test/discr/lp_discr_test.py b/odl/test/discr/lp_discr_test.py index b97a63b57df..69cf1917358 100644 --- a/odl/test/discr/lp_discr_test.py +++ b/odl/test/discr/lp_discr_test.py @@ -7,6 +7,7 @@ # obtain one at https://mozilla.org/MPL/2.0/. from __future__ import division + import numpy as np from packaging.version import parse as parse_version import pytest @@ -17,8 +18,7 @@ from odl.space.npy_tensors import NumpyTensor from odl.space.weighting import ConstWeighting from odl.util.testutils import ( - all_equal, all_almost_equal, noise_elements, simple_fixture) - + all_almost_equal, all_equal, noise_elements, simple_fixture) USE_ARRAY_UFUNCS_INTERFACE = ( parse_version(np.__version__) >= parse_version('1.13')) diff --git a/odl/test/discr/partition_test.py b/odl/test/discr/partition_test.py index 333c272fb15..77022d3ac4c 100644 --- a/odl/test/discr/partition_test.py +++ b/odl/test/discr/partition_test.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -7,12 +7,12 @@ # obtain one at https://mozilla.org/MPL/2.0/. from __future__ import division -import pytest + import numpy as np +import pytest import odl -from odl.util.testutils import all_equal, all_almost_equal - +from odl.util.testutils import all_almost_equal, all_equal # ---- RectPartition ---- # diff --git a/odl/test/operator/operator_test.py b/odl/test/operator/operator_test.py index 80c89f9b5bf..914e4e92149 100644 --- a/odl/test/operator/operator_test.py +++ b/odl/test/operator/operator_test.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -7,22 +7,23 @@ # obtain one at https://mozilla.org/MPL/2.0/. from __future__ import division -import pytest -import numpy as np + import sys +import numpy as np +import pytest + import odl -from odl import (Operator, OperatorSum, OperatorComp, - OperatorLeftScalarMult, OperatorRightScalarMult, - FunctionalLeftVectorMult, OperatorRightVectorMult, - MatrixOperator, OperatorLeftVectorMult, - OpTypeError, OpDomainError, OpRangeError) -from odl.operator.operator import _function_signature, _dispatch_call_args +from odl import ( + FunctionalLeftVectorMult, MatrixOperator, OpDomainError, Operator, + OperatorComp, OperatorLeftScalarMult, OperatorLeftVectorMult, + OperatorRightScalarMult, OperatorRightVectorMult, OperatorSum, + OpRangeError, OpTypeError) +from odl.operator.operator import _dispatch_call_args, _function_signature from odl.util.testutils import ( all_almost_equal, noise_element, noise_elements, simple_fixture) from odl.util.utility import getargspec - # --- Fixtures --- # diff --git a/odl/test/operator/oputils_test.py b/odl/test/operator/oputils_test.py index 54b0d63f0d5..8c7071cfa62 100644 --- a/odl/test/operator/oputils_test.py +++ b/odl/test/operator/oputils_test.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -7,6 +7,7 @@ # obtain one at https://mozilla.org/MPL/2.0/. from __future__ import division + import numpy as np import pytest @@ -107,6 +108,7 @@ def _call(self, x): def test_matrix_representation_wrong_domain(): """Verify that the matrix representation function gives correct error""" + class MyOp(odl.Operator): """Small test operator.""" def __init__(self): @@ -125,6 +127,7 @@ def _call(self, x, out): def test_matrix_representation_wrong_range(): """Verify that the matrix representation function gives correct error""" + class MyOp(odl.Operator): """Small test operator.""" def __init__(self): diff --git a/odl/test/operator/pspace_ops_test.py b/odl/test/operator/pspace_ops_test.py index 7f7e1f572b7..34cc1aa8c98 100644 --- a/odl/test/operator/pspace_ops_test.py +++ b/odl/test/operator/pspace_ops_test.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -7,12 +7,12 @@ # obtain one at https://mozilla.org/MPL/2.0/. from __future__ import division + import pytest import odl from odl.util.testutils import all_almost_equal, simple_fixture - base_op = simple_fixture( 'base_op', [odl.IdentityOperator(odl.rn(3)), diff --git a/odl/test/operator/tensor_ops_test.py b/odl/test/operator/tensor_ops_test.py index 57f7d199478..9ef70b87c6b 100644 --- a/odl/test/operator/tensor_ops_test.py +++ b/odl/test/operator/tensor_ops_test.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2018 The ODL contributors # # This file is part of ODL. # @@ -9,18 +9,18 @@ """Unit tests for `tensor_ops`.""" from __future__ import division -import pytest + import numpy as np +import pytest import scipy import odl from odl.operator.tensor_ops import ( - PointwiseNorm, PointwiseInner, PointwiseSum, MatrixOperator) + MatrixOperator, PointwiseInner, PointwiseNorm, PointwiseSum) from odl.space.pspace import ProductSpace from odl.util import moveaxis from odl.util.testutils import ( - all_almost_equal, all_equal, simple_fixture, noise_element, noise_elements) - + all_almost_equal, all_equal, noise_element, noise_elements, simple_fixture) matrix_dtype = simple_fixture( name='matrix_dtype', diff --git a/odl/util/utility.py b/odl/util/utility.py index 828ee0ed5d1..02f96b71e64 100644 --- a/odl/util/utility.py +++ b/odl/util/utility.py @@ -1307,7 +1307,8 @@ def attribute_repr_string(inst_str, attr_str): def method_repr_string(inst_str, meth_str, arg_strs=None, - allow_mixed_seps=True): + allow_mixed_seps=True, + force_manual_arg_fmt=False): r"""Return a repr string for a method that respects line width. This function is useful to generate a ``repr`` string for a derived @@ -1336,6 +1337,10 @@ class that is created through a method, for instance :: In case some of the ``arg_strs`` span multiple lines, it is usually advisable to set ``allow_mixed_seps`` to ``False`` since the result tends to be more readable that way. + force_manual_arg_fmt : bool, optional + If ``True``, force the function to not use newlines around the + argument part and not to indent. This can be useful when the arguments + already contain appropriate newlines and are already indented. Returns ------- @@ -1411,7 +1416,9 @@ class that is created through a method, for instance :: # Method call part arg_str_oneline = ', '.join(arg_strs) - if meth_line_start_len + 1 + len(arg_str_oneline) + 1 <= linewidth: + if ('\n' not in arg_str_oneline and + meth_line_start_len + 1 + len(arg_str_oneline) + 1 <= linewidth + ): meth_call_str = '(' + arg_str_oneline + ')' elif not arg_str_oneline: meth_call_str = '(\n)' @@ -1425,7 +1432,10 @@ class that is created through a method, for instance :: for arg_str, sep in zip_longest(arg_strs, arg_seps, fillvalue=''): full_arg_str += arg_str + sep - meth_call_str = '(\n' + indent(full_arg_str) + '\n)' + if force_manual_arg_fmt: + meth_call_str = '(' + full_arg_str + ')' + else: + meth_call_str = '(\n' + indent(full_arg_str) + '\n)' return '.'.join(init_parts) + meth_call_str From 25cedfbc39f4babd731e9d115c24f06f605e2500 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Mon, 25 Jun 2018 22:43:44 +0200 Subject: [PATCH 3/9] STY: pep8 fixes --- odl/operator/default_ops.py | 2 ++ odl/operator/tensor_ops.py | 3 ++- odl/util/utility.py | 3 ++- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/odl/operator/default_ops.py b/odl/operator/default_ops.py index 2be2611b88e..f8210591f35 100644 --- a/odl/operator/default_ops.py +++ b/odl/operator/default_ops.py @@ -1872,6 +1872,7 @@ def __repr__(self): inner_parts = signature_string_parts(posargs, optargs) return repr_string(self.__class__.__name__, inner_parts) + class ComplexModulusSquared(Operator): """Operator that computes the squared complex modulus (absolute value).""" @@ -2090,6 +2091,7 @@ def __repr__(self): inner_parts = signature_string_parts(posargs, optargs) return repr_string(self.__class__.__name__, inner_parts) + if __name__ == '__main__': from odl.util.testutils import run_doctests run_doctests() diff --git a/odl/operator/tensor_ops.py b/odl/operator/tensor_ops.py index 5879b0187e4..06b79e33d7c 100644 --- a/odl/operator/tensor_ops.py +++ b/odl/operator/tensor_ops.py @@ -570,7 +570,8 @@ class PointwiseInnerAdjoint(PointwiseTensorFieldOperator): def _call(self, f, out): """Implement ``self(vf, out)``.""" for vfi, oi, ran_wi, dom_wi in zip( - op.vecfield, out, adj_ran_weights, op.weights): + op.vecfield, out, adj_ran_weights, op.weights + ): vfi.multiply(f, out=oi) if not np.isclose(ran_wi, dom_wi): oi *= dom_wi / ran_wi diff --git a/odl/util/utility.py b/odl/util/utility.py index 02f96b71e64..52f455d0a02 100644 --- a/odl/util/utility.py +++ b/odl/util/utility.py @@ -1416,7 +1416,8 @@ class that is created through a method, for instance :: # Method call part arg_str_oneline = ', '.join(arg_strs) - if ('\n' not in arg_str_oneline and + if ( + '\n' not in arg_str_oneline and meth_line_start_len + 1 + len(arg_str_oneline) + 1 <= linewidth ): meth_call_str = '(' + arg_str_oneline + ')' From 9efe5ac3bcfbcd14ef7a2f0351e74e2d4b8a016e Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Mon, 25 Jun 2018 23:25:32 +0200 Subject: [PATCH 4/9] DOC/TST: fix various doctest failures --- odl/solvers/util/callback.py | 2 +- odl/tomo/geometry/conebeam.py | 10 ++-------- odl/tomo/geometry/parallel.py | 10 ++-------- 3 files changed, 5 insertions(+), 17 deletions(-) diff --git a/odl/solvers/util/callback.py b/odl/solvers/util/callback.py index 31f333b296d..e18cccd40c2 100644 --- a/odl/solvers/util/callback.py +++ b/odl/solvers/util/callback.py @@ -173,7 +173,7 @@ def __repr__(self): >>> callback = odl.solvers.CallbackPrint() >>> operator = odl.ScalingOperator(r3, 2.0) >>> callback * operator - CallbackPrint() * ScalingOperator(rn(3), 2.0) + CallbackPrint() * ScalingOperator(rn(3), scalar=2.0) """ return '{!r} * {!r}'.format(self.callback, self.operator) diff --git a/odl/tomo/geometry/conebeam.py b/odl/tomo/geometry/conebeam.py index 31578f60c08..93092b870d1 100644 --- a/odl/tomo/geometry/conebeam.py +++ b/odl/tomo/geometry/conebeam.py @@ -533,10 +533,7 @@ def __getitem__(self, indices): >>> geom[::2, :] FanFlatGeometry( - nonuniform_partition( - [ 0.5, 2.5], - min_pt=0.0, max_pt=4.0 - ), + nonuniform_partition([ 0.5, 2.5], min_pt=0.0, max_pt=4.0), uniform_partition(-1.0, 1.0, 20), src_radius=50.0, det_radius=100.0 @@ -1183,10 +1180,7 @@ def __getitem__(self, indices): >>> geom[::2] ConeFlatGeometry( - nonuniform_partition( - [ 0.5, 2.5], - min_pt=0.0, max_pt=4.0 - ), + nonuniform_partition([ 0.5, 2.5], min_pt=0.0, max_pt=4.0), uniform_partition([-1., -1.], [ 1., 1.], (20, 20)), src_radius=50.0, det_radius=100.0, diff --git a/odl/tomo/geometry/parallel.py b/odl/tomo/geometry/parallel.py index 84e45a732ad..a078e39c803 100644 --- a/odl/tomo/geometry/parallel.py +++ b/odl/tomo/geometry/parallel.py @@ -682,10 +682,7 @@ def __getitem__(self, indices): >>> geom[::2] Parallel2dGeometry( - nonuniform_partition( - [ 0.5, 2.5], - min_pt=0.0, max_pt=4.0 - ), + nonuniform_partition([ 0.5, 2.5], min_pt=0.0, max_pt=4.0), uniform_partition(-1.0, 1.0, 20) ) """ @@ -1446,10 +1443,7 @@ def __getitem__(self, indices): >>> geom[::2] Parallel3dAxisGeometry( - nonuniform_partition( - [ 0.5, 2.5], - min_pt=0.0, max_pt=4.0 - ), + nonuniform_partition([ 0.5, 2.5], min_pt=0.0, max_pt=4.0), uniform_partition([-1., -1.], [ 1., 1.], (20, 20)) ) """ From c80607c4e61112d291a0e6fc84beb4ca5d368598 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Mon, 27 Aug 2018 00:05:16 +0200 Subject: [PATCH 5/9] ENH: make repr precision more flexible --- odl/discr/diff_ops.py | 10 +++++----- odl/discr/discr_ops.py | 6 +++--- odl/discr/grid.py | 6 +++--- odl/discr/lp_discr.py | 12 ++++++------ odl/discr/partition.py | 15 +++++++-------- odl/operator/default_ops.py | 18 +++++++++--------- odl/operator/operator.py | 14 +++++++------- odl/operator/tensor_ops.py | 6 +++--- odl/util/utility.py | 24 ++++++++++++++++++++++-- 9 files changed, 65 insertions(+), 46 deletions(-) diff --git a/odl/discr/diff_ops.py b/odl/discr/diff_ops.py index 17aab220562..6eecb9bccda 100644 --- a/odl/discr/diff_ops.py +++ b/odl/discr/diff_ops.py @@ -16,7 +16,7 @@ from odl.operator.tensor_ops import PointwiseTensorFieldOperator from odl.space import ProductSpace from odl.util import ( - REPR_PRECISION, npy_printoptions, repr_string, signature_string_parts, + npy_printoptions, repr_precision, repr_string, signature_string_parts, writable_array) __all__ = ('PartialDerivative', 'Gradient', 'Divergence', 'Laplacian') @@ -193,7 +193,7 @@ def __repr__(self): ('pad_mode', self.pad_mode, 'constant'), ('pad_const', self.pad_const, 0)] optmod = ['', '!r', '', '', ''] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs, mod=['!r', optmod]) return repr_string(self.__class__.__name__, inner_parts, @@ -414,7 +414,7 @@ def __repr__(self): ('pad_mode', self.pad_mode, 'constant'), ('pad_const', self.pad_const, 0)] optmod = ['!r', '', '', ''] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs, mod=['!r', optmod]) return repr_string(self.__class__.__name__, inner_parts, @@ -624,7 +624,7 @@ def __repr__(self): ('pad_mode', self.pad_mode, 'constant'), ('pad_const', self.pad_const, 0)] optmod = ['!r', '', '', ''] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs, mod=['!r', optmod]) return repr_string(self.__class__.__name__, inner_parts, @@ -792,7 +792,7 @@ def __repr__(self): ('pad_mode', self.pad_mode, 'constant'), ('pad_const', self.pad_const, 0)] optmod = ['!r', '', ''] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs, mod=['!r', optmod]) return repr_string(self.__class__.__name__, inner_parts, diff --git a/odl/discr/discr_ops.py b/odl/discr/discr_ops.py index 6ec5a280e78..a92303664ad 100644 --- a/odl/discr/discr_ops.py +++ b/odl/discr/discr_ops.py @@ -17,8 +17,8 @@ from odl.set import IntervalProd from odl.space import FunctionSpace, tensor_space from odl.util import ( - REPR_PRECISION, attribute_repr_string, normalized_scalar_param_list, - npy_printoptions, repr_string, resize_array, safe_int_conv, + attribute_repr_string, normalized_scalar_param_list, npy_printoptions, + repr_precision, repr_string, resize_array, safe_int_conv, signature_string_parts, writable_array) from odl.util.numerics import _SUPPORTED_RESIZE_PAD_MODES @@ -470,7 +470,7 @@ def __repr__(self): posargs = [self.domain, self.range] optargs = [('pad_mode', self.pad_mode, 'constant'), ('pad_const', self.pad_const, 0)] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs, mod=['!r', '']) return repr_string(self.__class__.__name__, inner_parts, diff --git a/odl/discr/grid.py b/odl/discr/grid.py index 0795e54eb10..e9f3a77b40c 100644 --- a/odl/discr/grid.py +++ b/odl/discr/grid.py @@ -18,8 +18,8 @@ from odl.set import IntervalProd, Set from odl.util import ( - REPR_PRECISION, array_str, normalized_index_expression, - normalized_scalar_param_list, npy_printoptions, repr_string, safe_int_conv, + array_str, normalized_index_expression, normalized_scalar_param_list, + npy_printoptions, repr_precision, repr_string, safe_int_conv, signature_string_parts) __all__ = ('RectGrid', 'uniform_grid', 'uniform_grid_fromintv') @@ -997,7 +997,7 @@ def __repr__(self): posargs = self.coord_vectors posmod = array_str - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, [], mod=[posmod, '']) return repr_string(ctor, inner_parts) diff --git a/odl/discr/lp_discr.py b/odl/discr/lp_discr.py index 4b13c061763..b40933365fd 100644 --- a/odl/discr/lp_discr.py +++ b/odl/discr/lp_discr.py @@ -26,11 +26,11 @@ from odl.space.entry_points import tensor_space_impl from odl.space.weighting import ConstWeighting from odl.util import ( - REPR_PRECISION, apply_on_boundary, array_str, attribute_repr_string, - dtype_str, is_complex_floating_dtype, is_floating_dtype, is_numeric_dtype, + 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) + normalized_scalar_param_list, npy_printoptions, repr_precision, + repr_string, safe_int_conv, signature_string_parts) __all__ = ('DiscreteLp', 'DiscreteLpElement', 'uniform_discr_frompartition', 'uniform_discr_fromspace', @@ -602,7 +602,7 @@ def __repr__(self): if self.dtype in (float, complex, int, bool): optmod[3] = '!s' - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs, mod=[posmod, optmod]) @@ -611,7 +611,7 @@ def __repr__(self): posargs = [self.fspace, self.partition, self.tspace] optargs = [('interp', self.interp, 'nearest')] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs) return repr_string(ctor, inner_parts, allow_mixed_seps=True) diff --git a/odl/discr/partition.py b/odl/discr/partition.py index 831cc5b7c8e..89ed083e8f4 100644 --- a/odl/discr/partition.py +++ b/odl/discr/partition.py @@ -23,10 +23,9 @@ 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, npy_printoptions, repr_string, safe_int_conv, - signature_string_parts) + array_str, attribute_repr_string, normalized_index_expression, + normalized_nodes_on_bdry, normalized_scalar_param_list, npy_printoptions, + repr_precision, repr_string, safe_int_conv, signature_string_parts) __all__ = ('RectPartition', 'uniform_partition_fromintv', 'uniform_partition_fromgrid', 'uniform_partition', @@ -898,7 +897,7 @@ def __repr__(self): optargs = [('nodes_on_bdry', self.nodes_on_bdry, False)] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs, mod=[posmod, '']) else: @@ -930,7 +929,7 @@ def __repr__(self): optargs.append(('min_pt', self.min_pt[0], None)) optmod.append('') else: - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): optargs.append( ('min_pt', array_str(self.min_pt), '')) optmod.append('!s') @@ -940,12 +939,12 @@ def __repr__(self): optargs.append(('max_pt', self.max_pt[0], None)) optmod.append('') else: - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): optargs.append( ('max_pt', array_str(self.max_pt), '')) optmod.append('!s') - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs, mod=[posmod, optmod]) diff --git a/odl/operator/default_ops.py b/odl/operator/default_ops.py index f8210591f35..8f9be6a9d54 100644 --- a/odl/operator/default_ops.py +++ b/odl/operator/default_ops.py @@ -20,8 +20,8 @@ from odl.set import ComplexNumbers, Field, LinearSpace, RealNumbers from odl.space import ProductSpace from odl.util import ( - REPR_PRECISION, attribute_repr_string, method_repr_string, - npy_printoptions, repr_string, signature_string_parts) + attribute_repr_string, method_repr_string, npy_printoptions, + repr_precision, repr_string, signature_string_parts) __all__ = ('ScalingOperator', 'ZeroOperator', 'IdentityOperator', 'LinCombOperator', 'MultiplyOperator', 'PowerOperator', @@ -224,7 +224,7 @@ def __repr__(self): posargs = [self.domain] optargs = [('scalar', self.scalar, None), ('range', self.range, self.domain)] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs) return repr_string(self.__class__.__name__, inner_parts) @@ -481,7 +481,7 @@ def __repr__(self): optargs = [('domain', self.domain, getattr(self.multiplicand, 'space', None)), ('range', self.range, self.domain)] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs) return repr_string(self.__class__.__name__, inner_parts) @@ -601,7 +601,7 @@ def __repr__(self): posargs = [self.domain] optargs = [('exponent', self.exponent, None), ('range', self.range, self.domain)] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs) return repr_string(self.__class__.__name__, inner_parts) @@ -742,7 +742,7 @@ def __repr__(self): optargs = [('domain', self.domain, getattr(self.vector, 'space', None)), ('range', self.range, self.domain.field)] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs) return repr_string(self.__class__.__name__, inner_parts) @@ -975,7 +975,7 @@ def __repr__(self): optargs = [('domain', self.domain, getattr(self.vector, 'space', None)), ('range', self.range, RealNumbers())] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs) return repr_string(self.__class__.__name__, inner_parts) @@ -1067,7 +1067,7 @@ def __repr__(self): optargs = [('domain', self.domain, self.range), ('range', self.range, getattr(self.constant, 'space', None))] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs) return repr_string(self.__class__.__name__, inner_parts) @@ -1636,7 +1636,7 @@ def __repr__(self): posargs = [self.domain] optargs = [('range', self.range, self.domain.complex_space), ('scalar', self.scalar, 1.0)] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs) return repr_string(self.__class__.__name__, inner_parts) diff --git a/odl/operator/operator.py b/odl/operator/operator.py index 3f8c4f0c567..3de08d136ac 100644 --- a/odl/operator/operator.py +++ b/odl/operator/operator.py @@ -20,7 +20,7 @@ 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, + cache_arguments, npy_printoptions, repr_precision, repr_string, signature_string_parts) __all__ = ('Operator', 'OperatorComp', 'OperatorSum', 'OperatorVectorSum', @@ -1311,7 +1311,7 @@ def derivative(self, point): def __repr__(self): """Return ``repr(self)``.""" posargs = [self.operator, self.vector] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, []) return repr_string(self.__class__.__name__, inner_parts, allow_mixed_seps=False) @@ -1687,7 +1687,7 @@ def adjoint(self): def __repr__(self): """Return ``repr(self)``.""" posargs = [self.operator, self.scalar] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, []) return repr_string(self.__class__.__name__, inner_parts, allow_mixed_seps=False) @@ -1871,7 +1871,7 @@ def adjoint(self): def __repr__(self): """Return ``repr(self)``.""" posargs = [self.operator, self.scalar] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, []) return repr_string(self.__class__.__name__, inner_parts, allow_mixed_seps=False) @@ -1990,7 +1990,7 @@ def adjoint(self): def __repr__(self): """Return ``repr(self)``.""" posargs = [self.functional, self.vector] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, []) return repr_string(self.__class__.__name__, inner_parts, allow_mixed_seps=False) @@ -2106,7 +2106,7 @@ def adjoint(self): def __repr__(self): """Return ``repr(self)``.""" posargs = [self.operator, self.vector] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, []) return repr_string(self.__class__.__name__, inner_parts, allow_mixed_seps=False) @@ -2229,7 +2229,7 @@ def adjoint(self): def __repr__(self): """Return ``repr(self)``.""" posargs = [self.operator, self.vector] - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, []) return repr_string(self.__class__.__name__, inner_parts, allow_mixed_seps=False) diff --git a/odl/operator/tensor_ops.py b/odl/operator/tensor_ops.py index 06b79e33d7c..8368a7288b5 100644 --- a/odl/operator/tensor_ops.py +++ b/odl/operator/tensor_ops.py @@ -21,8 +21,8 @@ from odl.space.base_tensors import TensorSpace from odl.space.weighting import ArrayWeighting from odl.util import ( - REPR_PRECISION, array_str, attribute_repr_string, dtype_repr, moveaxis, - npy_printoptions, repr_string, signature_string_parts, writable_array) + array_str, attribute_repr_string, dtype_repr, moveaxis, npy_printoptions, + repr_precision, repr_string, signature_string_parts, writable_array) __all__ = ('PointwiseNorm', 'PointwiseInner', 'PointwiseSum', 'MatrixOperator', 'SamplingOperator', 'WeightedSumSamplingOperator', @@ -386,7 +386,7 @@ def __repr__(self): if self.is_weighted: optargs.append(('weighting', array_str(self.weights), '')) optmod.append('!s') - with npy_printoptions(precision=REPR_PRECISION): + with npy_printoptions(precision=repr_precision()): inner_parts = signature_string_parts(posargs, optargs, mod=['!r', optmod]) return repr_string(self.__class__.__name__, inner_parts, diff --git a/odl/util/utility.py b/odl/util/utility.py index 52f455d0a02..5879e31b376 100644 --- a/odl/util/utility.py +++ b/odl/util/utility.py @@ -9,13 +9,13 @@ """Utilities mainly for internal use.""" from __future__ import absolute_import, division, print_function +from future.moves.itertools import zip_longest import inspect import sys from builtins import object from collections import OrderedDict from functools import wraps -from future.moves.itertools import zip_longest from itertools import product import numpy as np @@ -29,10 +29,11 @@ 'real_dtype', 'complex_dtype', 'is_string', 'nd_iterator', 'conj_exponent', 'writable_array', 'run_from_ipython', 'NumpyRandomSeed', 'cache_arguments', 'unique', - 'REPR_PRECISION') + 'repr_precision') REPR_PRECISION = 4 # For printing scalars and array entries +NPY_DEFAULT_PRECISION = np.get_printoptions()['precision'] TYPE_MAP_R2C = {np.dtype(dtype): np.result_type(dtype, 1j) for dtype in np.sctypes['float']} @@ -1441,6 +1442,25 @@ class that is created through a method, for instance :: return '.'.join(init_parts) + meth_call_str +def repr_precision(): + """Return the print precision for floating point numbers. + + By default, the value of the global variable ``REPR_PRECISION`` + is returned, unless the user changed the value of the ``'precision'`` + entry in `numpy.get_printoptions`, e.g., via :: + + with npy_printoptions(precision=10): + ... + + """ + cur_npy_prec = np.get_printoptions()['precision'] + if cur_npy_prec == NPY_DEFAULT_PRECISION: + # No change by the user + return REPR_PRECISION + else: + return cur_npy_prec + + def run_from_ipython(): """If the process is run from IPython.""" return '__IPYTHON__' in globals() From 958e0a0593898fd58f51580da8817c6ac6cf34a1 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Wed, 29 Aug 2018 22:32:58 +0200 Subject: [PATCH 6/9] MAINT: remove operator scaling and lico --- odl/operator/default_ops.py | 104 +++++++++++++++--------------------- 1 file changed, 42 insertions(+), 62 deletions(-) diff --git a/odl/operator/default_ops.py b/odl/operator/default_ops.py index 8f9be6a9d54..fe2d2faa275 100644 --- a/odl/operator/default_ops.py +++ b/odl/operator/default_ops.py @@ -30,45 +30,6 @@ 'ComplexModulus', 'ComplexModulusSquared') -def _scale_op(operator, scalar): - """Scale an operator, optimizing for ``scalar=0`` and ``scalar=1``.""" - if scalar == 0: - return ZeroOperator(operator.domain, operator.range) - elif scalar == 1: - return operator - else: - return scalar * operator - - -def _lico_ops(a, op1, b, op2): - """Linear combination of operators, optimizing trivial cases.""" - if op1.domain != op2.domain or op1.range != op2.range: - raise ValueError('domain/range mismatch between {!r} and {!r}' - .format(op1, op2)) - dom, ran = op1.domain, op1.range - if a == 0: - if b == 0: - return ZeroOperator(dom, ran) - elif b == 1: - return op2 - else: - return b * op2 - elif a == 1: - if b == 0: - return op1 - elif b == 1: - return op1 + op2 - else: - return op1 + b * op2 - else: - if b == 0: - return a * op1 - elif b == 1: - return a * op1 + op2 - else: - return a * op1 + b * op2 - - class ScalingOperator(Operator): """Operator of multiplication with a scalar. @@ -1538,24 +1499,33 @@ def inverse(self): return ZeroOperator(self.range, self.domain) if self.domain.is_real: - # Real domain - # Optimizations for simple cases. + # Real domain, with optimizations for simple cases if self.scalar.real == self.scalar: - return _scale_op(RealPart(self.range, self.domain), - 1 / self.scalar.real) + # embedding x -> (a + i*0) * x + op = RealPart(self.range, self.domain) + if self.scalar.real != 1: + op = (1 / self.scalar.real) * op + return op elif 1j * self.scalar.imag == self.scalar: - return _scale_op(ImagPart(self.range, self.domain), - 1 / self.scalar.imag) + # embedding x -> (0 + i*b) * x + op = ImagPart(self.range, self.domain) + if self.scalar.imag != 1: + op = (1 / self.scalar.imag) * op + return op else: - # General case - inv_scalar = (1 / self.scalar).conjugate() - return _lico_ops( - inv_scalar.real, RealPart(self.range, self.domain), - inv_scalar.imag, ImagPart(self.range, self.domain)) + # embedding x -> (a + i*b) * x + inv_scalar = 1 / self.scalar + re_op = RealPart(self.range, self.domain) + if inv_scalar.real != 1: + re_op = inv_scalar.real * re_op + im_op = ImagPart(self.range, self.domain) + if inv_scalar.imag != 1: + im_op = inv_scalar.imag * im_op + + return re_op + im_op else: # Complex domain - return ComplexEmbedding(self.range, self.domain, - self.scalar.conjugate()) + return ComplexEmbedding(self.range, self.domain, 1 / self.scalar) @property def adjoint(self): @@ -1602,19 +1572,29 @@ def adjoint(self): return ZeroOperator(self.range, self.domain) if self.domain.is_real: - # Real domain - # Optimizations for simple cases. + # Real domain, with optimizations for simple cases if self.scalar.real == self.scalar: - return _scale_op(self.scalar.real, - ComplexEmbedding(self.range, self.domain)) + # embedding x -> (a + i*0) * x + op = RealPart(self.range, self.domain) + if self.scalar.real != 1: + op = self.scalar.real * op + return op elif 1j * self.scalar.imag == self.scalar: - return _scale_op(self.scalar.imag, - ImagPart(self.range, self.domain)) + # embedding x -> (0 + i*b) * x + op = ImagPart(self.range, self.domain) + if self.scalar.imag != 1: + op = self.scalar.imag * op + return op else: - # General case - return _lico_ops( - self.scalar.real, RealPart(self.range, self.domain), - self.scalar.imag, ImagPart(self.range, self.domain)) + # embedding x -> (a + i*b) * x + re_op = RealPart(self.range, self.domain) + if self.scalar.real != 1: + re_op = self.scalar.real * re_op + im_op = ImagPart(self.range, self.domain) + if self.scalar.imag != 1: + im_op = self.scalar.imag * im_op + + return re_op + im_op else: # Complex domain return ComplexEmbedding(self.range, self.domain, From f30969464273eaf8b6487dd136eb508d8c8df304 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Wed, 29 Aug 2018 23:36:39 +0200 Subject: [PATCH 7/9] MAINT: improve some docs and fix minor issues --- odl/discr/diff_ops.py | 1 - odl/discr/discr_mappings.py | 35 +++++++++++++++-------------------- odl/discr/discr_ops.py | 4 ++-- odl/discr/grid.py | 1 - odl/discr/lp_discr.py | 6 ++++-- odl/discr/partition.py | 5 +++-- odl/operator/default_ops.py | 9 +++------ odl/operator/pspace_ops.py | 16 ++++++++++------ 8 files changed, 37 insertions(+), 40 deletions(-) diff --git a/odl/discr/diff_ops.py b/odl/discr/diff_ops.py index 6eecb9bccda..183d3a8a4e4 100644 --- a/odl/discr/diff_ops.py +++ b/odl/discr/diff_ops.py @@ -498,7 +498,6 @@ def __init__(self, domain=None, range=None, method='forward', Verify adjoint: >>> g = div.range.element(data ** 2) - >>> adj_div_g = div.adjoint(g) >>> div(f).inner(g) / f.inner(div.adjoint(g)) 1.0 """ diff --git a/odl/discr/discr_mappings.py b/odl/discr/discr_mappings.py index 69a55541288..5725e1d3c0c 100644 --- a/odl/discr/discr_mappings.py +++ b/odl/discr/discr_mappings.py @@ -417,16 +417,14 @@ def __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) + >>> space = odl.uniform_discr([0, 0], [1, 1], shape=(4, 2)) + >>> interp_op = odl.NearestInterpolation( + ... space.fspace, space.partition, space.tspace) >>> interp_op NearestInterpolation( FunctionSpace(IntervalProd([ 0., 0.], [ 1., 1.])), uniform_partition([ 0., 0.], [ 1., 1.], (4, 2)), - rn((4, 2)) + rn((4, 2), weighting=0.125) ) """ posargs = [self.range, self.partition, self.domain] @@ -485,16 +483,14 @@ def __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) + >>> space = odl.uniform_discr([0, 0], [1, 1], shape=(4, 2)) + >>> interp_op = odl.LinearInterpolation( + ... space.fspace, space.partition, space.tspace) >>> interp_op LinearInterpolation( FunctionSpace(IntervalProd([ 0., 0.], [ 1., 1.])), uniform_partition([ 0., 0.], [ 1., 1.], (4, 2)), - rn((4, 2)) + rn((4, 2), weighting=0.125) ) """ posargs = [self.range, self.partition, self.domain] @@ -636,17 +632,15 @@ def __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']) + >>> space = odl.uniform_discr([0, 0], [1, 1], shape=(4, 2)) + >>> interp_op = odl.PerAxisInterpolation( + ... space.fspace, space.partition, space.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)), + rn((4, 2), weighting=0.125), schemes=('linear', 'nearest') ) """ @@ -965,7 +959,8 @@ def _evaluate(self, indices, norm_distances, out=None): """Evaluate linear interpolation. Modified for in-place evaluation and treatment of out-of-bounds - points by implicitly assuming 0 at the next node.""" + points by implicitly assuming 0 at the next node. + """ # slice for broadcasting over trailing dimensions in self.values vslice = (slice(None),) + (None,) * (self.values.ndim - len(indices)) diff --git a/odl/discr/discr_ops.py b/odl/discr/discr_ops.py index a92303664ad..e3b04a071d4 100644 --- a/odl/discr/discr_ops.py +++ b/odl/discr/discr_ops.py @@ -111,7 +111,7 @@ def inverse(self): >>> resampling = odl.Resampling(coarse_discr, fine_discr) >>> resampling_inv = resampling.inverse - The inverse is proper left inverse if the resampling goes from a + The inverse is a proper left inverse if the resampling maps from a coarser to a finer sampling: >>> resampling_inv(resampling([0.0, 1.0, 0.0])) @@ -171,7 +171,7 @@ class ResizingOperator(Operator): ("padded") according to a provided parameter ``pad_mode``. All resizing operator variants are linear, except constant padding - with constant != 0. + with constant != 0, which is affine. See `the online documentation `_ diff --git a/odl/discr/grid.py b/odl/discr/grid.py index e9f3a77b40c..f88caaf0c37 100644 --- a/odl/discr/grid.py +++ b/odl/discr/grid.py @@ -508,7 +508,6 @@ def __eq__(self, other): def __hash__(self): """Return ``hash(self)``.""" - # TODO(doable, #841): update (what exactly) coord_vec_str = tuple(cv.tobytes() for cv in self.coord_vectors) return hash((type(self), coord_vec_str)) diff --git a/odl/discr/lp_discr.py b/odl/discr/lp_discr.py index b40933365fd..23cab9c9ed1 100644 --- a/odl/discr/lp_discr.py +++ b/odl/discr/lp_discr.py @@ -521,7 +521,8 @@ def __repr__(self): Examples -------- - Uniform: + Uniformly discretized space, resulting in ``repr`` using the + `uniform_discr` factory function: >>> space = odl.uniform_discr(0, 1, 4, dtype='float32') >>> space @@ -533,7 +534,8 @@ def __repr__(self): [ 0., 0., 0.], [ 1., 1., 1.], (2, 4, 8), dtype=complex ) - Non-uniform: + For non-uniform discretizations, the output is more verbose and + uses the class constructor: >>> rect = odl.IntervalProd([0, 0], [1, 1]) >>> fspace = odl.FunctionSpace(rect) diff --git a/odl/discr/partition.py b/odl/discr/partition.py index 89ed083e8f4..e20615eb96c 100644 --- a/odl/discr/partition.py +++ b/odl/discr/partition.py @@ -497,13 +497,14 @@ def __getitem__(self, indices): Examples -------- - Take every second grid point. Note that is is in general non-uniform: + Take every second, starting at index 2 (note that this generally + results in a non-uniform partition): >>> partition = odl.uniform_partition(0, 10, 10) >>> 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: + A more advanced example: >>> intvp = odl.IntervalProd([-1, 1, 4, 2], [3, 6, 5, 7]) >>> grid = odl.RectGrid([-1, 0, 3], [2, 4], [5], [2, 4, 7]) diff --git a/odl/operator/default_ops.py b/odl/operator/default_ops.py index fe2d2faa275..4b13048a608 100644 --- a/odl/operator/default_ops.py +++ b/odl/operator/default_ops.py @@ -333,9 +333,9 @@ def __init__(self, multiplicand, domain=None, range=None): rn(3).element([ 1., 4., 9.]) For a scalar or `element-like` multiplicand, ``domain`` (and - ``range``) should be given: + optionally ``range``) should be given: - >>> op = odl.MultiplyOperator(x, domain=r3.field, range=r3) + >>> op = odl.MultiplyOperator(x, domain=r3.field) >>> op(3) rn(3).element([ 3., 6., 9.]) >>> out = r3.element() @@ -1668,10 +1668,7 @@ def __init__(self, domain, range=None): def _call(self, x): """Return ``self(x)``.""" squared_mod = x.real ** 2 + x.imag ** 2 - if hasattr(squared_mod, 'ufuncs'): - return squared_mod.ufuncs.sqrt() - else: - return np.sqrt(squared_mod) + return squared_mod.ufuncs.sqrt() def derivative(self, x): r"""Return the derivative operator in the "C = R^2" sense. diff --git a/odl/operator/pspace_ops.py b/odl/operator/pspace_ops.py index 6a3eff749aa..638bdecf305 100644 --- a/odl/operator/pspace_ops.py +++ b/odl/operator/pspace_ops.py @@ -519,8 +519,8 @@ def __repr__(self): Examples -------- - All rows and columns with an operator, i.e., domain and range - are unambiguous: + All rows and columns contain at least one `Operator` instance, + so both domain and range can be inferred: >>> r3 = odl.rn(3) >>> pspace = odl.ProductSpace(r3, r3) @@ -533,10 +533,12 @@ def __repr__(self): [0, IdentityOperator(rn(3))]] ) - Domain not completely determined (column with all zeros): + Second column all zeros, this the domain is not determined and + has to be specified: >>> prod_op = odl.ProductSpaceOperator([[I, 0], - ... [I, 0]], domain=pspace) + ... [I, 0]], + ... domain=pspace) >>> prod_op ProductSpaceOperator( [[IdentityOperator(rn(3)), 0], @@ -544,10 +546,12 @@ def __repr__(self): domain=ProductSpace(rn(3), 2) ) - Range not completely determined (row with all zeros): + Second row all zeros, thus the range is undetermined and has to + be given: >>> prod_op = odl.ProductSpaceOperator([[I, I], - ... [0, 0]], range=pspace) + ... [0, 0]], + ... range=pspace) >>> prod_op ProductSpaceOperator( [[IdentityOperator(rn(3)), IdentityOperator(rn(3))], From 7aa2b31357ebef7a6c18776744202f497bf7c8cf Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Wed, 29 Aug 2018 23:37:16 +0200 Subject: [PATCH 8/9] MAINT: make standalone ComponentEmbedding class from ComponentProjection.adjoint --- odl/operator/pspace_ops.py | 128 +++++++++++++++++++++++++++---------- 1 file changed, 96 insertions(+), 32 deletions(-) diff --git a/odl/operator/pspace_ops.py b/odl/operator/pspace_ops.py index 638bdecf305..e4180574407 100644 --- a/odl/operator/pspace_ops.py +++ b/odl/operator/pspace_ops.py @@ -17,10 +17,9 @@ 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) +from odl.util import array_str, repr_string, signature_string_parts -__all__ = ('ProductSpaceOperator', 'ComponentProjection', +__all__ = ('ProductSpaceOperator', 'ComponentProjection', 'ComponentEmbedding', 'BroadcastOperator', 'ReductionOperator', 'DiagonalOperator') @@ -694,50 +693,115 @@ def adjoint(self): [ 4., 5., 6.] ]) """ - op = self + return ComponentEmbedding(self.domain, self.index) - class ComponentProjectionAdjoint(Operator): + 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) + """ + posargs = [self.domain, self.index] + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts) - """Adjoint operator to `ComponentProjection`.""" - def _call(self, x, out=None): - """Extend ``x`` to the full space.""" - if out is None: - out = self.range.zero() - else: - out.set_zero() +class ComponentEmbedding(Operator): + + r"""Embedding of a space into a product space at a given index. + + For a product space :math:`X = X_1 \times \dots \times X_n`, the + component embedding is defined as + + .. math:: + E_i: X_i \to X,\quad E_i(x_i) = (0, \dots, 0, x_i, 0, \dots, 0). + + More generally, for an index set :math:`I \subset \{1, \dots, n\}`, + the embedding operator :math:`E_I` inserts the components + of a given vector into the larger vector at the indices :math:`i \in I`, + and sets the other components to zero. + + This is the (unweighted) adjoint operator to `ComponentProjection`. + """ + + def __init__(self, space, index): + """Initialize a new instance. + + Parameters + ---------- + space : `ProductSpace` + Space to embed into. + index : int, slice, or list + Indices defining the subspace. If ``index`` is not an integer, + the `Operator.domain` of this operator is also a `ProductSpace`. + + Examples + -------- + Embedding with a single integer component: + + >>> pspace = odl.ProductSpace(odl.rn(1), odl.rn(2), odl.rn(3)) + >>> emb = odl.ComponentEmbedding(pspace, 0) + >>> x0 = emb.domain.element([1]) + >>> emb(x0) + ProductSpace(rn(1), rn(2), rn(3)).element([ + [ 1.], + [ 0., 0.], + [ 0., 0., 0.] + ]) - out[op.index] = x - return out + Embedding with multiple indices: + + >>> emb = odl.ComponentEmbedding(pspace, [0, 2]) + >>> x02 = emb.domain.element([[1], [2, 3, 4]]) + >>> emb(x02) + ProductSpace(rn(1), rn(2), rn(3)).element([ + [ 1.], + [ 0., 0.], + [ 2., 3., 4.] + ]) + """ + if np.isscalar(index): + self.__index = int(index) + elif isinstance(index, slice): + self.__index = index + else: + self.__index = list(index) + super(ComponentEmbedding, self).__init__( + space[index], space, linear=True) - @property - def adjoint(self): - """Adjoint of the adjoint, the original operator.""" - return op + @property + def index(self): + """Index, indices or slice defining the subspace.""" + return self.__index - def __repr__(self): - """Return ``repr(self)``. + def _call(self, x, out=None): + """Extend ``x`` to the full space.""" + if out is None: + out = self.range.zero() + else: + out.set_zero() - 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') + out[self.index] = x + return out - return ComponentProjectionAdjoint(self.range, self.domain, linear=True) + @property + def adjoint(self): + """Adjoint operator, the projection onto the index set.""" + return ComponentProjection(self.range, 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) + >>> pspace = odl.rn(1) * odl.rn(2) + >>> odl.ComponentEmbedding(pspace, 0) + ComponentEmbedding(ProductSpace(rn(1), rn(2)), 0) """ - posargs = [self.domain, self.index] + posargs = [self.range, self.index] inner_parts = signature_string_parts(posargs, []) return repr_string(self.__class__.__name__, inner_parts) From e9c807c2a3a5e5af7ab30dd5fad179618300bd2d Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Thu, 30 Aug 2018 23:50:53 +0200 Subject: [PATCH 9/9] BUG: fix wrong range of MultiplyOperator --- odl/operator/default_ops.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/odl/operator/default_ops.py b/odl/operator/default_ops.py index 4b13048a608..93b1790f099 100644 --- a/odl/operator/default_ops.py +++ b/odl/operator/default_ops.py @@ -346,7 +346,7 @@ def __init__(self, multiplicand, domain=None, range=None): domain = multiplicand.space if range is None: - range = domain + range = multiplicand.space super(MultiplyOperator, self).__init__(domain, range, linear=True)