Hi fellow Julians,
I would like to correctly pre-allocate the space for storing the result of some custom linear operators. In particular, the eltype of the resulting collection should depend on the type of the coefficients (so to speak) of the operator, and on the eltype of the input: given all different types T <: Real and the combinations of T vs Complex{T}, I thought there might be some “best” way to do this using promotions.
Before thinking, I asked myself: how is this done in LinearAlgebra, for example to define (*) for Matrix objects? There, promote_op is used on a helper function matprod to infer the eltype of the array which will store the result:
https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/matmul.jl#L48
However, if I look for promote_op, I find a scary warning message attached to it:
https://github.com/JuliaLang/julia/blob/1b92f51ba3544edfcdb70a2f43ac1bb3a7bb5543/base/promotion.jl#L389
Can anyone clarify on this? Is it simply the case that promote_op is safe for linear operations?
Let me just share this, from Martin Holters:
The Masters office. Behind his desk, The Master is deeply concentrated, working in front of his computer.
A knock on door.
Master: Yes!
A young programming apprentice enters.
Apprentice: Excuse me, Master. I was just implementing a function which needs to allocate an array for its return values, and wanted to use promote_op to decide the element type. But now I’ve learned that might be a bad idea. Say, Master, what should I do now?
Master: I see you do read the docs. That is good. You will progress quickly. I have high hopes with you. Now to your problem: Fear not! Just use return_type directly and all will be good.
Apprentice: Thank you, Master, for sharing your wisdom.
The apprentice leaves. The Master returns to his work.
Briefly thereafter, it knocks again.
Master: Yes!
The young programming apprentice enters.
Apprentice: Excuse me, Master, I’ve followed your advice and use return_type, but often, the type will be too wide. Sometimes I even get Any. How can I get a narrower type?
Master: Ah, that my young apprentice, is simple: Just sprinkle enough @pures in your code.
Apprentice: Thank you, Master, for sharing your wisdom.
The apprentice leaves.
The Master rips off his mask and laughs a hollow, frightening laugh.
2 Likes
Thank you — now I know my question makes sense 
Any hint on where to look exactly for return_type? It doesn’t seem to appear anywhere in the documentation…
Edit: there’s in fact a bit of context in the PR you linked, thanks again.
Martin’s koan is highly sarcastic. You don’t want to be using @pure or return_type, either.
This is indeed a very hard problem, and it’s a pain to do correctly. I’d love to see this improve — I’ve spent just a bit of time brainstorming solutions but haven’t seen any silver bullets yet.
The key is that you want the element type to be entirely determined by the values. So the “easy” solution is to just perform your operation on the first element(s), get the resulting value, and then use its type to allocate the container. There are two gotchas that make this hard though:
- What do you do if the next element returns something different?
- What do you do in the empty case?
If you get a new element type that doesn’t fit, the standard protocol is to split your function such that you can re-allocate a wider element type and recurse. That is, you want an inner function that can handle picking up at an arbitrary place in your iteration through the elements so that you can “swap out” the destination array to the wider one. Something like this: base/array.jl#L625-L645.
Of course, this is all predicated on having values. What do you do about the empty case? The ideal case is if you can just throw an error — but again, that’s not always possible. So this is where we reach in and ask inference for a little help. This is where it’s appropriate to use Base.return_type.
5 Likes
Thanks for clarifying.
My case may be simpler than what you described: I know that all elements in the result will have the same type. Just like for matrix-vector multiplication, each operand has all elements of the same type, and linear combinations of them should result in something according to how types are promoted when doing scalar + and * operations. That’s why the promote_op seemed like one way to go – until I read what people have to say about it.
Related topic: Alternatives to `Base.promote_op(op, ::Type...)` for `op::Type`?
This is such a general idiom in Julia that it would be great if it could end up in the manual (in due course).
2 Likes