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
#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.
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.
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.
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[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
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[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)
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