Нарезка и трансляция многомерных массивов в Джулии: пример meshgrid

Question

Недавно я начал изучать Юлию, кодируя простую реализацию самоорганизующихся карт. Я хочу, чтобы размер и размеры карты указывались пользователем, а это значит, что я не могу использовать циклы для работы с массивами карт, потому что я не знаю заранее, сколько слоев циклов мне понадобится. Поэтому мне абсолютно необходимы функции широковещания и нарезки, которые работают на массивах произвольных измерений.

Прямо сейчас мне нужно построить массив индексов карты. Скажем, моя карта определена массивом размера mapsize = (5, 10, 15), мне нужно построить массив indices размера (3, 5, 10, 15), где indices[:, a, b, c] должен вернуть [a, b, c].

Я родом из фона Python/NumPy, в котором решение уже задается определенной «функции», MGRID:

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

Я не ожидал, что Джулия иметь такую ​​функцию на ГЭТ -go, поэтому я обратился к вещанию. В NumPy трансляция основана на наборе правил, которые я нахожу вполне ясными и логичными. Вы можете использовать массивы разных размеров в одном выражении, пока размеры в каждом матче измерения или один из него является 1:

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

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

Чтобы помочь с этим, вы можете также использовать numpy.newaxis или Нет, чтобы легко добавлять новые размеры в массив:

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

Это помогает широковещательные массивы легко:

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) 

Используя это, создавая indices массив довольно 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 

Общий случай, конечно, немного сложнее, но можно обойти, используя список понимание и ломтиков:

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

Итак, вернемся к Джулии. Я попытался применить такое же обоснование и в итоге получил эквивалент списка arrays приведенного выше кода. Это закончилось тем, что было довольно проще, чем на партнерской благодаря NumPy выражению соединения синтаксиса:

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

Теперь я застрял здесь, как я не знаю, как применить вещание в мой список генерации массивов здесь ... Функции broadcast[!] требуют применения функции f, и у меня ее нет. Я попытался с помощью для цикла, чтобы попытаться заставить вещания:

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

Но это дает мне ошибку: ERROR: convert has no method matching convert(::Type{Int64}, ::Array{Int64,3})

я делаю это правильный путь? Упустил ли я что-то важное? Любая помощь приветствуется.

Answer 1

Радиовещание в Юлии было смоделировано в значительной степени для вещания в NumPy, поэтому вы должны надеяться, что оно будет подчиняться более или менее таким же простым правилам (не уверен, что путь к размерам прокладки, когда не все входы имеют одинаковое количество размеры одинаковы, так как массивы Julia являются столбцами).

Однако ряд полезных вещей, таких как newaxis, индексирование и broadcast_arrays еще не реализованы. (Надеюсь, они это сделают.) Также обратите внимание, что индексация работает по-иному в Джулии по сравнению с NumPy: когда вы оставляете индексы для отстающих измерений в NumPy, остальные индексы делят на двоеточия. В Julia их можно было использовать вместо них вместо них.

Я не уверен, действительно ли вам нужна функция meshgrid, большинство вещей, которые вы хотели бы использовать для нее, можно было бы сделать, используя оригинальные записи вашего массива arrays с вещательными операциями. Основная причина, по которой meshgrid полезна в Matlab, заключается в том, что это ужасно при трансляции.

Но это довольно просто сделать то, что вы хотите сделать с помощью 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 

мне пришлось добавить начальное одноплодное измерение на записях arrays выстраиваться с осями indices (newaxis был полезен здесь ...). Затем я просматриваю каждую координату, создаю подмассиву (вид) в соответствующей части indices и заполняю ее. (Индексирование будет по умолчанию возвращать подмассивы в Julia 0.4, но пока мы должны явно использовать sub).

Звонок broadcast! просто оценивает функцию тождества identity(x)=x на входе arr=arrays[i], передает в форму вывода. Для этого нет никакой эффективности в использовании функции identity; broadcast! создает специализированную функцию, основанную на заданной функции, количестве аргументов и количестве измерений результата.

Answer 2

Я думаю, что это то же самое, что и функциональность MATLAB meshgrid. Я никогда не думал об обобщении более чем в двух измерениях, поэтому его немного сложнее обвести голову.

Во-первых, здесь моя совершенно общая версия, которая является своего рода сумасшедшим, но я не могу придумать лучшего способа сделать это без генерации кода для стандартных размеров (например, 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 

сейчас то, что я имею в виду генерации кода является то, что для 2D случая, вы можете просто сделать

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 

и для 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 

Испытание, например.

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

Я думаю mgrid засовывает их всех в один тензор, так что вы могли бы сделать это как этот

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

который еще немного отличается:

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

Если вы работаете Джулию 0,4, вы можете сделать это:

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

Было бы еще лучше, если бы можно было бы просто сказать

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

я смотреть на это :-).