Итерация по массивам
Объект итератора nditer, представленный в NumPy 1.6, предоставляет множество гибких способов посещения всех элементов одного или нескольких массивов систематическим образом. Эта страница знакомит с некоторыми базовыми способами использования объекта для вычислений над массивами в Python, а затем завершается описанием того, как можно ускорить внутренний цикл на Cython. Поскольку Python-интерфейс nditer является относительно простым отображением API итератора C-массивов, эти идеи также окажутся полезными при работе с итерацией по массивам из C или C++.
Итерация по одному массиву
Самая базовая задача, которую можно выполнить с помощью nditer, — посетить каждый элемент массива. Каждый элемент предоставляется по одному с использованием стандартного Python-интерфейса итератора.
Пример
>>> a = np.arange(6).reshape(2,3) >>> for x in np.nditer(a): ... print(x, end=' ') ... 0 1 2 3 4 5
Важным моментом при этой итерации является то, что порядок выбран для соответствия расположению элементов в памяти массива, а не с использованием стандартного порядка C или Fortran. Это делается для повышения эффективности доступа, отражая идею, что по умолчанию необходимо просто посетить каждый элемент без учета определенного порядка. Мы можем увидеть это, выполняя итерацию по транспонированному массиву, в сравнении с копированием этой транспонированной матрицы в порядке C.
Пример
>>> a = np.arange(6).reshape(2,3) >>> for x in np.nditer(a.T): ... print(x, end=' ') ... 0 1 2 3 4 5
>>> for x in np.nditer(a.T.copy(order='C')): ... print(x, end=' ') ... 0 3 1 4 2 5
Элементы как a , так и a.T проходятся в том же порядке, а именно в том порядке, в котором они хранятся в памяти, тогда как элементы a.T.copy(order=’C’) проходятся в другом порядке, поскольку они были помещены в другое расположение в памяти.
Управление порядком итерации
Иногда важно посетить элементы массива в определенном порядке, независимо от расположения элементов в памяти. Объект nditer предоставляет параметр order для управления этим аспектом итерации. Значение по умолчанию, имеющее описанное выше поведение, — order=’K’ (сохранить существующий порядок). Это может быть переопределено с order=’C’ для порядка C и order=’F’ для порядка Fortran.
Пример
>>> a = np.arange(6).reshape(2,3) >>> for x in np.nditer(a, order='F'): ... print(x, end=' ') ... 0 3 1 4 2 5 >>> for x in np.nditer(a.T, order='C'): ... print(x, end=' ') ... 0 3 1 4 2 5
Изменение значений массивов
По умолчанию nditer обрабатывает входной операнд как объект только для чтения. Для возможности изменения элементов массива необходимо указать режим чтения-записи или записи только используя флаги ‘readwrite’ или ‘writeonly’ для каждого операнда.
В таком случае nditer вернёт массивы буферов, которые можно модифицировать. Однако, поскольку nditer должен скопировать данные этого буфера обратно в исходный массив после завершения итерации, необходимо сигнализировать о завершении итерации одним из двух способов. Вы можете:
- использовать nditer как менеджер контекста, используя оператор
with, и временные данные будут записаны обратно при выходе из контекста; - вызвать метод
closeитератора после завершения итерации, что запустит запись обратно.
После вызова close или выхода из контекста nditer уже нельзя итерировать.
Пример
>>> a = np.arange(6).reshape(2,3)
>>> a
array([[0, 1, 2],
[3, 4, 5]])
>>> with np.nditer(a, op_flags=['readwrite']) as it:
... for x in it:
... x[...] = 2 * x
...
>>> a
array([[ 0, 2, 4],
[ 6, 8, 10]])
Использование внешнего цикла
Во всех предыдущих примерах элементы a предоставлялись итератором по одному за раз, потому что вся логика цикла находится внутри самого итератора. Хотя это просто и удобно, это не очень эффективно. Лучший подход — перенести одномерный внутренний цикл в свой код, снаружи итератора. Таким образом, векторизованные операции NumPy можно применять к более крупным кускам посещаемых элементов.
nditer попытается предоставить блоки, максимально возможные для внутреннего цикла. Принудительно задавая порядок ‘C’ и ‘F’, мы получаем различные размеры внешнего цикла. Этот режим включён путём указания флага итератора.
Обратите внимание, что при использовании по умолчанию, сохраняющем порядок в памяти, итератор может предоставить один одномерный блок, в то время как при принудительном использовании порядка Fortran, он должен предоставить три блока по два элемента каждый.
Пример
>>> a = np.arange(6).reshape(2,3) >>> for x in np.nditer(a, flags=['external_loop']): ... print(x, end=' ') ... [0 1 2 3 4 5]
>>> for x in np.nditer(a, flags=['external_loop'], order='F'): ... print(x, end=' ') ... [0 3] [1 4] [2 5]
Отслеживание индекса или многомерного индекса
Во время итерации вы можете захотеть использовать индекс текущего элемента в вычислении. Например, вы можете захотеть посетить элементы массива в порядке памяти, но использовать индекс в порядке C, Fortran или многомерный индекс для поиска значений в другом массиве.
Протокол Python-итератора не имеет естественного способа запросить эти дополнительные значения из итератора, поэтому мы вводим альтернативный синтаксис для итерации с nditer. Этот синтаксис явно работает с самим объектом итератора, поэтому его свойства легко доступны во время итерации. С этим конструктом цикла текущее значение доступно путем индексирования в итератор, а отслеживаемый индекс — это свойство index или multi_index, в зависимости от запрошенного.
К сожалению, Python-интерпретатор печатает значения выражений внутри цикла во время каждой итерации. Мы изменили вывод в примерах, использующих этот синтаксис цикла, чтобы он был более удобочитаем.
Пример
>>> a = np.arange(6).reshape(2,3)
>>> it = np.nditer(a, flags=['f_index'])
>>> while not it.finished:
... print("%d <%d>" % (it[0], it.index), end=' ')
... it.iternext()
...
0 <0> 1 <2> 2 <4> 3 <1> 4 <3> 5 <5>
>>> it = np.nditer(a, flags=['multi_index'])
>>> while not it.finished:
... print("%d <%s>" % (it[0], it.multi_index), end=' ')
... it.iternext()
...
0 <(0, 0)> 1 <(0, 1)> 2 <(0, 2)> 3 <(1, 0)> 4 <(1, 1)> 5 <(1, 2)>
>>> it = np.nditer(a, flags=['multi_index'], op_flags=['writeonly'])
>>> with it:
.... while not it.finished:
... it[0] = it.multi_index[1] - it.multi_index[0]
... it.iternext()
...
>>> a
array([[ 0, 1, 2],
[-1, 0, 1]])
Отслеживание индекса или многомерного индекса несовместимо с использованием внешнего цикла, потому что требуется разное значение индекса для каждого элемента. Если вы попытаетесь объединить эти флаги, объект nditer сгенерирует исключение.
Пример
>>> a = np.zeros((2,3)) >>> it = np.nditer(a, flags=['c_index', 'external_loop']) Traceback (most recent call last): File "<stdin>", line 1, in <module> ValueError: Iterator flag EXTERNAL_LOOP cannot be used if an index or multi-index is being tracked
Буферизация элементов массива
При принудительном указании порядка итерации мы наблюдали, что опция внешнего цикла может предоставлять элементы меньшими блоками, потому что элементы не могут быть посещены в нужном порядке с постоянным шагом. При написании кода на C это обычно нормально, однако в чистом Python-коде это может значительно снизить производительность.
Включив режим буферизации, блоки, предоставляемые итератором внутреннему циклу, могут быть больше, значительно снижая издержки интерпретатора Python. В примере с принудительным порядком итерации Fortran внутренний цикл получает доступ ко всем элементам сразу, когда буферизация включена.
Пример
>>> a = np.arange(6).reshape(2,3) >>> for x in np.nditer(a, flags=['external_loop'], order='F'): ... print(x, end=' ') ... [0 3] [1 4] [2 5]
>>> for x in np.nditer(a, flags=['external_loop','buffered'], order='F'): ... print(x, end=' ') ... [0 3 1 4 2 5]
Итерация как конкретный тип данных
Иногда необходимо рассматривать массив как тип данных, отличный от того, в котором он хранится. Например, можно выполнить все вычисления с плавающей точкой 64-битной точности, даже если обрабатываемые массивы представляют собой числа с плавающей точкой 32-битной точности. Кроме случаев написания низкоуровневого C-кода, лучше позволить итератору обрабатывать копирование или буферизацию, а не преобразовывать тип данных в внутреннем цикле самостоятельно.
Существует два механизма, которые позволяют это сделать: временные копии и режим буферизации. С временными копиями создаётся копия всего массива с новым типом данных, затем итерация выполняется по копии. Доступ к записи разрешён с помощью режима, который обновляет исходный массив после завершения всех итераций. Основной недостаток временных копий заключается в том, что временная копия может занимать много памяти, особенно если тип итерационных данных имеет больший размер элемента, чем исходный.
Режим буферизации уменьшает проблему с объемом памяти и является более дружественным к кэшу, чем создание временных копий. За исключением особых случаев, когда весь массив необходим сразу вне итератора, буферизация рекомендуется вместо создания временных копий. В NumPy буферизация используется ufunc и другими функциями для поддержки гибких входных данных с минимальной нагрузкой на память.
В наших примерах мы будем обрабатывать входной массив с комплексным типом данных, чтобы мы могли извлекать квадратные корни из отрицательных чисел. Без включения режима копирования или буферизации итератор генерирует исключение, если тип данных не совпадает точно.
Пример
>>> a = np.arange(6).reshape(2,3) - 3 >>> for x in np.nditer(a, op_dtypes=['complex128']): ... print(np.sqrt(x), end=' ') ... Traceback (most recent call last): File "<stdin>", line 1, in <module> TypeError: Iterator operand required copying or buffering, but neither copying nor buffering was enabled
В режиме копирования флаг ‘copy’ задаётся для каждого операнда. Это делается для управления для каждого операнда. Режим буферизации задаётся как флаг итератора.
Пример
>>> a = np.arange(6).reshape(2,3) - 3 >>> for x in np.nditer(a, op_flags=['readonly','copy'], ... op_dtypes=['complex128']): ... print(np.sqrt(x), end=' ') ... 1.73205080757j 1.41421356237j 1j 0j (1+0j) (1.41421356237+0j)
>>> for x in np.nditer(a, flags=['buffered'], op_dtypes=['complex128']): ... print(np.sqrt(x), end=' ') ... 1.73205080757j 1.41421356237j 1j 0j (1+0j) (1.41421356237+0j)
Итератор использует правила приведения типов NumPy, чтобы определить, разрешено ли конкретное преобразование. По умолчанию он применяет «безопасное» приведение. Это означает, например, что он выдаст исключение, если вы попытаетесь обработать массив 64-битных чисел с плавающей точкой как массив 32-битных чисел с плавающей точкой. Во многих случаях правило «same_kind» является наиболее разумным, так как оно позволит преобразование из 64-битных в 32-битные числа с плавающей точкой, но не из чисел с плавающей точкой в целые числа или из комплексных чисел в числа с плавающей точкой.
Пример
>>> a = np.arange(6.)
>>> for x in np.nditer(a, flags=['buffered'], op_dtypes=['float32']):
... print(x, end=' ')
...
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: Iterator operand 0 dtype could not be cast from dtype('float64') to dtype('float32') according to the rule 'safe'
>>> for x in np.nditer(a, flags=['buffered'], op_dtypes=['float32'], ... casting='same_kind'): ... print(x, end=' ') ... 0.0 1.0 2.0 3.0 4.0 5.0
>>> for x in np.nditer(a, flags=['buffered'], op_dtypes=['int32'], casting='same_kind'):
... print(x, end=' ')
...
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: Iterator operand 0 dtype could not be cast from dtype('float64') to dtype('int32') according to the rule 'same_kind'
Обратите внимание на преобразование обратно к исходному типу данных при использовании операнда чтения-записи или записи только. Типичный случай — реализация внутреннего цикла в терминах чисел с плавающей точкой 64-битной точности и использование «same_kind» приведения, чтобы обработать и другие типы чисел с плавающей точкой. В режиме только чтения массив целых чисел может быть предоставлен, но режим чтения-записи вызовет исключение, поскольку преобразование обратно в массив нарушит правило приведения.
Пример
>>> a = np.arange(6)
>>> for x in np.nditer(a, flags=['buffered'], op_flags=['readwrite'],
... op_dtypes=['float64'], casting='same_kind'):
... x[...] = x / 2.0
...
Traceback (most recent call last):
File "<stdin>", line 2, in <module>
TypeError: Iterator requested dtype could not be cast from dtype('float64') to dtype('int64'), the operand 0 dtype, according to the rule 'same_kind'
Итерация массивов с векторизованным вычислением
NumPy имеет набор правил для работы с массивами различных форм, которые применяются всякий раз, когда функции принимают несколько операндов для поэлементного объединения. Это называется векторизованным вычислением. Объект nditer может применять эти правила за вас, когда вам нужно написать такую функцию.
В качестве примера мы выведем результат векторизованного вычисления одномерного и двумерного массивов.
Пример
>>> a = np.arange(3)
>>> b = np.arange(6).reshape(2,3)
>>> for x, y in np.nditer([a,b]):
... print("%d:%d" % (x,y), end=' ')
...
0:0 1:1 2:2 0:3 1:4 2:5
При возникновении ошибки векторизованного вычисления итератор выводит исключение, которое включает в себя формы входных данных для облегчения диагностики проблемы.
Пример
>>> a = np.arange(2)
>>> b = np.arange(6).reshape(2,3)
>>> for x, y in np.nditer([a,b]):
... print("%d:%d" % (x,y), end=' ')
...
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ValueError: operands could not be broadcast together with shapes (2) (2,3)
Итератор, выделяющий массивы вывода
В функциях NumPy часто встречается случай, когда выводы выделяются на основе векторизованного вычисления входных данных, а также имеется необязательный параметр «out», где результат будет помещён, если он предоставлен. Объект nditer предоставляет удобный idiom, который очень легко поддерживает этот механизм.
Мы покажем, как это работает, создав функцию square, которая возводит в квадрат свой вход. Начнём с минимальной функции определения без поддержки параметра «out».
Пример
>>> def square(a): ... with np.nditer([a, None]) as it: ... for x, y in it: ... y[...] = x*x ... return it.operands[1] ... >>> square([1,2,3]) array([1, 4, 9])
По умолчанию, nditer использует флаги «allocate» и «writeonly» для операндов, которые передаются как None. Это означает, что мы смогли предоставить только два операнда итератору, а он справился с остальным.
При добавлении параметра «out» нам нужно явно указать эти флаги, потому что если кто-то передает массив как «out», итератор по умолчанию установит «readonly», и наш внутренний цикл завершится ошибкой. Причина, по которой «readonly» является значением по умолчанию для входных массивов, заключается в предотвращении случайного запуска операции сокращения. Если значение по умолчанию было бы «readwrite», любая операция вещания также запускала бы операцию сокращения, тема, которая рассматривается позже в этом документе.
Пока мы здесь, давайте также представим флаг «no_broadcast», который предотвратит вещание результата. Это важно, потому что мы хотим получить только одно входное значение для каждого результата. Агрегирование более одного входного значения — это операция сокращения, которая требует специальной обработки. Она уже вызвала бы ошибку, потому что сокращения должны быть явно включены во флаге итератора, но сообщение об ошибке, которое получается при отключении вещания, намного понятнее для конечных пользователей. Чтобы узнать, как обобщить функцию квадратного корня до сокращения, посмотрите на функцию суммы квадратов в разделе о Cython.
Для полноты мы также добавим флаги «external_loop» и «buffered», так как они обычно необходимы для повышения производительности.
Пример
>>> def square(a, out=None): ... it = np.nditer([a, out], ... flags = ['external_loop', 'buffered'], ... op_flags = [['readonly'], ... ['writeonly', 'allocate', 'no_broadcast']]) ... with it: ... for x, y in it: ... y[...] = x*x ... return it.operands[1] ...
>>> square([1,2,3]) array([1, 4, 9])
>>> b = np.zeros((3,)) >>> square([1,2,3], out=b) array([ 1., 4., 9.]) >>> b array([ 1., 4., 9.])
>>> square(np.arange(6).reshape(2,3), out=b) Traceback (most recent call last): File "<stdin>", line 1, in <module> File "<stdin>", line 4, in square ValueError: non-broadcastable output operand with shape (3) doesn't match the broadcast shape (2,3)
Итерация внешнего произведения
Любая бинарная операция может быть расширена до операции над массивом в виде внешнего произведения, как в outer, а объект nditer предоставляет способ достижения этого путем явного сопоставления осей операндов. Это также возможно с индексированием newaxis, но мы покажем вам, как напрямую использовать параметр nditer op_axes для достижения этого без промежуточных представлений.
Мы выполним простое внешнее произведение, разместив размеры первого операнда перед размерами второго операнда. Параметр op_axes требует одного списка осей для каждого операнда и предоставляет отображение осей итератора на оси операнда.
Предположим, что первый операнд является одномерным, а второй операнд — двумерным. Итератор будет иметь три измерения, поэтому op_axes будет содержать два списка из 3 элементов. Первый список выбирает одну ось первого операнда и имеет значение -1 для остальных осей итератора, в результате чего получается [0, -1, -1]. Второй список выбирает две оси второго операнда, но не должен перекрываться с осями, выбранными для первого операнда. Его список равен [-1, 0, 1]. Операнд результата отображается на оси итератора стандартным образом, поэтому мы можем указать None вместо создания другого списка.
Операция во внутреннем цикле — это простое умножение. Все, что связано с внешним произведением, обрабатывается настройкой итератора.
Пример
>>> a = np.arange(3)
>>> b = np.arange(8).reshape(2,4)
>>> it = np.nditer([a, b, None], flags=['external_loop'],
... op_axes=[[0, -1, -1], [-1, 0, 1], None])
>>> with it:
... for x, y, z in it:
... z[...] = x*y
... result = it.operands[2] # same as z
...
>>> result
array([[[ 0, 0, 0, 0],
[ 0, 0, 0, 0]],
[[ 0, 1, 2, 3],
[ 4, 5, 6, 7]],
[[ 0, 2, 4, 6],
[ 8, 10, 12, 14]]])
Обратите внимание, что после закрытия итератора мы не можем получить доступ к operands и должны использовать ссылку, созданную внутри менеджера контекста.
Итерация сокращения
Всякий раз, когда у записываемого операнда меньше элементов, чем весь пространство итерации, этот операнд подвергается сокращению. Объект nditer требует, чтобы любой операнд сокращения был помечен как читаемо-записываемый, и допускает сокращения только при предоставлении «reduce_ok» в качестве флага итератора.
В качестве простого примера рассмотрим суммирование всех элементов в массиве.
Пример
>>> a = np.arange(24).reshape(2,3,4) >>> b = np.array(0) >>> with np.nditer([a, b], flags=['reduce_ok', 'external_loop'], ... op_flags=[['readonly'], ['readwrite']]) as it: ... for x,y in it: ... y[...] += x ... >>> b array(276) >>> np.sum(a) 276
Дело немного сложнее, когда сочетаются сокращение и выделенные операнды. Перед началом итерации любой операнд сокращения должен быть инициализирован своими начальными значениями. Вот как мы можем это сделать, суммируя по последней оси a.
Пример
>>> a = np.arange(24).reshape(2,3,4)
>>> it = np.nditer([a, None], flags=['reduce_ok', 'external_loop'],
... op_flags=[['readonly'], ['readwrite', 'allocate']],
... op_axes=[None, [0,1,-1]])
>>> with it:
... it.operands[1][...] = 0
... for x, y in it:
... y[...] += x
... result = it.operands[1]
...
>>> result
array([[ 6, 22, 38],
[54, 70, 86]])
>>> np.sum(a, axis=2)
array([[ 6, 22, 38],
[54, 70, 86]])
Для буферизированного сокращения требуется еще одна корректировка во время настройки. Обычно построение итератора включает копирование первого буфера данных из читаемых массивов в буфер. Любой операнд сокращения является читаемым, поэтому он может быть прочитан в буфер. К сожалению, инициализация операнда после завершения этой операции буферизации не отразится в буфере, с которого начинается итерация, и будут получены некорректные результаты.
Флаг итератора «delay_bufalloc» предназначен для совместного существования итератора, выделяемого итератором, с буферизацией. Когда этот флаг установлен, итератор оставит свои буферы неинициализированными до получения сброса, после чего он будет готов к обычной итерации. Вот как предыдущий пример выглядит, если мы также включим буферизацию.
Пример
>>> a = np.arange(24).reshape(2,3,4)
>>> it = np.nditer([a, None], flags=['reduce_ok', 'external_loop',
... 'buffered', 'delay_bufalloc'],
... op_flags=[['readonly'], ['readwrite', 'allocate']],
... op_axes=[None, [0,1,-1]])
>>> with it:
... it.operands[1][...] = 0
... it.reset()
... for x, y in it:
... y[...] += x
... result = it.operands[1]
...
>>> result
array([[ 6, 22, 38],
[54, 70, 86]])
Размещение внутреннего цикла в Cython
Те, кто хочет получить действительно высокую производительность от своих операций низкого уровня, должны серьезно рассмотреть прямой доступ к API итераций, предоставляемому на C, но для тех, кто не знаком с C или C++, Cython является хорошей промежуточной платформой с разумными компромиссами в производительности. Для объекта nditer это означает, что итератор позаботится о вещании, преобразовании типов и буферизации, а внутренний цикл передается Cython.
В качестве примера мы создадим функцию суммы квадратов. Для начала давайте реализуем эту функцию простым способом на Python. Мы хотим поддержать параметр «axis», аналогичный функции numpy sum, поэтому нам нужно будет создать список для параметра op_axes. Вот как это выглядит.
Пример
>>> def axis_to_axeslist(axis, ndim): ... if axis is None: ... return [-1] * ndim ... else: ... if type(axis) is not tuple: ... axis = (axis,) ... axeslist = [1] * ndim ... for i in axis: ... axeslist[i] = -1 ... ax = 0 ... for i in range(ndim): ... if axeslist[i] != -1: ... axeslist[i] = ax ... ax += 1 ... return axeslist ... >>> def sum_squares_py(arr, axis=None, out=None): ... axeslist = axis_to_axeslist(axis, arr.ndim) ... it = np.nditer([arr, out], flags=['reduce_ok', 'external_loop', ... 'buffered', 'delay_bufalloc'], ... op_flags=[['readonly'], ['readwrite', 'allocate']], ... op_axes=[None, axeslist], ... op_dtypes=['float64', 'float64']) ... with it: ... it.operands[1][...] = 0 ... it.reset() ... for x, y in it: ... y[...] += x*x ... return it.operands[1] ... >>> a = np.arange(6).reshape(2,3) >>> sum_squares_py(a) array(55.0) >>> sum_squares_py(a, axis=-1) array([ 5., 50.])
Чтобы преобразовать эту функцию в Cython, мы заменяем внутренний цикл (y[…] += x*x) кодом Cython, специализированным для типа данных float64. С включенным флагом «external_loop» массивы, предоставляемые внутреннему циклу, всегда будут одномерными, поэтому практически не требуется проверки.
Вот список sum_squares.pyx:
import numpy as np
cimport numpy as np
cimport cython
def axis_to_axeslist(axis, ndim):
if axis is None:
return [-1] * ndim
else:
if type(axis) is not tuple:
axis = (axis,)
axeslist = [1] * ndim
for i in axis:
axeslist[i] = -1
ax = 0
for i in range(ndim):
if axeslist[i] != -1:
axeslist[i] = ax
ax += 1
return axeslist
@cython.boundscheck(False)
def sum_squares_cy(arr, axis=None, out=None):
cdef np.ndarray[double] x
cdef np.ndarray[double] y
cdef int size
cdef double value
axeslist = axis_to_axeslist(axis, arr.ndim)
it = np.nditer([arr, out], flags=['reduce_ok', 'external_loop',
'buffered', 'delay_bufalloc'],
op_flags=[['readonly'], ['readwrite', 'allocate']],
op_axes=[None, axeslist],
op_dtypes=['float64', 'float64'])
with it:
it.operands[1][...] = 0
it.reset()
for xarr, yarr in it:
x = xarr
y = yarr
size = x.shape[0]
for i in range(size):
value = x[i]
y[i] = y[i] + value * value
return it.operands[1]
На этом компьютере создание модуля из файла .pyx выглядело так, но, возможно, вам придется найти некоторые руководства по Cython, чтобы узнать подробности для вашей конфигурации системы.
$ cython sum_squares.pyx $ gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -I/usr/include/python2.7 -fno-strict-aliasing -o sum_squares.so sum_squares.c
Запуск этого из интерпретатора Python даёт те же результаты, что и наш исходный Python/NumPy код.
Пример
>>> from sum_squares import sum_squares_cy >>> a = np.arange(6).reshape(2,3) >>> sum_squares_cy(a) array(55.0) >>> sum_squares_cy(a, axis=-1) array([ 5., 50.])
Небольшой замер времени в IPython показывает, что уменьшенные накладные расходы и выделение памяти внутреннего цикла Cython обеспечивают значительное ускорение по сравнению с простым кодом Python и выражением с использованием встроенной функции sum NumPy.
>>> a = np.random.rand(1000,1000) >>> timeit sum_squares_py(a, axis=-1) 10 loops, best of 3: 37.1 ms per loop >>> timeit np.sum(a*a, axis=-1) 10 loops, best of 3: 20.9 ms per loop >>> timeit sum_squares_cy(a, axis=-1) 100 loops, best of 3: 11.8 ms per loop >>> np.all(sum_squares_cy(a, axis=-1) == np.sum(a*a, axis=-1)) True >>> np.all(sum_squares_py(a, axis=-1) == np.sum(a*a, axis=-1)) True
© 2005–2019 NumPy Developers
Licensed under the 3-clause BSD License.
https://docs.scipy.org/doc/numpy-1.16.1/reference/arrays.nditer.html
