Krojenie i rozgłaszanie wielowymiarowych tablic w Julii: meshgrid example

Question

Niedawno zacząłem uczyć się Julii, kodując prostą implementację samoorganizujących się map. Chcę, aby rozmiar i wymiary mapy były określone przez użytkownika, co oznacza, że ​​naprawdę nie mogę używać pętli do pracy na macierzach map, ponieważ nie wiem z góry, ile warstw pętli będę potrzebować. Tak więc absolutnie potrzebuję funkcji rozgłaszania i cięcia, które działają na tablicach o dowolnych wymiarach.

Teraz muszę zbudować tablicę wskaźników mapy. Powiedzmy, że moja mapa jest zdefiniowana przez tablicę o rozmiarze mapsize = (5, 10, 15), potrzebuję skonstruować tablicę indices o rozmiarze (3, 5, 10, 15), gdzie indices[:, a, b, c] powinien zwrócić .

Pochodzę z Python/Numpy tle, w których rozwiązanie jest już podany przez konkretnego „funkcje”, mgrid:

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

Nie spodziewałem Julia mieć taką funkcję get -go, więc zwróciłem się do nadawania. W NumPy emisja jest oparta na zbiorze reguł, które uważam za dość jasne i logiczne. Można używać tablic o różnych wymiarach w tym samym wyrażeniu tak długo, jak rozmiary w każdym meczu wymiar lub jeden z nich jest 1:

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

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

Aby w tym pomóc, można również użyć numpy.newaxis albo żaden z łatwością dodaj nowy wymiar do swojej tablicy:

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

pomaga tablice broadcast prościej:

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) 

Używanie tego, tworząc tablicę indices jest raczej 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 

ogólnym przypadku jest oczywiście nieco bardziej skomplikowane, ale można obejść stosując listowych i plasterki:

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

Wracając do Julii. Próbowałem zastosować ten sam rodzaj uzasadnienia i ostatecznie osiągnąłem ekwiwalent listy kodów powyżej arrays. To okazało się być raczej prostsze niż NumPy CPF dzięki składni wyrażeń złożonych:

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

Teraz utknąłem tutaj, a ja naprawdę nie wiem, jak zastosować nadawanie na mojej liście generowania tablic tutaj ... Funkcje broadcast[!] wymagają podania funkcji f, a ja nie mam żadnych. Próbowałem za pomocą pętli for spróbować zmuszając nadawanie:

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

Ale to daje mi błąd: ERROR: convert has no method matching convert(::Type{Int64}, ::Array{Int64,3})

robię to we właściwy sposób? Czy przeoczyłem coś ważnego? Każda pomoc jest doceniana.

Answer 1

Radiofonia w Julii została prawie w całości zaprojektowana na nadawanie w NumPy, więc powinieneś mieć nadzieję, że będzie ona przestrzegać mniej więcej tych samych prostych reguł (nie wiesz, czy sposób na dopełnienie wymiarów, gdy nie wszystkie wejścia mają taką samą liczbę wymiary są takie same, ponieważ tablice Julii mają kolumnę główną).

Jednak wiele przydatnych rzeczy, takich jak indeksowanie newaxis i broadcast_arrays, nie zostało jeszcze zaimplementowanych (jeszcze). (Mam nadzieję, że to zrobią). Zwróć też uwagę, że indeksowanie działa nieco inaczej w Julii w porównaniu do NumPy: kiedy zostawiasz indeksy dla końcowych wymiarów w NumPy, pozostałe indeksy domyślnie przyjmują dwukropki. W Julii można powiedzieć, że domyślnie zamiast nich.

Nie jestem pewien, czy faktycznie potrzebujesz funkcji siatki, większość rzeczy, które chciałbyś jej użyć, może być wykonane przy użyciu oryginalnych wpisów z tablicy arrays z operacjami rozgłaszania. Głównym powodem, dla którego meshgrid jest przydatny w programie Matlab, jest to, że jest straszny w przypadku nadawania.

Ale to jest dość łatwe do osiągnięcia tego, co chcesz zrobić za pomocą broadcast! funkcję:

# 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 

musiałem dodać początkowy wymiar singleton na wpisów arrays linii z osiami indices (newaxis było przydatne tutaj ...). Następnie przechodzę przez każdą współrzędną, tworząc podtablicę (widok) na odpowiedniej części indices i wypełniając ją. (Indeksowanie domyślnie powróci do podrzędnych podarunków w Julia 0.4, ale na razie musimy wyraźnie użyć sub).

Wezwanie do broadcast! tylko ocenia funkcję tożsamości identity(x)=x na wejściu arr=arrays[i], transmisje do kształtu wyjściowego. Nie ma utraty wydajności w korzystaniu z funkcji identity; broadcast! generuje wyspecjalizowaną funkcję opartą na podanej funkcji, liczbie argumentów i liczbie wymiarów wyniku.

Answer 2

Domyślam się, że to jest to samo, co funkcja MATLAB meshgrid. Nigdy tak naprawdę nie myślałem o uogólnieniu do więcej niż dwóch wymiarów, więc nieco trudniej mi się rozejrzeć.

pierwsze, tutaj jest moja wersja całkowicie ogólnego, który jest trochę szalony, ale nie mogę myśleć o lepszym sposobem, aby to zrobić bez generowania kodu dla typowych wymiarach (np 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 

teraz , co mam na myśli generowania kodu jest to, że dla przypadku 2D, można po prostu zrobić

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 

i dla przypadku 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 

Przetestuj za pomocą np.

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

Chyba mgrid je wszystkie do jednego pchnięcia tensora, więc można zrobić jak ten

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

która nadal jest nieco inna:

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

Jeśli używasz julia 0.4, możesz to zrobić:

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

To byłoby jeszcze ładniej, jeśli można tak powiedzieć

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

Zajrzę do tego :-).