Skip to content

Commit

Permalink
WIP: implement discretized tensor-valued function spaces
Browse files Browse the repository at this point in the history
  • Loading branch information
Holger Kohr committed Dec 1, 2017
1 parent 348095a commit aec99e2
Show file tree
Hide file tree
Showing 6 changed files with 306 additions and 96 deletions.
13 changes: 10 additions & 3 deletions odl/discr/discretization.py
Original file line number Diff line number Diff line change
Expand Up @@ -375,6 +375,13 @@ def __eq__(self, other):
def __getitem__(self, indices):
"""Return ``self[indices]``.
.. note::
This is a basic implementation that does little more than
propagating the indexing to the `tensor` attribute. In
particular, the result will **not** be wrapped as a
`DiscretizedSpaceElement` again. Subclasses should take
care of that by overriding ``__getitem__``.
Parameters
----------
indices : int or `slice`
Expand Down Expand Up @@ -417,9 +424,9 @@ def sampling(self, ufunc, **kwargs):
Parameters
----------
ufunc : ``self.space.fspace`` element
The continuous function that should be samplingicted.
The continuous function that should be sampled.
kwargs :
Additional arugments for the sampling operator implementation
Additional arguments for the sampling operator implementation.
Examples
--------
Expand All @@ -445,7 +452,7 @@ def interpolation(self):
Returns
-------
interpolation_op : `FunctionSpaceMapping`
Operatior representing a continuous interpolation of this
Operator representing a continuous interpolation of this
element.
Examples
Expand Down
112 changes: 92 additions & 20 deletions odl/discr/lp_discr.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,8 @@
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)
normalized_scalar_param_list, safe_int_conv, normalized_nodes_on_bdry,
normalized_index_expression, simulate_slicing)

__all__ = ('DiscreteLp', 'DiscreteLpElement',
'uniform_discr_frompartition', 'uniform_discr_fromspace',
Expand Down Expand Up @@ -100,8 +101,8 @@ def __init__(self, fspace, partition, tspace, interp='nearest', **kwargs):
else:
# Got sequence of strings
if len(interp) != partition.ndim:
raise ValueError('expected {} (ndim) entries in `interp`, '
'got {}'.format(partition.ndim, len(interp)))
raise ValueError('expected {} entries in `interp`, got {}'
''.format(partition.ndim, len(interp)))

self.__interp_byaxis = tuple(str(s).lower() for s in interp)
if any(s not in _SUPPORTED_INTERP for s in self.interp_byaxis):
Expand Down Expand Up @@ -491,33 +492,104 @@ def __getitem__(self, indices):
A single integer slices along the first axis (index does not
matter as long as it lies within the bounds):
>>> rn = odl.rn((3, 4, 5, 6))
>>> rn[0]
rn((4, 5, 6))
>>> rn[2]
rn((4, 5, 6))
>>> space = odl.uniform_discr([0, 0, 0, 0], [1, 2, 3, 4],
... shape=(2, 4, 6, 8))
>>> space[0]
uniform_discr([ 0., 0., 0.], [ 2., 3., 4.], (4, 6, 8))
>>> space[1]
uniform_discr([ 0., 0., 0.], [ 2., 3., 4.], (4, 6, 8))
Multiple indices slice into corresponding axes from the left:
>>> rn[0, 0]
rn((5, 6))
>>> rn[0, 1, 1]
rn(6)
>>> space[0, 0]
uniform_discr([ 0., 0.], [ 3., 4.], (6, 8))
>>> space[0, 1, 1]
uniform_discr(0.0, 4.0, 8)
Ellipsis (``...``) and ``slice(None)`` (``:``) can be used to keep
one or several axes intact:
>>> rn[0, :, 1, :]
rn((4, 6))
>>> rn[..., 0, 0]
rn((3, 4))
>>> space[0, :, 1, :]
uniform_discr([ 0., 0.], [ 2., 4.], (4, 8))
>>> space[..., 0, 0]
uniform_discr([ 0., 0.], [ 1., 2.], (2, 4))
With slices, parts of axes can be selected:
>>> rn[0, :3, 1:4, ::2]
rn((3, 3, 3))
>>> space[0, :3, 1:4, :]
uniform_discr([0. , 0.5, 0. ], [ 1.5, 2.5, 4. ], (3, 3, 8))
Indexing is also supported for vector- or tensor-valued function
spaces. In that case, leftmost indices are associated with the
output dimensions, rightmost indices with the input (domain)
dimensions:
>>> space = odl.uniform_discr(0, 1, 5, dtype=(float, (2, 3)))
>>> space.shape
(2, 3, 5)
>>> space[0] # (2, 3)-valued -> (3,)-valued (first row)
uniform_discr(0.0, 1.0, 5, dtype='('float64', (3,))')
>>> space[:, 0] # (2, 3)-valued -> (2,)-valued (first column)
uniform_discr(0.0, 1.0, 5, dtype='('float64', (2,))')
>>> space[..., :2] # last axis = domain
uniform_discr(0.0, 0.4, 2, dtype='('float64', (2, 3))')
"""
pass
# Unwrap tensor for passing on to underlying indexing methods
if isinstance(indices, self.element_type):
indices = indices.tensor

if (isinstance(indices, (self.tspace.element_type,
self.tspace.array_type)) and
indices.dtype == bool):
# Raising here as signal for `DiscreteLpElement.__getitem__`
# This is an optimization that avoids counting the number of
# `True` entries here *and* in the actual element indexing
# operation. We do it lazily instead.
raise TypeError('cannot index space with a boolean element or '
'array')

# The overall logic of the slicing is as follows:
# - indices are normalized to be a tuple of length `ndim`
# - the `tspace` is indexed with the whole index tuple
# - the `partition` is indexed with the "in" part of the tuple
# only, i.e., the last part excluding `n_out = len(shape_out)`
# axes
# - the `out_shape` for the function space is determined using the
# first `n_out` index tuple entries
# - the new function space is constructed from `new_partition.set`,
# `fspace.scalar_out_dtype` and the new `out_shape`

# Normalize the index expression based on the full shape, and
# split into "in" and "out" parts
indices = normalized_index_expression(indices, self.shape)
indices_out = indices[:len(self.shape_out)]
indices_in = indices[len(self.shape_out):]

# Index tensor space
res_tspace = self.tspace[indices]

# Index partition, manually removing collapsed axes
res_part = self.partition[indices_in]
_, collapsed_axes_in, _ = simulate_slicing(self.shape_in, indices_in)
remaining_axes_in = [i for i in range(len(self.shape_in))
if i not in collapsed_axes_in]
res_part = res_part.byaxis[remaining_axes_in]

# Determine new fspace
sliced_shape_out, _, _ = simulate_slicing(self.shape_out, indices_out)
res_fspace = FunctionSpace(
res_part.set,
out_dtype=(self.fspace.scalar_out_dtype, sliced_shape_out))

# Further attributes for the new space
res_interp = [self.interp_byaxis[i] for i in range(len(self.shape_in))
if i in remaining_axes_in]
res_labels = [self.axis_labels[i] for i in range(len(self.shape_in))
if i in remaining_axes_in]

# Create new space
return DiscreteLp(res_fspace, res_part, res_tspace, interp=res_interp,
axis_labels=res_labels)

@property
def byaxis_in(self):
Expand Down Expand Up @@ -668,7 +740,7 @@ def __repr__(self):
sep=[',\n', ', ', ',\n'],
mod=['!r', '!s'])

return '{}({})'.format(ctor, indent(inner_str))
return '{}(\n{}\n)'.format(ctor, indent(inner_str))

def __str__(self):
"""Return ``str(self)``."""
Expand Down
78 changes: 19 additions & 59 deletions odl/space/base_tensors.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,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,
normalized_index_expression,
simulate_slicing,
array_str, dtype_str, signature_string, indent, writable_array)
from odl.util.ufuncs import TensorSpaceUfuncs
from odl.util.utility import TYPE_MAP_R2C, TYPE_MAP_C2R
Expand Down Expand Up @@ -358,15 +358,21 @@ def __contains__(self, other):
def __getitem__(self, indices):
"""Return ``self[indices]``.
For most cases, indexing is implemented such that for an element
``x in space``, the statement ::
For all supported cases, indexing is implemented such that for an
element ``x in space``, the statement ::
x[indices] in space[indices]
is ``True``.
Space indexing does not work with index arrays, boolean arrays and
"fancy indexing".
Space indexing works with
- integers,
- `slice` objects,
- index arrays, i.e., 1D array-like objects containing integers,
and combinations of the above. It does not work with boolean arrays
or more advanced "fancy indexing".
.. note::
This method is a default implementation that propagates only
Expand Down Expand Up @@ -404,66 +410,20 @@ def __getitem__(self, indices):
>>> rn[0, :3, 1:4, ::2]
rn((3, 3, 3))
Lists with integers will repeat parts in an axis (this is not very
useful for spaces, more so for elements):
Array-like objects (must all have the same 1D shape) of integers are
treated as follows: if their common length is ``n``, a new axis
of length ``n`` is created at the position of the leftmost index
array, and all index array axes are collapsed (Note: this is
not so useful for spaces, more so for elements):
>>> rn[0, 0, [0, 1, 0, 2], :]
rn((4, 6))
An exception is indexing with lists of integers of equal length
in each axis. In this case a flat space with size equal to the
list lengths is returned (again, this is more useful for arrays):
>>> rn[:, [1, 1], [0, 2], :]
rn((3, 2, 6))
>>> rn[[2, 0], [3, 3], [0, 1], [5, 2]]
rn(2)
"""
if isinstance(indices, tuple) and len(indices) == 1:
indices = indices[0]

if (isinstance(indices, (Tensor, self.array_type)) and
indices.dtype == bool):
# Raising here as signal for `Tensor.__getitem__`
raise TypeError('cannot index space with an element or array')

indices = normalized_index_expression(indices, self.shape)

# List-of-lists-of-integers indexing handled separately here
if all(isinstance(idx, list) for idx in indices):
# TODO: not good to hard-code Numpy array here.
# Better solution? Use len() only and let fail if in doubt?
indices = [np.asarray(idx, dtype='int64') for idx in indices]
if not all(idx.ndim == 1 and idx.size == indices[0].size
for idx in indices):
# Not entirely correct error type, but it's supposed to
# signal invalid indexing to `Tensor.__getitem__`
raise TypeError(
'list-of-lists indexing requires the sublists to '
'be one-dimensional and of the same size')

new_shape = [indices[0].size]

else:
new_shape = []
for i in range(len(indices)):
if isinstance(indices[i], Integral):
# Skip axis in shape, it's being squeezed
pass
elif isinstance(indices[i], slice):
start, stop, step = indices[i].indices(self.shape[i])
new_shape.append(int(np.ceil((stop - start) / step)))
elif isinstance(indices[i], list):
if not all(isinstance(j, Integral) for j in indices[i]):
raise TypeError(
'index list {!r} contains invalid entries'
''.format(indices[i]))
new_shape.append(len(indices[i]))
else:
raise TypeError('got invalid element {!r} in `indices`'
''.format(indices[i]))

# Add implicit `[:]` for "too short" indices
new_shape.extend(self.shape[len(indices):])

new_shape, _, _ = simulate_slicing(self.shape, indices)
return type(self)(shape=new_shape, dtype=self.dtype)

def __eq__(self, other):
Expand Down
6 changes: 4 additions & 2 deletions odl/space/npy_tensors.py
Original file line number Diff line number Diff line change
Expand Up @@ -763,6 +763,8 @@ def __hash__(self):
return hash((super(NumpyTensorSpace, self).__hash__(),
self.weighting))

# TODO: add __getitem__ with weight propagation when #1001 is fixed

@property
def byaxis(self):
"""Return the subspace defined along one or several dimensions.
Expand All @@ -784,7 +786,7 @@ def byaxis(self):
"""
space = self

class NpyTensorSpacebyaxis(object):
class NpyTensorSpaceByAxis(object):

"""Helper class for indexing by axis."""

Expand All @@ -810,7 +812,7 @@ def __repr__(self):
"""Return ``repr(self)``."""
return repr(space) + '.byaxis'

return NpyTensorSpacebyaxis()
return NpyTensorSpaceByAxis()

def __repr__(self):
"""Return ``repr(self)``."""
Expand Down
28 changes: 18 additions & 10 deletions odl/util/normalize.py
Original file line number Diff line number Diff line change
Expand Up @@ -224,15 +224,22 @@ def normalized_index_expression(indices, shape, int_to_slice=False):
if len(indices) < ndim and Ellipsis not in indices:
indices.append(Ellipsis)

# Turn Ellipsis into the correct number of slice(None)
if Ellipsis in indices:
if indices.count(Ellipsis) > 1:
raise ValueError('cannot use more than one Ellipsis.')
# Turn Ellipsis into the correct number of slice(None). We need a slightly
# more elaborate loop because we need to avoid comparing arrays to
# something using `==`, i.e., we cannot use `indices.count()` or
# `... in indices`.
ell_idx = None
for i, idx in enumerate(indices):
if idx is Ellipsis:
if ell_idx is not None:
raise ValueError('cannot use more than one Ellipsis.')
else:
ell_idx = i

eidx = indices.index(Ellipsis)
if ell_idx is not None:
extra_dims = ndim - len(indices) + 1
indices = (indices[:eidx] + [slice(None)] * extra_dims +
indices[eidx + 1:])
indices = (indices[:ell_idx] + [slice(None)] * extra_dims +
indices[ell_idx + 1:])

# Turn single indices into length-1 slices if desired
for (i, idx), n in zip(enumerate(indices), shape):
Expand All @@ -241,7 +248,7 @@ def normalized_index_expression(indices, shape, int_to_slice=False):
idx += n

if idx >= n:
raise IndexError('Index {} is out of bounds for axis '
raise IndexError('index {} is out of bounds for axis '
'{} with size {}.'
''.format(idx, i, n))
if int_to_slice:
Expand All @@ -251,8 +258,9 @@ def normalized_index_expression(indices, shape, int_to_slice=False):
if any(s.start == s.stop and s.start is not None or
s.start == n
for s, n in zip(indices, shape) if isinstance(s, slice)):
raise ValueError('Slices with empty axes not allowed.')
if None in indices:
raise ValueError('slices with empty axes not allowed.')
# Avoid array comparison with `==` also here
if any(idx is None for idx in indices):
raise ValueError('creating new axes is not supported.')
if len(indices) > ndim:
raise IndexError('too may indices: {} > {}.'
Expand Down
Loading

0 comments on commit aec99e2

Please sign in to comment.