Итерация по массивам
Объект итератора 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после завершения итерации, что вызовет запись обратно.
nditer больше нельзя использовать для итерации после вызова close или выхода из его контекста.
Пример
>>> 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, чтобы определить, разрешено ли конкретное преобразование. По умолчанию он применяет правило «safe» преобразования. Это означает, например, что он вызовет исключение, если вы попытаетесь рассматривать массив с плавающей точкой 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 будет содержать два списка из трёх элементов. Первый список выбирает одну ось первого операнда, а остальные оси итератора имеют значение -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.15.4/reference/arrays.nditer.html
