5

Let's say I have a function named trig, which returns two outputs:

function trig(x)
    return(sin(x), cos(x))
end

If I want to evaluate trig over many values, I can use the map function:

julia> out = map(trig, (0:(pi/12):(pi/2)))

out is a 7 element array and in each element, there is a tuple containing two elements:

julia> out
7-element Array{Tuple{Float64,Float64},1}:
(0.0,1.0)
(0.258819,0.965926)
(0.5,0.866025)
(0.707107,0.707107)
(0.866025,0.5)
(0.965926,0.258819)
(1.0,6.12323e-17)

My question is: What is the best way to disentangle my sines and cosines so that I have two arrays with 7 elements each? Is it possible to broadcast trig without creating the superfluous array of tuples and instead directly create the two arrays that I am actually interested in?

For the time being, I am invoking map again in order to extract values from out in order to populate the arrays I want, but I don't think it is the best way to do this:

sines = map(x->x[1], out)
cosines = map(x->x[2], out)

For the purpose of this question, assume trig is a computationally expensive function. So, please do not give me an answer that requires trig to be evaluated more than once.

  • I don't think there is any way to do this sensibly without re-allocation, since, e.g. `0.0` and `0.258819` are not next to each other in memory (the `1.0` is in-between). To see this, use `reinterpret(Float64, out, (2, 7))`, which has no overhead, but note that `sines` is the first row of the output matrix, and `cosines` is the second row, and since arrays are column-major order, extracting `sines` as a `Vector` *must* trigger re-allocation. In other words, I suspect `map` will be as good a way as any of solving this... – Colin T Bowers May 29 '17 at 04:31
  • 1
    Alternatively, just provide your own method definition for `trig` for the `AbstractVector` input case and allocate the output of calls to the `trig` method in the question to two vectors explicitly. There's probably a clever one-liner to do this but I can't see it off the top of my head (my version would have a pre-allocation line and an explicit loop) – Colin T Bowers May 29 '17 at 04:43
  • 1
    See also https://stackoverflow.com/q/44010033 and https://stackoverflow.com/q/41183481/6172490. – tim May 29 '17 at 08:00
  • 1
    Note that in the specific example you can just do `x, y = trig(0:(pi/12):(pi/2))`, because there are array versions of `sin` and `cos` – Michael K. Borregaard May 29 '17 at 10:19

1 Answers1

1

Thank you tim for sharing your answer to an earlier question that I must have overlooked in my search before asking this one. Before today, I had never heard of the getIndex function, yet it seems that getindex is the function I want, provided that I vectorize it by putting a dot in front:

julia> @time sine_map = map(x->x[1], out)
  0.051494 seconds (13.32 k allocations: 602.941 KB)
7-element Array{Float64,1}:
 0.0
 0.258819
 0.5
 0.707107
 0.866025
 0.965926
 1.0

julia> @time sine_geti = getindex.(out, 1)
  0.029560 seconds (9.24 k allocations: 416.910 KB)
7-element Array{Float64,1}:
 0.0
 0.258819
 0.5
 0.707107
 0.866025
 0.965926
 1.0

julia> @time cosine_map = map(x->x[2], out)
  0.037328 seconds (13.32 k allocations: 602.941 KB)
7-element Array{Float64,1}:
 1.0
 0.965926
 0.866025
 0.707107
 0.5
 0.258819
 6.12323e-17

julia> @time cosey_geti = getindex.(out, 2)
  0.024785 seconds (9.24 k allocations: 416.910 KB)
7-element Array{Float64,1}:
 1.0
 0.965926
 0.866025
 0.707107
 0.5
 0.258819
 6.12323e-17

Reducing the number of allocations by 30% is nothing to sneeze at. Thank you.

I think I am safe making this even more concise:

@time sines, cosines = map(x->getindex.(out, x), 1:2)
  0.062047 seconds (20.81 k allocations: 956.831 KB)
2-element Array{Array{Float64,1},1}:
 [0.0,0.258819,0.5,0.707107,0.866025,0.965926,1.0]
 [1.0,0.965926,0.866025,0.707107,0.5,0.258819,6.12323e-17]

And thank you Colin T Bowers for suggesting that I can define a custom method for trig. I will definitely consider doing this if getindex fails to deliver the performance I want.

  • You can try `sines, cosines = map(collect, zip(out...))` too. May not be very efficiency though. – 张实唯 May 30 '17 at 16:16
  • 2
    Don't benchmark in global scope. Wrap the code in a function or use `BenchmarkTools` and `@btime getindex.($out, 1)` instead. – tim May 30 '17 at 16:57