Perculiar behavior with two dimensional array of vectors in Julia

Question

I want to create a two dimensional array of 2D vectors (to represent a vector field).

My code is something like this

N=10
dx=1/(N-1)
dy=1/(N-1)

#initial data
U=fill(zeros(2), (N, N))

for i=1:1:N
    for j=1:1:N
        U[i,j][1]=(i-1)*dx
        U[i,j][2]=(j-1)*dy
    end
end

print(U[5, 7])

The result is [1.0, 1.0], which is not what I want. I have no idea why. However, if I change the code to something like this

N=10
dx=1/(N-1)
dy=1/(N-1)

#initial data
U=fill(zeros(2), (N, N))

for i=1:1:N
    for j=1:1:N
        U[i,j]=[(i-1)*dx, (i-1)*dx]
    end
end

print(U[5, 7])

Then it print out the correct result, which is [0.4444444444444444, 0.6666666666666666]. So, what going on?

Answer 1

This behaviour is expected. Note the following:

julia> x = fill(zeros(1), 2)
2-element Array{Array{Float64,1},1}:
 [0.0]
 [0.0]

julia> x[1][1] = 5.0
5.0

julia> x[2][1]
5.0

I managed to change x[2][1] just by changing x[1][1]. You can probably guess the problem at this point. You are populating all the elements of your matrix with the same vector. Therefore when you mutate one, you are mutating all.

To get the behaviour you want, you could build your initial matrix like this:

x = [ zeros(2) for n in 1:N, m in 1:N ]

The key point here is to consider whether the first argument to your fill call is, or contains, a mutable. If it does not, then it'll work like you expected, e.g. fill(0.0, 2). But if it does contain a mutable, then the output of fill will contain pointers to the single mutable object, and you'll get the behaviour you've encountered above.

Note, my use of the word "contains" here is important, since an immutable that contains a mutable will still result in a pointer to the single mutable object and hence the behaviour you have encountered. So for example:

struct T1 ; x::Vector{Float64} ; end
T1() = T1(zeros(1))
x = fill(T1(), 2)
x[1].x[1] = 5.0
x[2].x[1]

still mutates the second element of x.