# Speeding up force calculations and mutable structs

**URL:** <https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228>\
**Category:** Performance\
**Created:** [August 4, 2020, 3:41am UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228 "2020-08-04T03:41:39Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Michael\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_wang/32/8535_2.png) [@Michael\_Wang](https://discourse.julialang.org/u/Michael_Wang)\
**Post date:** [August 4, 2020, 3:41am UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/1 "2020-08-04T03:41:39Z")

</div>

I put together a simple package to simulate many interacting Brownian particles. I’m trying to speed up the force calculation (for Lennard-Jones). Below is the force calculation using cell lists to find neighbors. The part I want to bring your attention to is towards the end (although the rest of the code could be relevant). If it helps, I could link the package or attach more code. I apologize if this is difficult to follow. I’m not sure what code to link here since it’s all interconnected.

```julia
function compute_pair_interaction!(lj::LennardJones; period_x::Float64 = -1.0, period_y::Float64 = -1.0)
    @use_threads lj.multithreaded for particle in lj.particles
        x, y = particle.x, particle.y
        i = trunc(Int64, x / lj.cell_list.cell_spacing_x)
        j = trunc(Int64, y / lj.cell_list.cell_spacing_y)
        @inbounds for di = -1 : 1, dj = -1 : 1
            idi = mod(i + di, lj.cell_list.num_cells_x) + 1
            jdj = mod(j + dj, lj.cell_list.num_cells_y) + 1
            pid = lj.cell_list.start_pid[idi, jdj]
            while pid > 0
                neighbor = lj.cell_list.particles[pid]
                Δx = wrap_displacement(x - neighbor.x; period = period_x)
                Δy = wrap_displacement(y - neighbor.y; period = period_y)
                Δr² = Δx^2 + Δy^2

                σ² = (lj.σ < 0.0 ? particle.R + neighbor.R : lj.σ)^2
                cutoff² = (lj.cutoff < 0.0 ? LJ_CONST * σ² : lj.cutoff^2)
                if 0.0 < Δr² < cutoff²
                    coef = lj_coef(lj.ϵ, σ², Δr²)
                    f_x, f_y = coef * Δx, coef * Δy

                    #THIS PART SEEMS TO BE SLOW
                    particle.f_x += f_x
                    particle.f_y += f_y
                    if lj.use_newton_3rd
                        neighbor.f_x -= f_x
                        neighbor.f_y -= f_y
                    end
                end
                pid = lj.cell_list.next_pid[pid]
            end
        end
    end
end

```

The chunk of code

```julia
particle.f_x += f_x
particle.f_y += f_y
if lj.use_newton_3rd
    neighbor.f_x -= f_x
    neighbor.f_y -= f_y
end

```

seems to take up about half the computation time (using `@benchmark`, the whole function takes about 1 ms). `particle` and `neighbor` are a custom type `Particle` containing all the information for a particle:

```julia
mutable struct Particle
    ptype::Symbol # particle type

    x::Float64 # x position
    y::Float64 # y position
    θ::Float64 # orientation

    f_x::Float64 # x force
    f_y::Float64 # y force
    t_θ::Float64 # torque in z

    R::Float64 # particle radius
    γ_trans::Float64 # translational friction
    γ_rot::Float64 # rotational friction
    D_trans::Float64 # translational diffusivity
    D_rot::Float64 # rotational diffusivity

    active_force::Union{AbstractActiveForce, Nothing} # active force/self-propulsion
end

```

I did some tests by commenting out different parts of the calculation.

**TEST 1:** I commented out

```julia
particle.f_x += f_x
particle.f_y += f_y

```

**TEST 2:** I commented out

```julia
if lj.use_newton_3rd
    neighbor.f_x -= f_x
    neighbor.f_y -= f_y
end

```

**TEST 3:** I commented out both chunks.

Note that I have also tested this without multithreading and get similar speed-up. Now, TEST 1 and TEST 2 both give ~1 ms for `@benchmark`. TEST 3, however, gives ~500 \mu\textrm{s} (half the time) for `@benchmark`. It seems like they individually don’t contribute a lot to the time but collectively they do. _ **Is this part slow because of my use of mutable structs?** _ There are a few other places where I use a similar `struct.field +=` and those don’t seem to contribute nearly as much time. I know that most other codes store particle positions in an array, so if it comes to it, I’ll probably have to rewrite everything, but for the moment I like the organization of `Particle`.

Another thing that is a little strange is that I am testing with `lj.use_newton_3rd` set as `false`. Logically, this means that TEST 1 should be equivalent to TEST 3, but it still takes ~1 ms instead of ~500 \mu\textrm{s}. It only goes down to ~500 \mu\textrm{s} if I replace `if lj.use_newton_3rd` with `if false`. I’m guessing this has to do with there being a lot of background code being generated with `if lj.use_newton_3rd`?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [August 4, 2020, 3:54am UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/2 "2020-08-04T03:54:18Z")

</div>

My guess is that what’s happening is that the compiler realizes that if you don’t modify either `particle` or `neighbor` that your function doesn’t have any effects. Thus the compiler is free to remove most of the code of your function.

---

<div class="post-metadata">

**Author:** ![Michael\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_wang/32/8535_2.png) [@Michael\_Wang](https://discourse.julialang.org/u/Michael_Wang)\
**Post date:** [August 4, 2020, 4:08am UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/3 "2020-08-04T04:08:47Z")

</div>

Thanks for your response! Do you mean that the compiler doesn’t know the outcome of the if statements until runtime so it must generate code for both `true` and `false`? If so, is there anyway around this? I guess a minimal example is the following

```julia
using BenchmarkTools

mutable struct mystruct
    x::Float64
    y::Float64

    b::Bool

    u::Float64
    v::Float64
end

function myfunction(arr::Array{mystruct, 1})
    for a in arr
        for _ = 1 : 10
            val = cos(1.0)
            
            a.x += val
            a.y += val
            if a.b
                a.u -= val
                a.v -= val
            end
        end
    end
end

test = [mystruct(0.0, 0.0, false, 0.0, 0.0) for _ = 1 : 100000]
@benchmark myfunction($test)

```

**TEST 1:** comment out

```julia
a.x += val
a.y += val

```

**TEST 2:** comment out

```julia
if a.b
    a.u -= val
    a.v -= val
end

```

**TEST 3:** comment out both.

The whole function takes ~1.4 ms. TEST 1 and TEST 2 still give ~1.4 ms. TEST3 gives ~400 \mu\textrm{s}. I think this captures the issue.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [August 4, 2020, 7:05am UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/4 "2020-08-04T07:05:38Z")

</div>

There is no real point benchmarking against a function with the assignments removed, because the compiler might just cut everything else as well. Instead benchmark alternate methods that actually do what you want.

But from reading it, writing to mutable structs neighbors and particle is not the fastest way to do this. Just use a variable `neighbor_f_x` etc inside the for loop, not `neighbor.f_x`. Then copy it over outside the loop if you have to.

---

<div class="post-metadata">

**Author:** ![bernhard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bernhard/32/2619_2.png) [@bernhard](https://discourse.julialang.org/u/bernhard)\
**Post date:** [August 4, 2020, 8:03am UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/5 "2020-08-04T08:03:04Z")

</div>

Have you considered profiling?  
The packages pprof.jl and profileview.jl may help you understand which lines/function calls take up a lot of time.

---

<div class="post-metadata">

**Author:** ![Duane\_Wilson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/duane_wilson/32/9506_2.png) [@Duane\_Wilson](https://discourse.julialang.org/u/Duane_Wilson)\
**Post date:** [August 4, 2020, 11:26am UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/6 "2020-08-04T11:26:51Z")

</div>

There was a good talk at Juliacon about this. Maybe look at:

> **[GitHub - jw3126/Setfield.jl: Update deeply nested immutable structs.](https://github.com/jw3126/Setfield.jl)**
>
> Update deeply nested immutable structs. Contribute to jw3126/Setfield.jl development by creating an account on GitHub.

This will let you write immutable structs but wrap them in lens to simulate mutability.

---

<div class="post-metadata">

**Author:** ![Michael\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_wang/32/8535_2.png) [@Michael\_Wang](https://discourse.julialang.org/u/Michael_Wang)\
**Post date:** [August 4, 2020, 2:00pm UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/7 "2020-08-04T14:00:21Z")

</div>

> [@Raf](#):
>
> There is no real point benchmarking against a function with the assignments removed, because the compiler might just cut everything else as well. Instead benchmark alternate methods that actually do what you want.
> 
> But from reading it, writing to mutable structs neighbors and particle is not the fastest way to do this. Just use a variable `neighbor_f_x` etc inside the for loop, not `neighbor.f_x` . Then copy it over outside the loop if you have to.

Ya, I realized that that is not the proper way to test things. I tried an alternative where I created an array `f_thread = [[0.0, 0.0] for _ = 1 : Threads.nthread()]` so that each thread can store the force with `f_thread[Threads.threadid()][1] += f_x` and similarly for `f_y`. This way I can move `particle.f_x += ...` outside the neighbor finding while loop. I still get pretty much the exact same performance. I also tried it without threads by using a stored `tempf = [[0.0, 0.0] for _ = 1 : N_particles]` and updated that inside the force calculation with `tempf[n][1] += f_x`. Then outside the entire force calculation I updated `particle.f_x`. This was actually slower (by about x1.5-x2), presumably because I now have two loops over the particles: one for calculating force and one for updating force.

---

<div class="post-metadata">

**Author:** ![Michael\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_wang/32/8535_2.png) [@Michael\_Wang](https://discourse.julialang.org/u/Michael_Wang)\
**Post date:** [August 4, 2020, 2:28pm UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/8 "2020-08-04T14:28:56Z")

</div>

> [@bernhard](#):
>
> Have you considered profiling?  
> The packages pprof.jl and profileview.jl may help you understand which lines/function calls take up a lot of time.

I tried profileview.jl. I ran `@profview` on my minimal example. The results varied but quite frequently there were a lot of calls to **abstractinterpretation.jl** and **typeinfer.jl**. I saw a similar thing with my `compute_pair_interactions!` when I used the built in julia profiler. Is this due to my custom structs? Here is an example. The only thing I see related to my function is at the end.

```julia
=========================================================
 ╎8 @Base\task.jl:356; (::Atom.var"#184#185")()
 ╎ 8 @Atom\src\eval.jl:41; macro expansion
 ╎ 8 @Base\essentials.jl:709; invokelatest(::Any, ::Any, ::Vararg{Any,N} where N)
 ╎ 8 @Base\essentials.jl:710; invokelatest(::Any, ::Any, ::Vararg{Any,N} where N; kwargs::Base.Iterators.Pairs{Union{},Union...
 ╎ 8 @Atom\src\eval.jl:113; eval(::String, ::Int64, ::String, ::String, ::Bool)
 ╎ 8 C:\Users\Michael Wang\.julia\packages\Media\ItEPc\src\dynamic.jl:24; macro expansion
 ╎ ╎ 8 @Atom\src\eval.jl:116; macro expansion
 ╎ ╎ 8 @Atom\src\repl.jl:127; hideprompt(::Atom.var"#200#205"{String,Int64,String,Bool})
 ╎ ╎ 8 @Atom\src\eval.jl:117; #200
 ╎ ╎ 8 @Base\logging.jl:514; with_logger
 ╎ ╎ 8 @Base\logging.jl:408; with_logstate(::Function, ::Any)
 ╎ ╎ ╎ 8 @Atom\src\eval.jl:118; (::Atom.var"#201#206"{String,Int64,String,Bool})()
 ╎ ╎ ╎ 8 @Atom\src\eval.jl:9; withpath(::Function, ::String)
 ╎ ╎ ╎ 8 @CodeTools\src\utils.jl:30; withpath(::Atom.var"#202#207"{String,Int64,String,Bool}, ::String)
 ╎ ╎ ╎ 8 @Atom\src\eval.jl:119; (::Atom.var"#202#207"{String,Int64,String,Bool})()
 ╎ ╎ ╎ 8 @CodeTools\src\eval.jl:30; include_string(::Module, ::String, ::String, ::Int64)
 ╎ ╎ ╎ ╎ 8 @Base\loading.jl:1096; include_string
4╎ ╎ ╎ ╎ 8 @Base\loading.jl:1088; include_string(::Function, ::Module, ::String, ::String)
 ╎ ╎ ╎ ╎ 4 @Base\compiler\typeinfer.jl:601; typeinf_ext(::Core.MethodInstance, ::UInt64)
 ╎ ╎ ╎ ╎ 4 @Base\compiler\typeinfer.jl:570; typeinf_ext(::Core.MethodInstance, ::Core.Compiler.Params)
 ╎ ╎ ╎ ╎ 3 @Base\compiler\typeinfer.jl:12; typeinf(::Core.Compiler.InferenceState)
 ╎ ╎ ╎ ╎ ╎ 3 @Base\compiler\abstractinterpretation.jl:1326; typeinf_nocycle(::Core.Compiler.InferenceState)
 ╎ ╎ ╎ ╎ ╎ 3 @Base\compiler\abstractinterpretation.jl:1270; typeinf_local(::Core.Compiler.InferenceState)
 ╎ ╎ ╎ ╎ ╎ 3 @Base\compiler\abstractinterpretation.jl:1005; abstract_eval(::Any, ::Array{Any,1}, ::Core.Compiler.InferenceState)
 ╎ ╎ ╎ ╎ ╎ 3 @Base\compiler\abstractinterpretation.jl:911; abstract_call(::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Core.Compiler.Infe...
 ╎ ╎ ╎ ╎ ╎ 3 @Base\compiler\abstractinterpretation.jl:926; abstract_call(::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Core.Compiler.Inf...
 ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\abstractinterpretation.jl:903; abstract_call_known(::Any, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Core...
 ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\typeutils.jl:46; argtypes_to_type
 ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\utilities.jl:39; anymap(::typeof(Core.Compiler.widenconst), ::Array{Any,1})
1╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\typelattice.jl:206; widenconst(::Core.Compiler.Const)
 ╎ ╎ ╎ ╎ ╎ ╎ 2 @Base\compiler\abstractinterpretation.jl:904; abstract_call_known(::Any, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Core...
 ╎ ╎ ╎ ╎ ╎ ╎ 2 @Base\compiler\abstractinterpretation.jl:134; abstract_call_gf_by_type(::Any, ::Array{Any,1}, ::Any, ::Core.Compiler.Inference...
 ╎ ╎ ╎ ╎ ╎ ╎ 2 @Base\compiler\abstractinterpretation.jl:266; abstract_call_method_with_const_args(::Any, ::Any, ::Array{Any,1}, ::Core.Simpl...
 ╎ ╎ ╎ ╎ ╎ ╎ 2 @Base\compiler\typeinfer.jl:33; typeinf(::Core.Compiler.InferenceState)
 ╎ ╎ ╎ ╎ ╎ ╎ 2 @Base\compiler\optimize.jl:174; optimize(::Core.Compiler.OptimizationState, ::Any)
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\driver.jl:134; run_passes(::Core.CodeInfo, ::Int64, ::Core.Compiler.OptimizationState)
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\driver.jl:127; slot2reg
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\slot2ssa.jl:626; construct_ssa!(::Core.CodeInfo, ::Core.Compiler.IRCode, ::Core.Compiler.DomT...
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\slot2ssa.jl:519; compute_live_ins(::Core.Compiler.CFG, ::Core.Compiler.SlotInfo)
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\sort.jl:780; sort##kw
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\sort.jl:780; #sort#8
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\abstractarray.jl:971; copymutable
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\array.jl:350; copyto!
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\array.jl:326; copyto!
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\array.jl:335; _copyto_impl!(::Array{Tuple{Int64,Int64,Bool},1}, ::Int64, ::Array{Tuple...
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\array.jl:293; unsafe_copyto!
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\driver.jl:138; run_passes(::Core.CodeInfo, ::Int64, ::Core.Compiler.OptimizationState)
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\inlining.jl:74; ssa_inlining_pass!
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\inlining.jl:999; assemble_inline_todo!(::Core.Compiler.IRCode, ::Core.Compiler.OptimizationSt...
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\inlining.jl:982; process_simple!(::Core.Compiler.IRCode, ::Int64, ::Core.Compiler.Params)
 ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\inlining.jl:20; with_atype
1╎ ╎ ╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\inlining.jl:18; Core.Compiler.Signature(::Any, ::Any, ::Any, ::Any)
 ╎ ╎ ╎ ╎ 1 @Base\compiler\typeinfer.jl:33; typeinf(::Core.Compiler.InferenceState)
 ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\optimize.jl:174; optimize(::Core.Compiler.OptimizationState, ::Any)
 ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\driver.jl:138; run_passes(::Core.CodeInfo, ::Int64, ::Core.Compiler.OptimizationState)
 ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\inlining.jl:74; ssa_inlining_pass!
 ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\inlining.jl:1046; assemble_inline_todo!(::Core.Compiler.IRCode, ::Core.Compiler.OptimizationState)
 ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\ssair\inlining.jl:703; analyze_method!(::Int64, ::Core.Compiler.Signature, ::Any, ::Core.SimpleVector, :...
1╎ ╎ ╎ ╎ ╎ ╎ 1 @Base\compiler\utilities.jl:131; specialize_method
 ╎1 @Base\threadingconstructs.jl:48; (::SoftSquishyMatter.var"#315#threadsfor_fun#68"{Float64,Float64,LennardJones,Array{Particle,1}})()
 ╎ 1 @Base\threadingconstructs.jl:81; (::SoftSquishyMatter.var"#315#threadsfor_fun#68"{Float64,Float64,LennardJones,Array{Particle,1}}...
1╎ 1 @SoftSquishyMatter\src\pairinteractions.jl:42; macro expansion
Total snapshots: 9

```

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [August 4, 2020, 2:36pm UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/9 "2020-08-04T14:36:17Z")

</div>

No, that doesn’t look normal. You’re seeing a profile of the Julia compiler itself, not of your code. To avoid any compilation steps showing up in profiled code, I typically do:

```julia
my_function() # ensure everything is compiled
Profile.clear() # ensure the profiler is cleared
@profile my_function() # actually collect profiling data

```

The ProfileView.jl package also makes visualizing profiler results much easier.

---

<div class="post-metadata">

**Author:** ![Michael\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_wang/32/8535_2.png) [@Michael\_Wang](https://discourse.julialang.org/u/Michael_Wang)\
**Post date:** [August 4, 2020, 3:23pm UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/10 "2020-08-04T15:23:23Z")

</div>

I did

```julia
test = [mystruct(0.0, 0.0, false, 0.0, 0.0) for _ = 1 : 100000]
myfunction(test)
Profile.clear()
@profile myfunction(test)
Profile.print()

```

for my minimal example above. It took a few times to collect samples, but I finally got the following

```julia
Overhead ╎ [+additional indent] Count File:Line; Function
=========================================================
 ╎1 @Base\task.jl:356; (::Atom.var"#184#185")()
 ╎ 1 @Atom\src\eval.jl:41; macro expansion
 ╎ 1 @Base\essentials.jl:709; invokelatest(::Any, ::Any, ::Vararg{Any,N} where N)
 ╎ 1 @Base\essentials.jl:710; invokelatest(::Any, ::Any, ::Vararg{Any,N} where N; kwargs::Base.Iterators.Pairs{Union{}...
 ╎ 1 @Atom\src\eval.jl:113; eval(::String, ::Int64, ::String, ::String, ::Bool)
 ╎ 1 ...s\Michael Wang\.julia\packages\Media\ItEPc\src\dynamic.jl:24; macro expansion
 ╎ ╎ 1 @Atom\src\eval.jl:116; macro expansion
 ╎ ╎ 1 @Atom\src\repl.jl:127; hideprompt(::Atom.var"#200#205"{String,Int64,String,Bool})
 ╎ ╎ 1 @Atom\src\eval.jl:117; #200
 ╎ ╎ 1 @Base\logging.jl:514; with_logger
 ╎ ╎ 1 @Base\logging.jl:408; with_logstate(::Function, ::Any)
 ╎ ╎ ╎ 1 @Atom\src\eval.jl:118; (::Atom.var"#201#206"{String,Int64,String,Bool})()
 ╎ ╎ ╎ 1 @Atom\src\eval.jl:9; withpath(::Function, ::String)
 ╎ ╎ ╎ 1 @CodeTools\src\utils.jl:30; withpath(::Atom.var"#202#207"{String,Int64,String,Bool}, ::String)
 ╎ ╎ ╎ 1 @Atom\src\eval.jl:119; (::Atom.var"#202#207"{String,Int64,String,Bool})()
 ╎ ╎ ╎ 1 @CodeTools\src\eval.jl:30; include_string(::Module, ::String, ::String, ::Int64)
 ╎ ╎ ╎ ╎ 1 @Base\loading.jl:1096; include_string
 ╎ ╎ ╎ ╎ 1 @Base\loading.jl:1088; include_string(::Function, ::Module, ::String, ::String)
 ╎ ╎ ╎ ╎ 1 ...uments\NYU\Research\Julia Simulation Package\TEST.jl:41; myfunction(::Array{mystruct,1})
 ╎ ╎ ╎ ╎ 1 @Base\Base.jl:34; setproperty!
Total snapshots: 1

```

Line 41 is `a.y += val` which is the last thing since `a.b` is `false`. ProfileView just shows a bunch of equal length colored lines stacked on top of each other with no other features. It seems like there isn’t much of a way to speed up things outside of rewriting everything and using a different algorithm/container (e.g. storing everything in arrays).

---

<div class="post-metadata">

**Author:** ![Michael\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_wang/32/8535_2.png) [@Michael\_Wang](https://discourse.julialang.org/u/Michael_Wang)\
**Post date:** [August 4, 2020, 3:26pm UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/11 "2020-08-04T15:26:24Z")

</div>

> [@Duane\_Wilson](#):
>
> There was a good talk at Juliacon about this. Maybe look at:
> 
> ### [jw3126/Setfield.jl](https://github.com/jw3126/Setfield.jl)
> 
> Update deeply nested immutable structs. Contribute to jw3126/Setfield.jl development by creating an account on GitHub.
> 
> This will let you write immutable structs but wrap them in lens to simulate mutability.

Thanks I’ll take a look at this!

---

<div class="post-metadata">

**Author:** ![mikkoku](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikkoku/32/16274_2.png) [@mikkoku](https://discourse.julialang.org/u/mikkoku)\
**Post date:** [August 4, 2020, 3:39pm UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/12 "2020-08-04T15:39:51Z")

</div>

It looks like your code is too fast and profiling is only able to collect 1 sample, which in statistical terms is a really small sample and won’t give representative results. In general the profiled code should take a couple of seconds to run and produce at least a thousand samples.

If you are using Juno there is also `@profiler` macro.

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [August 4, 2020, 3:57pm UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/13 "2020-08-04T15:57:50Z")

</div>

Furthermore, I think profiling `myfunction()` is unlikely to be very informative in this case because the function is so simple. There’s no actual cost to getting or setting a field of a mutable struct (assuming that field has a concrete type in the struct definition) beyond the memory access itself, so the code produced by the compiler is just going to be a sequence of loads, additions, and stores. You’re not going to see a bunch of overhead from calling the `setfield` function or the `+` function because, by the time the compiler has optimized the code, those function calls will already have been inlined or eliminated.

If your original question is:

> [@Michael\_Wang](#):
>
> Is this part slow because of my use of mutable structs?

I’m pretty confident in saying “no”. If you are frequently _creating_ new instances of some mutable struct in a loop then yes, that can cause memory allocation and slow down your code. But just reading or writing the fields of a mutable struct is very cheap.

The only reason this might not be true is cache locality. Immutable structs in Julia can be stored inline in a `Vector`, meaning that adjacent particles in your array will have their data stored adjacently in memory. Mutable structs are not stored inline, so adjacent particles will not necessarily have their data stored nearby in memory. That additional indirection can cause your CPU to spend more time moving data in and out of its low-level caches.

Before going down the road of thinking about cache locality, however, I would suggest first trying all of the easy things:

- Is your function type-stable?
- Is `@use_threads` doing the right thing? That is, how does single-threaded performance (without any threading macro usage) compare?
- Does your function allocate memory in its inner loop?

---

<div class="post-metadata">

**Author:** ![Michael\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_wang/32/8535_2.png) [@Michael\_Wang](https://discourse.julialang.org/u/Michael_Wang)\
**Post date:** [August 4, 2020, 5:08pm UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/14 "2020-08-04T17:08:57Z")

</div>

- I used `@code_warntype` on `compute_pair_interaction!` and got no type issues. I also did the same for a simplified version where I removed all `for` loops and `if` statements (basically chose one `particle` and one `neighbor` and computed the force instead of looping over all of them). `@code_warntype` didn’t output anything bad. I’m not entirely sure this means my code is type-stable because there was an incident in the past where `@code_warntype` came back fine but it turned out that there was a statement like `one_of_many_structs_grouped_under_abstract_type.shared_field` that needed to be rewritten as `....shared_field::Float64` in order to have performant code.

- `@use_threads` is a macro that becomes `Threads.@threads` if `lj.multithreaded` is `true`; otherwise it becomes a normal `for` loop without threads (see [post](https://discourse.julialang.org/t/putting-threads-threads-or-any-macro-in-an-if-statement/41406/9)). With three threads, this gives about a x2 boost to speed compare to a single thread. Note that if I replace `@use_threads lj.multithreaded` with `Threads.@threads`, the performance is nearly identical.

- There are about 17 allocations (2.7 KiB) per function call. 16 of those allocations comes from the threading regardless if I use `@use_threads` or `Threads.@threads`. I haven’t quite found where that remaining allocation comes from, though it seems insignificant. Is it possible for `@time` or `@btime` to miss some allocations in a function?

Edit: I will try Duane\_Wilson’s suggestion with Setfield and see if switching all my structs to immutable has any effect.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [August 5, 2020, 4:22am UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/15 "2020-08-05T04:22:54Z")

</div>

I think you are missing my point because I used the wrong example… In your original example you are writing to `particle.x` many times inside the loop when you don’t have to. Just use a local variable and write that to `particle.x` outside the loop. I find that can improve performance sometimes.

And I hadn’t realised `neighbor` was being loaded from the cell list inside the loop - so you obviously need to update that inside the loop then too my mistake. Setfield may help, but writing to the mutable struct may actually be faster as well, it’s hard to say. The memory layout of your struct is going to matter, possibly StructsOfArrays.jl will help there.

---

<div class="post-metadata">

**Author:** ![bernhard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bernhard/32/2619_2.png) [@bernhard](https://discourse.julialang.org/u/bernhard)\
**Post date:** [August 5, 2020, 7:14am UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/16 "2020-08-05T07:14:59Z")

</div>

Have you considered the possibility that the code is already reasonably close to an optimal solution?  
Do you have timings from an alternative implementation?  
It seems that you are mainly looping over pids. How many pids are there for a typicall call?

---

<div class="post-metadata">

**Author:** ![Michael\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_wang/32/8535_2.png) [@Michael\_Wang](https://discourse.julialang.org/u/Michael_Wang)\
**Post date:** [August 5, 2020, 3:09pm UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/17 "2020-08-05T15:09:18Z")

</div>

If I understand you correctly, do you mean take `particle.x` outside the `while pid > 0` loop? If so, I tried that by using a local `f_x,f_y`, updating those, and then setting `particle.f_x,particle.f_y` outside the while loop. The performance is the essentially the same (I think probably because for this test system, the loop only works once. Actually your suggestion is probably how I should do it if the cells are larger and contain more than one particle). If I define a local variable like `f_x,f_y` inside the loop that is threaded, is it shared by all the threads or does each thread get its own copy?

---

<div class="post-metadata">

**Author:** ![Michael\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_wang/32/8535_2.png) [@Michael\_Wang](https://discourse.julialang.org/u/Michael_Wang)\
**Post date:** [August 5, 2020, 3:18pm UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/18 "2020-08-05T15:18:59Z")

</div>

Yes that is a good point. For non-interacting particles, my simulations runs at around 30k steps/s while with interactions it goes down to around 1k-2k, so I hoping to get somewhere in between those two numbers! But I guess there might not be much more I can do.

For very short ranged interactions like hard-ish disks in 2D (for which I am focused on now), the cell sizes are chosen so that you can essentially have only one particle per cell (possibly in very very rare instances two particles might somehow cram into one cell). However, in other types of simulations (e.g. different sized particles), there may be many pids per cell. There are other types of cell lists that I can use, but so far I’ve only implemented the simplest.

Edit: I do not know of other completely different implementations, though I am looking. I think this approach with cell lists is the most common used in other packages (e.g. hoomd-blue). There are definitely some improvements I can do to how I organize my data.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [August 5, 2020, 3:32pm UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/19 "2020-08-05T15:32:37Z")

</div>

In your original code it’s inside a nested loop, writing multiple times to a mutable struct. It should be faster to do that once at the end and use local variables the rest of the time, especially when you say that is the most expensive part of the function.

---

<div class="post-metadata">

**Author:** ![Michael\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/michael_wang/32/8535_2.png) [@Michael\_Wang](https://discourse.julialang.org/u/Michael_Wang)\
**Post date:** [August 5, 2020, 3:39pm UTC](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228/20 "2020-08-05T15:39:49Z")

</div>

Oh my bad, I was looking at the wrong loop. Yes, I can definitely lift it outside the nested loops `di,dj` loops. Just to check, if I define `f_x, f_y = 0.0, 0.0` inside the outer loop (the one that is threaded), are these local to each thread i.e. each thread updates it’s own copy of `f_x, f_y`? Thanks!

Edit: I did a simply test and it seems to be the case (they are at least local to the loop, but not sure if threads might mess with them)

[Next page](https://discourse.julialang.org/t/speeding-up-force-calculations-and-mutable-structs/44228.md?page=2)