def linear_indx(r, x, y, c, i_max, j_max):
    d0, d1, d2, d3 = r.shape
np.put(r, np.arange(i_max)[:, (None)] * d1 * d2 * d3 + np.arange(j_max) *
    d2 * d3 + x * d3 + y, c)
return r

def org_app(r, x, y, c, i_max, j_max):
    for i in range(i_max):
        for j in range(j_max):
            r[i, j, x[i, j], y[i, j]] = c[i, j]
return r
Оптимизация операции индексирования numpy ndarray