Can you use `.=`

across an array of SVectors to reduce the number of allocations?

Yes.

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.

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
#Pkg.add("StaticArrays")
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[1] <= boxSize && 0 <= vec[2] <= 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[1]) > boxSize/2 && abs(diffVec[2]) > boxSize/2
modifiedPos = positions[i] + boxSize * SVector(sign(diffVec[1]) , sign(diffVec[2]))
elseif abs(diffVec[1]) > boxSize/2
modifiedPos = positions[i] + boxSize * SVector(sign(diffVec[1]) , 0)
elseif abs(diffVec[2]) > boxSize/2
modifiedPos = positions[i] + boxSize * SVector(0 , sign(diffVec[2]))
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[1],mypositions),map(x -> x[2],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.

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]
```

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.

Concerning your original question: while you can use broadcasting with static arrays, that is about convenience, not about saving allocations:

```
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.

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[1] # 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

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[1] = 10
10
```

But this last `x[1] = 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)

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`

.

Exactly.