Perculiar behavior with two dimensional array of vectors in Julia

user_13472777 0
Viewed 106

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?

julia
2 Answers

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.

user_16678950

The issue is that fill(zeros(2), (N, N)) points to a single 2-element vector. So if you change the contents of one, they all change:

julia> U = fill(ones(2), 2)
2-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [1.0, 1.0]

julia> U[1][1] = 2
2

julia> U
2-element Vector{Vector{Float64}}:
 [2.0, 1.0]
 [2.0, 1.0]

Thus, in your first snippet, you only change the contents of the unique 2-element vector. Check your U, it should be filled with [1.0, 1.0].

In your second snippet however, you allocate new vectors in each entry of U, so you don't have this problem.

A solution is to preallocate differently. You could try, e.g.,

U = Array{Vector{Float64}}(undef, (N, N))
user_56021490
Related