julia
julia copied to clipboard
Broadcasting tuple functions
After reading several recent questions about divrem
and sincos
(not in Base) broadcasting on arrays, I think the current situation could be improved. To summarize, the issue arises when trying to compute sin.(A), cos.(A)
(where A
is a vector) in an efficient way. (minmax
, extrema
, divrem
, fldmod
, fldmod1
, reim
, are similar)
The approach with tmp = sincos.(A); first.(tmp), last.(tmp)
is unsatisfactory because it requires an intermediate allocation.
The approach with sin.(A), cos.(A)
is unsatisfactory because there may be efficiency gains from computing sin
and cos
at once. I suspect some of the desire for an unzip function is caused by this problem, which has no convenient solution at the moment, aside from writing out the loop. Thinking about this problem, I have a few proposals.
Proposal A
It may make sense to allow
(sin, cos).(A)
or
broadcast((sin, cos), A)
which returns sin.(A), cos.(A)
but perhaps computed more efficiently. (In other words, an "unzipped broadcast".)
However this raises the question about whether sincos(x)
and divrem(x, y)
themselves might better be called (sin, cos)(x)
...
Proposal B
If this syntax is too revolutionary, an alternative is to offer unzip
(#13942), along with a way to make the inner broadcast lazy. That is, perhaps we could have
unzip(sincos@.(A))
which expands to
unzip(Broadcast(sincos, A))
where Broadcast
is a minimal iterable object which calls broadcast
upon collect
, but unzip
on it can be specialized to avoid the intermediate allocation.
Here I'm proposing @.
as an idiom for a lazy broadcast; this has a reasonably nice parallel to @[...]
for lazy indexing (i.e. viewing). This approach has the advantage that lazy broadcasts have themselves been a commonly requested feature.
(Now, unzip
might itself be reasonably expected to be lazy, so we might wish to require collect(unzip(sincos@.(A)))
.)
Proposal C
Ideally sin.(A), cos.(A)
would be as efficient as computing sincos.(A)
. It may be possible to intercept this pattern in lowering, and lower the construction of tuples that are broadcast over the same symbol arguments to Base._tuple_broadcast((sin, cos), A)
which by default falls back to broadcast(sin, A), broadcast(cos, A)
, but allows specialization. I don't particularly like this kind of solution, because it does break referential transparency.
Personally I am a fan of proposal B, not least because the other features it would introduce (unzip
and lazy Broadcast
) are ones that I would want independently, and it is quite nice that they work well in tandem.
Currently you can do
julia> (t::Tuple{F,G}){F,G}(x) = (t[1](x), t[2](x))
julia> (sin, cos)(0)
(0.0,1.0)
julia> (sin, cos).([-1, 0, 1]) # broadcast((sin, cos), [-1,0,1])
3-element Array{Tuple{Float64,Float64},1}:
(-0.841471,0.540302)
(0.0,1.0)
(0.841471,0.540302)
Right, but that's not the result you'd want. The point of this is to return a tuple of arrays instead of an array of tuples.
Oh, right, I missed the intention of the OP. So the example I gave above is probably a reason not to choose the Proposal A approach.
I have to confess I kind of like the idea of tuples being callable by calling each of their components on the given arguments. But I suspect others may not care for that behavior and would prefer that an operator have the effect of "tupling" functions, i.e.
funprod(f::Function, g::Function) = (args...)->(f(args...), g(args...)
Part of the reason this is appealing is that there doesn't seem to be anything else that calling a tuple could sensibly mean and the function which returns the tuple of the result of calling each element of the tuple on a given set of arguments is the correct product in the category of functions.
I agree there should be an efficient and convenient solution for this problem. Some function return multiple values.
I have come up with a solution that is faster than first.(v)
and friends. It does not avoid the allocation of an intermediary array of tuples (this is beyond my skills), but makes a single pass.
See: https://github.com/spalato/Destruct.jl
Feel free to steal.
I like Proposal A, (sin, cos).(A)
, especially because it requires no modifications to the Julia language, and in fact it could be fully implemented in an external package ("TupleCast"?).
It would be great to have a syntax which can "unzip" tuples from a single function, eg like sincos
in the example. Some calculations return tuples because the results calculated share computation.
I'm not convinced you need a "syntax" (i.e. language changes) for this.
For example, unzip(sincos).(x)
can work today with appropriate library support (unzip(f)
would return a special functor object with its own broadcast
method), and it is a lot clearer to me than unzip(sincos@.(A))
. Similarly for (sin, cos).(A)
.
It would be great if someone worked on a "TupleCast" package for this sort of thing. It's a fair amount of work only because broadcast
is probably the most complicated function in Julia.
Of course, if "tuplecasting" is sufficiently successful, it might get inducted into Base, but it's better to experiment with this sort of thing outside of Base first where it won't get set in stone too quickly.
I would argue against this for Tuples
.
Tuples
have kinda solid semantic meaning already.
They are collections of positional arguments to functions; and they are multiple return values from, functions.
They are not callable, therefore it would be weird if they were broadcast-callable.
I also think this is just kinda niche. I agree it comes up, but not so much that it is needful to have a special syntax for it.
If one wanted to broadcast, defining the scalar operation of calling things, then broadcastinmg that would do fine, no?
callf(f, xs...) = f(xs...)
callf.((sin, cos), A)
It would be kind of pleasant though, if (sin, cos)(x)
worked, which would immediately imply (sin, cos).(v)
as a corollary.
funprod(f::Function, g::Function) = (args...)->(f(args...), g(args...)
I love this. I intuitively tried the tuple version, found it didn’t work, and found this thread. Here’s a more generic implementation:
> funprod(functions...) = (args...)->map(x->x(args...), functions)
> funprod(sin,cos,tan)(0)
(0.0, 1.0, 0.0)
I like the idea of this being equivalent to (sin,cos,tan)(0)
.
funprod(f::Function, g::Function) = (args...)->(f(args...), g(args...)
I love this. I intuitively tried the tuple version, found it didn’t work, and found this thread. Here’s a more generic implementation:
> funprod(functions...) = (args...)->map(x->x(args...), functions) > funprod(sin,cos,tan)(0) (0.0, 1.0, 0.0)
I like the idea of this being equivalent to
(sin,cos,tan)(0)
.
this is broadcasting function tuples, not tuple functions
I woulda thought [sin, cos].(1)
would be like [sin(1), cos(1)]
, broadcasting the array of functions over the argument.
is this the place where pipe finally shines?
julia> 2 .|> (sin, cos)
(0.9092974268256817, -0.4161468365471424)
edit, thanks @giordano . So we want (sin,cos).(A)
to be (sin.(A), cos.(A))
not (sin(A), cos(A))
, that feels like one to many level of .
for me personally
is this the place where pipe finally shines?
I don't think so, since the point of the feature requested is to broadcast tuples/arrays in tuple functions:
julia> (2, 3) .|> (sin, cos)
(0.9092974268256817, -0.9899924966004454)
which is hardly what people want here
It's not obvious to me what (2, 3) .|> (sin, cos)
is supposed to do.
It seems obvious to me that it would do (sin(2), cos(3))
.
IMHO, as it is now, the behavior of Broadcast with several outputs is quite useless. The generated array of tuples breaks somehow the logic and always needs allocations.
Consider a simple elemental function with the signature:
a, b = f(a, b)
If I want to evaluate this function in-place in a vectorized manner for two arrays va
, vb
I would intuitively try:
va, vb .= f.(va, vb)
.
However, f.(va, vb)
yields, for some reason unknown to me, an array of tuples (AoT) which totally breaks the performance of the code and is probably useless in many cases as output. For instance, to my knowledge, code which works on an AoT can hardly be vectorized (see discussion about AoS vs. SoA) and an AoT result can not be fed again recursevly in the same way in the same function. Thus, at the moment, I am forced to write my own for loop to get this done. Writing a for loop isn't a problem in general, but it destroys the generic approach that Broadcast offers. E.g., when having a function which may mash up scalars and arrays I have to consider many special cases manually which is a lot of pain considering Broadcast could fundamentally solve this for me...
The fact that
( x .|> f ) !== ( f.(x) )
when f
is a tuple of functions is jarring to me. If one wanted a way to broadcast M
functions over N
arguments I would suggest using arrays (which currently also doesn't work, because Array
s are also not callable):
[sin cos].([1,2,3,4,5])
I would also prefer (sin, cos).(1)
over (sin, cos)(1)
for the "piping symmetry" reason. Perhaps (sin, cos)(1)
could be for the very special cases (as mentioned above) of functions that can more efficiently be computed in tandem, like sin
and cos
, while (f, g).(x)
must be used for the general case of arbitrary functions f, g
.
I'm in complete agreement with https://github.com/JuliaLang/julia/issues/22129#issuecomment-1209696536.
This fact that you can't return a tuple of arrays causes many headaches for functions that compute on tensors. It's a major weak point in the language and I hope that it gets addressed. In general, I try to avoid returning tuples because of the way broadcasting on functions returning tuples currently works. But it often results in awkward library design. The alternative is to return the array of tuples (AoT), reallocate, and pay the performance penalty.
Is there some intermediate, best-practice workaround for this?
It would be great if someone worked on a "TupleCast" package for this sort of thing.
There is StructArrays, with which @piever showed me this:
using StructArrays
function unzip_broadcast(f, args...)
bc = Broadcast.instantiate(Broadcast.broadcasted(f, args...))
StructArrays.components(StructArray(bc))
end
unzip_broadcast((x, y) -> (x, x/y), [1,2,3], [4 5])
# ([1 1; 2 2; 3 3], [0.25 0.2; 0.5 0.4; 0.75 0.6])
so broadcasting lazystack over unzip_broadcast achieves what is desired: I can generate a tuple of arrays from broadcasting over a function with multiple return values and then stack each array of the tuple. However, this is completely non-obvious and these helper packages are sufficiently sophisticated such that I don't understand the performance implications. But the code is less verbose. I also feel there is a gap in Julia where manipulating tensors should be more straightforward. Thanks again for the help!
@prittjam I had made this a while back (before 1.0). https://github.com/spalato/Destruct.jl
It worked in 1.0. It presumably still works. It's not optimal as it does work from the array of tuples. However, it is reasonably effective, preserves array shape, works with mixed types. Has no dependencies outside Base. The function has 10 lines of code.
I can't say I love the idea of making Tuple
s callable.
It seems like most of the desire for this functionality would be satisfied by a variant of map
that returns a tuple of containers rather than container of tuples for Tuple
-output functions. Eventually we'd want to have a broadcast
variant as well, but a map
-like can probably cover most use cases so is an easier start and won't require a debate over if/what syntax enables this for broadcasting.
Any updates on how best to make broadcasting on a function that returns a Tuple to give a Tuple of Arrays instead an Array of Tuples?
there are a lot of ways. one might be
A = foo.(arg) # returns an Array of Tuples
ntuple(i -> getindex.(A, i), length(first(A))) # extract as Tuple of Arrays
@bclyons12: I am not aware of any solution that avoids interim allocations (and results in separate objects, not a view).
I think that this should be implemented first in a l package though, following the techniques of collect
, or maybe more generic as in BangBang.jl Then if it proves useful it could end up in Base with dedicated syntax.
Probably the most straightforward solution is to use StructArrays. It's a very popular and stable package.
julia> A = StructArray(x=[1,2], y=[3,4])
julia> f(a) = (u=a.x+a.y, v=a.x-a.y)
julia> B = f.(A)
2-element StructArray ...
# actual arrays, not views
julia> B.u
2-element Vector{Int64}:
4
6