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