Slicing e broadcasting di array multidimensionali in Julia: esempio meshgrid

Question

Recentemente ho iniziato a studiare Julia codificando una semplice implementazione di Self Organizing Maps. Voglio che le dimensioni e le dimensioni della mappa siano specificate dall'utente, il che significa che non posso davvero usare i loop per lavorare sugli array di mappe perché non so in anticipo quanti strati di loop avrò bisogno. Quindi ho assolutamente bisogno di funzioni di broadcasting e slicing che funzionano su array di dimensioni arbitrarie.

Al momento, ho bisogno di costruire una serie di indici della mappa. La mia mappa è definita da una serie di dimensioni mapsize = (5, 10, 15), ho bisogno di costruire un array indices della dimensione (3, 5, 10, 15) dove indices[:, a, b, c] deve restituire [a, b, c].

io vengo da un background Python/NumPy, in cui la soluzione è già data da una "funzione" specifica, MGrid:

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

non mi aspettavo Julia di avere una tale funzione sul get -solo, così ho iniziato a trasmettere. In NumPy, la trasmissione è basata su un insieme di regole che trovo abbastanza chiare e logiche. È possibile utilizzare array di dimensioni diverse nella stessa espressione fino a quando le dimensioni in ogni partita di quota o una parte di esso è 1:

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

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

per aiutare con questo, è anche possibile utilizzare numpy.newaxis o None per facilmente aggiungere nuove dimensioni per l'array:

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

questo aiuta le matrici di trasmissione facilmente:

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) 

Usando questo, creando la matrice indices è piuttosto straightfo rward:

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

Il caso generale è, naturalmente, un po 'più complicato, ma può essere lavorato in giro con la lista di comprensione e fette:

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)) 

Ma torniamo a Julia. Ho provato ad applicare lo stesso tipo di motivazione e ho finito per ottenere l'equivalente della lista arrays del codice sopra. Questo ha finito per essere piuttosto semplice rispetto alle NumPy controparte grazie alla sintassi delle espressioni composto:

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

Ora mi sono bloccato qui, come io non so davvero come applicare la trasmissione alla mia lista di generare matrici qui ... Le funzioni broadcast[!] richiedono una funzione f da applicare e io non ne ho. Ho provato ad utilizzare un ciclo for per provare a forzare il broadcasting:

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

Ma questo mi dà un errore: ERROR: convert has no method matching convert(::Type{Int64}, ::Array{Int64,3})

sto facendo questo il modo giusto? Ho trascurato qualcosa di importante? Qualsiasi aiuto è apprezzato.

Answer 1

La trasmissione in Julia è stata modellata praticamente sulla trasmissione in NumPy, quindi dovresti scoprire che obbedisce più o meno alle stesse semplici regole (non sono sicuro se il modo di riempire le dimensioni quando non tutti gli input hanno lo stesso numero di le dimensioni sono le stesse però, dal momento che gli array Julia sono colonne maggiori).

Un numero di cose utili come l'indicizzazione newaxis e broadcast_arrays non sono state ancora implementate (ancora). (Spero lo faranno.) Si noti inoltre che l'indicizzazione funziona in modo leggermente diverso in Julia rispetto a NumPy: quando si omettono gli indici per le dimensioni finali in NumPy, gli indici rimanenti vengono impostati automaticamente su due punti. In Julia si potrebbe dire che per impostazione predefinita sono invece quelli.

Non sono sicuro che se si ha effettivamente bisogno di una funzione meshgrid, la maggior parte delle cose che si vorrebbe usarlo potrebbe essere eseguita utilizzando le voci originali dell'array arrays con le operazioni di trasmissione. Il motivo principale per cui meshgrid è utile in MATLAB è perché è terribile alla trasmissione.

Ma è abbastanza semplice da realizzare ciò che si vuole fare usando la funzione broadcast!:

# assume mapsize is a vector with the desired shape, e.g. 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 

ho dovuto aggiungere una dimensione iniziale di Singleton sulle voci di arrays per allinearsi con gli assi di indices (newaxis era stato utile qui ...). Quindi passo attraverso ogni coordinata, creo un sottotitolo (una vista) sulla parte rilevante di indices e lo compilo. (L'indicizzazione sarà predefinita per restituire i sottotitoli in Julia 0.4, ma per ora dobbiamo usare esplicitamente lo sub).

La chiamata a broadcast! valuta la funzione di identità identity(x)=x sull'ingresso arr=arrays[i], trasmette la forma dell'output. Non c'è perdita di efficienza nell'uso della funzione identity per questo; broadcast! genera una funzione specializzata basata sulla funzione data, numero di argomenti e numero di dimensioni del risultato.

Answer 2

Credo che questo è lo stesso della funzionalità MATLAB meshgrid. Non ho mai veramente pensato alla generalizzazione a più di due dimensioni, quindi è un po 'più difficile capirlo.

In primo luogo, ecco la mia versione del tutto generale, che è un pò pazzo, ma non riesco a pensare ad un modo migliore per farlo senza generare codice per dimensioni comuni (ad esempio, 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 

Ora , ciò che intendo per la generazione del codice è che, per il caso 2D, si può semplicemente fare

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 

e per il caso 3D:

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 con, ad es.

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

immagino mgrid infila tutti in un unico tensore, così si potrebbe fare che in questo modo

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

che è ancora un po 'diversa:

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 

Answer 3

Se stai usando julia 0.4, è possibile effettuare questa operazione:

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

Sarebbe ancora più bello se si potesse solo dire

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

Guarderò in quella :-).