diff --git a/odl/__init__.py b/odl/__init__.py index 835f4027e86..eea1aa6298e 100644 --- a/odl/__init__.py +++ b/odl/__init__.py @@ -65,8 +65,10 @@ from . import solvers from . import tomo from . import trafos -from . import ufunc_ops from . import util +if False: + # Broken due to lack of `ufuncs` + from . import ufunc_ops # Add `test` function to global namespace so users can run `odl.test()` from .util import test diff --git a/odl/space/base_tensors.py b/odl/space/base_tensors.py index cc94efd347e..c50de7d23f4 100644 --- a/odl/space/base_tensors.py +++ b/odl/space/base_tensors.py @@ -1,4 +1,4 @@ -# Copyright 2014-2018 The ODL contributors +# Copyright 2014-2019 The ODL contributors # # This file is part of ODL. # @@ -9,7 +9,6 @@ """Base classes for implementations of tensor spaces.""" from __future__ import print_function, division, absolute_import -from builtins import object from numbers import Integral import numpy as np @@ -18,8 +17,7 @@ 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) -from odl.util.ufuncs import TensorSpaceUfuncs + dtype_str, signature_string) from odl.util.utility import TYPE_MAP_R2C, TYPE_MAP_C2R @@ -507,403 +505,11 @@ class Tensor(LinearSpaceElement): """Abstract class for representation of `TensorSpace` elements.""" - def asarray(self, out=None): - """Extract the data of this tensor as a Numpy array. - - This method should be overridden by subclasses. - - Parameters - ---------- - out : `numpy.ndarray`, optional - Array to write the result to. - - Returns - ------- - asarray : `numpy.ndarray` - Numpy array of the same data type and shape as the space. - If ``out`` was given, the returned object is a reference - to it. - """ - raise NotImplementedError('abstract method') - - def __getitem__(self, indices): - """Return ``self[indices]``. - - This method should be overridden by subclasses. - - Parameters - ---------- - indices : index expression - Integer, slice or sequence of these, defining the positions - of the data array which should be accessed. - - Returns - ------- - values : `TensorSpace.dtype` or `Tensor` - The value(s) at the given indices. Note that depending on - the implementation, the returned object may be a (writable) - view into the original array. - """ - raise NotImplementedError('abstract method') - - def __setitem__(self, indices, values): - """Implement ``self[indices] = values``. - - This method should be overridden by subclasses. - - Parameters - ---------- - indices : index expression - Integer, slice or sequence of these, defining the positions - of the data array which should be written to. - values : scalar, `array-like` or `Tensor` - The value(s) that are to be assigned. - - If ``index`` is an integer, ``value`` must be a scalar. - - If ``index`` is a slice or a sequence of slices, ``value`` - must be broadcastable to the shape of the slice. - """ - raise NotImplementedError('abstract method') - @property def impl(self): """Name of the implementation back-end of this tensor.""" return self.space.impl - @property - def shape(self): - """Number of elements per axis.""" - return self.space.shape - - @property - def dtype(self): - """Data type of each entry.""" - return self.space.dtype - - @property - def size(self): - """Total number of entries.""" - return self.space.size - - @property - def ndim(self): - """Number of axes (=dimensions) of this tensor.""" - return self.space.ndim - - def __len__(self): - """Return ``len(self)``. - - The length is equal to the number of entries along axis 0. - """ - return len(self.space) - - @property - def itemsize(self): - """Size in bytes of one tensor entry.""" - return self.space.itemsize - - @property - def nbytes(self): - """Total number of bytes in memory occupied by this tensor.""" - return self.space.nbytes - - def astype(self, dtype): - """Return a copy of this element with new ``dtype``. - - Parameters - ---------- - dtype : - Scalar data type of the returned space. Can be provided - in any way the `numpy.dtype` constructor understands, e.g. - as built-in type or as a string. Data types with non-trivial - shapes are not allowed. - - Returns - ------- - newelem : `Tensor` - Version of this element with given data type. - """ - raise NotImplementedError('abstract method') - - def __repr__(self): - """Return ``repr(self)``.""" - maxsize_full_print = 2 * np.get_printoptions()['edgeitems'] - self_str = array_str(self, nprint=maxsize_full_print) - if self.ndim == 1 and self.size <= maxsize_full_print: - return '{!r}.element({})'.format(self.space, self_str) - else: - return '{!r}.element(\n{}\n)'.format(self.space, indent(self_str)) - - def __str__(self): - """Return ``str(self)``.""" - return array_str(self) - - def __bool__(self): - """Return ``bool(self)``.""" - if self.size > 1: - raise ValueError('The truth value of an array with more than one ' - 'element is ambiguous. ' - 'Use np.any(a) or np.all(a)') - else: - return bool(self.asarray()) - - def __array__(self, dtype=None): - """Return a Numpy array from this tensor. - - Parameters - ---------- - dtype : - Specifier for the data type of the output array. - - Returns - ------- - array : `numpy.ndarray` - """ - if dtype is None: - return self.asarray() - else: - return self.asarray().astype(dtype, copy=False) - - def __array_wrap__(self, array): - """Return a new tensor wrapping the ``array``. - - Parameters - ---------- - array : `numpy.ndarray` - Array to be wrapped. - - Returns - ------- - wrapper : `Tensor` - Tensor wrapping ``array``. - """ - if array.ndim == 0: - return self.space.field.element(array) - else: - return self.space.element(array) - - def __array_ufunc__(self, ufunc, method, *inputs, **kwargs): - """Interface to Numpy's ufunc machinery. - - This method is called by Numpy version 1.13 and higher as a single - point for the ufunc dispatch logic. An object implementing - ``__array_ufunc__`` takes over control when a `numpy.ufunc` is - called on it, allowing it to use custom implementations and - output types. - - This includes handling of in-place arithmetic like - ``npy_array += custom_obj``. In this case, the custom object's - ``__array_ufunc__`` takes precedence over the baseline - `numpy.ndarray` implementation. It will be called with - ``npy_array`` as ``out`` argument, which ensures that the - returned object is a Numpy array. For this to work properly, - ``__array_ufunc__`` has to accept Numpy arrays as ``out`` arguments. - - See the `corresponding NEP`_ and the `interface documentation`_ - for further details. See also the `general documentation on - Numpy ufuncs`_. - - .. note:: - This basic implementation casts inputs and - outputs to Numpy arrays and evaluates ``ufunc`` on those. - For `numpy.ndarray` based data storage, this incurs no - significant overhead compared to direct usage of Numpy arrays. - - For other (in particular non-local) implementations, e.g., - GPU arrays or distributed memory, overhead is significant due - to copies to CPU main memory. In those classes, the - ``__array_ufunc__`` mechanism should be overridden in favor of - a native implementations if possible. - - .. note:: - If no ``out`` parameter is provided, this implementation - just returns the raw array and does not attempt to wrap the - result in any kind of space. - - Parameters - ---------- - ufunc : `numpy.ufunc` - Ufunc that should be called on ``self``. - method : str - Method on ``ufunc`` that should be called on ``self``. - Possible values: - - ``'__call__'``, ``'accumulate'``, ``'at'``, ``'outer'``, - ``'reduce'``, ``'reduceat'`` - - input1, ..., inputN: - Positional arguments to ``ufunc.method``. - kwargs: - Keyword arguments to ``ufunc.method``. - - Returns - ------- - ufunc_result : `Tensor`, `numpy.ndarray` or tuple - Result of the ufunc evaluation. If no ``out`` keyword argument - was given, the result is a `Tensor` or a tuple - of such, depending on the number of outputs of ``ufunc``. - If ``out`` was provided, the returned object or tuple entries - refer(s) to ``out``. - - References - ---------- - .. _corresponding NEP: - https://docs.scipy.org/doc/numpy/neps/ufunc-overrides.html - - .. _interface documentation: - https://docs.scipy.org/doc/numpy/reference/arrays.classes.html\ -#numpy.class.__array_ufunc__ - - .. _general documentation on Numpy ufuncs: - https://docs.scipy.org/doc/numpy/reference/ufuncs.html - - .. _reduceat documentation: - https://docs.scipy.org/doc/numpy/reference/generated/\ -numpy.ufunc.reduceat.html - """ - # --- Process `out` --- # - - # Unwrap out if provided. The output parameters are all wrapped - # in one tuple, even if there is only one. - out_tuple = kwargs.pop('out', ()) - - # Check number of `out` args, depending on `method` - if method == '__call__' and len(out_tuple) not in (0, ufunc.nout): - raise ValueError( - "ufunc {}: need 0 or {} `out` arguments for " - "`method='__call__'`, got {}" - ''.format(ufunc.__name__, ufunc.nout, len(out_tuple))) - elif method != '__call__' and len(out_tuple) not in (0, 1): - raise ValueError( - 'ufunc {}: need 0 or 1 `out` arguments for `method={!r}`, ' - 'got {}'.format(ufunc.__name__, method, len(out_tuple))) - - # We allow our own tensors, the data container type and - # `numpy.ndarray` objects as `out` (see docs for reason for the - # latter) - valid_types = (type(self), type(self.data), np.ndarray) - if not all(isinstance(o, valid_types) or o is None - for o in out_tuple): - return NotImplemented - - # Assign to `out` or `out1` and `out2`, respectively - out = out1 = out2 = None - if len(out_tuple) == 1: - out = out_tuple[0] - elif len(out_tuple) == 2: - out1 = out_tuple[0] - out2 = out_tuple[1] - - # --- Process `inputs` --- # - - # Convert inputs that are ODL tensors or their data containers to - # Numpy arrays so that the native Numpy ufunc is called later - inputs = tuple( - np.asarray(inp) if isinstance(inp, (type(self), type(self.data))) - else inp - for inp in inputs) - - # --- Get some parameters for later --- # - - # Arguments for `writable_array` and/or space constructors - out_dtype = kwargs.get('dtype', None) - if out_dtype is None: - array_kwargs = {} - else: - array_kwargs = {'dtype': out_dtype} - - # --- Evaluate ufunc --- # - - # Trivial context used to create a single code path for the ufunc - # evaluation. For `None` output parameter(s), this is used instead of - # `writable_array`. - class CtxNone(object): - """Trivial context manager class. - - When used as :: - - with CtxNone() as obj: - # do stuff with `obj` - - the returned ``obj`` is ``None``. - """ - __enter__ = __exit__ = lambda *_: None - - if method == '__call__': - if ufunc.nout == 1: - # Make context for output (trivial one returns `None`) - if out is None: - out_ctx = CtxNone() - else: - out_ctx = writable_array(out, **array_kwargs) - - # Evaluate ufunc - with out_ctx as out_arr: - kwargs['out'] = out_arr - res = ufunc(*inputs, **kwargs) - - # Return result (may be a raw array or a space element) - return res - - elif ufunc.nout == 2: - # Make contexts for outputs (trivial ones return `None`) - if out1 is not None: - out1_ctx = writable_array(out1, **array_kwargs) - else: - out1_ctx = CtxNone() - if out2 is not None: - out2_ctx = writable_array(out2, **array_kwargs) - else: - out2_ctx = CtxNone() - - # Evaluate ufunc - with out1_ctx as out1_arr, out2_ctx as out2_arr: - kwargs['out'] = (out1_arr, out2_arr) - res1, res2 = ufunc(*inputs, **kwargs) - - # Return results (may be raw arrays or space elements) - return res1, res2 - - else: - raise NotImplementedError('nout = {} not supported' - ''.format(ufunc.nout)) - - else: # method != '__call__' - # Make context for output (trivial one returns `None`) - if out is None: - out_ctx = CtxNone() - else: - out_ctx = writable_array(out, **array_kwargs) - - # Evaluate ufunc method - if method == 'at': - with writable_array(inputs[0]) as inp_arr: - res = ufunc.at(inp_arr, *inputs[1:], **kwargs) - else: - with out_ctx as out_arr: - kwargs['out'] = out_arr - res = getattr(ufunc, method)(*inputs, **kwargs) - - # Return result (may be scalar, raw array or space element) - return res - - # Old ufuncs interface, will be deprecated when Numpy 1.13 becomes minimum - - @property - def ufuncs(self): - """Access to Numpy style universal functions. - - These default ufuncs are always available, but may or may not be - optimized for the specific space in use. - - .. note:: - This interface is will be deprecated when Numpy 1.13 becomes - the minimum required version. Use Numpy ufuncs directly, e.g., - ``np.sqrt(x)`` instead of ``x.ufuncs.sqrt()``. - """ - return TensorSpaceUfuncs(self) - def show(self, title=None, method='', indices=None, force_show=False, fig=None, **kwargs): """Display the function graphically. diff --git a/odl/space/npy_tensors.py b/odl/space/npy_tensors.py index 56da33ad4ed..ca4c395c168 100644 --- a/odl/space/npy_tensors.py +++ b/odl/space/npy_tensors.py @@ -404,7 +404,7 @@ def element(self, inp=None, data_ptr=None, order=None): else: arr = np.empty(self.shape, dtype=self.dtype, order=order) - return self.element_type(self, arr) + return self.element_type(arr) elif inp is None and data_ptr is not None: if order is None: @@ -416,7 +416,7 @@ def element(self, inp=None, data_ptr=None, order=None): as_numpy_array = np.ctypeslib.as_array(as_ctype_array) arr = as_numpy_array.view(dtype=self.dtype) arr = arr.reshape(self.shape, order=order) - return self.element_type(self, arr) + return self.element_type(arr) elif inp is not None and data_ptr is None: if inp in self and order is None: @@ -425,8 +425,9 @@ def element(self, inp=None, data_ptr=None, order=None): # Try to not copy but require dtype and order if given # (`order=None` is ok as np.array argument) - arr = np.array(inp, copy=False, dtype=self.dtype, ndmin=self.ndim, - order=order) + arr = np.array( + inp, copy=False, dtype=self.dtype, ndmin=self.ndim, order=order + ) # Make sure the result is writeable, if not make copy. # This happens for e.g. results of `np.broadcast_to()`. if not arr.flags.writeable: @@ -434,7 +435,7 @@ def element(self, inp=None, data_ptr=None, order=None): if arr.shape != self.shape: raise ValueError('shape of `inp` not equal to space shape: ' '{} != {}'.format(arr.shape, self.shape)) - return self.element_type(self, arr) + return self.element_type(arr) else: raise TypeError('cannot provide both `inp` and `data_ptr`') @@ -846,924 +847,28 @@ def __repr__(self): @property def element_type(self): """Type of elements in this space: `NumpyTensor`.""" - return NumpyTensor - - -class NumpyTensor(Tensor): - - """Representation of a `NumpyTensorSpace` element.""" - - def __init__(self, space, data): - """Initialize a new instance.""" - Tensor.__init__(self, space) - self.__data = data - - @property - def data(self): - """The `numpy.ndarray` representing the data of ``self``.""" - return self.__data - - def asarray(self, out=None): - """Extract the data of this array as a ``numpy.ndarray``. - - This method is invoked when calling `numpy.asarray` on this - tensor. - - Parameters - ---------- - out : `numpy.ndarray`, optional - Array in which the result should be written in-place. - Has to be contiguous and of the correct dtype. - - Returns - ------- - asarray : `numpy.ndarray` - Numpy array with the same data type as ``self``. If - ``out`` was given, the returned object is a reference - to it. - - Examples - -------- - >>> space = odl.rn(3, dtype='float32') - >>> x = space.element([1, 2, 3]) - >>> x.asarray() - array([ 1., 2., 3.], dtype=float32) - >>> np.asarray(x) is x.asarray() - True - >>> out = np.empty(3, dtype='float32') - >>> result = x.asarray(out=out) - >>> out - array([ 1., 2., 3.], dtype=float32) - >>> result is out - True - >>> space = odl.rn((2, 3)) - >>> space.one().asarray() - array([[ 1., 1., 1.], - [ 1., 1., 1.]]) - """ - if out is None: - return self.data - else: - out[:] = self.data - return out - - def astype(self, dtype): - """Return a copy of this element with new ``dtype``. - - Parameters - ---------- - dtype : - Scalar data type of the returned space. Can be provided - in any way the `numpy.dtype` constructor understands, e.g. - as built-in type or as a string. Data types with non-trivial - shapes are not allowed. - - Returns - ------- - newelem : `NumpyTensor` - Version of this element with given data type. - """ - return self.space.astype(dtype).element(self.data.astype(dtype)) - - @property - def data_ptr(self): - """A raw pointer to the data container of ``self``. - - Examples - -------- - >>> import ctypes - >>> space = odl.tensor_space(3, dtype='uint16') - >>> x = space.element([1, 2, 3]) - >>> arr_type = ctypes.c_uint16 * 3 # C type "array of 3 uint16" - >>> buffer = arr_type.from_address(x.data_ptr) - >>> arr = np.frombuffer(buffer, dtype='uint16') - >>> arr - array([1, 2, 3], dtype=uint16) - - In-place modification via pointer: - - >>> arr[0] = 42 - >>> x - tensor_space(3, dtype='uint16').element([42, 2, 3]) - """ - return self.data.ctypes.data - - def __eq__(self, other): - """Return ``self == other``. - - Returns - ------- - equals : bool - True if all entries of ``other`` are equal to this - the entries of ``self``, False otherwise. - - Examples - -------- - >>> space = odl.rn(3) - >>> x = space.element([1, 2, 3]) - >>> y = space.element([1, 2, 3]) - >>> x == y - True - - >>> y = space.element([-1, 2, 3]) - >>> x == y - False - >>> x == object - False - - Space membership matters: - - >>> space2 = odl.tensor_space(3, dtype='int64') - >>> y = space2.element([1, 2, 3]) - >>> x == y or y == x - False - """ - if other is self: - return True - elif other not in self.space: - return False - else: - return np.array_equal(self.data, other.data) - - def copy(self): - """Return an identical (deep) copy of this tensor. - - Parameters - ---------- - None - - Returns - ------- - copy : `NumpyTensor` - The deep copy - - Examples - -------- - >>> space = odl.rn(3) - >>> x = space.element([1, 2, 3]) - >>> y = x.copy() - >>> y == x - True - >>> y is x - False - """ - return self.space.element(self.data.copy()) - - def __copy__(self): - """Return ``copy(self)``. - - This implements the (shallow) copy interface of the ``copy`` - module of the Python standard library. - - See Also - -------- - copy - - Examples - -------- - >>> from copy import copy - >>> space = odl.rn(3) - >>> x = space.element([1, 2, 3]) - >>> y = copy(x) - >>> y == x - True - >>> y is x - False - """ - return self.copy() - - def __getitem__(self, indices): - """Return ``self[indices]``. - - Parameters - ---------- - indices : index expression - Integer, slice or sequence of these, defining the positions - of the data array which should be accessed. - - Returns - ------- - values : `NumpyTensorSpace.dtype` or `NumpyTensor` - The value(s) at the given indices. Note that the returned - object is a writable view into the original tensor, except - for the case when ``indices`` is a list. - - Examples - -------- - For one-dimensional spaces, indexing is as in linear arrays: + space = self - >>> space = odl.rn(3) - >>> x = space.element([1, 2, 3]) - >>> x[0] - 1.0 - >>> x[1:] - rn(2).element([ 2., 3.]) - - In higher dimensions, the i-th index expression accesses the - i-th axis: - - >>> space = odl.rn((2, 3)) - >>> x = space.element([[1, 2, 3], - ... [4, 5, 6]]) - >>> x[0, 1] - 2.0 - >>> x[:, 1:] - rn((2, 2)).element( - [[ 2., 3.], - [ 5., 6.]] - ) + class NumpyTensor(Tensor, np.ndarray): - Slices can be assigned to, except if lists are used for indexing: + """Representation of a `NumpyTensorSpace` element.""" - >>> y = x[:, ::2] # view into x - >>> y[:] = -9 - >>> x - rn((2, 3)).element( - [[-9., 2., -9.], - [-9., 5., -9.]] - ) - >>> y = x[[0, 1], [1, 2]] # not a view, won't modify x - >>> y - rn(2).element([ 2., -9.]) - >>> y[:] = 0 - >>> x - rn((2, 3)).element( - [[-9., 2., -9.], - [-9., 5., -9.]] - ) - """ - # Lazy implementation: index the array and deal with it - if isinstance(indices, NumpyTensor): - indices = indices.data - arr = self.data[indices] - - if np.isscalar(arr): - if self.space.field is not None: - return self.space.field.element(arr) - else: + def __new__(cls, data): + arr = np.asarray(data).view(cls) + arr.__space = space return arr - else: - if is_numeric_dtype(self.dtype): - weighting = self.space.weighting - else: - weighting = None - space = type(self.space)( - arr.shape, dtype=self.dtype, exponent=self.space.exponent, - weighting=weighting) - return space.element(arr) - - def __setitem__(self, indices, values): - """Implement ``self[indices] = values``. - - Parameters - ---------- - indices : index expression - Integer, slice or sequence of these, defining the positions - of the data array which should be written to. - values : scalar, array-like or `NumpyTensor` - The value(s) that are to be assigned. - - If ``index`` is an integer, ``value`` must be a scalar. - - If ``index`` is a slice or a sequence of slices, ``value`` - must be broadcastable to the shape of the slice. - - Examples - -------- - For 1d spaces, entries can be set with scalars or sequences of - correct shape: - - >>> space = odl.rn(3) - >>> x = space.element([1, 2, 3]) - >>> x[0] = -1 - >>> x[1:] = (0, 1) - >>> x - rn(3).element([-1., 0., 1.]) - - It is also possible to use tensors of other spaces for - casting and assignment: - - >>> space = odl.rn((2, 3)) - >>> x = space.element([[1, 2, 3], - ... [4, 5, 6]]) - >>> x[0, 1] = -1 - >>> x - rn((2, 3)).element( - [[ 1., -1., 3.], - [ 4., 5., 6.]] - ) - >>> short_space = odl.tensor_space((2, 2), dtype='short') - >>> y = short_space.element([[-1, 2], - ... [0, 0]]) - >>> x[:, :2] = y - >>> x - rn((2, 3)).element( - [[-1., 2., 3.], - [ 0., 0., 6.]] - ) - - The Numpy assignment and broadcasting rules apply: - - >>> x[:] = np.array([[0, 0, 0], - ... [1, 1, 1]]) - >>> x - rn((2, 3)).element( - [[ 0., 0., 0.], - [ 1., 1., 1.]] - ) - >>> x[:, 1:] = [7, 8] - >>> x - rn((2, 3)).element( - [[ 0., 7., 8.], - [ 1., 7., 8.]] - ) - >>> x[:, ::2] = -2. - >>> x - rn((2, 3)).element( - [[-2., 7., -2.], - [-2., 7., -2.]] - ) - """ - if isinstance(indices, type(self)): - indices = indices.data - if isinstance(values, type(self)): - values = values.data - - self.data[indices] = values - - @property - def real(self): - """Real part of ``self``. - - Returns - ------- - real : `NumpyTensor` - Real part of this element as a member of a - `NumpyTensorSpace` with corresponding real data type. - - Examples - -------- - Get the real part: - - >>> space = odl.cn(3) - >>> x = space.element([1 + 1j, 2, 3 - 3j]) - >>> x.real - rn(3).element([ 1., 2., 3.]) - - Set the real part: - - >>> space = odl.cn(3) - >>> x = space.element([1 + 1j, 2, 3 - 3j]) - >>> zero = odl.rn(3).zero() - >>> x.real = zero - >>> x - cn(3).element([ 0.+1.j, 0.+0.j, 0.-3.j]) - - Other array-like types and broadcasting: - - >>> x.real = 1.0 - >>> x - cn(3).element([ 1.+1.j, 1.+0.j, 1.-3.j]) - >>> x.real = [2, 3, 4] - >>> x - cn(3).element([ 2.+1.j, 3.+0.j, 4.-3.j]) - """ - if self.space.is_real: - return self - elif self.space.is_complex: - real_space = self.space.astype(self.space.real_dtype) - return real_space.element(self.data.real) - else: - raise NotImplementedError('`real` not defined for non-numeric ' - 'dtype {}'.format(self.dtype)) - - @real.setter - def real(self, newreal): - """Setter for the real part. - - This method is invoked by ``x.real = other``. - Parameters - ---------- - newreal : array-like or scalar - Values to be assigned to the real part of this element. - """ - self.real.data[:] = newreal + def __init__(self, data): + pass - @property - def imag(self): - """Imaginary part of ``self``. + def __array_finalize__(self, obj): + self.__space = NumpyTensorSpace(self.shape, self.dtype) - Returns - ------- - imag : `NumpyTensor` - Imaginary part this element as an element of a - `NumpyTensorSpace` with real data type. - - Examples - -------- - Get the imaginary part: + @property + def space(self): + return self.__space - >>> space = odl.cn(3) - >>> x = space.element([1 + 1j, 2, 3 - 3j]) - >>> x.imag - rn(3).element([ 1., 0., -3.]) - - Set the imaginary part: - - >>> space = odl.cn(3) - >>> x = space.element([1 + 1j, 2, 3 - 3j]) - >>> zero = odl.rn(3).zero() - >>> x.imag = zero - >>> x - cn(3).element([ 1.+0.j, 2.+0.j, 3.+0.j]) - - Other array-like types and broadcasting: - - >>> x.imag = 1.0 - >>> x - cn(3).element([ 1.+1.j, 2.+1.j, 3.+1.j]) - >>> x.imag = [2, 3, 4] - >>> x - cn(3).element([ 1.+2.j, 2.+3.j, 3.+4.j]) - """ - if self.space.is_real: - return self.space.zero() - elif self.space.is_complex: - real_space = self.space.astype(self.space.real_dtype) - return real_space.element(self.data.imag) - else: - raise NotImplementedError('`imag` not defined for non-numeric ' - 'dtype {}'.format(self.dtype)) - - @imag.setter - def imag(self, newimag): - """Setter for the imaginary part. - - This method is invoked by ``x.imag = other``. - - Parameters - ---------- - newimag : array-like or scalar - Values to be assigned to the imaginary part of this element. - - Raises - ------ - ValueError - If the space is real, i.e., no imagninary part can be set. - """ - if self.space.is_real: - raise ValueError('cannot set imaginary part in real spaces') - self.imag.data[:] = newimag - - def conj(self, out=None): - """Return the complex conjugate of ``self``. - - Parameters - ---------- - out : `NumpyTensor`, optional - Element to which the complex conjugate is written. - Must be an element of ``self.space``. - - Returns - ------- - out : `NumpyTensor` - The complex conjugate element. If ``out`` was provided, - the returned object is a reference to it. - - Examples - -------- - >>> space = odl.cn(3) - >>> x = space.element([1 + 1j, 2, 3 - 3j]) - >>> x.conj() - cn(3).element([ 1.-1.j, 2.-0.j, 3.+3.j]) - >>> out = space.element() - >>> result = x.conj(out=out) - >>> result - cn(3).element([ 1.-1.j, 2.-0.j, 3.+3.j]) - >>> result is out - True - - In-place conjugation: - - >>> result = x.conj(out=x) - >>> x - cn(3).element([ 1.-1.j, 2.-0.j, 3.+3.j]) - >>> result is x - True - """ - if self.space.is_real: - if out is None: - return self - else: - out[:] = self - return out - - if not is_numeric_dtype(self.space.dtype): - raise NotImplementedError('`conj` not defined for non-numeric ' - 'dtype {}'.format(self.dtype)) - - if out is None: - return self.space.element(self.data.conj()) - else: - if out not in self.space: - raise LinearSpaceTypeError('`out` {!r} not in space {!r}' - ''.format(out, self.space)) - self.data.conj(out.data) - return out - - def __ipow__(self, other): - """Return ``self **= other``.""" - try: - if other == int(other): - return super(NumpyTensor, self).__ipow__(other) - except TypeError: - pass - - np.power(self.data, other, out=self.data) - return self - - def __int__(self): - """Return ``int(self)``.""" - return int(self.data) - - def __long__(self): - """Return ``long(self)``. - - This method is only useful in Python 2. - """ - return long(self.data) - - def __float__(self): - """Return ``float(self)``.""" - return float(self.data) - - def __complex__(self): - """Return ``complex(self)``.""" - if self.size != 1: - raise TypeError('only size-1 tensors can be converted to ' - 'Python scalars') - return complex(self.data.ravel()[0]) - - def __array_ufunc__(self, ufunc, method, *inputs, **kwargs): - """Interface to Numpy's ufunc machinery. - - This method is called by Numpy version 1.13 and higher as a single - point for the ufunc dispatch logic. An object implementing - ``__array_ufunc__`` takes over control when a `numpy.ufunc` is - called on it, allowing it to use custom implementations and - output types. - - This includes handling of in-place arithmetic like - ``npy_array += custom_obj``. In this case, the custom object's - ``__array_ufunc__`` takes precedence over the baseline - `numpy.ndarray` implementation. It will be called with - ``npy_array`` as ``out`` argument, which ensures that the - returned object is a Numpy array. For this to work properly, - ``__array_ufunc__`` has to accept Numpy arrays as ``out`` arguments. - - See the `corresponding NEP`_ and the `interface documentation`_ - for further details. See also the `general documentation on - Numpy ufuncs`_. - - .. note:: - This basic implementation casts inputs and - outputs to Numpy arrays and evaluates ``ufunc`` on those. - For `numpy.ndarray` based data storage, this incurs no - significant overhead compared to direct usage of Numpy arrays. - - For other (in particular non-local) implementations, e.g., - GPU arrays or distributed memory, overhead is significant due - to copies to CPU main memory. In those classes, the - ``__array_ufunc__`` mechanism should be overridden to use - native implementations if possible. - - .. note:: - When using operations that alter the shape (like ``reduce``), - or the data type (can be any of the methods), - the resulting array is wrapped in a space of the same - type as ``self.space``, propagating space properties like - `exponent` or `weighting` as closely as possible. - - Parameters - ---------- - ufunc : `numpy.ufunc` - Ufunc that should be called on ``self``. - method : str - Method on ``ufunc`` that should be called on ``self``. - Possible values: - - ``'__call__'``, ``'accumulate'``, ``'at'``, ``'outer'``, - ``'reduce'``, ``'reduceat'`` - - input1, ..., inputN : - Positional arguments to ``ufunc.method``. - kwargs : - Keyword arguments to ``ufunc.method``. - - Returns - ------- - ufunc_result : `Tensor`, `numpy.ndarray` or tuple - Result of the ufunc evaluation. If no ``out`` keyword argument - was given, the result is a `Tensor` or a tuple - of such, depending on the number of outputs of ``ufunc``. - If ``out`` was provided, the returned object or tuple entries - refer(s) to ``out``. - - Examples - -------- - We apply `numpy.add` to ODL tensors: - - >>> r3 = odl.rn(3) - >>> x = r3.element([1, 2, 3]) - >>> y = r3.element([-1, -2, -3]) - >>> x.__array_ufunc__(np.add, '__call__', x, y) - rn(3).element([ 0., 0., 0.]) - >>> np.add(x, y) # same mechanism for Numpy >= 1.13 - rn(3).element([ 0., 0., 0.]) - - As ``out``, a Numpy array or an ODL tensor can be given (wrapped - in a sequence): - - >>> out = r3.element() - >>> res = x.__array_ufunc__(np.add, '__call__', x, y, out=(out,)) - >>> out - rn(3).element([ 0., 0., 0.]) - >>> res is out - True - >>> out_arr = np.empty(3) - >>> res = x.__array_ufunc__(np.add, '__call__', x, y, out=(out_arr,)) - >>> out_arr - array([ 0., 0., 0.]) - >>> res is out_arr - True - - With multiple dimensions: - - >>> r23 = odl.rn((2, 3)) - >>> x = y = r23.one() - >>> x.__array_ufunc__(np.add, '__call__', x, y) - rn((2, 3)).element( - [[ 2., 2., 2.], - [ 2., 2., 2.]] - ) - - The ``ufunc.accumulate`` method retains the original `shape` and - `dtype`. The latter can be changed with the ``dtype`` parameter: - - >>> x = r3.element([1, 2, 3]) - >>> x.__array_ufunc__(np.add, 'accumulate', x) - rn(3).element([ 1., 3., 6.]) - >>> np.add.accumulate(x) # same mechanism for Numpy >= 1.13 - rn(3).element([ 1., 3., 6.]) - >>> x.__array_ufunc__(np.add, 'accumulate', x, dtype=complex) - cn(3).element([ 1.+0.j, 3.+0.j, 6.+0.j]) - - For multi-dimensional tensors, an optional ``axis`` parameter - can be provided: - - >>> z = r23.one() - >>> z.__array_ufunc__(np.add, 'accumulate', z, axis=1) - rn((2, 3)).element( - [[ 1., 2., 3.], - [ 1., 2., 3.]] - ) - - The ``ufunc.at`` method operates in-place. Here we add the second - operand ``[5, 10]`` to ``x`` at indices ``[0, 2]``: - - >>> x = r3.element([1, 2, 3]) - >>> x.__array_ufunc__(np.add, 'at', x, [0, 2], [5, 10]) - >>> x - rn(3).element([ 6., 2., 13.]) - - For outer-product-type operations, i.e., operations where the result - shape is the sum of the individual shapes, the ``ufunc.outer`` - method can be used: - - >>> x = odl.rn(2).element([0, 3]) - >>> y = odl.rn(3).element([1, 2, 3]) - >>> x.__array_ufunc__(np.add, 'outer', x, y) - rn((2, 3)).element( - [[ 1., 2., 3.], - [ 4., 5., 6.]] - ) - >>> y.__array_ufunc__(np.add, 'outer', y, x) - rn((3, 2)).element( - [[ 1., 4.], - [ 2., 5.], - [ 3., 6.]] - ) - - Using ``ufunc.reduce`` produces a scalar, which can be avoided with - ``keepdims=True``: - - >>> x = r3.element([1, 2, 3]) - >>> x.__array_ufunc__(np.add, 'reduce', x) - 6.0 - >>> x.__array_ufunc__(np.add, 'reduce', x, keepdims=True) - rn(1).element([ 6.]) - - In multiple dimensions, ``axis`` can be provided for reduction over - selected axes: - - >>> z = r23.element([[1, 2, 3], - ... [4, 5, 6]]) - >>> z.__array_ufunc__(np.add, 'reduce', z, axis=1) - rn(2).element([ 6., 15.]) - - Finally, ``add.reduceat`` is a combination of ``reduce`` and - ``at`` with rather flexible and complex semantics (see the - `reduceat documentation`_ for details): - - >>> x = r3.element([1, 2, 3]) - >>> x.__array_ufunc__(np.add, 'reduceat', x, [0, 1]) - rn(2).element([ 1., 5.]) - - References - ---------- - .. _corresponding NEP: - https://docs.scipy.org/doc/numpy/neps/ufunc-overrides.html - - .. _interface documentation: - https://docs.scipy.org/doc/numpy/reference/arrays.classes.html\ -#numpy.class.__array_ufunc__ - - .. _general documentation on Numpy ufuncs: - https://docs.scipy.org/doc/numpy/reference/ufuncs.html - - .. _reduceat documentation: - https://docs.scipy.org/doc/numpy/reference/generated/\ -numpy.ufunc.reduceat.html - """ - # Remark: this method differs from the parent implementation only - # in the propagation of additional space properties. - - # --- Process `out` --- # - - # Unwrap out if provided. The output parameters are all wrapped - # in one tuple, even if there is only one. - out_tuple = kwargs.pop('out', ()) - - # Check number of `out` args, depending on `method` - if method == '__call__' and len(out_tuple) not in (0, ufunc.nout): - raise ValueError( - "ufunc {}: need 0 or {} `out` arguments for " - "`method='__call__'`, got {}" - ''.format(ufunc.__name__, ufunc.nout, len(out_tuple))) - elif method != '__call__' and len(out_tuple) not in (0, 1): - raise ValueError( - 'ufunc {}: need 0 or 1 `out` arguments for `method={!r}`, ' - 'got {}'.format(ufunc.__name__, method, len(out_tuple))) - - # We allow our own tensors, the data container type and - # `numpy.ndarray` objects as `out` (see docs for reason for the - # latter) - valid_types = (type(self), type(self.data), np.ndarray) - if not all(isinstance(o, valid_types) or o is None - for o in out_tuple): - return NotImplemented - - # Assign to `out` or `out1` and `out2`, respectively - out = out1 = out2 = None - if len(out_tuple) == 1: - out = out_tuple[0] - elif len(out_tuple) == 2: - out1 = out_tuple[0] - out2 = out_tuple[1] - - # --- Process `inputs` --- # - - # Convert inputs that are ODL tensors to Numpy arrays so that the - # native Numpy ufunc is called later - inputs = tuple( - inp.asarray() if isinstance(inp, type(self)) else inp - for inp in inputs) - - # --- Get some parameters for later --- # - - # Arguments for `writable_array` and/or space constructors - out_dtype = kwargs.get('dtype', None) - if out_dtype is None: - array_kwargs = {} - else: - array_kwargs = {'dtype': out_dtype} - - exponent = self.space.exponent - weighting = self.space.weighting - - # --- Evaluate ufunc --- # - - # Trivial context used to create a single code path for the ufunc - # evaluation. For `None` output parameter(s), this is used instead of - # `writable_array`. - class CtxNone(object): - """Trivial context manager class. - - When used as :: - - with CtxNone() as obj: - # do stuff with `obj` - - the returned ``obj`` is ``None``. - """ - __enter__ = __exit__ = lambda *_: None - - if method == '__call__': - if ufunc.nout == 1: - # Make context for output (trivial one returns `None`) - if out is None: - out_ctx = CtxNone() - else: - out_ctx = writable_array(out, **array_kwargs) - - # Evaluate ufunc - with out_ctx as out_arr: - kwargs['out'] = out_arr - res = ufunc(*inputs, **kwargs) - - # Wrap result if necessary (lazily) - if out is None: - if is_floating_dtype(res.dtype): - # Weighting contains exponent - spc_kwargs = {'weighting': weighting} - else: - # No `exponent` or `weighting` applicable - spc_kwargs = {} - out_space = type(self.space)(self.shape, res.dtype, - **spc_kwargs) - out = out_space.element(res) - - return out - - elif ufunc.nout == 2: - # Make contexts for outputs (trivial ones return `None`) - if out1 is not None: - out1_ctx = writable_array(out1, **array_kwargs) - else: - out1_ctx = CtxNone() - if out2 is not None: - out2_ctx = writable_array(out2, **array_kwargs) - else: - out2_ctx = CtxNone() - - # Evaluate ufunc - with out1_ctx as out1_arr, out2_ctx as out2_arr: - kwargs['out'] = (out1_arr, out2_arr) - res1, res2 = ufunc(*inputs, **kwargs) - - # Wrap results if necessary (lazily) - # We don't use exponents or weightings since we don't know - # how to map them to the spaces - if out1 is None: - out1_space = type(self.space)(self.shape, res1.dtype) - out1 = out1_space.element(res1) - if out2 is None: - out2_space = type(self.space)(self.shape, res2.dtype) - out2 = out2_space.element(res2) - - return out1, out2 - - else: - raise NotImplementedError('nout = {} not supported' - ''.format(ufunc.nout)) - - else: # method != '__call__' - # Make context for output (trivial one returns `None`) - if out is None: - out_ctx = CtxNone() - else: - out_ctx = writable_array(out, **array_kwargs) - - # Evaluate ufunc method - with out_ctx as out_arr: - if method != 'at': - # No kwargs allowed for 'at' - kwargs['out'] = out_arr - res = getattr(ufunc, method)(*inputs, **kwargs) - - # Shortcut for scalar or no return value - if np.isscalar(res) or res is None: - # The first occurs for `reduce` with all axes, - # the second for in-place stuff (`at` currently) - return res - - # Wrap result if necessary (lazily) - if out is None: - if is_floating_dtype(res.dtype): - if res.shape != self.shape: - # Don't propagate weighting if shape changes - weighting = NumpyTensorSpaceConstWeighting(1.0, - exponent) - spc_kwargs = {'weighting': weighting} - else: - spc_kwargs = {} - - out_space = type(self.space)(res.shape, res.dtype, - **spc_kwargs) - out = out_space.element(res) - - return out + return NumpyTensor def _blas_is_applicable(*args): @@ -1808,11 +913,11 @@ def _lincomb_impl(a, x1, b, x2, out): if size < THRESHOLD_SMALL: # Faster for small arrays - out.data[:] = a * x1.data + b * x2.data + out.data[:] = a * x1 + b * x2 return elif (size < THRESHOLD_MEDIUM or - not _blas_is_applicable(x1.data, x2.data, out.data)): + not _blas_is_applicable(x1, x2, out)): def fallback_axpy(x1, x2, n, a): """Fallback axpy implementation avoiding copy.""" @@ -1833,23 +938,23 @@ def fallback_copy(x1, x2, n): return x2 axpy, scal, copy = (fallback_axpy, fallback_scal, fallback_copy) - x1_arr = x1.data - x2_arr = x2.data - out_arr = out.data + x1_arr = x1 + x2_arr = x2 + out_arr = out else: # Need flat data for BLAS, otherwise in-place does not work. # Raveling must happen in fixed order for non-contiguous out, # otherwise 'A' is applied to arrays, which makes the outcome # dependent on their respective contiguousness. - if out.data.flags.f_contiguous: + if out.flags.f_contiguous: ravel_order = 'F' else: ravel_order = 'C' - x1_arr = x1.data.ravel(order=ravel_order) - x2_arr = x2.data.ravel(order=ravel_order) - out_arr = out.data.ravel(order=ravel_order) + x1_arr = x1.ravel(order=ravel_order) + x2_arr = x2.ravel(order=ravel_order) + out_arr = out.ravel(order=ravel_order) axpy, scal, copy = scipy.linalg.blas.get_blas_funcs( ['axpy', 'scal', 'copy'], arrays=(x1_arr, x2_arr, out_arr))