Slicing and translating multidimensional arrays in Julia: an example of meshgrid

Question

Recently, I began to study Julia, coding a simple implementation of Self Organizing Maps. I want the size and dimensions of the map to be specified by the user, which means that I cannot use loops to work with map arrays, because I do not know in advance how many layers of loops I will need. Therefore, I absolutely need the broadcast and slicing functions that work on arrays of arbitrary sizes.

Now I need to build an array of map indices. Let's say my map is defined by an array of size mapsize = (5, 10, 15) , I need to build an array of indices size (3, 5, 10, 15) , where indices[:, a, b, c] should return [a, b, c] .

I come from the Python / NumPy background, in which the solution is already defined by a specific "function", mgrid:

 indices = numpy.mgrid[:5, :10, :15] print indices.shape # gives (3, 5, 10, 15) print indices[:, 1, 2, 3] gives [1, 2, 3] 

I did not expect that Julia would have such a function on the go, so I turned to broadcasting. In NumPy, translation is based on a set of rules that I find quite clear and logical. You can use arrays of different sizes in one expression if the sizes in each dimension are the same or one of them is 1:

 (5, 10, 15) broadcasts to (5, 10, 15) (10, 1) (5, 1, 15) also broadcasts to (5, 10, 15) (1, 10, 1) 

To help with this, you can also use numpy.newaxis or None to easily add new dimensions to your array:

 array = numpy.zeros((5, 15)) array[:,None,:] has shape (5, 1, 15) 

This makes it easy to cast arrays:

 A = numpy.arange(5) B = numpy.arange(10) C = numpy.arange(15) bA, bB, bC = numpy.broadcast_arrays(A[:,None,None], B[None,:,None], C[None,None,:]) bA.shape == bB.shape == bC.shape = (5, 10, 15) 

Using this, creating an indices array is quite simple:

 indices = numpy.array(numpy.broadcast_arrays(A[:,None,None], B[None,:,None], C[None,None,:])) (indices == numpy.mgrid[:5,:10,:15]).all() returns True 

The general case, of course, is a little more complicated, but it can be used using an understanding of lists and fragments:

 arrays = [ numpy.arange(i)[tuple([None if m!=n else slice(None) for m in range(len(mapsize))])] for n, i in enumerate(mapsize) ] indices = numpy.array(numpy.broadcast_arrays(*arrays)) 

So back to Julia. I tried to apply the same rationale and ended up with the equivalent arrays list of the above code. This turned out to be simpler than the NumPy counterpart due to the syntax of the compound expression:

 arrays = [ (idx = ones(Int, length(mapsize)); idx[n] = i;reshape([1:i], tuple(idx...))) for (n,i)=enumerate(mapsize) ] 

Now I'm stuck here, since I really don’t know how to apply translation in my list of generator arrays here ... The broadcast[!] Functions require the use of the f function, and I don’t have any. I tried using the for loop to try to force broadcasting:

 indices = Array(Int, tuple(unshift!([i for i=mapsize], length(mapsize))...)) for i=1:length(mapsize) A[i] = arrays[i] end 

But this gives me an error: ERROR: convert has no method matching convert(::Type{Int64}, ::Array{Int64,3})

Am I doing it right? Am I missing something important? Any help is appreciated.

Answer 1

The broadcast in Julia was modeled to a large extent for broadcasting in NumPy, so you hopefully find that it obeys more or less the same simple rules (not sure if the path to the dimensions of the strip, when not all inputs have the same number of measurements, same as Julia arrays are columns).

However, a number of useful functions, such as newaxis indexing and broadcast_arrays , are not yet implemented. (I hope they do.) Also note that indexing works differently in Julia compared to NumPy: when you leave indexes for lagging dimensions in NumPy, the rest of the indexes are colon-separated. In Julia, they could be used by default instead of them.

I'm not sure if you really need the meshgrid function, most of the things you would like to use for it could be done using the source records of your arrays array with broadcast operations. The main reason meshgrid is useful in Matlab is because it is terrible at translation.

But it's pretty simple to do what you want to do using the broadcast! function broadcast! :

 # assume mapsize is a vector with the desired shape, eg mapsize = [2,3,4] N = length(mapsize) # Your line to create arrays below, with an extra initial dimension on each array arrays = [ (idx = ones(Int, N+1); idx[n+1] = i;reshape([1:i], tuple(idx...))) for (n,i) in enumerate(mapsize) ] # Create indices and fill it one coordinate at a time indices = zeros(Int, tuple(N, mapsize...)) for (i,arr) in enumerate(arrays) dest = sub(indices, i, [Colon() for j=1:N]...) broadcast!(identity, dest, arr) end 

I had to add an initial singleton dimension in the arrays record to align with the indices axes ( newaxis was useful here ...). Then I look at each coordinate, create a subarray (view) in the corresponding part of indices and fill it. (Indexing will by default return subarrays in Julia 0.4, but for now we must explicitly use sub ).

Call broadcast! it simply evaluates the identity function identity(x)=x at the input arr=arrays[i] , is translated into the output form. There is no efficiency in using the identity function for this; broadcast! generates a specialized function based on this function, the number of arguments and the number of dimensions of the result.

Answer 2

If you are using julia 0.4 you can do this:

 julia> function mgrid(mapsize) T = typeof(CartesianIndex(mapsize)) indices = Array(T, mapsize) for I in eachindex(indices) indices[I] = I end indices end 

It would be even better if you could just say

 indices = [I for I in CartesianRange(CartesianIndex(mapsize))] 

I will consider this :-).

Answer 3

I assume this is the same as the MATLAB meshgrid functionality. I never thought about generalizing in more than two dimensions, so it's a little harder to circle your head.

Firstly, here is my fully general version, which is a little crazy, but I can't think of a better way to do this without writing code for general measurements (e.g. 2, 3)

 function numpy_mgridN(dims...) X = Any[zeros(Int,dims...) for d in 1:length(dims)] for d in 1:length(dims) base_idx = Any[1:nd for nd in dims] for i in 1:dims[d] cur_idx = copy(base_idx) cur_idx[d] = i X[d][cur_idx...] = i end end @show X end X = numpy_mgridN(3,4,5) @show X[1][1,2,3] # 1 @show X[2][1,2,3] # 2 @show X[3][1,2,3] # 3 

Now, what I mean by code generation is that for a 2D case, you can just do

 function numpy_mgrid(dim1,dim2) X = [i for i in 1:dim1, j in 1:dim2] Y = [j for i in 1:dim1, j in 1:dim2] return X,Y end 

and for the 3D case:

 function numpy_mgrid(dim1,dim2,dim3) X = [i for i in 1:dim1, j in 1:dim2, k in 1:dim3] Y = [j for i in 1:dim1, j in 1:dim2, k in 1:dim3] Z = [k for i in 1:dim1, j in 1:dim2, k in 1:dim3] return X,Y,Z end 

Test with for example

 X,Y,Z=numpy_mgrid(3,4,5) @show X @show Y @show Z 

I think mgrid drags them all into one tensor, so you can do it like this

 all = cat(4,X,Y,Z) 

which is still a little different:

 julia> all[1,2,3,:] 1x1x1x3 Array{Int64,4}: [:, :, 1, 1] = 1 [:, :, 1, 2] = 2 [:, :, 1, 3] = 3 julia> vec(all[1,2,3,:]) 3-element Array{Int64,1}: 1 2 3