Remove `Element` classes!?

[kohr-h]

This is a somewhat more radical extension of #1457 to remove complexity in ODL. I'd love to have some feedback on it.

---

With a bit of distance to the weeds of the code, I've started thinking a bit about decisions we made (implicitly or explicitly) early on, and how well they served us so far. All in all I think we made good choices in many places, and the usefulness, flexibility and expressive power of ODL are indications of that.

At the same time, however, there are a couple of places where our solutions are more complex than they have to be. Structures like Discretization were simply mapped directly from mathematics, and sometimes there wasn't so much consideration put into the question "do we need this abstraction?" Rather, it was easy to do at the time and corresponded nicely to the math. I did that a lot.

Some of these rather innocent decisions have come back to hit us really hard. To me, the single most important case of this was the decision to introduce ...Element classes that wrap data containers. I can speak from the experience of implementing __array_ufunc__ for the Tensor class that it was a huge pain to write and test the code for all corner cases, just to make sure that Tensor is almost indistinguishable from NumPy arrays (but still different enough to be annoying quite regularly, see below).
Another think I've grown increasingly annoyed with is the necessity to use .asarray() all the time, just to be sure that every array functionality works as expected and we're not subject by arbitrary limitations of Tensor, either because it's ambiguous (some cases with indexing), or just not yet implemented.
And finally, all the tedious work with the CuPy backend, just to wrap it into a Tensor class and make it behave nicely, seems kind of pointless work when the array itself already works as expected.

So the big question is: Why don't we work with arrays directly?

Why don't we scrap all those element classes and use arrays directly? That would probably remove 1/4 of the whole code while touching almost no functionality. What exactly would we lose? I can't think of anything more important than these:

• x.inner(y) would have to become space.inner(x, y), same with norm and dist.
• We couldn't ask for x.space. So what? In the vast majority we have an operator that we can ask.
• Anything else??

On the other hand, I think the benefits would be massive:

• Given the right abstractions in TensorSpace and the likes, support for new array implementations right away. In the worst case, a wrapper module with name aliases would be needed, to provide a common interface (similar to the get_array_module function in CuPy)
• A huge reduction in complexity by removing all the wrapping code that needs to be maintained and probably has additional bugs in it.
• No more asarray casting, no more worrying "is this supported by the Tensor class?"
• No more decision required regarding "should this method support array-like input or not?" For instance, operator calls do, but inner and similar methods don't, which is confusing.
• One thing less to worry about for users who know NumPy already, but not ODL.

Of course, it would be quite a big change. But in my opinion, we'd benefit from it big time.

[adler-j]

Now this is a ballsy suggestion. I like it.

The main issue, by far, is of course the lack of x.space as you mention. But perhaps that is not too bad.

My suggestion here would be that you make a small trial branch for this and simply try it out. Then we'll see where and how vector.space is used.

[ozanoktem]

This sounds like a large change to the core API for ODL. Wouldn't this break a lot of code? It obviously would for parts that query the underlying space of a vector space element, right?

[adler-j]

It is a major change (to say the least) but it would probably make ODL much more accessible, especially to more engineering-y audiences.

[kohr-h]

It would be a rather large change, but my gut feeling is that it would make a lot of code simpler, while only requiring moderate adaptions to existing code.

Of course, gut feeling isn't proof. I'll have a go at it and see where the caveats lie. (In fact, I have a branch that attempts this, but at some point I lost focus and changed a lot of other things, so it isn't representative any more.)

One interesting aspect will be product spaces, in particular product space elements. Here my go-to solution would be np.array(..., dtype=object) for heterogeneous elements and an array of same type with extra dimensions for homogeneous elements.

[adler-j]

One interesting aspect will be product spaces, in particular product space elements. Here my go-to solution would be np.array(..., dtype=object) for heterogeneous elements and an array of same type with extra dimensions for homogeneous elements.

It's not obvious that this is preferable to e.g. lists. We have to consider that.

[ozanoktem]

From a conceptual point of view, it is advantageous to let the vector space element carry information about its space. This e.g. makes it easy to define a noise model that takes values on data space, we can just ask data which space it belongs to and use that when defining the noise model. If I have not misunderstood things, with the suggested changes one suddenly needs to remember the name of the data space that was once created.

[kohr-h]

It's not obvious that this is preferable to e.g. lists. We have to consider that.

True, although one advantage is support for multi-indexing, which we wouldn't have for nested lists.

[kohr-h]

From a conceptual point of view, it is advantageous to let the vector space element carry information about its space. This e.g. makes it easy to define a noise model that takes values on data space, we can just ask data which space it belongs to and use that when defining the noise model. If I have not misunderstood things, with the suggested changes one suddenly needs to remember the name of the data space that was once created.

A valid point. However, as of now, we always generate noise in a given space, not from a given element. See https://github.com/odlgroup/odl/blob/master/odl/phantom/noise.py
Only exception is Poisson noise, where you might need an extra space argument, although it would also work without.

[adler-j]

The places where the space is implicitly passed along with a vector are e.g. show, where we use the space to give context to the values. Fixing this would not be trivial (well, except for simply adding a space.show(array)) but It might be doable.

[kohr-h]

Right, show would have to go into the space.

[olivierverdier]

I agree with the proposal of @kohr-h . I have the same gut feeling that it is a good, simplifying change, and that it will probably not break that much code. x.space was nice, but its price (in term of complexity) is probably too high.

[JevgenijaAksjonova]

I think it is generally bad to break backwards compatibility. Using numpy arrays directly is nice, but isn't it possible to add such possibility without removing old classes?