Broadcasting with vectors of static vectors

General Usage
1

Any advice on how to broadcast vectors of SVectors in Julia 0.7?

eg:

julia> using StaticArrays

julia> const Vec = SVector{3, Float64}
SArray{Tuple{3},Float64,1,3}

julia> a=[Vec(i,i,i) for i in 1:5]
5-element Array{SArray{Tuple{3},Float64,1,3},1}:
 [1.0, 1.0, 1.0]
 [2.0, 2.0, 2.0]
 [3.0, 3.0, 3.0]
 [4.0, 4.0, 4.0]
 [5.0, 5.0, 5.0]

julia> a .+ Vec(2,2,2)
ERROR: DimensionMismatch("arrays could not be broadcast to a common size")
Stacktrace:
 [1] _bcs1 at ./broadcast.jl:444 [inlined]
 [2] _bcs at ./broadcast.jl:438 [inlined]
 [3] broadcast_shape at ./broadcast.jl:432 [inlined]
 [4] combine_axes at ./broadcast.jl:427 [inlined]
 [5] instantiate at ./broadcast.jl:266 [inlined]
 [6] materialize(::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(+),Tuple{Array{SArray{Tuple{3},Float64,1,3},1},SArray{Tuple{3},Float64,1,3}}}) at ./broadcast.jl:729
 [7] top-level scope at none:0

Or is this just some crazy dream?

2 Likes
2

Broadcasting works only on the outermost shape of each and every argument. As you’ve written it, you’re attempting to broadcast a 5-element vector and a 3-element vector.

Now, what you want is for the 3-element vector to behave like a “scalar” — that is, broadcast should repeat it for each element of the 5-element vector. To do that, you can “protect” it inside a 1-element container like a tuple:

julia> a .+ (Vec(2,2,2),)
5-element Array{SArray{Tuple{3},Float64,1,3},1}:
 [3.0, 3.0, 3.0]
 [4.0, 4.0, 4.0]
 [5.0, 5.0, 5.0]
 [6.0, 6.0, 6.0]
 [7.0, 7.0, 7.0]

Note that this doesn’t really have anything to do with StaticArrays — except for the fact that you might be slightly more likely to think of a static vector as a “scalar-like” thing unto itself.

8 Likes
3

@mbauman, won’t this apply + rather than .+? This should yield the same result, but + under-performs .+ for speed (IIRC).

The poor-man’s solution that I was going to suggest was to write your own addition function that iterates over vectors and then applies .+ to each vector. (Or, better yet, that goes element-by-element.)

4

Yes, it will apply + elementwise — broadcasting only works at the outermost level. It should have the same performance in this case, but it’s true that it won’t fuse if there are other operations.

5

How about
[v .+ Vec(2, 2, 2) for v in a]
?

6

In case you ever wonder how to deal with something thats nested very deep, I think I found an alright way of making a function that recursively broadcasts a single argument function. Here’s your example:

using StaticArrays
const Vec = SVector{3, Float64}

isbroadcastscalar(x) = Base.Broadcast.BroadcastStyle(typeof(x)) == Base.Broadcast.DefaultArrayStyle{0}()
function recursive_broadcast(f::Function)
    function (args...)
        if isbroadcastscalar(arr)
            f(arr)
        else
            recursive_broadcast(f).(arr)
        end
    end
end

a=[Vec(i,i,i) for i in 1:5]
b = Vec(2,2,2)

julia> recursive_broadcast(x -> x + b)(a)
5-element Array{SArray{Tuple{3},Any,1,3},1}:
 [[3.0, 3.0, 3.0], [3.0, 3.0, 3.0], [3.0, 3.0, 3.0]]
 [[4.0, 4.0, 4.0], [4.0, 4.0, 4.0], [4.0, 4.0, 4.0]]
 [[5.0, 5.0, 5.0], [5.0, 5.0, 5.0], [5.0, 5.0, 5.0]]
 [[6.0, 6.0, 6.0], [6.0, 6.0, 6.0], [6.0, 6.0, 6.0]]
 [[7.0, 7.0, 7.0], [7.0, 7.0, 7.0], [7.0, 7.0, 7.0]]
7

Thanks everyone! I didn’t know wrapping with a tuple would work to force something to act as a scalar. Good to know!

8

It is documented in the Arrays chapter:

treats any argument that is not an array, tuple or Ref (except for Ptr) as a “scalar”

but perhaps that could be emphasized further.

9

When you say that “it” is documented - are you referring to the fact that my initial example will attempt to broadcast a 5-element vector and a 3-element vector, or that wrapping in another container negates this?

I understood why my initial example behaved as it did - that is obvious/documented enough IMO. I wasn’t clear on how to change the expression to make it do what I wanted. I haven’t seen it documented that wrapping an argument in a one-element container makes that argument act as a scalar. I would take a stab at helping document it, but I don’t fully understand why it works like that.

10

I meant that the ways of “escaping” broadcast are documented.

11

Perhaps the language of “protecting” or “escaping” makes this feel like a special case — but it’s not at all. It’s just that a one-element container will repeat that one element such that it fills the final shape of the broadcast operation. That is, it just falls out of how broadcast works on the outer shapes of containers.

We could probably do a better job of explaining this in the manual, and an example like this would definitely be a nice addition.

4 Likes
12

This! This is what I was missing in my understanding. Thanks!

I’ll try to find some time to contribute to the manual (the least I can do given how much value I have got out of julia.)

2 Likes