# Array cumulation

Hello,

I am trying to, quickly, accumulate an array without over-riding previous allocations. In particular consider the following example:

``````B = zeros(2,2);
A = [0.1 0.3; 0.2 0.5];
idx = [1 1; 1 1];
``````

I would like to add to B, according to `idx` all the elements in `A`. Since `idx` all maps to the first element, I would hope to get:

``````B = [1.1 0; 0 0];
``````

but with the following command I get:

``````@views B[idx] .+= A
B = [0.5 0; 0 0];
``````

I.e. only the last allocation is recorded, overwriting everything else. Is there a way out of this?

Thank you!

You may want to opt for a for-loop:

``````for (i, Aᵢ) in zip(idx, A)
B[i] += Aᵢ
end
``````
``````julia> B
2×2 Array{Float64,2}:
1.1  0.0
0.0  0.0
``````
1 Like

I was trying to avoid a loop because I have to do this many many times. I was under the impression there was a quick way of doing this, but maybe not?

For-loops are very fast in Julia. If you need to do it many times, you can just write a small function like so:

``````function cumulate!(B, A, idx)
for (i, Aᵢ) in zip(idx, A)
B[i] += Aᵢ
end
end
``````
``````julia> cumulate!(B, A, idx)

julia> B
2×2 Array{Float64,2}:
1.1  0.0
0.0  0.0
``````
4 Likes

I think Julia is specifically preventing you from doing this. It wants `B .+= A` to mean `B .= B .+ A` with the right hand side evaluated first, and since this contains views of the same data as on the left, it makes copies to be safe. I’m not sure I can always predict what it’s going to do in such cases, but I do think that anything which relies on broadcasting being done sequentially (or worse, in a particular order) is fishy. In this case, one way to check is by timing it – the 15MB allocation tells you it’s in fact making two copies, although I’m not sure what the second one is:

``````julia> N = 10^3; A = rand(N,N); B = zeros(N,N); idx = ones(Int, N,N);

julia> @time B[idx] .= A;
0.003550 seconds (6 allocations: 208 bytes)

julia> @time B[idx] .+= A;  # has B[idx] on the right
0.007897 seconds (7 allocations: 7.630 MiB)

julia> @time @views B[idx] .+= A;  # has aliased view on the right
0.007404 seconds (12 allocations: 15.259 MiB)

julia> @time B[idx];
0.004902 seconds (2 allocations: 7.629 MiB)

julia> @time @view B[idx];
0.000582 seconds (3 allocations: 128 bytes)
``````

As @stillyslalom says, if you do want it done sequentially, then you want a loop. There is no reason this need be slower – broadcasting is itself ultimately just loops.

2 Likes