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
#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

image

Using .= produced this:

image

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