odlgroup / odl

Operator Discretization Library https://odlgroup.github.io/odl/
Mozilla Public License 2.0
368 stars 105 forks source link

Indexing of multi-dimensional spaces #965

Open kohr-h opened 7 years ago

kohr-h commented 7 years ago

When indexing spaces with ndim axes, it would be natural to read the axes in left-to-right order. What I mean by that is this:

>>> space.shape == shape
True
>>> space[0].shape == shape[1:]
True
>>> space[0, 0].shape == space[0][0].shape == shape[2:]
True

and so on. Now the top one and the right part of the bottom one are as in ProductSpace, except that it doesn't implement shape, but conceptually that's how it is.

Question: Should we extend ProductSpace.__getitem__ to support tuple indexing by just doing this:

def __getitem__(self, indices):
    ...
    if isinstance(indices, tuple):
        if not indices:
            return self
        idx = indices[0]
        if isinstance(idx, Number):
            return self.spaces[idx][indices[1:]]
        elif isinstance(idx, slice):
            spaces = self.spaces[idx]
            return ProductSpace(*(space[indices[1:]]
                                  for space in spaces))
        else:
            raise TypeError('index tuple can only contain'
                            'ints or slices')
kohr-h commented 7 years ago

Of course, the question is where to stop. Probably it would be sensible to raise an error when trying to use a remaining index for a non-productspace.

adler-j commented 7 years ago

Your example makes perfect sense to me, and I agree that we want that.

With that said, we really want #157 now that we have true tensor spaces, and this could be part of that.

kohr-h commented 7 years ago

Agreed. This is now done, and I'll mirror this behavior for the indexing of TensorSpace.

kohr-h commented 7 years ago

I've commented in #861 about this, here's an update:

When indexing spaces there's always the ambiguity of indexing the "imaginary" array that the resulting space should wrap vs. indexing by axis.

The former is of limited use in TensorSpace type spaces apart from a somewhat convenient flow for indexing of tensors (index space and array the same way, then wrap the resulting array with the indexed space). In contrast, for DiscreteLp that indexing is interesting since the partition carries information on each space point that indexing will transport. Another interesting case is FunctionSpace with tensor-valued functions. Indexing should logically affect the output dimensions, but any indexing except by axis doesn't really make sense.

So anyway, I think that due to a not entirely clear-cut situation, the default should be to put indexing by axis on a byaxis property and decide later what to do in __getitem__.

For FunctionSpace, and by consequence also later for DiscreteLp, indexing can happen by input or output dimensions, so there will be byaxis_in and byaxis_out.