ENH: Tensor-valued DiscreteLp #1238

kohr-h · 2017-11-14T22:19:01Z

TODOs:

Rebase on master and resolve conflicts

Settle on a reasonable behavior of weights when removing/adding/slicing axes. Currently we have this (I need to check that this is correct):

Weighting	add axes	remove axes	slice
`ConstWeighting`	keep	keep	keep
`PerAxisWeighting`	add 1.0 factors	subset of factors	slice arrays / keep constants
`ArrayWeighting`	fallback no weighting	fallback no weighting	slice array (1)
others	fallback no weighing	fallback no weighing	fallback no weighing

(1) only applies if no axes are added/removed.

Add byaxis_out in DiscreteLp
Add support for dtype with shape in TensorSpace.astype
Use Weighting.__getitem__ to slice, rather than an external function. ~~Having my doubts, see below.~~ Maybe I'll do this after all.
Add support for None (new axis) in various __getitem__ methods.
- TensorSpace and its elements: no restriction
- Weighting: array -> fallback to None, others should work w/o problems
- DiscreteLp and its elements: new axes only in the "output" part of the dimensions, if left of everything they turn e.g. scalar functions into (1, ..., 1) tensor-valued ones.
Add methods newaxis and dropaxis for fine-grained control w.r.t. weighting and domain min/max
Change behavior of astype with shaped dtype as discussed
Migrate to default PerAxisWeighting everywhere
Fix large-scale tests
Decide 'threshold' vs. 'edgeitems' for printing, or make a separate issue.
Rename astra_cuda_bp_scaling_factor to correction_factor or so
Clarify meaning of "simulate" in simulate_slicing
More tests, in particular for the last-minute added functionality
- All kinds of slicing of spaces and elements
- __getitem__ and __setitem__ with more complex inputs than in the doctests
- Weight propagation (done for reduce and outer in ufuncs)
- ProductSpace.astype

adler-j

Some minor comments. Much needed functionality that should go in soon if possible! Looking great to me

adler-j · 2017-11-15T08:42:30Z

odl/discr/lp_discr.py

+
+            x[indices] in space[indices]
+
+        Space indexing does not work with index arrays, boolean arrays and


Why not index arrays?

Index arrays produce flat arrays when used with Numpy arrays, and I don't know the right space for such a flat array.

Works in some cases anyway:

arr =np.ones([5, 5]) arr[[1, 2], :] Out[3]: array([[ 1., 1., 1., 1., 1.], [ 1., 1., 1., 1., 1.]])

That works here, too. What I mean by index arrays (seems to need clarification) is

>>> arr[[0, 1, 4], [3, 0, 0]] array([ 1., 1., 1.])

I don't know what the name of that one either :/

The more I think about it, that should actually (however weird it feels) be a flat TensorSpace with the correct dtype. I.e. we should exactly mirror the numpy behaviour.

Edit: This is not a requirement for this PR.

Hm, it should actually be an easy change. Decent suggestion anyway.

The only headache with this is that DiscreteLp does not necessarily know how to construct a TensorSpace. We can, of course, assume that shape and dtype should always be enough, but in the case of CuPy spaces (see #1231) we would also match the device, or for pyTorch spaces (#1229) usage of pinned versus virtual memory should be the same.

So I don't think threre's a better way to do it than to drop down to the tensors and let them handle the space wrapping.

adler-j · 2017-11-24T13:52:42Z

odl/discr/lp_discr.py

-            self.tensor[indices] = values.tensor
+        if isinstance(indices, type(self)):
+            indices = indices.tensor.data
+        res_tens = self.tensor.__getitem__(indices)


why getitem instead of actual indexing?

Well yeah 🙄

adler-j · 2017-11-24T13:53:59Z

odl/discr/lp_discr.py

+        try:
+            res_space = self.space[indices]
+        except (ValueError, TypeError):
+            # TODO: adapt error type


We should be quite careful here, i guess this is OK for now, but we need to document it (similar to numpy arrays returning flat arrays in some cases, which is confusing). We also expect stuff like

vec[complicated_indices] = 2 * vec[complicated_indices]

to work.

That code would call __setitem__ anyway and not be affected by this. But yeah, catching to too many exceptions is a bit dangerous.

Do we have any tests that verifies that the above code actually works?

I think there are unit tests with "complicated" indices for __setitem__, but I'll check again.

adler-j · 2017-11-24T13:55:00Z

odl/space/base_tensors.py

+        rn(2)
+        """
+        new_shape, _, _ = simulate_slicing(self.shape, indices)
+        return type(self)(shape=new_shape, dtype=self.dtype)


Not quite sure we want a default implementation here, this could be very implementation specific (see e.g. weighting)

Fair enough, I'll remove the default implementation.

adler-j · 2017-11-24T13:55:18Z

odl/space/base_tensors.py

@@ -496,6 +568,11 @@ def available_dtypes():
        """
        raise NotImplementedError('abstract method')

+    @property
+    def array_type(self):


Is this needed for this PR?

Not really, I used it in DiscreteLp.__getitem__, but that can be done in a different way. Removed.

adler-j · 2017-11-24T13:56:03Z

odl/space/npy_tensors.py

@@ -2159,6 +2190,73 @@ def norm(self, x):
            return float(_pnorm_diagweight(x, self.exponent, self.array))


+class PerAxisWeighting(Weighting):


Leave out partial implementation from this PR if possible.

I thought I do it here already, but no, I'll make a separate PR.

odl/test/discr/lp_discr_test.py

adler-j · 2017-11-24T13:56:49Z

odl/test/discr/lp_discr_test.py

+    assert discr.axis_labels == ()
+    assert discr.tangent_bundle == odl.ProductSpace(field=odl.RealNumbers())
+    assert discr.complex_space == odl.uniform_discr([], [], (), dtype=complex)
+    hash(discr)


Might as well check repeatability assert hash(discr) == hash(discr)

adler-j · 2017-11-24T13:58:39Z

odl/util/numerics.py

+    >>> # arr[[2, 0], [3, 3], [0, 1], [5, 2]]
+    >>> simulate_slicing(shape, ([2, 0], [3, 3], [0, 1], [5, 2]))
+    ((2,), (1, 2, 3), 0)
+    """


Show usage with np.s_ which makes stuff like np.s_[..., 2] much easier to read

Bump on this

I thought I'd done this ...

pep8speaks · 2017-12-11T18:55:51Z

Checking updated PR...

No PEP8 issues.

Comment last updated on July 02, 2018 at 23:13 Hours UTC

kohr-h · 2017-12-16T11:16:53Z

Bump

adler-j

Partial review. More to come.

adler-j · 2017-12-17T00:12:57Z

odl/discr/lp_discr.py

+        For spaces of discretized vector- or tensor-valued functions,
+        this includes the output components as ::
+
+            shape = fspace.out_shape + partition.shape


Why not make this an Examples?

odl/discr/lp_discr.py

adler-j · 2017-12-17T00:17:55Z

odl/discr/lp_discr.py

+        indices = normalized_index_expression(indices, self.shape)
+        # Avoid array comparison with `==` if `indices.contains()` is used
+        if any(idx is None for idx in indices):
+            raise ValueError('creating new axes is not supported.')


Perhaps this should be a feature for the future? Perhaps we could create something like min_pt=[None, 0, 0], max_pt=[None, 1, 1] for this case and call it a degenerate axis?

Not sure about having "invalid" entries in min_pt and max_pt. As an alternative, we could have a method newaxis that takes an index (or a bunch of indices) and min and max.

I've updated the TODOs regarding this. IMO we can allow adding axes to the output part of the dimensions, i.e. users can do

>>> space = odl.uniform_discr(0, 1, 10) >>> space[None, ...] uniform_discr(0.0, 1.0, 10, dtype=(float, (1,)))

or

>>> space = odl.uniform_discr(0, 1, 10, dtype=(float, (2, 3))) >>> space[None, ...] uniform_discr(0.0, 1.0, 10, dtype=(float, (1, 2, 3))) >>> space[:, None, None, ...] uniform_discr(0.0, 1.0, 10, dtype=(float, (2, 1, 1, 3)))

For new axes in the "in" part of the dimensions, I'd prefer an explicit methods that takes min and max in the new axis (axes), as I wrote above.

adler-j · 2017-12-17T00:18:43Z

odl/discr/lp_discr.py

+        if any(idx is None for idx in indices):
+            raise ValueError('creating new axes is not supported.')
+
+        indices_out = indices[:len(self.shape_out)]


adler-j · 2017-12-17T13:35:43Z

odl/discr/lp_discr.py

-        if values in self.space:
-            self.tensor[indices] = values.tensor
+        if isinstance(indices, type(self)):
+            indices = indices.tensor.data


Here we're truly starting to assume quite a bit about the backend. Is there no way to push this behavior down into the Tensor backend?

I'll see what I can do.

My suggestion here is to make data abstract in Tensor and therefore safe to rely on.

adler-j · 2017-12-17T13:38:18Z

odl/discr/lp_discr.py

+        try:
+            iter(indices)
+        except TypeError:
+            if indices is None:


This could be done before the try?

adler-j · 2017-12-17T13:41:49Z

odl/discr/lp_discr.py

+        try:
+            res_space = self.space[indices]
+        except (ValueError, TypeError):
+            # TODO: adapt error type


Do we have any tests that verifies that the above code actually works?

adler-j · 2017-12-17T20:10:21Z

odl/discr/lp_discr.py

-                     weighting=weighting)
+    if (weighting is None and
+            is_numeric_dtype(dtype) and
+            exponent != float('inf')):


what about -inf?

Right, that should be covered as well.

odl/space/base_tensors.py

adler-j

Some stuff left to do, but looks great!

odl/space/npy_tensors.py

adler-j · 2017-12-18T10:16:01Z

odl/space/npy_tensors.py

+
+            x[indices] in space[indices]
+
+        is ``True``.


Perhaps too "legal" of a formulation, why not

For all supported cases, indexing is implemented such that for an element ``x in space`` :: x[indices] in space[indices]

adler-j · 2017-12-18T10:22:45Z

odl/space/npy_tensors.py

-# broadcast, i.e. between scalar and full-blown in dimensionality?
+
+def slice_weighting(weighting, space_shape, indices):
+    """Return a weighting for a space after indexing.


Shouldn't this be implemented via indexing of the weighting? It seems this style is not scalable nor extendable

I had this in mind as well. That would mean putting __getitem__ and byaxis on the weighting classes, and to raise an exception where it doesn't make sense. I'd be okay with that, it would also be more extensible.

I feel it should be done that way, much more future proof

Looking at the code again regarding this one, I'm having doubts about __getitem__. The issue is that in quite a few cases, the result after slicing is not the same class as before. That would seem okay for such simple things like (array of 1 element -> scalar) transformations, but here we make much bigger changes, like ArrayWeighting -> None. This would feel wrong in a __getitem__.

adler-j · 2017-12-18T10:23:51Z

odl/space/npy_tensors.py

+          conjugate and ":math:`\odot`" for pointwise product.
+
+        - For other exponents, only norm and dist are defined. In the
+          case of exponent :math:`\\infty`, the weighted norm is defined


Do we allow -inf?

I'll add that, it should behave the same.

adler-j · 2017-12-18T10:24:27Z

odl/space/npy_tensors.py

+          .. math::
+              \| a \|_{v, \\infty} := \| a \|_{\\infty},
+
+          otherwise it is


Change the order of this docstring, e.g. default first then special case

odl/test/trafos/fourier_test.py

adler-j · 2017-12-18T11:58:42Z

odl/tomo/backends/astra_cuda.py

-            # Copy data to GPU memory
+            # Get adjoint weighting functions, using an extra correction
+            # factor accounting for inconsistencies in certain ASTRA versions
+            extra_factor = astra_cuda_bp_scaling_factor(


prefer calling it "astra_correction" or similar

odl/tomo/operators/ray_trafo.py

adler-j · 2017-12-18T12:02:33Z

odl/util/numerics.py

+def simulate_slicing(shape, indices):
+    """Simulate slicing into a Numpy array with given indices.
+
+    This function is intended to simulate indexing of a Numpy array


What does "simulate" mean here?

Something like "pretending to do it, but really only looking at the resulting shape and stuff".

adler-j · 2017-12-18T12:03:06Z

odl/util/numerics.py

+    >>> # arr[[2, 0], [3, 3], [0, 1], [5, 2]]
+    >>> simulate_slicing(shape, ([2, 0], [3, 3], [0, 1], [5, 2]))
+    ((2,), (1, 2, 3), 0)
+    """


Bump on this

adler-j · 2018-06-29T08:41:45Z

So, this would be a lovely addition, whats the status?

kohr-h · 2018-06-29T11:04:29Z

So, this would be a lovely addition, whats the status?

I guess I'm also gonna have my rebasing adventure... I'll try to bring it a bit up to speed this weekend.

kohr-h · 2018-06-30T23:57:21Z

We're back on track here again. Rebased and loads of issues fixed.

Major points are the weighting propagation and the tests yet to be written as mentioned in the TODOs.

I also kept the adjoint_weightings helper in there, since I think that in this form it's a much less intrusive addition than #1177, which still feels a bit shaky conceptually. Here we just return two functions that can used by an adjoint implementation. It's really useful to implement the correct weighting for RayTransform, so I'd like to keep it in.

kohr-h · 2018-07-02T23:27:06Z

Here's a question: If we call space.astype((float, (2,))) (i.e. with a shaped dtype), we have the following behavior currently:

rn(3).astype((float, (2,)))) -> rn((2, 3))
uniform_discr(0, 1, 3).astype((float, (2,))) -> uniform_discr(0, 1, 3, dtype=(float, (2,))
uniform_discr(0, 1, 3, dtype=(float, (4,))).astype((float, (2,))) -> uniform_discr(0, 1, 3, dtype=(float, (2,))

So for rn we add dtype.shape to the left as usual, but for uniform_discr we replace the current shape of the dtype with the new one. I'm wondering if this is reasonable or not -- the alternative would be to always add to the left.

adler-j · 2018-07-03T05:43:35Z

Excellent question. To me, astype implies replacing the type, not the shape. Hence I'd say:

rn(3).astype((float, (2,)))) should not work

But these are all good:

uniform_discr(0, 1, 3).astype((float, (2,))) -> uniform_discr(0, 1, 3, dtype=(float, (2,))
uniform_discr(0, 1, 3, dtype=(float, (4,))).astype((float, (2,))) -> uniform_discr(0, 1, 3, dtype=(float, (2,))

kohr-h · 2018-07-03T08:17:11Z

Thanks for the answer. I think that's the most reasonable suggestions, too. In Numpy, you can't do np.ones(2).astype((float, (1,))) because the result doesn't broadcast.

So I'll go with what you write above. I was already about to implement that, but realized that removing that capability would break DiscreteLp.astype with shaped dtype. So I think I'll need to add something like TensorSpace.resize to have a method for resizing. The alternative would be to create the new space in the DiscreteLp method, which wouldn't be great (would assume knowledge of the tensor impl).

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

- Fix issue with labels for new axes - Add shape and ndim properties to DiscreteLpElement - Add examples to a couple of new methods - Fix typos

kohr-h added this to Work in progress in PR Status Nov 14, 2017

kohr-h mentioned this pull request Nov 15, 2017

Representing vector- and tensor-valued functions #908

Open

kohr-h mentioned this pull request Nov 23, 2017

ENH: add cupy tensor space #1231

Open

adler-j mentioned this pull request Nov 24, 2017

Behaviour of slicing for zero and one? #1247

Closed

adler-j reviewed Nov 24, 2017

View reviewed changes

kohr-h force-pushed the issue-908__discr_tens_valued branch 2 times, most recently from 7af758f to 6e1c24e Compare December 7, 2017 17:57

kohr-h force-pushed the issue-908__discr_tens_valued branch 2 times, most recently from 53d02c4 to bb31079 Compare December 11, 2017 19:46

kohr-h changed the title ~~WIP: Tensor-valued DiscreteLp~~ ENH: Tensor-valued DiscreteLp Dec 16, 2017

kohr-h moved this from Work in progress to Needs review in PR Status Dec 16, 2017

kohr-h mentioned this pull request Dec 16, 2017

ENH: add decorator to fix adjoint weightings #1177

Open

9 tasks

adler-j reviewed Dec 17, 2017

View reviewed changes

kohr-h moved this from Needs review to In revision in PR Status Dec 18, 2017

adler-j requested changes Dec 18, 2017

View reviewed changes

kohr-h force-pushed the issue-908__discr_tens_valued branch from a3be5df to f3d8880 Compare January 31, 2018 12:37

kohr-h force-pushed the issue-908__discr_tens_valued branch 3 times, most recently from a3a03f7 to 90671c1 Compare June 30, 2018 23:52

kohr-h added this to In Progress in Release 1.0.0 Sep 11, 2018

Holger Kohr and others added 2 commits September 12, 2018 21:42

ENH: support adding new axes to spaces

a165e9c

kohr-h added 23 commits September 12, 2018 21:43

ENH: speed up bool array indexing and other indexing ops

d6b96e3

ENH: make show work with tensor-valued functions

ac8be7f

ENH: make DiscreteLp.astype work properly with shaped dtype

08359d3

BUG: fix issue with axis_labels in slicing and show

b1f51f7

ENH: complete and improve byaxis in fspace and DiscreteLp

7eb0592

MAINT: harmonize names shape_out and dtype_out

a4ce850

BUG: remove const from per-axis weighting due to issues with __eq__

8d7de2f

BUG: fix byaxis of fspace and lp_discr

3fe5515

ENH: add npy_erroroptions contex mgr and silence test warnings

e3c6267

ENH: add ProductSpaceElement.astype, closes odlgroup#1397

38c3431

MAINT: fix/improve various tests

78622e5

BUG: Fix is_int to not accept floats like 2.0

1ba4320

ENH: add array_hash utility

587f4a3

MAINT: indexing and other fixes

013be00

- Fix issue with labels for new axes - Add shape and ndim properties to DiscreteLpElement - Add examples to a couple of new methods - Fix typos

ENH: support exponent -inf and harmonize norm defs

b62d023

DOC: fix typo in lp_discr doctest

f4a6fd5

BUG: fix is_int utility

d9f0874

TST: remove weighing in expected inf-norm

4829f1d

DOC: minor doctest fix in numerics.py

6e6187d

BUG: fix string failure in is_int

723d0bf

DOC: improve doc of is_int

abfdb0d

ENH: support astype with shaped dtype in DiscreteLp

688c0d5

BUG: fix weight propagation in default case in NumpyTensorSpace

0a93baf

kohr-h force-pushed the issue-908__discr_tens_valued branch from 07443c0 to 0a93baf Compare September 12, 2018 21:39

kohr-h added 4 commits September 12, 2018 23:44

MAINT: avoid ComplexWarning in c2r FT

e094e86

API: make data an abstract method of Tensor

527bf45

ENH: use np precision instead of hardcoded one in npy tensor repr

5270eb4

ENH: add newaxis to TensorSpace and use it in DiscreteLp.astype

b9d03b0

kohr-h force-pushed the issue-908__discr_tens_valued branch from 8b99801 to b9d03b0 Compare September 13, 2018 22:32

kohr-h mentioned this pull request Feb 23, 2019

WIP: Remove element classes #1475

Draft


		x[indices] in space[indices]

		Space indexing does not work with index arrays, boolean arrays and

		@@ -2159,6 +2190,73 @@ def norm(self, x):
		return float(_pnorm_diagweight(x, self.exponent, self.array))


		class PerAxisWeighting(Weighting):

ENH: Tensor-valued DiscreteLp #1238

Are you sure you want to change the base?

ENH: Tensor-valued DiscreteLp #1238

Conversation

kohr-h commented Nov 14, 2017 • edited

adler-j left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

adler-j Nov 24, 2017 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kohr-h Dec 7, 2017 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pep8speaks commented Dec 11, 2017 • edited

Comment last updated on July 02, 2018 at 23:13 Hours UTC

kohr-h commented Dec 16, 2017

adler-j left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kohr-h Jan 29, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

adler-j left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

adler-j commented Jun 29, 2018

kohr-h commented Jun 29, 2018

kohr-h commented Jun 30, 2018

kohr-h commented Jul 2, 2018

adler-j commented Jul 3, 2018

kohr-h commented Jul 3, 2018

kohr-h commented Nov 14, 2017 •

edited

adler-j Nov 24, 2017 •

edited

kohr-h Dec 7, 2017 •

edited

pep8speaks commented Dec 11, 2017 •

edited

kohr-h Jan 29, 2018 •

edited