Loop fusion and multiple return values

New to Julia
1

Assuming a function takes 2 arguments and returns 2 values, for example

f(x,y)=x^2,y^2

now x and y are large vectors of size N which cannot be splatted.
The ouptut of

f.(x,y)

is an array of tuples but I would like to have a tuple of size 2 (the number of return values) of 1-d arrays of size N

I found some tips using zip and splatting, but it turned out to be really slow and impractical especially for larger arrays.

For now I preallocated the output arrays and then did a manual for loop to fill in the values using

out_x=Array{Float64}(N)
out_y=Array{Float64}(N)
for i=1:N 
    out_x[i], out_y[i] = f(x[i],y[i])
end

However, I thought there must be a solution somehow using loop fusion or at least something more simple here.
I’m on Julia 0.6.3. Any advice?

edit: fixed wording and missing i-index in input parameters

2

Cf
https://github.com/JuliaLang/julia/issues/22129

3

There’s nothing wrong with the code you posted above (although I agree it’s a bit verbose). Your implementation is efficient and the code is easy to understand, so you might want to consider just writing a function to encapsulate that behavior and call it a day (after all, the best thing about Julia is that your hand-written function will be just as fast as any built-in solution).

But if you want a different approach, you could consider:

  1. Store the result tuples and then just grab the relevant element from each:
out = [f(x[i], y[i]) for i in 1:N]  # or f.(x, y)
out_x = first.(out)  # or getindex.(out, 1)
out_y = last.(out) # or getindex.(out, 2)
  1. Create a struct to hold your two results and use something like https://github.com/simonster/StructsOfArrays.jl to store them
  2. Consider whether it makes more sense to treat the pair (x, y) as a scalar in your downstream code. If you’re computing those results together, perhaps they should actually stay together? For example, rather than storing x coordinates and y coordinates, consider storing a vector of xy points. If you do that, then your original problem just goes away. This is not applicable in all cases, but it’s something I’ve found myself doing more as I’ve gotten more comfortable with Julia.
4 Likes
4

@rdeits was faster than me, but since I already wrote up the StructsOfArray solution, here is how to do it:

N = 100
xin, yin = rand(N), rand(N)
out = StructOfArrays(Tuple{Float64, Float64}, N)
out .= f.(xin, yin)
julia> out.arrays[1]
100-element Array{Float64,1}:
...

I think StructsOfArrays could use a usability update and offer something like:

soa(x::Vector, y::Vector)::StructOfArrays{Tuple{Float64, Float64}}
4 Likes
5

ZippedArrays has the unzip function with should be useful and works with arbitrary iterators.

6

Thanks! That is very interesting, the proposed unzip for Base the linked issue Base.unzip() · Issue #13942 · JuliaLang/julia · GitHub would be exactly whats needed here to keep the solution in one line like

xout, yout = unzip(f.(x,y))

unfortunately it does not seem to come in the 1.0 release

However, in the issue above there was https://github.com/spalato/Destruct.jl mentioned which does exactly the same as the yet missing unzip.

1 Like
7

Hello,

I realize this thread is older, but I would like to add a question anyway. Is there a way to do basically this if you have array slices on the left side, e.g. something along the line

function f(x)
  a=...
  b=...
  return a,b
end

myarray[:,1:2]=f.(x)

where x is a vector or iterator and a and b are vectors due to broadcasting ? This does not appear to work with Destruct.jl. If that would work one could also do fancy things like, assuming f.(x) returns four vectors,

a,myarray[:,1:2],c=f.(x)

you get the idea.

Of course this behaviour can be achieved with using an explicit loop which is only a small modification of the loop in the original question or (after googling this) fancy constructs like

function ttest(t)
         x=3.0*sin(t)
         y=2.0*cos(t)
         return x,y
       end
t=range(0.0,step=0.1,length=2)
delta=zeros(2,3)
map(x->delta[:,x].=collect(ttest.(t))[x], 1:size(collect(ttest.(t)),1))
map(x->delta[:,x].=[ttest.(t)...][x], 1:size(collect(ttest.(t)),1))

but it might be that in these cases the function is called twice as many times as needed compared to storing the result in a temporary vector, which in turn would be unneeded if the return values could be assigned to array slices directly.