# Trying to understand how to write good code (example)

**URL:** <https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101>\
**Category:** New to Julia\
**Created:** [February 14, 2017, 4:28pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101 "2017-02-14T16:28:09Z")\
**Posts on this page:** 15\
**Page:** 2

<div class="post-metadata">

**Author:** ![johnmyleswhite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnmyleswhite/32/31_2.png) [@johnmyleswhite](https://discourse.julialang.org/u/johnmyleswhite)\
**Post date:** [February 14, 2017, 11:46pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/21 "2017-02-14T23:46:25Z")

</div>

> [@DNF](#):
>
> In fact, I would expect not mutating to be faster; since the views do not allocate, you’re just applying the label rj to a view.

I’m a little lost at this point, but It seems like you’re still confusing variable names acting as mutable bindings with the mutability of the data structures those variable names are bound to.

A simple mental model for how variable names work is to assume that all variable names live in a dictionary. That dictionary is mutable. When you type `rj = x`, you are mutating that dictionary. The type of `x` is totally irrelevant and therefore the mutability of `x` is also irrelevant. The important thing is that you’re doing something like `mutate_name!("rj", x)`, whereas an expression like `rj[:] = x` is doing something like `mutate_value!(rj, :, x)`.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 15, 2017, 12:10am UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/22 "2017-02-15T00:10:33Z")

</div>

> [@johnmyleswhite](#):
>
> I’m a little lost at this point, but It seems like you’re still confusing variable names acting as mutable bindings with the mutability of the data structures those variable names are bound to.

I’m not the OP, btw, but also not completely confident of having the right mental model.

My expectation would be that mutating the name is cheap, while mutating values is a bit more expensive, while (at least in the current example) allocating a new array is quite expensive.

So:

```julia
rj = Wr[i,:,iw]

```

allocates a new matrix on the right, and binds the variable `rj` to it.

```julia
rj .= Wr[i,:,iw]

```

allocates a new array on the right, _and_ copies the data into `rj`. According to my benchmarks this is more expensive.

```julia
rj = @view Wr[i,:,iw]

```

makes a view, which just barely allocates a small amount of memory, and binds `rj` to it, while

```julia
rj .= @view Wr[i,:,iw]

```

makes a view on the right, and then, on top of that, copies the data from the view into `rj`. According to the benchmarks I’m looking at, the non-dotted version is cheaper in each case.

It _feels_ like I’m really wrong here, but I don’t see how.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 15, 2017, 12:22am UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/23 "2017-02-15T00:22:03Z")

</div>

> [@DNF](#):
>
> It feels like I’m really wrong here, but I don’t see how.

No, you’re right. It depends on what you’re trying to do though. If `rj` can, in the real application, be linked to the exact same memory as `Wr`, it’s best to just bind a new view to `rj`. However, if you plan on changing `rj` but don’t want `Wr` to change (or vice versa), then you need to keep the values separate and mutate.

As always, it’s about how you’re applying it that determines what the right method is.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 15, 2017, 3:41am UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/24 "2017-02-15T03:41:17Z")

</div>

> [@ChrisRackauckas](#):
>
> As always, it’s about how you’re applying it that determines what the right method is.

(One of the things that makes the present discussion difficult is that the code that was posted doesn’t actually _do_ anything. It just creates a bunch of matrix views or copies and then discards them. In that case, the optimal thing is just to delete all the code. 😉 )

---

<div class="post-metadata">

**Author:** ![mazzanti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mazzanti/32/32948_2.png) [@mazzanti](https://discourse.julialang.org/u/mazzanti)\
**Post date:** [February 15, 2017, 1:49pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/25 "2017-02-15T13:49:40Z")

</div>

OK guys, I keep working on this since I have to decide whether I shall invest more on Julia or surrender and move back to fortran. Since deleting all the code is not what I am looking for 🙂, I can put now an almost-entire version of the function

```julia
function Drift_Force(W::Walker)
    Np,dim,NW = size(W.r);
    F_drift = zeros(Np,dim,NW);
    ri = zeros(dim,1);
    rj = zeros(dim,1);
    rij = zeros(dim,1);
    aux = zeros(dim,1);
    for iw in 1:NW
        for i in 1:Np-1
            for k in 1:dim 
                ri[k] = W.r[i,k,iw];
            end;
            for j in i+1:Np
                for k in 1:dim 
                    rij[k] = W.r[i,k,iw]-ri[k]
                end;
                norm = sqrt(dot(rij,rij));
                aux = rij/norm;
                F_drift[i,:,iw] += aux;
                F_drift[j,:,iw] -= aux;                
            end;        
        end;    
    end;
end; 

```

which performs like this:  
0.070826 seconds (504.01 k allocations: 20.010 MB, 7.97% gc time)

while I see that most of the time/allocations go to the final steps of F\_drift updating, since if I remove the last two F\_drift lines

```julia
function Drift_Force(W::Walker)
    Np,dim,NW = size(W.r);
    F_drift = zeros(Np,dim,NW);
    ri = zeros(dim,1);
    rj = zeros(dim,1);
    rij = zeros(dim,1);
    aux = zeros(dim,1);
    for iw in 1:NW
        for i in 1:Np-1
            for k in 1:dim 
                ri[k] = W.r[i,k,iw];
            end;
            for j in i+1:Np
                for k in 1:dim 
                    rij[k] = W.r[i,k,iw]-ri[k]
                end;
                norm = sqrt(dot(rij,rij));
                aux = rij/norm;
            end;        
        end;    
    end;
end; 

```

then what I get is  
0.003131 seconds (20.17 k allocations: 2.169 MB)

that is, much less time, much less allocations etc etc. So the problem is now, filling the results into the F\_drift array (which, btw, is what the function must return in the end).  
Do you know how shall I write everything here for performance, or is it in its best way already?

Thanks for the inputs again,

Feran.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 15, 2017, 2:29pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/27 "2017-02-15T14:29:06Z")

</div>

Looks like you should do something like:

```julia
function Drift_Force(W::Walker)
    Np,dim,NW = size(W.r);
    F_drift = zeros(Np,dim,NW);
    ri = zeros(dim);
    rij = zeros(dim);
    @inbounds for iw in 1:NW
        for i in 1:Np-1
            for k in 1:dim 
                ri[k] = W.r[i,k,iw];
            end;
            for j in i+1:Np
                @simd for k in 1:dim 
                    rij[k] = W.r[i,k,iw]-ri[k]
                end
                normalize!(rij)
                @simd for k in 1:dim 
                    F_drift[i,k,iw] += rij[k]
                    F_drift[j,k,iw] -= rij[k]
                end             
            end
        end 
    end
    return F_drift # presumably you want to return the result!
end

```

(Note that you don’t need all of the semicolons, and there is no need for `ri` and `rij` to be 2d arrays.)

Note also that you might want to re-order your dimensions so that `F_drift = zeros(dim,Np,NW)` and similarly for `W.r`, since that way the loops over `k` will be contiguous in memory.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 15, 2017, 2:32pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/28 "2017-02-15T14:32:02Z")

</div>

> [@mazzanti](#):
>
> `rij[k] = W.r[i,k,iw]-ri[k]`

Is this supposed to be `rij[k] = W.r[j,k,iw]-ri[k]`? Otherwise it is always zero.

---

<div class="post-metadata">

**Author:** ![Ferran\_Mazzanti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ferran_mazzanti/32/7458_2.png) [@Ferran\_Mazzanti](https://discourse.julialang.org/u/Ferran_Mazzanti)\
**Post date:** [February 15, 2017, 7:45pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/29 "2017-02-15T19:45:55Z")

</div>

Yes indeed, it is suppossed to be rij[k] = W.r[j,k,iw]-ri[k]… my fault ☹

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 15, 2017, 8:45pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/30 "2017-02-15T20:45:31Z")

</div>

You have some unnecessary cruft in your code. For example, there seems to be no point to defining `ri`, you can just write

```julia
rij[k] = W.r[j, k, iw] - W.r[i, k, iw]

```

without wasting time on `ri`. Also I don’t quite understand your code, so there might be algorithmic speedups possible, beyond the pure implemetation speedups.

Anyway, here are three versions of your code. The first one is basically your own version. The second is a cleaned-up version. The third one assumes that `dim` is always equal to `3`, and uses `StaticArrays`:

```julia
using StaticArrays, BenchmarkTools

function Drift_Force1(W::Array{Float64, 3})
    Np,dim,NW = size(W);
    F_drift = zeros(Np,dim,NW);
    ri = zeros(dim,1);
    rj = zeros(dim,1);
    rij = zeros(dim,1);
    aux = zeros(dim,1);
    for iw in 1:NW
        for i in 1:Np-1
            for k in 1:dim 
                ri[k] = W[i,k,iw];
            end
            for j in i+1:Np
                for k in 1:dim 
                    rij[k] = W[j,k,iw]-ri[k]
                end
                norm = sqrt(dot(rij,rij));
                aux = rij/norm;
                F_drift[i,:,iw] += aux;
                F_drift[j,:,iw] -= aux;                
            end        
        end    
    end
    return F_drift
end

function Drift_Force2(W::Array{Float64, 3})
    dim = size(W, 2)
    F_drift = zeros(W)
    rij = zeros(dim)
    @inbounds for iw in 1:size(W, 3)
        for i in 1:size(W, 1)-1
            for j in i+1:size(W, 1)
                for k in 1:dim
                    rij[k] = W[j, k, iw] - W[i, k, iw]
                end
                normalize!(rij)
                for k in 1:dim
                    F_drift[i, k, iw] += rij[k]
                    F_drift[j, k, iw] -= rij[k]
                end
            end
        end    
    end
    return F_drift
end

function Drift_Force3(W::Array{Float64, 3})
    dim = 3
    F_drift = zeros(W)
    @inbounds for iw in 1:size(W, 3)
        for i in 1:size(W, 1)-1
            for j in i+1:size(W, 1)
                rij = SVector(W[j, 1, iw] - W[i, 1, iw],
                              W[j, 2, iw] - W[i, 2, iw],
                              W[j, 3, iw] - W[i, 3, iw])
                rij = normalize(rij)
                for k in 1:dim
                    F_drift[i, k, iw] += rij[k]
                    F_drift[j, k, iw] -= rij[k]
                end
            end
        end    
    end
    return F_drift
end

WW = rand(64,3,10)
maximum(Drift_Force1(WW) .- Drift_Force2(WW))
maximum(Drift_Force1(WW) .- Drift_Force3(WW))
@benchmark Drift_Force1(WW)
@benchmark Drift_Force2(WW)
@benchmark Drift_Force3(WW)

```

On my computer version 1 runs in 26ms, version 2 in 2.2ms, and version 3 in 190μs, so that’s well over 100x speed-up. If you re-arrange the input matrix so that the size is 3 along the _first_ dimension, there’s an extra ~10% to gain.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [February 15, 2017, 9:40pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/31 "2017-02-15T21:40:39Z")

</div>

You don’t seem to have used `StaticArrays` in the code you posted.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 15, 2017, 9:42pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/32 "2017-02-15T21:42:33Z")

</div>

Try to scroll down a bit 😉

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [February 15, 2017, 9:43pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/33 "2017-02-15T21:43:16Z")

</div>

Oops, thanks.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [February 15, 2017, 9:44pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/34 "2017-02-15T21:44:11Z")

</div>

Seems like it shouldn’t be hard to parametrize this on the dimension.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 15, 2017, 9:49pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/35 "2017-02-15T21:49:48Z")

</div>

Perhaps. I don’t see it. Some `Val{3}` kind of magic? Do you mean the number of dimensions or the size of the dimensions?

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [February 16, 2017, 8:36pm UTC](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101/36 "2017-02-16T20:36:52Z")

</div>

> <https://stackoverflow.com/questions/33002572/creating-copies-in-julia-with-operator/33003055#33003055>

I should probably update that answer to include `.=` which is a shorthand for calling `broadcast!` which in turn, by default, does a bunch of `[]=` indexed assignments (i.e. calls to `setindex!`.

[Previous page](https://discourse.julialang.org/t/trying-to-understand-how-to-write-good-code-example/2101.md?page=1)