Сложная проблема нарезки на Python / numpy
https://stackoverflow.com/questions/1587367
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
У меня проблема с некоторыми дурацкими вещами.Мне нужен массив numpy, который вел бы себя необычным образом, возвращая фрагмент как представление данных, которые я вырезал, а не копию.Итак, вот пример того, что я хочу сделать:
Допустим, у нас есть простой массив, подобный этому:
a = array([1, 0, 0, 0])
Я хотел бы обновить последовательные записи в массиве (перемещаясь слева направо) предыдущей записью из массива, используя синтаксис, подобный этому:
a[1:] = a[0:3]
Это привело бы к следующему результату:
a = array([1, 1, 1, 1])
Или что-то вроде этого:
a[1:] = 2*a[:3]
# a = [1,2,4,8]
Чтобы проиллюстрировать далее, я хочу, чтобы поведение было следующего типа:
for i in range(len(a)):
if i == 0 or i+1 == len(a): continue
a[i+1] = a[i]
За исключением того, что мне нужна скорость numpy.
Поведение numpy по умолчанию заключается в том, чтобы скопировать фрагмент, так что на самом деле я получаю следующее:
a = array([1, 1, 0, 0])
У меня уже есть этот массив как подкласс ndarray, поэтому я могу вносить в него дополнительные изменения, если потребуется, мне просто нужно, чтобы фрагмент с правой стороны постоянно обновлялся по мере обновления фрагмента с левой стороны.
Мне снится или это волшебство возможно?
Обновить:Это все потому, что я пытаюсь использовать итерацию Гаусса-Зайделя для решения задачи линейной алгебры, более или менее.Это особый случай, связанный с гармоническими функциями, я пытался не вдаваться в подробности, потому что в этом действительно нет необходимости и, вероятно, это еще больше запутает ситуацию, но вот что получилось.
Алгоритм таков:
while not converged:
for i in range(len(u[:,0])):
for j in range(len(u[0,:])):
# skip over boundary entries, i,j == 0 or len(u)
u[i,j] = 0.25*(u[i-1,j] + u[i+1,j] + u[i, j-1] + u[i,j+1])
Верно?Но вы можете сделать это двумя способами, Jacobi включает обновление каждого элемента с его соседями без учета обновлений, которые вы уже сделали, до тех пор, пока цикл while не зациклится, чтобы сделать это в циклах, вы скопировали бы массив, а затем обновили один массив из скопированного массива.Однако Gauss-Seidel использует информацию, которую вы уже обновили для каждой из записей i-1 и j-1, поэтому нет необходимости в копировании, цикл должен по существу "знать", поскольку массив был переоценен после каждого обновления отдельного элемента.То есть, каждый раз, когда мы вызываем запись типа u [i-1,j] или u[i, j-1], там будет информация, вычисленная в предыдущем цикле.
Я хочу заменить эту медленную и уродливую ситуацию с вложенным циклом одной хорошей чистой строкой кода, используя numpy slicing:
u[1:-1,1:-1] = 0.25(u[:-2,1:-1] + u[2:,1:-1] + u[1:-1,:-2] + u[1:-1,2:])
Но результатом является итерация Якоби, потому что когда вы берете фрагмент:u[:,-2,1:-1] вы копируете данные, таким образом, фрагмент не знает о каких-либо сделанных обновлениях.Теперь numpy все еще зацикливается, верно?Это не параллельный, это просто более быстрый способ выполнения цикла, который выглядит как параллельная операция в python.Я хочу использовать это поведение, своего рода взломав numpy, чтобы возвращать указатель вместо копии, когда я беру фрагмент.Верно?Затем каждый раз, когда numpy выполняет цикл, этот фрагмент "обновляется" или на самом деле просто реплицирует все, что произошло в обновлении.Чтобы сделать это, мне нужны фрагменты с обеих сторон массива в качестве указателей.
В любом случае, если есть какой-то действительно очень умный человек, который настолько крут, но я в значительной степени смирился с мыслью, что единственный ответ - это цикл на C.
Решение
Поздний ответ, но это появилось в Google, поэтому я, вероятно, укажу на документ, который хотел ОП.Ваша проблема ясна:при использовании срезов NumPy создаются временные объекты.Оберните свой код в быстрый вызов weave.blitz, чтобы избавиться от временных объектов и добиться желаемого поведения.
Прочтите раздел weave.blitz на сайте Руководство по производительностиPython для получения полной информации.
Другие советы
накапливать предназначен для того, чтобы делать то, что вы хотите;то есть, чтобы распространить операцию на массив.Вот пример:
from numpy import *
a = array([1,0,0,0])
a[1:] = add.accumulate(a[0:3])
# a = [1, 1, 1, 1]
b = array([1,1,1,1])
b[1:] = multiply.accumulate(2*b[0:3])
# b = [1 2 4 8]
Другой способ сделать это — явно указать массив результатов в качестве входного массива.Вот пример:
c = array([2,0,0,0])
multiply(c[:3], c[:3], c[1:])
# c = [ 2 4 16 256]
Просто используйте цикл.Я не могу сразу придумать какой-либо способ заставить оператор среза вести себя так, как вы говорите, что хотите, за исключением может быть путем создания подклассов numpy array и переопределить соответствующий метод с помощью какого-то вуду Python...но что более важно, идея о том, что a[1:] = a[0:3] должен скопировать первое значение a в следующие три слота мне кажется совершенно бессмысленным.Я полагаю, что это может легко сбить с толку любого, кто посмотрит на ваш код (по крайней мере, первые несколько раз).
Это не правильная логика.Я попробую объяснить это с помощью букв.
Изображение array = abcd с элементами a,b,c,d.
Сейчас, array[1:] означает от элемента в позиции 1 (начиная с 0) на.
В этом случае:bcd и array[0:3] значит от персонажа в позиции 0 до третьего символа (того, что в позиции 3-1) в этом случае: 'abc'.
Пишу что-то вроде:
array[1:] = array[0:3]
означает:заменять bcd с abc
Чтобы получить желаемый результат, теперь в Python вы должны использовать что-то вроде:
a[1:] = a[0]
Должно быть, это как-то связано с назначением среза.Однако операторы, как вы, возможно, уже знаете, следуют вашему ожидаемому поведению:
>>> a = numpy.array([1,0,0,0])
>>> a[1:]+=a[:3]
>>> a
array([1, 1, 1, 1])
Если в вашей реальной задаче, где есть ваш пример, уже есть нули, то это решает ее.В противном случае, за дополнительную плату, установите их на ноль либо путем умножения на ноль, либо путем присвоения нуля (в зависимости от того, что быстрее).
редактировать:У меня была еще одна мысль.Вы можете предпочесть это:
numpy.put(a,[1,2,3],a[:3])
Numpy должен проверять, совпадает ли целевой массив с входным массивом при выполнении вызова setkey.К счастью, есть способы обойти это.Сначала я попробовал использовать numpy.put вместо
In [46]: a = numpy.array([1,0,0,0])
In [47]: numpy.put(a,[1,2,3],a[0:3])
In [48]: a
Out[48]: array([1, 1, 1, 1])
А затем, судя по документации, я попробовал использовать Flatiters (a.flat)
In [49]: a = numpy.array([1,0,0,0])
In [50]: a.flat[1:] = a[0:3]
In [51]: a
Out[51]: array([1, 1, 1, 1])
Но это не решает проблему, которую вы имели в виду
In [55]: a = np.array([1,0,0,0])
In [56]: a.flat[1:] = 2*a[0:3]
In [57]: a
Out[57]: array([1, 2, 0, 0])
Это не удается, поскольку умножение выполняется до присваивания, а не параллельно, как хотелось бы.
Numpy предназначен для многократного применения одной и той же операции параллельно в массиве.Чтобы сделать что-то более сложное, если только вы не сможете найти разложение его на такие функции, как numpy.cumsum и numpy.cumprod, вам придется прибегнуть к чему-то вроде scipy.weave или написать функцию на C.(См. ПроизводительностьPython страницу для более подробной информации.) (Кроме того, я никогда не использовал плетение, поэтому не могу гарантировать, что оно сделает то, что вы хотите.)
Вы могли бы взглянуть на np.lib.stride_tricks .
На этих замечательных слайдах есть кое-какая информация:http://mentat.za.net/numpy/numpy_advanced_slides/
с помощью stride_tricks, начинающихся со слайда 29.
Я не совсем понимаю этот вопрос, поэтому не могу предложить ничего более конкретного - хотя я бы, вероятно, сделал это на cython или fortran с помощью f2py или с помощью weave.На данный момент мне больше нравится fortran, потому что к тому времени, когда вы добавляете все необходимые аннотации типов в cython, я думаю, что в конечном итоге он выглядит менее понятным, чем fortran.
Здесь приведено сравнение этих подходов:
www.сципи.организация/ PerformancePython
(не могу опубликовать больше ссылок, так как я новый пользователь) с примером, который похож на ваш случай.
В конце концов я столкнулся с той же проблемой, что и вы.Мне пришлось прибегнуть к использованию итерации Якоби и ткача:
while (iter_n < max_time_steps):
expr = "field[1:-1, 1:-1] = (field[2:, 1:-1] "\
"+ field[:-2, 1:-1]+"\
"field[1:-1, 2:] +"\
"field[1:-1, :-2] )/4."
weave.blitz(expr, check_size=0)
#Toroidal conditions
field[:,0] = field[:,self.flow.n_x - 2]
field[:,self.flow.n_x -1] = field[:,1]
iter_n = iter_n + 1
Он работает и работает быстро, но это не метод Гаусса-Зейделя, поэтому сходимость может быть немного сложной.Единственный вариант выполнения Гаусса-Зейделя в виде традиционного цикла с индексами.
я бы предложил Cython вместо цикла в c.там мощь быть каким-то причудливым способом заставить ваш пример работать, используя множество промежуточных шагов...но поскольку вы уже знаете, как писать это на языке C, просто напишите этот небольшой фрагмент как функцию cython, и позвольте магии cython облегчить вам остальную работу.
