Как мы уже говорили, Python — язык для быстрой разработки. Однако чистый Python не предназначен для написания быстрых программ. Это интерпретируемый язык, поэтому программы на Python выполняются медленнее аналогов на C, C++ или Fortran. С другой стороны, математики, физики, биологи и инженеры часто применяют Python для решения вычислительных задач.
Нет ли тут противоречия? Как интерпретируемый язык может быть эффективен в вычислительной математике?
Оказывается, все дело в библиотеках. Python отлично подходит на роль промежуточной среды, оболочки, «клея» между библиотеками, написанными на разных языках.
В этом уроке мы поговорим о наиболее фундаментальной библиотеке для работы с вычислительной математикой — Numpy.
Многие другие пакеты для работы с данными и вычислениями используют базовые объекты этой библиотеки. В числе таких пакетов OpenCV — открытая библиотека для работы с распознаванием образов — и Pandas — библиотека, ориентированная на анализ данных.
Измерение скорости
Замеры времени (а в более общем случае, профилирование) — не такая простая процедура, как может показаться. Обычно проще всего замерить астрономическое время между началом и концом выполнения задачи, усредняя результаты нескольких опытов. Однако почти все ОС многозадачны, поэтому процессор (и его вычислительные ядра) редко отдаются одной задаче в единоличное пользование. Как следствие, в измерениях лучше использовать именно процессорное время. Механизмы замера времени сами вмешиваются в процесс и немного влияют на результат — как и в любом физическом эксперименте.
В Python для замера времени работы кода служит библиотека timeit.
Вы можете познакомиться с возможностями модуля в соответствующем разделе документации.
Например, мы можем замерить три разных способа заполнить список из миллиона квадратных корней.
import timeit
print(timeit.timeit("[sqrt(x) for x in range(1000000)]",
"from math import sqrt", number=1))
print(timeit.timeit("for i in range(1000000): a.append(sqrt(i))",
"from math import sqrt; a=[]", number=1))
print(timeit.timeit("list(map(sqrt, range(1000000)))",
"from math import sqrt; a=[]", number=1))
Как видим, в этой версии интерпретатора (3.7) предпочтительно использовать map.
Самый медленный способ — это, конечно же, динамическое расширение существующего списка (append). Причем, чем больше список — тем медленнее он меняет свой размер. Это вызвано необходимостью иногда переносить данные из одного места списка в другое.
Несмотря на относительную быстроту (0,15 секунд на извлечение 1 000 000 квадратных корней), скорость можно увеличить еще примерно в 10 раз. Давайте посмотрим как.
Массивы в Numpy
Основной объект в Numpy — многомерный массив.
Массивы — одна из базовых структур данных, которая позволяет моделировать многие объекты, относящиеся как к науке, так и к обычной жизни: список покупок, результаты наблюдения температуры, матрицы и векторы, изображения, видео и т. д. Напомним, что в чистом Python нет типа данных с именем массив, и нам приходится моделировать его с помощью списков.
Другое дело — Numpy. За тип массива здесь отвечает объект array.
Как же создать массив?
Во-первых, массив можно сделать из обычного списка:
В Numpy элементы одного массива должны быть однородны (одного типа). Это самое важное идеологическое отличие массивов от списков, в которых можно хранить объекты разной природы.
Numpy создаст массив из юникод-строк длиной 32. За тип элементов в большинстве случаев отвечает параметр dtype(data type). Обратите внимание на тип данных элементов массива.
Посмотрите так же на использование параметра dtype:
a = np.array([1, 3, 8])
a # => array([1, 3, 8])
a.dtype
# => dtype('int64')
type(a[0])
# => numpy.int64
a = np.array([1, 3, 8], dtype=np.float64)
type(a[0])
# => numpy.float64
Указание типов и работа с ними нужны, поскольку языки, на которых написана эта библиотека, строго типизованы. Вдобавок это увеличивает скорость обработки данных.
Размерность массива
Размерность массива можно в любой момент изменить операцией reshape().
Узнать размерность можно атрибутом shape.
Вообще говоря, размерность — всего лишь «синтаксический сахар». В памяти все может храниться как одномерный массив с пересчетом координат элемента. Таким образом, операция reshape() — просто изменение коэффициентов в алгоритме, а не перераспределение и копирование данных.
Например, фильм можно представить в виде 4-мерного массива кадров. Кадр — картинка, то есть трехмерный массив. Его можно представить и как двумерный массив пикселей, где каждый пиксель — одномерный массив из трех элементов: R, G, B.
Самое главное — при использовании функции reshape() произведение ее параметров должно быть равно количеству элементов в массиве.
Иначе мы получим ошибку:
a.reshape(2, 3, 4)
ValueError Traceback (most recent call last)
<ipython-input-31-a907d0800243> in <module>()
----> 1 a.reshape(2, 3, 4)
ValueError: cannot reshape array of size 100 into shape (2, 3, 4)
Кроме того, доступ можно организовать через списки с индексами:
a[[1], [4, 4, 7, 8]] # => array([14, 14, 17, 18])
Но самой удобной альтернативой обычному способу является тот, в котором в качестве «адреса» элемента используется кортеж координат:
a[(7, 9)] # => 79
Скобки, конечно же, можно опустить:
a[7, 9] # => 79
Массовые операции
В них вся суть. Ускорение, по сравнению с чистым Python, достигается за счет замены цикла обработки в Python вызовом одной функции, внутри которой находится скомпилированный машинный цикл. Чтобы упростить код и ускорить его выполнение, нужно убрать весь поток логики с «пробеганием» по массивам, заменив его вызовами базовых агрегатных функций или операторов.
Инициализация:
# заполняем единицами
np.ones(10)
array([ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
# заполняем единицами целого типа
np.ones(10, dtype=np.int32)
В этом случае сортировка происходит по последнему измерению. Обходя матрицу, мы сначала выбираем строку, а потом идем по этой строке, поэтому последнее измерение в данном случае — это строка. В результате мы отсортировали независимо каждую строку.
# а теперь укажем конкретное измерение
print(np.sort(a, axis=0))
Но самое интересное в том, что если в качестве значения параметра axis указать None, то матрица перед сортировкой будет линеаризована, то есть превращена в одномерный массив.
Обратите внимание: подобное поведение характерно не только для функции sort(), но и для многих других функций: min(), sum() и т. д.
Но об этом вы можете почитать самостоятельно на странице с документацией.
Вспомним PIL
Работая с библиотекой PIL, тоже можно использовать средства Numpy. Например, если мы хотим сделать изображение темнее оригинала, можем просто поделить его составляющие, например, на 10:
from PIL import Image
import numpy as np
# получим массив numpy из картинки, которую откроем из файла.
image = np.asarray(Image.open('images/Риана.jpg'))
# поделим все элементы массива на 10, приведем к типу uint8 (один байт
# без знака) преобразуем в изображение и сохраним в файл
Image.fromarray(np.uint8(image // 10)).save('r2.jpg')
Игра «Жизнь»
Несколько десятилетий назад Джон Конуэй придумал один из самых известных клеточных автоматов, который назвал игрой «Жизнь». Простота правил сочетается в ней с богатством результатов. Многие компьютерные инженеры хоть раз обращались к программированию и исследованию этой игры, которая послужила интересной моделью для многих отраслей науки.
Клеточный автомат — модель однородного пространства с некоторыми клетками. Каждая клетка может находиться в одном из нескольких состояний и иметь некоторое количество соседей. Задаются правила перехода из одного состояния в другое в зависимости от текущего состояния клетки и ее соседей.
Пространство «Жизни» — бесконечное поле клеток.
Каждая клетка имеет 8 соседей (сверху, снизу, справа, слева и по диагонали). Клетка может иметь два состояния — живое (на клетке стоит фишка) и мертвое (фишки нет).
Правила изменения следующие:
Если клетка была живой, она выживет, если у нее 2 или 3 соседа. Если соседей 4, 5, 6, 7 или 8, она умирает от перенаселенности, а если 0 или 1 — от одиночества
Новая клетка рождается в поле, у которого есть ровно 3 соседа
Время в этой игре дискретно и считается поколениями. Все начинается с начальной расстановки фишек (0 поколение), в дальнейшем рассматривается эволюция клеточного пространства в 1, 2, 3 и т. д. поколении. Процессы смерти и рождения происходят одновременно, после чего строится следующее поколение.
Давайте попробуем написать игру «Жизнь», используя библиотеку Numpy. Пусть у нас будет поле 10×10, в центр которого поместим конструкцию, известную как глайдер. Мы скоро выясним, почему она так называется.
Поле имеет тип uint8, чтобы оно занимало меньше памяти. Каждый его элемент занимает ровно 1 байт (8 бит) и является целым беззнаковым (unsigned) числом в диапазоне от 0 до 255.
Живые клетки обозначаются единицей, а мертвые — нулем. Нужно решить, что делать на границах поля. Мы не можем обеспечить бесконечность в обоих направлениях, поэтому замкнем поле само на себя. Если выйти за нижнюю границу, окажемся наверху, а если за правую — появимся слева, и наоборот. Получается что-то вроде бублика. Такая фигура называется тор.
Для начала познакомимся с операцией roll, доступной для массивов. Она сдвигает исходный массив вдоль одного из измерений (в данном случае — строки или столбца).
Выполним на матрице следующую операцию: «если у клетки 3 соседа, то в следующем поколении на этом месте будет клетка; а если 2 — клетка будет при условии, что она была „жива“ в текущем поколении».
Для этого воспользуемся операторами | (или) и &**** (и).
# выделим клетки, у которых ровно три соседа
neighbors == 3
Объединить матрицы с логическими и целочисленными элементами можно, поскольку они в данном случае могут быть сведены друг к другу: 0 — False, 1 — True.
Проследим эволюцию глайдера на протяжении четырех поколений. Для этого создадим функцию next_population().
С визуализацией у нас не очень здорово, но видно, что глайдер «летит»: каждые четыре поколения он сдвигается вниз и вправо. Иными словами, он движется в правый нижний угол.
Почему глайдер не исчезает с поля, а странным образом сначала разбивается на части, а потом материализуется в левом верхнем углу?
Это связано с тем, что фактически мы строим нашу «жизнь» не на плоскости, а на торе. Поэтому, «улетая» вправо, клетка прилетает слева, улетая вверх, прилетает снизу, а углы нашего поля — это одна и та же ячейка.
