odl
odl copied to clipboard
odlgroup
Remove `Element` classes!?
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 becomespace.inner(x, y), same withnormanddist. - 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
TensorSpaceand 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 theget_array_modulefunction 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
asarraycasting, no more worrying "is this supported by theTensorclass?" - No more decision required regarding "should this method support array-like input or not?" For instance, operator calls do, but
innerand 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.
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.
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?
It is a major change (to say the least) but it would probably make ODL much more accessible, especially to more engineering-y audiences.
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.
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.
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.
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.
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.
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.
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.
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?
I spent some time thinking about this @JevgenijaAksjonova, but frankly I cannot come up with a solution. Of course we could try to keep the old classes, but we'd have ambiguity in what we return. E.g. what should operator(array) return? Right now it returns an element, but in this case it should probably return an array.
The operator(array) returning an element in its domain is very natural, but having it return a NumPy-array is not really good. Maybe I am misunderstanding something.
@JevgenijaAksjonova I totally agree with your general sentiment. There needs to be a good justification for a big change like this. I think once there is a working implementation of the proposal, we can see and judge, and then decide if we like it or not.
In general, I see this as a long-term investment into making ODL more maintainable and easy to take up for newcomers. We would move out of business that we should never have dealt with. It's part of the "let data just be data" philosophy from functional programming, as opposed to the "make data active by adding behavior" philosophy of object-oriented programming.
@ozanoktem It's up to us to decide what "element in its domain" means. We can (and would) consider arrays of correct shape and data type to be elements of a certain tensor space. There's nothing unnatural about that, and in that scenario, operators returning arrays would be perfectly okay.
@kohr-h , one of the very early design goals in ODL was to make sure we kept some information on what the array represented. Just because two arrays have the same size and same datatype does not mean they necessarily represent the same vector space element. The main reason was to ensure one does not make logical errors just because they are permitted from a pure array-algebraic point of view. I suggest we do not abandon that design principle unless there is very strong reasons to do so,
What I'm concerned about is that not having an abstract Element class makes extending ODL harder when it comes to other discretisations which do not allow the representation as a single array (like regions of interest in CT, where you have certain areas with a finer grid or Galerkin methods). This would lead to several potentially incompatible interfaces and would make it harder to write general solvers which work together with all of these interfaces. And here, I have not even mentioned things which happen in non-linear spaces like Lie groups or manifolds.
And it would also mean that any change of the implementation (like the one of ProductSpace elements) would potentially introduce backwards incompatibility.
For these reasons, I think that this change should also be seen together with a longer-term perspective of what ODL should become (Operator Discretisation Library might not only include the type of discretisations we have now).
Just because two arrays have the same size and same datatype does not mean they necessarily represent the same vector space element.
Why does an object that represents data have to represent more than that? Why does it have to know about spaces at all? In my view, it's sufficient that spaces know about what types of objects they can work with. The opposite direction is not necessary.
Your position reflects to a high degree how things are handled in statically typed languages. You say once and upfront what type your data is, and that's that. If you want to interpret it differently, you're forced to cast explicitly. Contrast that to the dynamic view of languages like Python. It is determined at runtime what the type of your data is, and then decided if an operation is allowed or not. That's the principle that I'm applying here as well.
If you look at how things evolved over time, we started with a very rigid type system where the boundaries between different spaces were very clearly drawn. The big downside of that is usability: if you get errors thrown in your face just because of tiny differences, you find yourself casting everywhere, which is not only potentially inefficient, but also defeats the purpose of the type system in the first place. So we started tearing down the boundaries, and now you can give an array to an operator directly, for instance (in other places it's not possible, we're not very consistent there). That was a big usability plus.
Long story short, I think the principle of "this element belongs to this space and no other space" is impractical, and we should abandon it.
The main reason was to ensure one does not make logical errors just because they are permitted from a pure array-algebraic point of view.
In my experience, the vast majority of logical errors can be avoided based on the shape of the data alone, and the next (much smaller) portion of errors by data type. Any further "help" from a type system is mainly a source of frustration because it forces manual casting.
I suggest we do not abandon that design principle unless there is very strong reasons to do so,
I have one: Wrapping data in a custom class and trying hard to make it look almost like the wrapped thing is a huge maintenance burden, and quite error prone. I can tell, since I implemented a big portion of it.
@sbanert I'm not suggesting to completely disallow custom element types. If there is a case where existing types don't do the job, then sure, make a custom one.
What I mean is that none of the currently implemented scenarios require such a custom type.
Regarding backwards compatibility, sure, code that does x.inner(y) or x.space or x.asarray() will have to be ported. But those are easy changes, and we can write a migration guide if necessary. I still consider the benefit to be worth the effort.
Why does an object that represents data have to represent more than that? Why does it have to know about spaces at all? In my view, it's sufficient that spaces know about what types of objects they can work with. The opposite direction is not necessary.
Because certain mathematical operations do not make sense if the "vector" does not represent the same type of element. As an example, addition two 10^6-vectors where one is the Taylor series coefficients and the other is some Wavelet coefficients doesn't make sense.
Your position reflects to a high degree how things are handled in statically typed languages. You say once and upfront what type your data is, and that's that. If you want to interpret it differently, you're forced to cast explicitly. Contrast that to the dynamic view of languages like Python. It is determined at runtime what the type of your data is, and then decided if an operation is allowed or not. That's the principle that I'm applying here as well.
My position reflects the mathematical view on things, i.e., to use "typing" to ensure that mathematical operations make sense. A "dynamic view" without any data types does not care about mathematical consistency, or how would it decide whether an operation is allowed or not?
Now, mathematicians also frequently make use of the fact that a vector space element belongs to multiple spaces. As an example, one may view a continuous function with compact support as an element in L^2 at one stage, then at some other stage view it as an element in C. Certain arguments may use the L^2 topology, other arguments may be based on other topologies related to C. Obviously we cannot bookkeep such stuff, but also forget about what an array represents is in my opinion not a good choice.
My final remark concerns the community we primarily target with ODL, which in my opinion is the math community and not the engineering community. In such case, I am not convinced that the suggested changes put forward by @kohr-h are ideal.
Because certain mathematical operations do not make sense if the "vector" does not represent the same type of element. As an example, addition two 10^6-vectors where one is the Taylor series coefficients and the other is some Wavelet coefficients doesn't make sense.
That's a too broad statement IMO. If the two vectors have the same size and the entry type allows addition, then addition makes perfect sense. You may not want to add them because the result is not useful, but we're not here to decide for users what they ought to do. We can and should only catch obvious errors.
Apart from that, there is no element in ODL that would even try to catch that type of usage error. If both operators map to rn(1e6), addition of their results is (and should be) allowed.
My position reflects the mathematical view on things, i.e., to use "typing" to ensure that mathematical operations make sense. A "dynamic view" without any data types does not care about mathematical consistency, or how would it decide whether an operation is allowed or not?
I think you're conflating "formally correct" with "meaningful". Just saying "makes sense" or "makes no sense" can be understood in both ways, and I think we need to disentangle the two. ODL never made an attempt to decide for users what operations are meaningful. We only ever checked "formal correctness", and I think it's the only thing we should be doing. Of course, one can start to push things from the "not meaningful" into the "formally incorrect" category by introducing element types. But in my view, the past has shown that this is a poor way of dealing with the issue, since it introduces a usability and maintainability nightmare.
My final remark concerns the community we primarily target with ODL, which in my opinion is the math community and not the engineering community. In such case, I am not convinced that the suggested changes put forward by @kohr-h are ideal.
I'll give an example so we're just imagining what the consequences would be. Consider the CG method. Here's the current version (without distracting details):
def conjugate_gradient(op, x, rhs, niter):
if op.domain != op.range:
raise ValueError('operator needs to be self-adjoint')
if x not in op.domain:
raise TypeError('`x` {!r} is not in the domain of `op` {!r}'
''.format(x, op.domain))
r = op(x)
r.lincomb(1, rhs, -1, r) # r = rhs - A x
p = r.copy()
d = op.domain.element() # Extra storage for storing A x
sqnorm_r_old = r.norm() ** 2 # Only recalculate norm after update
if sqnorm_r_old == 0: # Return if no step forward
return
for _ in range(niter):
op(p, out=d) # d = A p
inner_p_d = p.inner(d)
if inner_p_d == 0.0: # Return if step is 0
return
alpha = sqnorm_r_old / inner_p_d
x.lincomb(1, x, alpha, p) # x = x + alpha*p
r.lincomb(1, r, -alpha, d) # r = r - alpha*d
sqnorm_r_new = r.norm() ** 2
beta = sqnorm_r_new / sqnorm_r_old
sqnorm_r_old = sqnorm_r_new
p.lincomb(1, r, beta, p) # p = s + b * p
This is what it would look like afterwards:
def conjugate_gradient(op, x, rhs, niter):
if op.domain != op.range:
raise ValueError('operator needs to be self-adjoint')
if x not in op.domain:
raise TypeError('`x` {!r} is not in the domain of `op` {!r}'
''.format(x, op.domain))
X = op.domain
r = op(x)
X.lincomb(1, rhs, -1, r, out=r) # was r.lincomb(1, rhs, -1, r)
p = r.copy()
d = op.domain.element()
sqnorm_r_old = X.norm(r) ** 2 # was r.norm() ** 2
if sqnorm_r_old == 0:
return
for _ in range(niter):
op(p, out=d)
inner_p_d = X.inner(p, d) # was p.inner(d)
if inner_p_d == 0.0:
return
alpha = sqnorm_r_old / inner_p_d
X.lincomb(1, x, alpha, p, out=x) # was x.lincomb(1, x, alpha, p)
X.lincomb(1, r, -alpha, d, out=r) # was r.lincomb(1, r, -alpha, d)
sqnorm_r_new = X.norm(r) ** 2 # was r.norm() ** 2
beta = sqnorm_r_new / sqnorm_r_old
sqnorm_r_old = sqnorm_r_new
X.lincomb(1, r, beta, p, out=p) # was p.lincomb(1, r, beta, p)
I think this code makes it actually more clear in what spaces we're operating. Here, domain and range are identical, but when they're not, the code much more clearly shows which operation happens in which space.
I spent some time thinking about this @JevgenijaAksjonova, but frankly I cannot come up with a solution. Of course we could try to keep the old classes, but we'd have ambiguity in what we return. E.g. what should
operator(array)return? Right now it returns an element, but in this case it should probably return an array.
What about returning an array if the array is given as input and space element in the opposite case?
What about returning an array if the array is given as input and space element in the opposite case?
The problem is that this would break current code, where we always return a space element, regardless of if the input is an array or space element.
What about returning an array if the array is given as input and space element in the opposite case?
The problem is that this would break current code, where we always return a space element, regardless of if the input is an array or space element.
Sure, there would be some damage, but smaller. So it still could be a good compromise. The biggest problem I see, is that it will complicate development in the future, as one will have to think about these 2 options. =/
Sure, there would be some damage, but smaller. So it still could be a good compromise. The biggest problem I see, is that it will complicate development in the future, as one will have to think about these 2 options. =/
To me, the main issue with this would be to guarantee that op(x) in op.range is always True. If we don't guarantee this, we'll be causing a lot of confusion.
But to guarantee this would mean to extend x in space to arrays x, and then we'd basically be very close to my proposal.
I am torn in this discussion: from one perspective it is good to reduce the amount of code that needs to be maintained. On the other hand, I tend to also agree with @ozanoktem @sbanert: I like the idea that the elements "know what they are". For example if one wants to implement a discretization.
You can kind of read it between the lines in some places, but I would like to know what the principle in the original design was? Why was all this mathematical abstraction implemented? What purpose was it actually supposed to fill?
From the write-ups by @kohr-h, it seems to currently fill very little purpose. Is that because we are not utilizing all the potential that it gives? Or simply because it does not make much of a difference if you have it or not?
In response to @aringh, we initially envisioned the vector spaces and elements to be quite lightweight, e.g. that they should basically satisfy the vector space axioms (scalar multiplication, addition, possibly an inner product).
However we quickly found that a significant part of implementations go beyond this. Users want to be able to work with the elements as arrays in several regards, e.g. pointwise multiplications, taking e.g. the pointwise exponential, slicing, etc. In the end we've spent A LOT of time to add these features, but meanwhile the objects have become more and more like their underlying implementations.
Frankly I do not doubt that the elements in themselves carry a lot of value, but the implementation is very cumbersome as is given that we basically need to write passthroughs to every single operation that our users want to do.
Perhaps we should look into subclassing instead, e.g. https://docs.scipy.org/doc/numpy-1.13.0/user/basics.subclassing.html
Frankly, the more I think about it the more I wonder. Why did we not go with sub-classing from the start? I seem to remember that we did some initial forages into this but stopped, although I cannot remember why.
In response to @aringh, there were two aims with ODL: (a) Enable users to write code that directly corresponds to the mathematical specification of an algorithm, and (b) offer means for users to share software components and/or implementations. A key aspect of code sharing was to stress compartmentalisation with associated documentation and unit tests, another was to separate out the discretisation from the continuum representation.
Another aspect was to offer consistency checks to prevent erroneous mathematical operations even though they are permitted from a computer science point of view. This is the source of the discussions right now concerning what information a vector space element is suppose to carry. The few times I have used ODL, it has been quite helpful to query an operator regarding its range to ensure it is properly set up. If this ability is lost, i.e., if the range is just some NymPy array, then I think we have lost important functionality.