How to think about indices and cartesian coordinates?

New to Julia
, , ,
1

Hello everyone!

I am new to Julia and I am writing a simple simulation that uses a cartesian grid and numerical derivatives. I would like to confirm if my thinking is correct.

I start by defining simple grids for X, Y and a scalar field f which in this case is just Y.

n = 3
X = ones(Float64, n) * collect(1:n)'
Y = collect(1:n) * ones(Float64, n)'
f = Y

f is what I would expect:

3×3 Matrix{Float64}:
 1.0  1.0  1.0
 2.0  2.0  2.0
 3.0  3.0  3.0

That is, f[3,1] #=> 3

I also define a simple gradient function:

function ∇(q::Matrix{Float64})
    ∇q_x = zeros(n,n)
    ∇q_y = zeros(n,n)

    for j in 2:n-1, i in 2:n-1
        ∇q_x[i,j] = (q[i,j+1] - q[i,j-1])/2
        ∇q_y[i,j] = (q[i+1,j] - q[i-1,j])/2
    end

    Dict(:x => ∇q_x, :y => ∇q_y)
end

Is it correct to think of the first index (i) as the y direction, and the second index (j) as the x direction? That feels odd to me.
I understand it may just be a matter of convention, but when I try to visualize it, it does appear that Julia prefers it that way.

using Plots; gr()
heatmap(f, c=:blues)
quiver!(X, Y, quiver=(∇(f)[:x], ∇(f)[:y]), c=:white, linewidth=5)

Screen Shot 2023-03-12 at 3.43.40 PM
Screen Shot 2023-03-12 at 3.43.40 PM1444×962 126 KB

Any thoughts are greatly appreciated!

2

I think this is more about how the plotting library does the heatmap (not julia itself). For example, in CairoMakie,

fig,ax = heatmap(ff, colormap=:blues)
arrows!(ax,X[1,:], Y[:,1], ∇(f)[:x], ∇(f)[:y], color=:white, linewidth=5)
fig

Produces this:

imagen
imagen1754×1340 107 KB

There is not a “correct” way for showing this, so that is why some libraries (even outside julia) have different approaches.

For Makie it has been discussed for example here:

heatmap permutes data · Issue #1791 · MakieOrg/Makie.jl · GitHub
Correct axes bug for heatmaps & images by IanButterworth · Pull Request #98 · JuliaPlots/AbstractPlotting.jl · GitHub
X / Y axes swapped in heat map · Issue #295 · MakieOrg/Makie.jl · GitHub
Resolved: Heatmap transposes matrices · Issue #205 · MakieOrg/Makie.jl · GitHub

As for your simulations, as long as you are consistent, there should not be a problem, I think!

(Edit: to show the CairoMakie also with arrows)

3 Likes
3

It is really neither correct nor incorrect, it is really up to you to choose a convention for yourself. Unfortunately.

The problem is that we are used to thinking of x as the first dimension and y as the second dimension, while in linear algebra the first dimension is ‘downwards’ and the second dimension is ‘across’.

Every time, I struggle with choosing a convention, and I personally believe everything would be much nicer if the first index always was the x-direction, and also indicated column (the opposite of today).

So, it’s up to you.

I urge you to not get into the habit of using collect. It is wasteful and pointless here. Forget you ever heard of that function. Remove it from all your codes, and don’t use until your code stops working without it.

5 Likes
4

To elaborate a little on @DNF’s good advice, in case it’s not clear:

julia> n = 3
3

julia> xs = 1:n
1:3

julia> ys = 1:n
1:3

julia> f = [y for y in ys, x in xs] # or any other function of (x,y)
3×3 Matrix{Int64}:
 1  1  1
 2  2  2
 3  3  3

You never need to explicitly form your X or Y matrices. You can either broadcast, or use a multidimensional array comprehension as above. Here’s an example of broadcasting:

julia> fun(x,y) = y
fun (generic function with 1 method)

julia> fun.(xs', ys)
3×3 Matrix{Int64}:
 1  1  1
 2  2  2
 3  3  3
2 Likes
5

Glad I’m not the only person who dislikes this function. I always had an uneasy feeling about it. IMO, most of what collect is used for should be functionality of the Array constructor.

To elaborate a little further: not only is using comprehensions more legible, but it’s also faster:

julia> using BenchmarkTools
       a = @btime ones(Float64, 3) * collect(1:3)'
       b = @btime Float64[j for i=1:3, j=1:3]
       c = @btime ((x,y)->Float64(x)).((1:3)',1:3)
       a == b == c
  60.652 ns (3 allocations: 288 bytes)
  16.917 ns (1 allocation: 128 bytes)
  16.817 ns (1 allocation: 128 bytes)
true
2 Likes
6

Great to know! That’s what I assumed but I wondered about the convention in the Julia community. Thank you for introducing me to CairoMakie.

I agree with you. For now I am inclined to adopt the convention where the first index is the y-coordinate since it is “downwards”, even if it feels a bit odd. In any case, it is great to know that I’m not alone in struggling to choose a convention.

Thank you for the great advice @DNF! And thank you @PeterSimon and @uniment for elaborating and providing examples. That was very helpful for better understand the problem with collect. I was able to shave a significant amount of time from my simulation by eliminating all the calls to collect!

1 Like