# Broadcasting over arrays of SVectors

Can you use `.=` across an array of SVectors to reduce the number of allocations?

Yes.

4 Likes

This is the kind of thing that it’s easy to just try yourself first, and then ask if you have difficulties or specific questions, backed up by actual code that shows precisely what the question is.

5 Likes

You’re right, I should have asked a better question. I didn’t want to show my entire code because I don’t think it’s fair to ask other people to debug my code. But, I’ll post it here and if someone can see the problem quickly, that’s great. If not, that’s OK too. Here is the code:

``````#using Pkg
using StaticArrays

struct MD
nParticles::Int64
boxSize::Float64
positions:: Array{SVector{2,Float64},1}
velocities:: Array{SVector{2,Float64},1}

end

function initializeParticles(nParticles::Int64,boxSize::Float64,lVecs::Array{Float64,2},random::Bool = false):: MD
positions = zeros(SVector{2,Float64},nParticles)
velocities = zeros(SVector{2,Float64},nParticles)
deltaR = 0.5
v0=1.

idx = 0
for i =-boxSize-10:boxSize+10
for j =-boxSize-10:boxSize+10
vec = i * lVecs[1,:] + j * lVecs[2,:] + [0.6, 0.6]
if random
rng = MersenneTwister(1234)
rVec = rand(rng,2)
vec +=  deltaR * rVec/norm(rVec)
end
if 0 <= vec <= boxSize && 0 <= vec <= boxSize && idx < nParticles
idx += 1
positions[idx] = vec

rng = MersenneTwister(1234)
vVec = rand(rng,2)
velocities[idx] = v0 * vVec/norm(vVec)

end
end
end
if idx < nParticles
println("Didn't find all of the requested particles")
println("Requested: ", nParticles)
println("Found: ", idx)
end
MD(nParticles,boxSize,positions,velocities)

end

function forceOnSingleParticle(positions::Array{SVector{2,Float64},1},particle::SVector{2,Float64},boxSize::Float64):: SVector{2,Float64}
fVec = SVector(0,0)
modifiedPos = SVector{2,Float64}
diffVec = SVector{2,Float64}
for i=1:size(positions,1)
diffVec = particle - positions[i]

if abs(diffVec) > boxSize/2 && abs(diffVec) > boxSize/2
modifiedPos = positions[i] + boxSize * SVector(sign(diffVec) , sign(diffVec))
elseif abs(diffVec) > boxSize/2
modifiedPos = positions[i] + boxSize * SVector(sign(diffVec) , 0)
elseif abs(diffVec) > boxSize/2
modifiedPos = positions[i] + boxSize * SVector(0 , sign(diffVec))
else
modifiedPos = positions[i]
end
diffVec = particle - modifiedPos
distance = norm(diffVec)

if distance > 0.5
fVec += 24 * (2/distance^13 - 1/distance^7) * 1/distance * diffVec
end
end
return fVec

end

mymod(x::SVector,y::Number)::SVector = mod.(x,y)

function increaseTemp(positionsP::Array{SVector{2,Float64},1},positionsPP::Array{SVector{2,Float64},1},R::Float64)::Array{SVector{2,Float64},1}

positionsPP = positionsP - R * (positionsP - positionsPP)

end

function simulate(myMD::MD,dt::Float64,tMax::Float64)
forces = zeros(SVector{2,Float64},size(myMD.positions,1))
check = zeros(SVector{2,Float64},size(myMD.positions,1))
T::Float64 = 0.4
positions = myMD.positions
positionsP = myMD.positions
positionsPP = positions - myMD.velocities * dt
velocities = myMD.velocities
R::Float64 = 0.95
N::Int = size(positions,1)
while T > 0.2
for time=0:dt:tMax
#if any(positions > myMD.boxSize)
#    exit("lost particles")
#end
for n=1:N
forces[n] = forceOnSingleParticle(positions,positions[n],myMD.boxSize)
end
positions .= mymod.(2 .* positionsP .- positionsPP .+ forces .* dt^2,myMD.boxSize)
#velocities = (positions - positionsPP)/(2 * dt)
positionsPP = positionsP
positionsP = positions
end
#positionsPP = increaseTemp(positionsP,positionsPP,R)  #Increase temperature by factor of R
break
#T -=0.3#*= R
end

positions
end
using Plots
plotly()
dt = .0001
tMax = 0.9
lvecs = 1.07 * [[1.  0.]; [0.5  sqrt(3)/2]]
myMD = initializeParticles(16,4.,lvecs)
@time mypositions = simulate(myMD,dt,tMax)
plot(map(x -> x,mypositions),map(x -> x,mypositions),seriestype = :scatter,aspect_ratio=:equal,xlim = (0,4),ylim = (0,4))
``````

The key line of code here is line 102:

``````positions .= mymod.(2 .* positionsP .- positionsPP .+ forces .* dt^2,myMD.boxSize)
``````

As it stands, the array doesn’t seem to update. But if I replace the `.=` with just `=`, then things behave as expected.

1 Like

It definitely should be updating. Minimal example:

``````julia> x = [@SVector(rand(4)) for _ ∈ 1:5]
5-element Vector{SVector{4, Float64}}:
[0.38783754617174626, 0.23479994102954094, 0.9131772241755032, 0.5158254036871766]
[0.5377088066052653, 0.052031912322721086, 0.03460283315171697, 0.3172115916348688]
[0.008984784765390552, 0.16111534758886692, 0.653866563808464, 0.853094814027981]
[0.3532881187190555, 0.7784351001503331, 0.007738316370001996, 0.952768794661861]
[0.7145584888646308, 0.6944186486326702, 0.8473666362930881, 0.7211799815251871]

julia> @. x = 2x
5-element Vector{SVector{4, Float64}}:
[0.7756750923434925, 0.4695998820590819, 1.8263544483510064, 1.031650807374353]
[1.0754176132105306, 0.10406382464544217, 0.06920566630343394, 0.6344231832697376]
[0.017969569530781104, 0.32223069517773384, 1.307733127616928, 1.706189628055962]
[0.706576237438111, 1.5568702003006663, 0.015476632740003993, 1.905537589323722]
[1.4291169777292616, 1.3888372972653404, 1.6947332725861761, 1.4423599630503743]

julia> x
5-element Vector{SVector{4, Float64}}:
[0.7756750923434925, 0.4695998820590819, 1.8263544483510064, 1.031650807374353]
[1.0754176132105306, 0.10406382464544217, 0.06920566630343394, 0.6344231832697376]
[0.017969569530781104, 0.32223069517773384, 1.307733127616928, 1.706189628055962]
[0.706576237438111, 1.5568702003006663, 0.015476632740003993, 1.905537589323722]
[1.4291169777292616, 1.3888372972653404, 1.6947332725861761, 1.4423599630503743]
``````
3 Likes

I’ve tried simple examples like that too and convinced myself that this should work. I wondered if it had something to do with my `mymod` function so I took that out but it behaved the same.

It is updating:

``````p = copy(positions)
positions .= mymod.(2 .* positionsP .- positionsPP .+ forces .* dt^2,myMD.boxSize)
return p - positions

``````
``````julia> simulate(myMD,dt,tMax)
16-element Array{SArray{Tuple{2},Float64,1,2},1}:
[-6.340483565309141e-5, -7.957531458613332e-5]
[-6.346290042008595e-5, -7.919910996889712e-5]
[-6.328366500352178e-5, -7.921290476664566e-5]
[-6.111917607287687e-5, -7.935480894660785e-5]
[-6.331523406521988e-5, -7.907629591796805e-5]
[-6.113201763979781e-5, -7.921480146100279e-5]
[-6.110647734414165e-5, -7.920980100584174e-5]
[-6.096459454729697e-5, -7.932846379254954e-5]
[-6.110741666365271e-5, -7.90963577439463e-5]
[-6.0965533867030075e-5, -7.92150205306541e-5]
[-6.093999357137392e-5, -7.921002007593714e-5]
[-5.8756777145951844e-5, -7.934852561852779e-5]
[-6.095283513829486e-5, -7.9070012589888e-5]
[-5.878834620753892e-5, -7.921191677007222e-5]
[-5.8609110791252306e-5, -7.922571156759872e-5]
[-5.866717555802481e-5, -7.884950695036252e-5]

``````

By chance weren’t you trying to test using something like:

``````p = positions
positions .= mymod.(2 .* positionsP .- positionsPP .+ forces .* dt^2,myMD.boxSize)
return p - positions

``````

In this case we don’t see any different between `p` and `positions`, but that is because `=` is only assigning a different name to the same vector.

``````function f!(positions,positionsP,positionsPP,forces,dt,boxsize)
positions .= mymod.(2 .* positionsP .- positionsPP .+ forces .* dt^2,boxsize)
return positions
end

function g!(positions,positionsP,positionsPP,forces,dt,boxsize)
for i in eachindex(positions)
positions[i] = 2 * positionsP[i] - positionsPP[i] + forces[i] * dt^2
positions[i] = mod.(positions[i],boxsize)
end
return positions
end

``````
``````julia> positions = rand(SVector{3,Float64},5);

julia> positionsP = rand(SVector{3,Float64},5);

julia> positionsPP = rand(SVector{3,Float64},5);

julia> forces = rand(SVector{3,Float64},5);

julia> dt = 0.1; boxsize = 0.2;

julia> f!(copy(positions),positionsP,positionsPP,forces,dt,boxsize) == g!(copy(positions),positionsP,positionsPP,forces,dt,boxsize)
true

julia> @btime f!(p,\$positionsP,\$positionsPP,\$forces,\$dt,\$boxsize) setup=(p=copy(positions))
163.913 ns (0 allocations: 0 bytes)

julia> @btime g!(p,\$positionsP,\$positionsPP,\$forces,\$dt,\$boxsize) setup=(p=copy(positions))
141.832 ns (0 allocations: 0 bytes)

``````

I was just visually inspecting the plot. I actually think that with `.=` there was some modification but it was way different (smaller) than what happened with `=`

What am I missing? I don’t think the line `@. x = 2x` worked? Edit: sorry, I guess it is the case “in the morning I get up and I wake up an hour later”.

Looks like each number got multiplied by 2 to me. ??

Using `=` produced this Using `.=` produced this: Which looks very close(if not equal) to the initial positions.

It is likely that the problem is in these lines:

``````    positions = myMD.positions
positionsP = myMD.positions

``````

Here, you are just assigning two different names to the same `myMD.positions` vector. Thus, that can cause you all sort of problems latter. You should use `positionsP = copy(positions)`.

And check your other assignments in the loop. You should use `positionsP .= positions` (with the dot), otherwise you will have the same problem. The `=` just assigns a new name to the variable, while `.=` is mutating the elements of the existing array, which is what you want in most cases.

2 Likes

Hmm. Great example. I was trying to eliminate all of the allocations too. I started with this:
`positions = 2 * positionsP - positionsPP + forces * dt^2`

and found that I could reduce allocations with use of `.`. I didn’t see that a loop over the array would have produced no allocations either. Gonna ponder on this a little. Thanks

Hmm. Thinking… Can you explain why those two lines didn’t cause problems when I used `=` in the update.

The thing is that in this line:

``````positions = mymod.(2 .* positionsP .- positionsPP .+ forces .* dt^2,myMD.boxSize)
``````

without the `.`, what you are doing is creating a new array on the right side, and naming it `positions`.
Therefore, in this case the `positions` vector got dissociated when updated from the `positionsP` vector (which is what you want).

Like this:

``````julia> x = [1,2] # create new array and name it x
2-element Array{Int64,1}:
1
2

julia> y = x # just to save the values, this only assings a new name to x
2-element Array{Int64,1}:
1
2

julia> x = x + [1,1] # This is creating a new vector on right side and naming it x
2-element Array{Int64,1}:
2
3

julia> y # continues to be the same as before, no mutation
1

julia> x .= x .+ [1,1]  # with the dot, you are mutating the `x` vector
2-element Array{Int64,1}:
3
4

``````

Reading suggestion: Assignment and mutation · JuliaNotes.jl

2 Likes

I see. Thanks for the help. I’ll keep studying.

That admittedly is confusing at the beginning. One gets used to it, though, and after that one notes that actually there is no way out from that if one wants to have a dynamically typed language.

The fact that the language is dynamic requires that we can do

``````julia> x = [1,2]
2-element Array{Int64,1}:
1
2

julia> x = π
π = 3.1415926535897...

``````

That is completely natural, but means that `=` is just the binding of a name to a value (which might be an array, or scalar, or whatever).

Thus, we need to be able to differentiate naming something from mutating something. Mutating is a function that acts on a mutable object. It is the `setindex!` function:

``````julia> x = [1,2]
2-element Array{Int64,1}:
1
2

julia> setindex!(x,10,1)
2-element Array{Int64,1}:
10
2
``````

Which, by convenience (obviously) can be called with the notation:

``````julia> x = 10
10

``````

But this last `x = 10` is a call to `setindex!`, not a name assignment as the other cases. (and a broadcasting of assignments, with `.=`, is just a loop calling `setindex!` for each element)

(ps, if you feel better, here I am asking the same thing two years ago)

1 Like

Can I try and say it in a different way? These lines:

``````positions = myMD.positions
positionsP = myMD.positions
``````

essentially created two names assigned to the same memory block. Then later when I did this:

``````positions .= mymod.(2 .* positionsP .- positionsPP .+ forces .* dt^2,myMD.boxSize)
``````

it mutated that memory block and hence changed the value of both `positions` and `positionsP` Later when I tried to shuffle past and present snapshots around

``````positionsPP = positionsP
positionsP = positions
``````

it made both arrays equal to the same block of memory. Did I say that right?

Take home message: Careful when doing `var1 = var2`.

3 Likes

Exactly.

1 Like