Итерация по массивам
Объект-итератор 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 для управления этим аспектом итерации. По умолчанию, с поведением, описанным выше, порядок ’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 печатает значения выражений внутри цикла while во время каждой итерации цикла. Мы изменили вывод в примерах, использующих эту конструкцию цикла, чтобы сделать его более читабельным.
Пример
>>> 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 буферизация используется функциями ufuncs и другими функциями для поддержки гибких входных данных с минимальными накладными расходами памяти.
В наших примерах мы будем обрабатывать входной массив с комплексным типом данных, чтобы мы могли извлекать квадратные корни отрицательных чисел. Без включения режимов копирования или буферизации итератор сгенерирует исключение, если тип данных не совпадает точно.
Пример
>>> 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 предоставляет удобный стиль, который очень упрощает поддержку этого механизма.
Покажем, как это работает, создав функцию 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, так и с выражением, использующим встроенную функцию NumPy sum.:
>>> 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.17.0/reference/arrays.nditer.html
