Программирование и научные вычисления на языке Python/§8

В третьем уроке мы узнали о списках, как об удобном способе хранения табулированных данных. Массив представляет собой объект, близкий к списку, но менее гибкий, а в вычислительном плане более эффективный. Когда мы используем компьютер для математических расчетов, мы часто сталкиваемся с огромным множеством чисел и связанных с ними арифметических операций. Хранение чисел в списках в таких случаях может привести к значительному снижению скорости работы программы, в то время как хранение в виде массивов чисел существенно ускоряет решение. Это может быть не очень важным для примеров этого курса, поскольку мы рассматриваем небольшие программы, работающие и с маленькими объемами данных, которые выдают результат в течение нескольких секунд. Тем не менее, более продвинутые приложения, особенно используемые для расчетов в промышленности и науке, прежде чем дать ответ, могут искать его недели и месяцы. Поэтому любая идея, уменьшающая время получения результата, всегда приветствуется. Однако, стоит сказать, что многие программисты изначально предъявляют слишком большое усердие в увеличении скорости, используя сложные конструкции, приводящие к тому, что программы становится дальше сложно поддерживать и совершенствовать. В первую очередь следует стремиться писать ясные, хорошо структурированные и легкие для понимания программы, а уже после этого, на следующем этапе вам будет гораздо проще разобраться как можно ускорить вычисления. В Python довольно часто самое ясное решение работает быстрее менее ясных.

Этот урок кратко знакомит нас с массивами — как они создаются и как могут использоваться. Работа с массивом обычно заканчивается большим количеством чисел, и довольно трудно понять, что они дают, если просто взглянуть на них. Поэтому такую информацию визуализируют в виде графиков кривых, о чем мы поговорим в следующем уроке. И там мы будем использовать массивы для хранения информации о координатах точек графика. То есть не только массивы требуют визуализации, но и графики требуют для себя массивов.

Векторы

Сейчас мы немного поговорим о векторах с тем предположением, что вы что-то слышали о векторах ранее. Это нам нужно как почва для того, чтобы начать работать с массивами и графиками.

Некоторые математические величины связаны с набором чисел. Например, точка на плоскости имеет две координаты, x и y, и точка ими и описывается как (x, y), где вместо символов можно подставить любые числа. То есть точка описывается в виде группы чисел, заключенных в скобки. Точка в трехмерном пространстве описывается схожим способом (x, y, z) или (x₁, x₂, x₃). Когда решаются n уравнений с n неизвестными, решение дает вам группу из n чисел (x₁, x₂, ... , x_n-1, x_n).

Такие величины, как (x, y), (x, y, z), (x₁, ... , x_n) могут быть представлены в виде векторов, идущих из начала координат в указанную точку. Например, вектор (x, y) идет из точки (0, 0) в точку (x, y), как и трехмерный вектор (x, y, z) идет из(0, 0, 0) в (x, y, z). Для последнего случая удобно ввести n-мерное пространство, где вектор идет из (0, ..., 0) в (x₁, ... , x_n). Векторы, как и массивы, можно визуализировать. На плоскости вектор можно представить в виде стрелки. Два вектора, имеющих одинаковое направление и длину, эквивалентны.

О векторе (x₁, ... , x_n) говорят, что он содержит n компонент. Каждое из чисел x₁, x₂, ... это компоненты вектора. Для того, чтобы записать вектор в Python, мы можем использовать списки или кортежи:

v1 = [x, y]
v2 = (-1, 2)
v3 = (x1, x2, x3)
from math import  exp
v4 = [exp(-i*0.1) for i in range(150)]

Здесь v1 и v2 — векторы на плоскости, v3 — вектор в трехмерном пространстве, а v4 — вектор в 150-мерном пространстве, состоящий из 150 значений экспоненциальной функции. Поскольку в Python (и многих других языках) индексация начинается с нуля, то более естественным записывать вектор вместо (x₁, x₂) как вектор (x₀, x₁). Это не общепринято в математике, но существенно сближает язык математики и язык программирования, что значительно облегчает понимание и уменьшает число потенциальных ошибок.

Невозможно представить как выглядит 150-мерное пространство. Переход от плоскости к пространству и тот бывает дается тяжело. Но представить как происходит переход к четырех-, пяти-, и-так-далее-мерному вектору в виде списка компонент не составляет труда.

Математические операции над векторами

С тех пор, как векторы были введены как массивы чисел имеющие длину и направление, они тут же оказались очень удобны в геометрии и физике. У скорости машины есть значение и направление, есть ускорение и позиция машины также есть точка, которую, как показано выше, можно представить в виде вектора. Грань треугольника также может быть рассмотрена как линия (стрелка), имеющая направление и длину.

В физике и геометрии, использующей векторы, очень важны применяемые математические операции. Давайте рассмотрим наиболее часто встречаемые операции и действующие математические правила. Для этого возьмем два вектора, (u₁, u₂) и (v₁, v₂) и для начала сложим их:

$~(u_{1},u_{2})+(v_{1},v_{2})=(u_{1}+v_{1},u_{2}+v_{2}).$

(8.1)

Для вычитания применяется такое же правило:

$~(u_{1},u_{2})-(v_{1},v_{2})=(u_{1}-v_{1},u_{2}-v_{2}).$

(8.2)

Вектор может быть умножен на число:

$~a\cdot (v_{1},v_{2})=(a\cdot v_{1},a\cdot v_{2})$

(8.3)

и скалярно на вектор, что даст число:

$~(u_{1},u_{2})\cdot (v_{1},v_{2})=u_{1}v_{1}+u_{2}v_{2}.$

(8.4)

Также возможно и векторное произведение, но рассматривать его здесь будет долго. Длина вектора определяется:

$~||(v_{1},v_{2})||={\sqrt {(v_{1},v_{2})\cdot (v_{1},v_{2})}}={\sqrt {{v_{1}}^{2}+{v_{2}}^{2}}}.$

(8.5)

Все эти операции можно по аналогии продолжить и на n-мерное пространство.

Векторные функции

Кроме операций, о которых мы напомнили себе выше, существуют и другие, играющие существенную роль в математических приложениях и особенно в таких средах как Matlab, Octave, Python и R. Эти операции вы вряд ли найдете в книгах посвященных математике, они относятся исключительно к потребностям, возникающим при программировании массивов. Для каждого элемента вектора, его компоненты, мы можем сопоставить функцию одной переменной f, тогда мы можем получить и некоторую векторную функцию, в которой компонентами служат функции компонент. Например, у нас есть вектор v = (v₀, ..., v_n-1). Тогда его векторная функция будет выглядеть как f(v) = (f(v₀), ..., f(v_n-1)). Например, синус от v будет записан: sin(v) = (sin(v₀), ..., sin(v_n-1)).

Векторное возведение в степень может означать: v^b = (v₀^b, ..., v_n-1^b). Особое векторное произведение ("asterix" multiplication) определяется как 'u * v' = (u₀v₀, u₁v₁, ... , u_n-1v_n-1). В компьютерных вычислениях возможна и операция прибавления скаляра к вектору — число прибавляется к каждому элементу вектора. Возможны и сложные выражения, с которыми мы столкнемся далее.

Снова отметим, что эти функции чаще всего мало имеют отношения к обычной математике векторов, в которой то же складывание вектора и скаляра невозможно, а возведение вектора в квадрат даст число, его длину в квадрате. Эти функции работают поочередно с каждым элементом и результатом функции является уже вектор таких элементов. Такие функции позволяют значительно ускорить работу с массивами, производя одновременно одни и те же действия над всеми элементами.

Использование списков

Представим, у нас есть функция f(x) и мы хотим применить ее к n числам x₁, x₂, ... , x_n-1, x_n. Мы можем составить n пар (x_i, f(x_i)), а можем создать два списка — один со значениями переменной, а другой с соответствующими значениями функции:

>>> def f(x):
...         return x**3
...
>>> n = 5 # number of points along the x axis
>>> dx = 1.0/(n-1) # spacing between x points  in [0,1]
>>> xlist = [i*dx for i in range(n)]
>>> ylist = [f(x) for x in xlist]
>>> pairs = [[x, y] for x, y in zip(xlist, ylist)]

Здесь для решения задачи мы использовали два приема: генерацию списков и двойной zip-проход по спискам. В списке pairs все элементы представляют собой списки из двух float чисел, в списках xlist и ylist все объекты float. Но список это довольно гибкий объект, и он может содержать объекты любых типов:

mylist = [2, 6.0, 'tmp.ps', [0,1]]

Также мы можем легко изменять, добавлять и удалять новые элементы из любого места списка. Эта гибкость списков делает их очень удобной для программистов, но в случае когда элементы однотипны и их число фиксировано, вместо списков используются массивы. Преимущества массивов в быстроте вычислений, меньшей занимаемой памяти и исключительно обширной математической поддержке таких данных. Поэтому массивы, как вы увидите в этом курсе, на практике (и в крупных математических пакетах) находят такое широкое применение. Списки отныне мы будем применять по назначению — когда нам будет нужно удалять и добавлять элементы и использовать в данных объекты различных типов.

Основы Numerical Python

Объект array может быть рассмотрен как вариант списка, но с учетом следующих допущений и возможностей:

• Все элементы массива представлены одним типом объектов, например целыми, действительными или комплексными числами, что делает их хранение и обработку наиболее эффективным.
• В тот момент, когда создается массив, число его элементов должно быть известно.
• Массивы не являются стандартной частью Python — они требуют специального дополнительного пакета, которым пользуются практически все, кто занимаются научными проектами на Python. Этот пакет называется Numerical Python или еще чаще NumPy, поскольку после его установки вызов осуществляется с помощью обычной инструкции импорта модуля: import numpy. Для того, чтобы установить NumPy, загрузите его с официального сайта проекта. На этой же странице вы обнаружите еще один пакет, который нам понадобится в дальнейшем — SciPy.
• С numpy широкий круг математических операций может быть решен непосредственно с помощью массивов, таким образом исключается потребность в циклах, проходящих по элементам массива. Это свойство носит названия векторизации (vectorization) или прорисовки.
• Массивы с одним индексом также называют векторами. Массивы с двумя индексами используются для создания матриц и представления табличной информации. Массив может содержать практически любое количество индексов, то есть быть n-мерным.

Как уже было сказано, после установки пакета, работа с модулем происходит обычным образом:

from numpy import *

Конвертирование списка r в массив a происходит привычным способом, но с помощью импортированной из numpy функции:

a = array(r)

Для того, чтобы создать массив из n нулевых элементов используем функцию zeros:

a = zeros(n)

Элементы по умолчанию являются float-объектами, второй аргумент функции позволяет изменить тип объектов, например, на int. Часто бывает нужно создать массив из элементов, равномерно распределенных в интервале [p, q]. Для этого в numpy есть функция linspace:

a = linspace(p, q, n)

Вообще говоря в numpy имеется огромное количество функций и внутренних модулей.

Доступ к элементу осуществляется так же как в списках, например, a[1]. Срезы тоже здесь работают, например срез a[1:-1], извлекает список всех элементов, кроме первого и последнего. Но в отличие от списка, здесь это не копия. Например:

b = a[1:-1]
b[2] = 0.1

изменится и массив a, его элемент a[3]=0.1.

К слову, о срезах

Срез в формате a[i:j:s] выбирает все элементы, начиная с i, заканчивая, но не включая, j с шагом s. Например, срез a[0:-1:2] выбирает каждый второй элемент, кроме последнего. Как и ранее, возможны пропуски аргументов, например a[::4] выберет каждый четвертый элемент. Можно взять и отрицательный шаг, тогда элементы будут идти в обратном порядке.

Задание координат и значений функций

Теперь, когда у нас есть эти простейшие операции, мы можем продолжить пример, в котором мы использовали списки:

>>> from numpy import  *
>>> x2 = array(xlist)
>>> y2 = array(ylist)
>>> x2
array([ 0.   ,   0.25,   0.5 ,   0.75,   1.   ])
>>> y2
array([ 0. ,   0.015625,  0.125    ,   0.421875,  1. ])

Вместо того, чтобы сначала создавать список, а потом конвертировать его в массив будет естественным сразу же создавать массив. Координаты, что мы задавали в xlist легко получить в виде массива с помощью функции linspace. Массив для значений мы создадим с помощью zeros, чтобы ему изначально была правильно отведена длина, в соответствии с числом элементов в xlist. Далее мы заполняем его с помощью цикла:

>>> from numpy import  *
>>> n = len(xlist)
>>> x2 = linspace(0, 1, n)
>>> y2 = zeros(n)
>>> for i in xrange(n):
...         y2[i] = f(x2[i])
...
>>> y2
array([ 0. ,    0.015625,  0.125    ,   0.421875,  1. ])

Заметьте, что в цикле мы используем вместо range другую функцию — xrange. Она является более предпочтительной для (обычно больших) массивов. Также отметим, что для y мы использовали генерацию списка, а для y2 — цикл for, поскольку массив это не список. Из положения можно выйти с помощью конвертирования:

>>> x2 = linspace(0, 1, n)
>>> y2 = array([f(xi) for xi in x2])

Тем не менее, есть лучший вариант, который объясняется далее.

Векторизация

Великолепным преимуществом массивов является то, что они могут обходиться без циклов и функция может применяться, как мы объясняли выше, к самому массиву и производить действия над всеми элементами:

>>> y2 = f(x2)
>>> y2
array([ 0. ,   0.015625,  0.125    ,   0.421875,  1. ])

И даже сложные составные выражения

r = sin(x)*cos(x)*exp(-x**2) + 2 + x**2

подвластны волшебству массивов:

r = zeros(len(x))
for i in xrange(len(x)):
    r[i] = sin(x[i])*cos(x[i])*exp(-x[i]**2) + 2 + x[i]**2

Это свойство и называется векторизацией. Существенный выигрыш в скорости по сравнению со списками происходит из-за того что в генерации списков используется относительно медленные циклы самого Python, в то время как векторизация их никак явно не использует, а задействует «быстрые циклы» внутри numpy. Кроме того, что векторизация существенно повышает скорость обработки, она делает код более понятным и ясным для чтения.

Но приведенный выше код не является "чистой векторизацией", в нем используется цикл for, без которого можно обойтись, если использовать тригонометрические функции из пакета numpy (которые поддаются векторизации):

from numpy import *
x=linspace(-pi,pi,11)

r = sin(x)*cos(x)*exp(-x**2) + 2 + x**2

Ссылки

• Проект научных вычислений SciPy — отсюда можно скачать пакеты NumPy и SciPy
• Краткое введение в NumPy
• Шпаргалка по массивам в SciPy (NumPy)
• Руководство по NumPy (англ.)