Yet another function/loop scope question

General Usage
, ,
1

Thanks for your patience with such a simple question, but I’m stumped. In a function, I want to:

  1. Declare a vector of vectors
  2. Iterate through two loops and update the original vector of vectors
  3. Return the vector of vectors

For some reason, the output of my function is set correctly throughout the loop, but then after the loop ends, these changes aren’t in scope. I don’t understand this next part either - for some reason, the 1st and 2nd row of the output are equal in the output of the function, but they were set to be different (as shown in the print statements in the code below).

How can I fix this scope issue?

Code:

using StaticArrays
using ComponentArrays
import LinearAlgebra: norm

r̅₁ = [0.0, 0.0, 0.0]
v̅₁ = [0.0, 0.0, 0.0]
m₁ = 5.972167867791379e24

r̅₂ = [0.0, 11681000.0, 0.0]
v̅₂ = [5134.0, 4226.0, 2787.0]
m₂ = 150.0

B = SVector{2}([ComponentArray((r̅=r̅₁, v̅=v̅₁)), ComponentArray((r̅=r̅₂, v̅=v̅₂))])
p = ComponentArray((G=6.6743e-11, m=[m₁, m₂]))


function gradient(B, p)
local ∂B = SVector{length(B)}(fill(ComponentArray((r̅=[0.0, 0.0, 0.0], v̅=[0.0, 0.0, 0.0])), size(B)))

for i = 1:length(B)

    ∂B[i].r̅ = B[i].v̅
    for j = 1:length(B)

        if i ≠ j
            ∂B[i].v̅ += ((p.m[i] * p.m[j]) / norm(B[j].r̅ .- B[i].r̅)^3 * (B[j].r̅ .- B[i].r̅))
        end    

    end
    ∂B[i].v̅ *= (p.G / p.m[i])

    println("∂B[", i, "] = ", ∂B[i])

end
println("∂B = ", ∂B)
end

gradient(B,p)

Output:

∂B[1] = (r̅ = [0.0, 0.0, 0.0], v̅ = [0.0, 7.337311123941769e-23, 0.0])
∂B[2] = (r̅ = [5134.0, 4226.0, 2787.0], v̅ = [0.0, -2.921310248692885, 0.0])
∂B = [(r̅ = [5134.0, 4226.0, 2787.0], v̅ = [0.0, -2.921310248692885, 0.0]), (r̅ = [5134.0, 4226.0, 2787.0], v̅ = [0.0, -2.921310248692885, 0.0])]
2

I am not sure I fully understand what you are trying to do, but as written, your function doesn’t return anything (technically it returns the result if println), which is nothing.

If you want to return the ∂B, just put that as the last expression.

Also, you should not need the local in the code above (but it’s harmless).

2 Likes
3

To simplify (for me) I dropped SVector. The problem is fill with your mutable ComponentArray
This works:

julia> function gradient(B, p)
           local ∂B = Vector()
           for i in 1:length(B)
               push!(∂B, ComponentArray((r̅=[0.0, 0.0, 0.0], v̅=[0.0, 0.0, 0.0])))
           end
           for i = 1:length(B)
               ∂B[i].r̅ = B[i].v̅
               for j = 1:length(B)
                   if i ≠ j
                       ∂B[i].v̅ += ((p.m[i] * p.m[j]) / norm(B[j].r̅ .- B[i].r̅)^3 * (B[j].r̅ .- B[i].r̅))
                   end    
               end
               ∂B[i].v̅ *= (p.G / p.m[i])

               println("∂B[", i, "] = ", ∂B[i])
           end
           println("∂B = ", ∂B)
       end
gradient (generic function with 1 method)

julia> gradient(B,p)
∂B[1] = (r̅ = [0.0, 0.0, 0.0], v̅ = [0.0, 7.33731e-23, 0.0])
∂B[2] = (r̅ = [5134.0, 4226.0, 2787.0], v̅ = [0.0, -2.92131, 0.0])
∂B = Any[(r̅ = [0.0, 0.0, 0.0], v̅ = [0.0, 7.33731e-23, 0.0]), (r̅ = [5134.0, 4226.0, 2787.0], v̅ = [0.0, -2.92131, 0.0])]

Example not working because of fill:

julia> A = fill([0,0], 2)
2-element Array{Array{Int64,1},1}:
 [0, 0]
 [0, 0]

julia> A[1][1] = 1
1

julia> A
2-element Array{Array{Int64,1},1}:
 [1, 0]
 [1, 0]
2 Likes
4 
help> fill

   [...]

   If x is an object reference, all elements will refer to the same object:

   julia> A = fill(zeros(2), 2);

   julia> A[1][1] = 42; # modifies both A[1][1] and A[2][1]

   julia> A
   2-element Array{Array{Float64,1},1}:
    [42.0, 0.0]
    [42.0, 0.0]

Note that in particular A[1] === A[2] in this example. The same holds for your original code.

2 Likes
5

Thanks for the quick answer! Yeah, when converting the code to a self-contained version, I accidentally left out the return statement. The local modifier was an attempt to fix my “scope” issue. As @remi-garcia pointed out, the real issue with my code wasn’t scope, it was an incorrect use of fill.

1 Like
6

Thanks for the quick answer & explanation, this fixed things!