ENH: implement tensor-valued DiscreteLp

Changes in detail:
- Add dtype with shape to DiscreteLp (mostly __repr__, factory
  functions and some downstream methods). As a consequence,
  `shape_[in,out]` and `ndim_[in,out]` are added for the
  different types of axes, as well as `scalar_dtype`.
- Add `PerAxisWeighting` and make it the default for
  `DiscreteLp` type spaces. Reason: this way the `tspace`
  knows how to deal with removed axes etc. This is important
  for a smooth experience with indexing and reductions over
  axes.
  Helpers for slicing weightings help structure this task.
- Implement `__getitem__` for `TensorSpace` and `DiscreteLp`,
  including (hopefully) reasonable propagation of weights.
  The new `simulate_slicing` utility function simplifies
  this task.
- Allow indexing with ODL tensors of boolean or integer dtype.
- Implement correct weighting for backprojections with
  non-uniform angles, using per-axis weighting and a new
  helper `adjoint_weightings` to apply the weightings in an
  efficient way.
  The correct weighting from the range of `RayTransform`
  is determined by the new `proj_space_weighting` helper.
- Change the space `_*_impl` methods to always expect and
  return Numpy arrays, and adapt the calling code.
- Change behavior of `norm` and `dist` to ignoring weights
  for `exponent=inf`, in accordance with math.
- Improve speed of `all_equal` for comparison of arrays.
- Account for `None` entries in indices in the
  `normalized_index_expression` helper, thus allowing
  creation of new axes.
- Remove `dicsr_sequence_space`, it was largely unused and
  just a maintenance burden. Use a regular `uniform-discr`
  from zero to `shape` instead.
- Remove `Weighting.equiv()` mehtods, never used and hard
  to maintain (n^2 possibilities).
- Remove the (largely useless) `_weighting` helper to create
  weighting instances since it would have been ambiguous
  with sequences of scalars (array or per axis?). Also remove
  the `npy_weighted_*` functions, they were useless, too.
- Remove some dead code from tomo/util.
- A bunch of minor fixes, as usual.

Closes: odlgroup#908
Closes: odlgroup#907
Closes: odlgroup#1113
Closes: odlgroup#965
Closes: odlgroup#286
Closes: odlgroup#267
Closes: odlgroup#1001

examples/tomo/ray_trafo_parallel_2d.py

@@ -28,7 +28,7 @@
 # projection data (or any element in the projection space).
 backproj = ray_trafo.adjoint(proj_data)
 
-# Shows a slice of the phantom, projections, and reconstruction
+# Show a slice of the phantom, projections, and reconstruction
 phantom.show(title='Phantom')
 proj_data.show(title='Projection data (sinogram)')
 backproj.show(title='Back-projected data', force_show=True)

odl/discr/discr_mappings.py

@@ -79,10 +79,11 @@ def __init__(self, map_type, fspace, partition, tspace, linear=False):
                              'of the function set {}'
                              ''.format(partition, fspace.domain, fspace))
 
-        if tspace.shape != partition.shape:
-            raise ValueError('`tspace.shape` not equal to `partition.shape`: '
-                             '{} != {}'
-                             ''.format(tspace.shape, partition.shape))
+        if tspace.shape != fspace.out_shape + partition.shape:
+            raise ValueError(
+                '`tspace.shape` not equal to '
+                '`fspace.out_shape +  partition.shape`: {} != {}'
+                ''.format(tspace.shape, partition.shape))
 
         domain = fspace if map_type == 'sampling' else tspace
         range = tspace if map_type == 'sampling' else fspace

odl/discr/discretization.py

@@ -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`
@@ -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
         --------
@@ -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