Python в научной работе. Введение

Здесь представлен учебный материал для начального ознакомления с Python, который я подготовил для своих коллег из Института структурной биологии в Гренобле. Вы можете использовать его для других целей, предлагать свои улучшения, но имейте в виду, что курс не предназначался для самообучения, т.е. информация в нем представляется не обязательно полностью.

Для примеров, приводимых в курсе, требуются модули NumPy и ScientificPython, которые должны быть установлены на вашем компьютере.


Конрад ХИНСЕН
hinsen@cnrs-orleans.fr

Перевод на русский язык: Ф.С.ЗАНЬКО
russianlutheran@gmail.com


Python для инженеров и исследователей

Интернет-адрес оригинального документа:
http://starship.python.net/crew/hinsen/
[Сама ссылка не работает. Осталась ее копия в архиве Интернета.- Прим. пер.]

Перевод на русский язык осуществлен с любезного разрешения автора. Разрешается свободное распространение и использование настоящего перевода для любых целей при условии сохранения текста перевода в неизменном виде.

О замеченных ошибках, неточностях, опечатках просьба сообщать по электронному адресу:
russianlutheran@gmail.com



СОДЕРЖАНИЕ



Часть 1: Python в роли калькулятора

Числа:

Целые 0, 1, 2, 3, -1, -2, -3
Вещественные 0., 3.1415926, -2.05e30, 1e-4
(обязательно должны содержать десятичную точку или экспоненту)
Мнимые/Комплексные 1j, -2.5j, 3+4j
(последнее является суммой)
Сложение 3+4, 42.+3, 1+0j
Вычитание 2-5, 3.-1, 3j-7.5
Умножение 4*3, 2*3.14, 1j*3j
Деление 1/3, 1./3., 5/3j
Возведение в степень 1.5**3, 2j**2, 2**-0.5

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

Предупреждение: 1/3 - целое число, следовательно, равно нулю. 1./3. - это, как и можно предположить, число вещественное.

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

Математические функции:

Стандартные математические функции (sqrt, log, log10, exp, sin, cos, tan, arcsin, arccos, arctan, sin, cosh плюс некоторые другие, которые будут упомянуты позже), а также константы pi и e не относятся к базовой части языка, а сосредоточены в отдельном модуле Numeric. Поэтому, прежде чем использовать, их нужно импортировать:

Импортировать функции можно по отдельности

from Numeric import sqrt
или группами
from Numeric import sin, cos, tan
или же возможно импортировать сразу все содержимое модуля:
from Numeric import *

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

import Numeric
и использовать "расширенные" названия функций: numpy.sqrt, numpy.exp и т.д.

Пример:

from Numeric import pi, sin
print sin(pi/3)

Строки:

Две формы представления: 'abc' or "abc"

Перенос строки обозначается как \n: "abc\ndef"

Объединение: "abc"+'def' дает "abcdef"

Повторение: 3*"ab" дает "ababab"

Векторы:

Векторы не относятся к основным типам данных в Python. Они определены в модуле Scientific.Geometry, откуда их и нужно импортировать:

from Scientific.Geometry import Vector

Обозначение: Vector(1,0,0)

Сложение, вычитание: Vector(1,0,0)+Vector(0,-1,3), Vector(,1,0)-Vector(1.5,4,0)

Умножение на скаляр: 3.5*Vector(1,1,0), Vector(0,0,1)*4.

Деление на скаляр: Vector(1,1,0)/2.

Скалярное произведение: Vector(1,2.5,0)*Vector(0,-1,3.1)

Векторное произведение: Vector(1,2.5,0).cross(Vector(0,-1,3.1))

Длина: Vector(2.5, 3.4, 1.).length()

Нормализация: Vector(1.,4.,2.).normal()

Угол между двумя векторами: Vector(1,2.5,0).angle(Vector(0,-1,3.1)) (результат в радианах)

Доступ к компонентам: Vector(1,0,3)[0], Vector(1,4.,3).normal()[1] (индекс меняет свое значение от 0 до 2)

Переменные:

Переменные используются для хранения и идентификации числовых значений:

x = 2.
sum = x + 25
greeting = "hello"
a_very_special_value = 42

Переменные могут иметь произвольную длину, содержать буквы и цифры. Но первым символом в названии переменной обязательно должна быть буква. Строчные и заглавные буквы считаются разными символами.


Ваша первая программа на Python

В программе задаются три точки, определяющие некий треугольник; на экран выводятся расстояния между всеми возможными парами этих точек.

from Scientific.Geometry import Vector

a = Vector(0, 1, 0)
b = Vector(4.3, -2.4, 0.005)
c = Vector(-3.2, 5.1, -3.)

print "a-b:", (a-b).length()
print "a-c:", (a-c).length()
print "b-c:", (b-c).length()

Упражнения


Часть 2: Работа с текстовыми файлами

Последовательности объектов

Последовательности объектов - списки, массивы, текстовые строки и т.п.,- включают в себя элементы, на которые можно ссылаться посредством индексов. Элементы любой последовательности объектов seq можно вызывать при помощи индексирования: seq[0], seq[1]. Индексирование можно проводить и с конца: seq[-1] - это последний элемент последовательности, seq[-2] - предпоследний и т.д.

Пример:

text = 'abc'
print text[1]
print text[-1]

Существует возможность вырезать подпоследовательности: seq[0:5] - это подпоследовательность, включающая в себя первые пять элементов с номерами 0, 1, 2, 3 и 4, но не элемент под номером 5. Допустимы и отрицательные индексы: выражение seq[1:-1] отбрасывает первый и последний элементы последовательности.

Пример:

text = 'Довольно длинная строчка.'
print text[2:10]
print text[-7:-1]

Длину последовательности можно определить с помощью функции len(seq).

Списки

Списки - это последовательности, которые могут включать в себя произвольные объекты (числа, строки, векторы, другие списки и т.п.):

some_prime_numbers = [2, 3, 5, 7, 11, 13]
names = ['Smith', 'Jones']
a_mixed_list = [3, [2, 'b'], Vector(1, 1, 0)]

Элементы и подпоследовательности списка могут быть изменены, если им присвоить новые значения:

names[1] = 'Python'
some_prime_numbers[3:] = [17, 19, 23, 29]

С помощью выражения list.append(new_element) к концу списка можно добавить новый элемент. Список также можно отсортировать (list.sort()) или изменить порядок следования элементов на противоположный: list.reverse().

Два списка могут быть объединены как и текстовые строки: [0, 1] + ['a', 'b'] дает [0, 1, 'a', 'b'].

Списки также можно повторять подобно текстовым строкам: 3*[0] дает [0, 0, 0].

Кортежи

Кортежи очень похожи на списки, за исключением того, что их никоим образом нельзя изменять. Будучи однажды создан, кортеж всегда будет включать в себя одни и те же элементы. Поэтому его можно применять в ситуациях, когда не имеет смысла использовать изменяемую последовательность (например, ключ в базе данных).

Пример:

tuple_1 = (1, 2)
tuple_2 = ('a', 'b')
combined_tuple = tuple_1 + 2*tuple_2

Циклическая обработка последовательности

Часто бывает необходимо произвести некоторую операцию над каждым элементом последовательности. Это называется циклической обработкой последовательности.

for prime_number in [2, 3, 5, 7]:
    square = prime_number**2
    print square

Замечание: Операции внутри цикла выделяются отступом.

Циклической обработке можно подвергнуть любую последовательность. Вот пример с текстовыми строками:

for vowel in 'aeiou':
    print 10*vowel

Циклическая обработка числового ряда - отдельный случай:

from Numeric import sqrt

for i in range(10):
    print i, sqrt(i)
Функция range(n) возвращает список, содержащий первые n целых чисел, т.е. от 0 до n-1. Выражение range(i, j) возвращает все целые от i до j-1, а выражение range(i, j, k) - числа i, i+k, i+2*k и т.д. вплоть до j-1.

Проверка условий

Чаще всего проверяются условия равенства и величины:
равно a == b
не равно a != b
больше a > b
меньше a < b
больше или равно a >= b
меньше или равно a <= b

Условия можно комбинировать с помощью операций and (И) и or (ИЛИ); операция отрицания вводится выражением not (НЕ). Результат проверки условия равен 1 для значения "истина" и 0 для значения "ложь".

Наиболее часто условия применяются для программирования выбора из нескольких возможных вариантов:

if a == b:
    print "equal"
elif a > b:
    print "greater"
else:
    print "less"
Может быть любое число ветвей elif (в том числе ни одной), а ветвь else является необязательной.

Кроме того, условия можно использовать для управления циклом:

x = 1.
while x > 1.e-2:
    print x
    x = x/2

Текстовые файлы

Текстовые файлы как объекты определены в модуле Scientific.IO.TextFile.

Чтение

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

from Scientific.IO.TextFile import TextFile

for line in TextFile('some_file'):
    print line[:-1]

Почему line[:-1]? Последний символ каждой строки - это символ перевода каретки, который мы печатать не хотим (он будет генерировать пустые строки между строчками нашего файла).

Текстовые файлы как объекты могут взаимодействовать и со сжатыми файлами. Если имя файла заканчивается на ".Z" или ".gz", значит оно ссылается на заархивированный файл. Программы, конечно же, получат разархивированную версию этих данных. Вы можете даже вместо имени файла использовать URL (Universal Resource Locator - унифицированный указатель информационного ресурса; адрес, используемый Web-браузерами) и, таким образом, считывать данные непосредственно из Internet!

Запись

Вместо чтения текстовые файлы можно открыть для записи:

from Scientific.IO.TextFile import TextFile

file = TextFile('a_compressed_file.gz', 'w')
file.write('The first line\n')
file.write('And the')
file.write(' second li')
file.write('ne')
file.write('\n')
file.close()

Файлы, открытые для записи, после окончания работы с ними следует закрывать, чтобы быть уверенным, что все данные действительно в него записаны. Конечно, по окончании работы программы все открытые файлы закрываются автоматически, но лучше на это не полагаться.

Обратите внимание, что для режима записи также доступна автоматическая архивация, но не адреса URL, поскольку большинство серверов в Internet не разрешают запись по очень веским причинам!

Некоторые полезные операции со строками

Модуль string содержит операции с обычными строками, которые определенно полезны при чтении и записи текстовых файлов. Здесь будут описаны только самые важные; полный список см. Python Library Reference.

Извлечение данных из строк

Функция strip(string) удаляет первый и последний пробелы, символы табуляции и конца строки из строчки. Функция split(string) возвращает список слов в строке; в их число входят и "слова", состоящие из пробелов. Разделитель слов может быть превращен в любую произвольную строку с помощью команды split(string, separator).

Извлечь числа из строк, можно с помощью функций atoi(string) (возвращает целое число) и atof(string) (возвращает вещественное число).

Для того чтобы найти какую-либо определенную строку, воспользуйтесь функцией find(string, text). Она возвращает порядковый номер первой строки string, в которой встречается текст text, или -1, если таких строк нет вообще.

Пример: Следующая программа считывает файл и распечатывает сумму всех чисел во второй колонке.

from Scientific.IO.TextFile import TextFile
import string

sum = 0.
for line in TextFile('data'):
    sum = sum + string.atof(string.split(line)[1])

print "The sum is: ", sum

Пример:Следующая программа выводит на экран все учетные записи пользователей на вашем компьютере:

from Scientific.IO.TextFile import TextFile
from string import split

for line in TextFile('/etc/passwd'):
    print split(line, ':')[0]

Преобразование данных в строку

Любой объект Python (числа, строки, векторы ...) может быть преобразован в строку. Для этого его достаточно записать между значками обратных апострофов:

from Scientific.Geometry import Vector

a = 42
b = 1./3.
c = Vector(0, 2, 1)

print `a` + ' ' + `b` + ' ' + `c`
Программа выведет на экран "42 0.333333333333 Vector(0,2,1)".

Другой путь преобразовать что-либо в строку - это использование функции str(data) (она является не частью модуля string, а входит в базовый язык). Оба рассмотренных метода не всегда дают один и тот же результат, хотя для чисел все получается одинаково. Вообще, str(data) дает более "приятное" представление заданной величины, тогда как обратные апострофы порождают представление, в котором отражается не только значение, но и тип данных. Например, если s - это строка, то str(s) - будет тоже s, тогда как `s` возвращает s, заключенное в апострофы, чтобы показать, что данные являются строкой. Практически, стоит попробовать оба варианта и использовать тот, который больше нравиться.

Функция join(words) обрабатывает список строк и возвращает их объединение; при этом слова разделяются пробелами. Поэтому последнюю строчку в приведенной выше программе можно было бы записать проще: print string.join(`a`, `b`, `c`). Разделитель слов, опять-таки, может быть заменен произвольной строкой.

Функции lower(string) и upper(string) переводят символы строки в нижний или верхний регистр, соответственно.

Функция ljust(string, width) возвращает строку string длиной по крайней мере width символов, выровненную по левому краю. Функции rjust и center работают сходным образом, только строка выравнивается по правому краю или по центру.

Некоторые не описанные здесь полезные функции

Python включает в себя очень большое количество разнообразных функций для работы с более или менее специализированными формами текста. Дать полное описание или даже просто перечислить их здесь невозможно. Всю необходимую информацию можно найти в Python Library Reference.

Прежде всего, много не упомянутых в этом документе функций содержится в модуле string.

Существует важный набор функций, называемых регулярными выражениями, специально предназначенных для поиска и изменения данных с помощью шаблонов. Эти функции расположены в модуле re. Они очень мощные, но синтаксис регулярных выражений (используемый также в таких инструментах Unics, как grep, и таких редакторах, как emacs) несколько сложен.

Модуль htmllib содержит функции для извлечения данных из HTML-файлов, которые традиционно используются в сети World-Wide Web. Модуль formatter предоставляет возможность генерировать HTML-файлы.

Файлы в формате языка Fortran

Программы на языке Fortran используют текстовые файлы, являющиеся эмуляциями пачки перфокарт; поэтому для них действуют иные соглашения о формате. Элементы данных идентифицируются скорее по их расположению в строке, чем по таким разделителям, как пробелы. Формат строки определяется спецификацией формата. Например, 2I4,3F12.5 означает два 4-символьных поля, содержащих целые числа, за которыми следуют три 12-символьных поля, содержащих вещественные числа с пятью знаками после запятой.

С помощью модуля Scientific.IO.FortranFormat возможно проводить преобразование между объектами данных Python и текстовыми строками в формате Fortran. Первый шаг - это создание объекта-формата, представляющего собой спецификацию формата Fortran. Второй шаг - создание из списка данных строки в формате Fortran (для вывода) или обратная операция (для ввода).

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

from Scientific.IO.TextFile import TextFile
from Scientific.IO.FortranFormat import FortranFormat, FortranLine
from Scientific.Geometry import Vector

generic_format = FortranFormat('A6')

atom_format = FortranFormat('A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,' +
                            '3X,3F8.3,2F6.2,7X,A4,2A2')
# Содержимое полей:
# тип записи, порядковый номер, название атома, индекс положения в
# периодической таблице (alternate location indicator)
# название остатка, порядковый номер цепочки, номер остаточной
# последовательности (residue sequence number), код замещения (insertion code),
# координата x, координата y, координата z, размещение (occupancy),
# температурный множитель,
# порядковый номер сегмента, обозначение элемента, заряд

for line in TextFile('protein.pdb'):
    record_type = FortranLine(line, generic_format)[0]
    if record_type == 'ATOM  ' or record_type == 'HETATM':
        data = FortranLine(line, atom_format)
        atom_name = data[2]
        position = Vector(data[8:11])
        print "Atom ", atom_name, " at position ", position

В следующем примере показано, как записывать данные в формате, совместимом с Fortran. В файле на выходе в первой колонке будет содержаться последовательность чисел, а во второй колонке - квадратные корни из тех же самых чисел.

from Scientific.IO.TextFile import TextFile
from Scientific.IO.FortranFormat import FortranFormat, FortranLine
from Numeric import sqrt

format = FortranFormat('2F12.5')

file = TextFile('sqrt.data', 'w')
for n in range(100):
    x = n/10.
    file.write(str(FortranLine([x, sqrt(x)], format)) + '\n')
file.close()

Упражнения


Часть 3: Более сложные вычисления

Отображения (mappings)

В отображениях хранятся связи между множеством ключей и множеством значений. Их можно рассматривать как обобщение последовательностей, поскольку последовательности подобны отображениям, с ключами в виде некоторого диапазона целых чисел. Отображения допускают ключи более общего характера, но при этом порядок элементов не задается. Кроме того, между последовательностями и отображениями имеются и другие существенные различия; так что их всегда следует отличать.

Доступ к отображениям осуществляется почти также как к последовательностям: mapping[key] возвращает значение, связанное с данным ключом, а mapping[key]=value изменяет его. Список всех ключей можно вывести с помощью команды mapping.keys(), а список всех значений - с помощью mapping.values(). При использовании выражения mapping.items() формируется список кортежей (key, value).

Словари

Наиболее часто используемые отображения (и единственные, рассматриваемые в этом курсе) - это словари. Ключом в них может быть любой неизменяемый объект, т.е. любой объект, который не изменяется в течение всего цикла своего существования. Это значит, что списки, словари и массивы не могут использоваться в качестве ключей, а целые числа, кортежи и строки (наряду со многими другими структурами данных) - могут. Значениями могут быть произвольные объекты.

В следующем примере создается словарь, содержащий атомные массы некоторых химических элементов. Затем добавляется новый пункт и, наконец, вычисляется масса молекулы воды

atomic_mass = {'H': 1., 'C': 12, 'S': 32}
atomic_mass['O'] = 16.
print atomic_mass['O'] + 2*atomic_mass['H']

Элемент словаря можно удалить командой del dictionary[key].

Массивы

Другой вид объектов-последовательностей, специально оптимизированный для использования в численных приложениях, - это массивы. Массивы и функции для работы с ними сосредоточены в модуле Numeric, куда входят и математические функции, представленные выше. Хорошая идея: начинать каждый интерактивный сеанс вычислений с команды from Numeric import *.

У массивов есть два свойства, отличающие их от других последовательностей: они могут быть многомерными, и их элементы все одного и того же типа, либо целые, либо вещественные, либо комплексные числа. (Существуют массивы символов, вещественных чисел с одинарной точностью и обычных объектов Python, но говорить о них в данном курсе мы не будем.) Размерность массива не может быть больше сорока.

Массивы создаются из других последовательностей с помощью функции Numeric.array. Многомерные массивы создаются из вложенных последовательностей. Необязательный второй аргумент обозначает тип массива (Integer, Float или Complex). Если тип не обозначен явно, подразумевается, что он совпадает с типом данных.

В примере, который приводится ниже, создается одномерный массив, состоящий из целых чисел, и трехмерный массив вещественных чисел:

from Numeric import *

integer_array = array([2, 3, 5, 7])
print integer_array.shape

real_array = array([ [ [0., 1.],
                       [2., 3.] ],
                     [ [4., 5.],
                       [6., 7.] ],
                     [ [8., 9.],
                       [10., 11.] ] ])
print real_array.shape
Шаблон массива, выводимый на экран в приведенном выше примере, представляет собой кортеж, состоящий из длин размерностей массива, (4,) для первого массива и (3, 2, 2) - для второго. Число измерений, т.е. длина шаблона, называется рангом массива. Стандартные числовые объекты (целые, вещественные и комплексные числа) могут рассматриваться как массивы нулевого ранга.

В модуле Numeric имеется несколько функций, с помощью которых можно создавать массивы особого рода. Функция zeros(shape, type) возвращает массив определенного шаблона (кортеж), все элементы которого равны нулю. Второй аргумент задает тип (Integer, Float или Complex) и является необязательным; по умолчанию - это Integer. Функция ones(shape,type) действует аналогичным образом и делает все элементы равными единице. Функция arrayrange(first, last, step) работает сходно с range, но возвращает массив (первого ранга) и допускает не только целочисленные, но и вещественные аргументы.

Массивы поддерживают более развитую форму индексирования, используемую для последовательностей. Кроме форм a[i] (один элемент) и a[i:j] (подмножество массива от i до j) существует форма a[i:j:k], расширяющая подмножество массива от i до j с шагом k. В случае многомерных массивов индексы/диапазоны отделяются друг от друга запятыми. Если число индексов превышает размерность, индексы присваиваются, начиная слева. Индексированным массивам можно присваивать значения, использовать их в списках.

Вот несколько примеров индексирования массивов:

import Numeric

# Создать массив с нулевыми элементами
a = Numeric.zeros((4, 2, 3), Numeric.Float)

# Присвоить заданному элементу значение 1
a[2, 1, 0] = 1.

# Присвоить значение 2.5 всем элементам с первым индексом 0
a[0] = 2.5

# Распечатать все элементы с первым индексом, равным 0, и
# последним индексом, равным 1
print a[0,:,1]

Существует два особых индекса, которые не указывают непосредственно на какие-либо элементы. Особый индекс ... (многоточие) "не обращает внимания" на размерность, так что присвоение индексов, следующих за ним, начинается справа. Например, выражение a[...,0] выделяет все элементы, последний индекс которых равен 0; размерность при этом может быть любой. Особый индекс NewAxis (определенный в модуле Numeric) добавляет новое измерение единичной длины. Например, если a имеет шаблон (2, 2), тогда a[:, Numeric.NewAxis, :] будет иметь шаблон (2, 1, 2) и тот же набор элементов, что и a.

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

Стандартные арифметические операции (сложение, умножение и т.д.), а также математические функции из модуля Numeric могут применяться к массивам, т.е. к каждому элементу данного массива. В случае бинарных операций (сложение и т.п.) шаблоны массивов должны соответствовать друг другу. Однако, "соответствие" не означает "равенство". Допустимо, например, умножить весь массив на какое-либо число или прибавить строку ко всем строкам матрицы. Строгая формулировка соответствия: сравните шаблоны двух массивов поэлементно, начиная справа, до тех пор, пока наименьший массив не закончится; если для всех размерностей длины равны, или одна из них равна единице, то массивы соответствуют друг другу. Менее формальное описание: размерности единичной длины повторяются в соответствии с размерностью другого массива.

Пример:

a = array([[1, 2]])               # шаблон (1, 2)
b = array([[10], [20], [30]])     # шаблон (3, 1)
После команды a+b a повторяется трижды вдоль первой оси, а b повторяется дважды вдоль второй. В итоге получается
a = array([[1, 2], [1, 2], [1, 2]])
b = array([[10, 10], [20, 20], [30, 30]])
Итак, теперь оба массива имеют одинаковый шаблон, и a+b будет равно array([[11, 12], [21, 22], [31, 32]]). Конечно же, массивы копируются не физически, так что опасности исчерпать память нет.

Бинарные операции существуют также в виде функций (в модуле Numeric): add(a, b) (сложить) эквивалентно a+b, также subtract (вычесть), multiply (умножить), divide (разделить) и power (возвести в степень) могут быть использованы вместо соответствующих бинарных операторов. Некоторые бинарные операции существуют только в виде функций:
maximum(a, b), minimum(a, b) наибольшее/наименьшее из a и b
equal(a, b), not_equal(a, b) проверка на равенство (возвращает 0/1)
greater(a, b), greater_equal(a, b) сравнение
less(a, b), less_equal(a, b) сравнение
logical_and(a, b), logical_or(a,b) логическое и/или
logical_not(a) логическое отрицание

Бинарные операции из этого списка могут применяться к комбинациям элементов одиночного массива. Например, add.reduce(a) вычисляет сумму всех элементов a вдоль первого измерения, а minimum.reduce(a) возвращает наименьший элемент массива a. Необязательный второй аргумент задает размерность явным образом. Он может быть положительным (0 = первое измерение) или отрицательным (-1 = последнее измерение). Другой вариант - функция accumulate: add.accumulate(a) возвращает массив, содержащий первый элемент массива а, сумму первых двух элементов, сумму первых трех элементов и т.д. Последний элемент поэтому равен add.reduce(a).

Еще один способ сформулировать операцию с помощью бинарных функций - использование outer. Функция add.outer(a, b) возвращает массив с комбинированной размерностью массивов a и b, элементы которого представляют все возможные комбинации сумм элементов a и b.

Дополнительные операции над массивами будут описаны ниже.

Функции

Вы уже использовали много функций Python и, возможно, захотели написать свою собственную. Функции дают возможность один раз запрограммировать (и протестировать) определенный код, а затем снова и снова его использовать. В настоящем курсе рассматриваются только функции Python; также возможно определять функции на С или каком-либо совместимом языке низкого уровня.

В следующем примере определяется функция distance, в которой вычисляется расстояние между двумя точками в пространстве. Аргументами функции являются два произвольных вектора:

from Scientific.Geometry import Vector

def distance(r1, r2):
    return (r1-r2).length()

print distance(Vector(1., 0., 0.), Vector(3., -2., 1.))

Легко заметить, что вызов функции состоит из трех этапов:

  1. Аргументы, указанные при вызове функции, присваиваются соответствующим переменным в описании функции.
  2. Исполняется код функции.
  3. Первоначальный вызов функции заменяется значением выражения return.

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

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

import Numeric

def cartesianToPolar(x, y):
    r = Numeric.sqrt(x**2+y**2)
    phi = Numeric.arccos(x/r)
    if y < 0.:
        phi = 2.*Numeric.pi-phi
    return r, phi

radius, angle = cartesianToPolar(1., 1.)
print radius
print angle

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

import Numeric

def printTable(function, first_value, last_value, step):
    for x in Numeric.arrayrange(first_value, last_value, step):
        print x, function(x)


def someFunction(x):
    return Numeric.sin(x**2)

printTable(someFunction, 0., 5., 0.1)

Модули

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

Для Python нет никакой разницы, с чем он имеет дело: со стандартной библиотекой или любым другим модулем. Просто стандартная библиотека поставляется вместе с компилятором. Механизмы импортирования всегда одни и те же, независимо от происхождения модулей.

По сути дела, модуль - это файл с описаниями (переменных, функций и т.д.). Название такого файла обязательно должно заканчиваться на .py; все остальное составляет название модуля. К примеру, команда import Numeric осуществляет поиск файла Numeric.py и запускает его на исполнение, делая содержащиеся в нем описания доступными импортирующей программе.

Однако, Python не ищет модули по всей файловой системе; это было бы неэффективно и опасно. Вместо этого используется стратегия, напоминающая метод поиска программ в системе Unix: перед запуском Python следует задать так называемую переменную окружения PYTHONPATH; ее значение представляет собой перечень названий каталогов, отделенных друг от друга двоеточием. Python сначала ищет модули в этих каталогах в указанном порядке; если модуль не найден, Python обратится к своей стандартной библиотеке.

Названия некоторых модулей, из тех, которые нам уже встречались, содержат одну или несколько точек. Это модули, являющиеся частью пакетов. Пакет - это модуль, включающий в себя другие модули. Пакеты используются для структурирования библиотек и для избежания конфликтов имен, когда библиотеки, написанные разными людьми, получают одинаковые названия. До сих пор единственным пакетом, с которым мы имели дело, был Scientific, включающий модули Geometry и IO (плюс остальные). Последние модули, в свою очередь, тоже являются пакетами; IO, например, содержит модули TextFile и FortranFormat, использовавшиеся в этом курсе. Полное описание пакета Scientific содержится в соответствующем руководстве.

Программы с параметрами

Иногда возникает необходимость написать на Python программу, входные параметры которой передаются через командную строку оболочки. Запуск такой программы происходит следующим образом:

python my_program.py a b c
Все, что нужно сделать в данном случае,- это определить способ доступа к параметрам командной строки. В модуле sys имеется переменная argv, значением которой является список всех параметров, начинающийся именем запускаемой Python-программы. Так, в приведенном выше примере sys.argv будет выглядеть как ['my_program.py', 'a', 'b', 'c'].

В Unix-системах программы на Python вообще можно сделать исполняемыми напрямую. Добиться этого можно в два приема: сперва введите в первой строке вашей Python-программы

#!/usr/local/bin/python
Возможно, путь придется изменить, если Python установлен в вашей системе где-то в другом месте. Затем сделайте этот файл исполняемым с помощью Unix-утилиты chmod.

Волшебную строчку, приведенную выше, Python проигнорирует, поскольку это комментарий (строка начинается с #). Но она распознается оболочками Unix, показывая, какую программу следует запустить для выполнения заданного файла.

Статистика

Модуль Scientific.Statistics содержит набор элементарных функций статистического анализа.

Функции average(data) (среднее), variance(data) (дисперсия) и standardDeviation(data) (среднеквадратическое отклонение) обрабатывают последовательность данных и возвращают статистическую величину, соответствующую их названию.

Подмодуль Histogram определяет объект-гистограмму. Histogram(данные, число_столбцов) обрабатывает последовательность данных и рассчитывает гистограмму, разделяя весь диапазон данных на заданное число столбцов; при этом вычисляется, какое количество данных попадает в каждый столбец. Необязательный третий аргумент - это кортеж, задающий диапазон, подлежащий анализу (нижнюю и верхнюю границу). Гистограмма - это массив второго ранга, представляющий собой последовательность столбцов, каждый из которых характеризуется координатой своей середины и высотой.

Построение графиков

Пока для Python не существует модуля для построения графиков с полным набором функций. На настоящий момент для этого применяются внешние программы. Есть два модуля, содержащие функцию рисования простых графиков: Scientific.Mathematica и Gnuplot. В них используются внешние программы, которые называются так же, как соответствующий модуль. В указанных модулях имеется функция plot, обе версии которой ведут себя схожим образом.

Функция plot может иметь дело с произвольным числом аргументов, каждый из которых задает некоторый массив данных. Все массивы данных будут выведены на одном графике. Массив данных представляет собой последовательность точек. Точка определяется двухэлементной последовательностью (координаты x и y) или одним числом, которое интерпретируется как координата y, в то время как координатой x будет порядковый номер точки в данной последовательности.

У функции plot есть необязательный аргумент в виде ключевого слова в форме file='something.ps'. Когда этот аргумент присутствует, данный график будет записан в определенный файл в формате PostScript.

Пример:

from Gnuplot import plot
from Scientific.Statistics.Histogram import Histogram
from Numeric import arrayrange, sqrt

plot(Histogram(sqrt(arrayrange(100.)), 10))

Упражнения


Часть 4: Высшая математика

Матрицы и линейная алгебра

Матрицы представляются в виде массивов второго ранга, в то время как векторы (в более общем смысле линейной алгебры, а не геометрические векторы) описываются массивами первого ранга. Модуль Numeric содержит несколько весьма полезных операций для работы с матрицами и векторами:

Модуль LinearAlgebra включает в себя набор обычных операций над матрицами. Матрицы представляются в виде двумерных массивов. Некоторые операции относятся только к квадратным матрицам; в этом случае производится соответствующая проверка входных данных. Все операции могут применяться как для вещественных, так и для комплексных матриц. Эти функции базируются на LAPACK-процедурах, тщательно тестировавшихся и эффективных.

Функция solve_linear_equations(a, b) решает систему линейных уравнений с квадратной невырожденной матрицей a и правосторонним вектором b. Несколько правосторонних векторов можно обрабатывать одновременно, сделав b двумерным массивом (т.е. последовательностью векторов). Функция inverse(a) вычисляет матрицу, обратную квадратной невырожденной матрице a, обращаясь к solve_linear_equations(a, b) с соответствующим b.

Функция eigenvalues(a) возвращает собственные значения квадратной матрицы a. Функция eigenvectors(a) возвращает и собственные значения, и собственные векторы, последние в виде двумерного массива (т.е. последовательности векторов).

Функция determinant(a) возвращает детерминант квадратной матрицы a.

Функция singular_value_decomposition(a) возвращает три массива V, S и WT, матричное произведение которых является первоначальной матрицей a. V и WT - это единичные матрицы (массивы второго ранга), в то время как S - вектор (массив первого ранга) диагональных элементов вырожденной матрицы (singular-value matrix). Эта функция применяется, главным образом, для проверки, является ли (и насколько) матрица плохо обусловленной.

Функция generalized_inverse(a) производит обобщенное обращение (также известное как псевдо-обращение или обращение Мура-Пенроуза) матрицы a. Область ее применения очень широка (линейные уравнения, метод наименьших квадратов).

Функция linear_least_squares(a, b) возвращает решение переопределенной системы линейных уравнений, полученное по методу наименьших квадратов. Необязательный третий аргумент фиксирует диапазон отбрасываемых несобственных значений (по умолчанию до 10-10). Возвращается четыре значения: само решение по методу наименьших квадратов, сумма квадратов разностей (т.е. минимизируемая в соответствии с методом наименьших квадратов величина), ранг матрицы a и несобственные значения a в порядке убывания.

Пример:

import Numeric; N = Numeric
import LinearAlgebra; LA = LinearAlgebra

# Создать матрицу
a = N.reshape(N.arrayrange(25.), (5, 5)) + N.identity(5)

# Вычислить обратную матрицу
inv_a = LA.inverse(a)

# Проверка операции обращения матрицы путем вывода на экран абсолютной
# величины наибольшего элемента a * a^{-1} - identity(5)
print "Inversion error:", \
      N.maximum.reduce(N.fabs(N.ravel(N.dot(a, inv_a)-N.identity(5))))

# Вывести собственные значения матрицы
print "Eigenvalues: ", LA.eigenvalues(a)
Обратите внимание, как можно создавать сокращенные обозначения для названий модулей. Модули - это объекты данных особого рода, поэтому их возможно использовать как значение переменной.

Преобразование Фурье

Модуль FFT содержит функции, с помощью которых можно выполнять быстрое преобразование Фурье (в различных модификациях) в одном или в двух измерениях с помощью пакета FFTPACK. Все функции могут применяться к массивам любого ранга; число размерностей, подвергающихся преобразованию, может быть определено и задано по умолчанию равным единице или двум.

Функция fft(a) производит стандартное комплексное одномерное быстрое преобразование Фурье, а функция inverse_fft(a) выполняет обратное преобразование. Обе функции имеют необязательный второй аргумент, определяющий длину последовательности данных (по умолчанию - вся длина массива); если его значение больше, чем сам массив, будет добавлено соответствующее число нулей. Необязательный третий аргумент задает размерность, по которой будет проводиться преобразование.

Функция fft2d(a) реализует двумерное быстрое преобразование Фурье. Необязательный второй аргумент определяет шаблон (кортеж длиной два) преобразуемого подмассива; по умолчанию преобразованию подвергается весь массив. Необязательный третий аргумент (также кортеж длиной два) задает две размерности преобразования.


Подбор параметров аппроксимирующей функции

Модуль LeastSquares содержит функцию leastSquaresFit(model, parameter, data), посредством которой можно проводить аппроксимацию по общему нелинейному методу наименьших квадратов с помощью алгоритма Левенберга-Марквардта.

Аргумент model задает аппроксимирующую функцию. Она будет вызываться с двумя аргументами: первый - это кортеж, содержащий параметры аппроксимации; второй представляет собой первый элемент данных (см. ниже). Значение, возвращаемое функцией, должно быть числом.

Аргумент parameter - это кортеж начальных значений параметров аппроксимации. Скорость сходимости метода зависит от качества начальных значений, поэтому их следует выбирать с осторожностью.

Аргумент data - это последовательность экспериментальных точек, для которых подбирается модель. Каждая точка данных - это последовательность длиной 2 или 3. Ее первый элемент определяет независимые переменные модели. Он передается функции модели в качестве ее первого параметра, но больше никак не используется. Второй элемент кортежа каждой точки данных - это предполагаемое значение (насколько возможно хорошо совпадающее с истинным), которое должно быть возвращено функцией модели. Третий элемент (по умолчанию равен 1.) - это статистическое отклонение экспериментальной точки, т.е. величина, обратная ее статистическому весу в процедуре аппроксимации.

Пример:

from LeastSquares import leastSquaresFit
import Numeric

# Математическая модель.
def exponential(parameters, x):
    a = parameters[0]
    b = parameters[1]
    return a*Numeric.exp(-b/x)

# Данные, которые необходимо аппроксимировать.
data = [(100, 4.999e-8),
        (200, 5.307e+2),
        (300, 1.289e+6),
        (400, 6.559e+7)]

fit = leastSquaresFit(exponential, (1e13, 4700), data)

print "Fitted parameters:", fit[0]
print "Fit error:", fit[1]

Сеточные функции

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

Функцию можно определить с помощью выражения InterpolatingFunction(grid,value), где grid - последовательность осей, длина которых задает размерность сетки. Аргумент values обязательно должен быть массивом с рангом равным или больше размерности сетки. Если он больше, то функция является многозначной (например векторное поле). Необязательный третий аргумент определяет значение функции за пределами сетки. Если он не задан, попытки оценивать значение функции за пределами сетки будут приводить к ошибке.

К интерполирующей функции можно обращаться как и к любой другой. Число аргументов должно строго соответствовать размерности сетки, причем сами аргументы должны быть целыми или вещественными числами.

Кроме простого определения значения функции интерполирующая функция предоставляет возможность использовать ряд полезных операций:

Пример:

from Scientific.Functions.Interpolation import InterpolatingFunction
import Numeric; N = Numeric

# Задать ось с 10 узлами
axis = N.arrayrange(0, 1.1, 0.1)

# Вычислить значения (здесь просто квадратный корень)
values = N.sqrt(axis)

# Создать интерполирующую функцию
f = InterpolatingFunction((axis,), values)

# Напечатать таблицу значений для сравнения
for x in N.arrayrange(0, 1., 0.13):
    print "%5.2f: %10.5f %10.5f" % (x, N.sqrt(x), f(x))

Визуализация (VRML)

При визуализации сложной информации, особенно в трех измерениях, не всегда возможно обойтись простыми графиками. Модуль Scientific.Visualization.VRML предоставляет возможность непосредственно конструировать трехмерные сцены, состоящие из стандартных объектов (кубы, сферы и т.д.). Такие сцены могут быть записаны в файлы формата VRML (Virtual Reality Modelling Language - язык моделирования виртуальной реальности), которые затем можно просматривать с помощью соответствующей программы просмотра VRML. Подобные программы просмотра VRML существуют в большинстве операционных систем. Более подробную информацию о VRML и средствах просмотра VRML можно почерпнуть в репозитории VRML Суперкомпьютерного центра в Сан-Диего.

Если, вместо того чтобы составлять файл и запускать средство просмотра VRML вручную,требуется вызывать VRML-сцены напрямую из ваших Python-программ, прежде всего, переменной окружения Unix VRMLVIEWER должно быть присвоено значение, соответствующее названию инструмента просмотра VRML, установленного в системе.

Замечание: в этом курсе рассматриваются лишь элементарные приемы использования модулей VRML. Более подробную информацию можно найти в исходном тексте программы.

Пример типичной последовательности действий:

  1. Создать сцену VRML с помощью команды Scene(objects), где objects - это список объектов VRML. В случае пустого списка создается пустая сцена.
  2. Создать объекты VRML (см. ниже) и разместить их на сцене посредством команды scene.addObject(object).
  3. Записать полученную сцену в файл с помощью инструкции scene.writeToFile(filename) или сразу просмотреть ее после команды scene.view().

Объекты VRML задаются данными геометрического характера, например, сфера определяется центром (вектор) и радиусом (число). Кроме того, объекты могут иметь одно или более свойств (properties), влияющих на их отображение. Самым важным свойством, обусловливающим цвет, прозрачность, отражательную способность и т.д., является material (материал). В настоящем курсе будут использоваться лишь диффузионные материалы (diffuse materials) заранее определенных цветов, т.е. непрозрачные неотражающие материалы. Обо всех имеющихся возможностях можно прочитать в любой книге по VRML.

Свойства всегда можно изменять, задаются они ключевыми словами. Материал, к примеру, задается параметром material=the_material. Объект, материал которого не определен, по умолчанию будет состоять из диффузионного материала белого цвета.

Материал создается командой DiffuseMaterial(color), где color - название цвета. Возможные значения: black (черный), white (белый), grey (серый), red (красный), green (зеленый), blue (синий), yellow (желтый), magenta (пурпурный) и cyan (голубой). Все они могут быть модифицированы путем добавления перед соответствующим названием 'light ' - "светлый" или 'dark ' - "темный" (обратите внимание на пробел), например 'light blue' - "светло-синий".

Доступные объекты VRML:

В приведенном ниже примере программа считывает PDB-файл и генерирует сцену, включающую в себя черные сферы для каждого атома С-альфа и красные прямые линии, соединяющие соседние атомы С-альфа.

from Scientific.IO.TextFile import TextFile
from Scientific.IO.FortranFormat import FortranFormat, FortranLine
from Scientific.Geometry import Vector
import string
from Scientific.Visualization import VRML

generic_format = FortranFormat('A6')
atom_format = FortranFormat('A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,' +
                            '3X,3F8.3,2F6.2,7X,A4,2A2')

# Считать PDB-файл и генерировать список всех положений С-альфа
c_alpha_positions = []
for line in TextFile('protein.pdb'):
    record_type = FortranLine(line, generic_format)[0]
    if record_type == 'ATOM  ' or record_type == 'HETATM':
        data = FortranLine(line, atom_format)
        atom_name = string.strip(data[2])
        position = Vector(data[8:11])
        if atom_name == 'CA':
            c_alpha_positions.append(position)

# Создать VRML-сцену, материалы и радиусы сфер, а также материал линий
scene = VRML.Scene([])
c_alpha_material = VRML.DiffuseMaterial('black')
c_alpha_radius = 0.2
line_material = VRML.DiffuseMaterial('red')

# Создать сферы и разместить их на сцене
for point in c_alpha_positions:
    scene.addObject(VRML.Sphere(point, c_alpha_radius,
                                material=c_alpha_material))

# Создать прямые линии и разместить их на сцене
for i in range(len(c_alpha_positions)-1):
    point1 = c_alpha_positions[i]
    point2 = c_alpha_positions[i+1]
    scene.addObject(VRML.Line(point1, point2, material=line_material))

# Просмотреть сцену
scene.view()

Другие научные модули

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


Упражнения


Часть 5: Создание новых типов данных

Что такое тип данных?

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

Python предоставляет возможность создавать новые типы данных, обладающие теми же возможностями, что и встроенные типы. По крайней мере с точки зрения пользователя здесь нет разницы. Действительно, многие типы данных, встречающиеся в настоящем курсе, не являются встроенными в Python: векторы, текстовые файлы, интерполирующие функции и т.д. Также возможно создавать новые типы данных на языке С; примером может служить тип данных array (массив) из модуля Numeric.

В Python поддерживается набор стандартных операций, которые могут использоваться с любым подходящим типом данных. Среди них можно выделить арифметические, индексирование, вызов функции и прочее. Помимо этого с типами данных могут быть связаны произвольные операции в форме методов. Методы похожи на функции, но зависят от заданного объекта данных. Мы уже видели несколько таких примеров: append (добавить) - это метод, определенный для списков, close (закрыть) - метод, определенный для файлов и т.п.

Есть множество причин, побуждающих создавать новые типы данных. Они помогают сохранять текст программы читаемым - намного проще написать (a+b)/2 для двух векторов a и b, чем что-нибудь вроде vector_scale(vector_add(a, b), 0.5), как требуется во многих языках программирования. Определение новых типов также содействует сокращению зависимости разных частей программы друг от друга. Программе, использующей векторы, не обязательно "знать", в каком виде хранится информация о векторах (для этого может быть использован список, кортеж, массив или три отдельные переменные). Таким образом, если метод хранения по каким-либо причинам придется изменить, это не окажет никакого влияния на остальные модули.

Объектно-ориентированное программирование

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

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

Классы

Определение нового типа данных называется классом. Класс определяет все операции над типом в форме методов (стандартные операции, например арифметические, входят в состав методов под особыми именами). Также в нем задаются начальные значения новых объектов.

В нижеследующем примере частично приводится определение класса Vector. Явно показаны только задание начального значения, сложение и вычисление длины, некоторые операции упрощены (за счет потери общности). Полное определение можно найти в исходном коде модуля Scientific.Geometry.

import Numeric

class Vector:

    def __init__(self, x, y, z):
        self.array = Numeric.array([x,y,z])

    def __add__(self, other):
	sum = self.array+other.array
	return Vector(sum[0], sum[1], sum[2])

    def length(self):
	return Numeric.sqrt(Numeric.add.reduce(self.array*self.array))

Методы определяются подобно функциям и ведут себя примерно также. Однако, их первый аргумент не совсем обычен: он обозначает сам объект, вызывающий данный метод. Этот аргумент принято называть self, но вместо этого вы можете использовать любое другое имя. При вызове метода v.length() (подразумевается, что v - это вектор), переменная self получает значение v.

Методы, чьи имена начинаются и заканчиваются символом двойного подчеркивания, имеют особое назначение; обычно их не вызывают явным образом (хотя это и возможно). Наиболее важный специальный метод - это __init__, который вызывается сразу же после создания объекта. Выражение Vector(1., 0., 1.) создает новый объект-вектор и затем вызывает его метод __init__, сохраняющий три координаты в массиве, который присваивается локальной переменной нового объекта.

Арифметические операции также можно применять как методы особого рода. Выражение a+b эквивалентно выражению a.__add__(b), другие операции также имеют свои эквиваленты; более подробно об этом см. Справочник по языку Python.

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

    def __repr__(self):
        return 'Vector(%s,%s,%s)' % (`self.array[0]`,
    				     `self.array[1]`,`self.array[2]`)

Атрибуты

Любой объект в Python может иметь произвольное количество атрибутов, во многом похожих на переменные, за исключением того, что они привязаны к определенному объекту, в то время как переменные привязываются к модулям или функциям. Фактически переменные, определенные в модулях, есть не что иное как атрибуты объектов-модулей. Для доступа к атрибутам всегда используется запись вида object.attribute. Названия методов также являются атрибутами, подобно тому как названия функций могут выступать в роли переменных.

В отличие от многих других объектно-ориентированных языков Python не пытается контролировать доступ к атрибутам. Любой код может использовать и даже изменять любой атрибут любого объекта. Например, вы можете запустить на выполнение следующие команды:

import Numeric
Numeric.sqrt = Numeric.exp
чтобы sqrt вел себя как exp. Очевидно, это не самая удачная идея, но Python не будет защищать вас от последствий собственной глупости. Конечно всегда есть некоторая вероятность случайного изменения атрибута, но на практике это не является проблемой.

"Все есть объект": вселенная Python

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

Есть несколько правил, определяющих, в каком пространстве имен существует то или иное имя. Определения в модулях содержатся в пространстве имен атрибутов объекта-модуля. Определения внутри функций и методов производятся во временном (во время исполнения программы) пространстве имен, код внутри функции также может использовать имена (но не присваивать им значения) модуля, содержащего эту функцию. Определения в классах попадают в пространство имен атрибутов данного класса. Наконец, объект, сконструированный из класса (экземпляр класса), конечно же, имеет свое пространство имен атрибутов, в котором и происходят все присвоения; но когда запрашиваемый атрибут не обнаруживается в этом пространстве имен, его поиск продолжается в пространстве имен класса. Так обычно происходит поиск методов.

Специализированные и расширяющиеся классы

Нередко несколько типов данных имеют что-то общее, например один может быть специализированной версией второго. Например, нормализованные векторы допустимо рассматривать как особый вид векторов. Или несколько типов данных могут совместно использовать многие опреации, но отличаться в конкретных деталях. Скажем, возможно определить типы данных, представляющие скалярные и векторные поля, с общими свойствами (например будут определенными на сетке в обоих случаях), но по ряду специфических операций, таких как "градиент" для скалярных полей и "дивергенция" для векторных полей, будут и отличаться. Оба типа данных можно реализовать в виде специализаций типа данных "поле", описывающего общие черты поведения типов и не предназначенного для непосредственного применения в программах (иногда такой случай называют абстрактным классом).

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

В приведенном ниже примере определяется класс, описывающий координатные оси в пространстве, т.е. векторы единичной длины, базирующийся на классе векторов, определенном выше:

class Direction(Vector):

    def __init__(self, x, y, z):
        Vector.__init__(self, x, y, z)
        self.array = self.array/self.length()
Единственный метод, который пришлось переопределять, - это задание начальных значений, где сейчас производится нормализация вектора. Обратите внимание, что метод задания начальных значений сначала вызывает метод инициализации класса Vector и только потом проводит нормализацию.

Класс Direction наследует все свои операции из класса Vector, таким образом, как будто их код был продублирован в новом классе. В частности, сумма двух направлений будет вектором, а не другим направлением. Для получения нормированного результата также следует переопределить метод __add__. Но поскольку сложение направлений не является столь уж полезной операцией, вряд ли стоит прилагать усилия в этом направлении.

Обработка ошибок

Если происходит ошибка, Python выводит на экран трассировку стека (список всех активных функций на момент возникновения ошибки) и делает останов. Часто это полезно, но не всегда. Возможно, вам хочется поработать над ошибками самостоятельно, например, организовать вывод предупреждения, запрос ввода пользователя или провести какие-то расчеты. Для любой программы Python предоставляет возможность обнаружить специфические условия, вызывающие возникновение ошибки, и работать с ними, как вы считаете нужным.

Для идентификации типа ошибки Python имеет несколько встроенных объектов-ошибок, например ValueError (обозначающий, что значение не подходит для выполнения операции, например отрицательные числа для извлечения квадратного корня) или TypeError (обозначающий несоответствие типа данных, например при нахождении логарифма от строки символов). Программа может фиксировать определенный объект-ошибку, определенную коллекцию или любую другую ошибку.

Общепринятый способ фиксации ошибки:

try:
    x = someFunction(a)
    anotherFunction(x)
except ValueError:
    print "Что-то не так"
else:
    print "Нет ошибок"
Сначала выполняется код, следующий после try:. Если встречается ValueError, то выполняется программный код, следующий за выражением except ValueError:, в противном случае - код после else:; последняя часть не обязательна. Для фиксации сразу нескольких типов ошибок примените кортеж (например, except (ValueError, TypeError)). Для выявления всех ошибок запишите выражение except: без аргументов. Чтобы по-разному обрабатывать разные типы ошибок, используйте команду except ...: code, состоящую из нескольких частей. Более подробно об этом и о других возможностях см. справочное руководство по языку.

Python-программы, безусловно, также могут генерировать ошибки, с помощью выражения raise ErrorObject или raise ErrorObject,"Explanation" (с выводом разъяснений для пользователя). В роли ErrorObject (объекта-ошибки) может выступить любой из стандартных объектов-ошибок, а также строка. Во многих модулях их собственные типы ошибок определяются в виде строк, причем другие модули могут их импортировать:

# Модуль A

AError = "Ошибка в модуле A"

def someFunction(x):
    raise AError
# Модуль B

import A

try:
    A.someFunction(0)
except A.AError:
    print "Что-то не так"

Упражнения

Часть 6: Графические интерфейсы пользователя (GUI)

GUI-программирование и Tk

Объекты GUI

Все современные комплекты инструментальных средств GUI, включая пакет Tk, используемый в Python, основаны на объектно-ориентированной модели пользовательского интерфейса. Наиболее распространенные типы объектов - это окна, поля ввода данных, кнопки, текстовые поля, графические поля, меню и т.д. Конструирование GUI подразумевает создание всех требуемых объектов, их правильное совместное расположение и задание связанных с ними действий.

Работа программы, управляемая событиями

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

Разработка интерфейса пользователя

В настоящем курсе рассматриваются только технические аспекты создания GUI. Столь же важен вопрос компоновки интерфейса пользователя, т.е. что, где и как размещать. По этой теме существует обширная литература, но самый полезный совет будет таким: изучайте хорошо написанные программы и пытайтесь им подражать. Наиболее распространенная ошибка начинающих - это чрезмерное стремление к новизне и оригинальности. Главная цель, однако,- не создание GUI, все функции которого сосредоточены в одном единственном окне, с использованием экзотических приемов, а проектирование интерфейса, в котором смогут быстро разобраться пользователи, основываясь на своем прошлом опыте работы с другими программами.

Tk и Python

Python не имеет своего собственного пакета для создания пользовательского интерфейса. Тем не менее, существует ряд интерфейсов для взаимодействия Python с несколькими такими пакетами. Самый популярный среди них - это Tk, изначально предназначенный для Tcl, очень простого языка сценариев. У Tk есть ряд преимуществ: легкость использования, завершенность и доступность для всех главных операционных систем. Однако, он довольно тесно привязан к Tcl, и несмотря на все усилия со стороны интерфейса Python скрыть этот факт, некоторые свойства Tk в сочетании с языком Python выглядят бессмысленными. Инструментарий GUI, созданный специально под Python, скорее всего, выглядел бы по-другому.

Взгляд на мир с точки зрения Tk

Tk знает о двух типах объектов: окнах и элементах управления (widgets - от "window" (окно) и "gadget" (приспособление)). Элементы управления - это объекты пользовательского интерфейса, "населяющие" окна: кнопки, текстовые поля и т.п. У элементов управления - иерархическая организация; у каждого элемента управления - свой родитель (другой элемент управления или, в случае элемента самого верхнего уровня,- окно). От родительского элемента зависит положение данного элемента управления на экране и многие его свойства. Окна также являются частью иерархии, но несколько иного рода: положение и содержимое окна совершенно не зависят от окна-родителя, но при закрытии последнего оно исчезает. Вершину иерархии окон занимает т.н. корневое окно, у которого родительского окна нет. При закрытии корневого окна основной цикл событий Tk прерывается.

Типичное приложение Tk начинает свою работу с конструирования элемента управления, присоединенного к корневому окну, и запуска основного цикла. Элементы управления задают функции, которые будут вызываться в ответ на те или иные события (действия пользователя). Вызов функции осуществляется основным циклом Tk после того, как произошло определенное событие. Эти функции могут вообще ничего не делать; какие-либо ограничения отсутствуют. Они могут открывать и закрывать окна и даже изменять сами элементы управления в существующих окнах. Кроме того, они могут завершить основной цикл Tk.

После создания элемента управления его создатель должен решить, где внутри родительского элемента он будет расположен. Tk предлагает на выбор две стратегии размещения элементов управления. Диспетчер компоновки place обеспечивает очень гибкую, но также и весьма сложную стратегию. С другой стороны, диспетчер pack (упаковщик) проще использовать, т.к. большую часть работы он проделывает сам. Это преимущество достигается за счет некоторой потери гибкости. В большинстве Tk-программ, как и во всех примерах данного курса, применяется только диспетчер компоновки pack.

Очень простой пример

Приводимая ниже программа проста настолько, насколько может быть простым приложение Tk. Она открывает окно, содержащее лишь одну кнопку. После нажатия кнопки на экран выводится сообщение. Программа прекращает свою работу после закрытия окна.

from Tkinter import *

def report():
    print "The button has been pressed"

object = Button(None, text = 'A trivial example', command = report)
object.pack()
object.mainloop()

Первое выражение импортирует все имена из модуля Tkinter. Из-за большого количества определений, содержащихся в этом модуле, импорт имен по отдельности крайне неудобен. Затем определяется функция реакции на событие, в роли которой может быть любая функция Python. В следующей строке создается элемент управления в виде кнопки. Его первый параметр - родительский элемент этого элемента управления; None означает корневое окно. Остальные параметры описывают свойства кнопки: текст на кнопке и команду, которая будет выполнена после нажатия на кнопку. Потом вызывается метод элемента управления pack(), означающий, что для размещения данного элемента управления будет использоваться диспетчер компоновки pack. Никаких указаний, как следует располагать объекты, нет, поскольку элемент управления один. Наконец, запускается основной цикл.

Расположение элементов управления

В окне Tk может располагаться только один элемент управления. Конечно, одноэлементные окна не слишком-то полезны, и существует способ сделать их более сложными. Речь идет о фреймах. Фрейм - это элемент управления, своеобразное дочернее окно. Обычно оно невидимо (хотя можно сделать видимой его внешнюю рамку), но содержит в себе другие элементы управления, как правило, видимые. Фрейм может включать в себя любое количество произвольных элементов управления.

Если в состав одного фрейма входит несколько элементов управления, должны быть заданы некие инструкции, определяющие их совместное расположение. Эта информация содержится в необязательных аргументах метода pack. Пример простейшего набора инструций: side=some_side с возможностью выбора TOP (верх), BOTTOM (низ), LEFT (слева) или RIGHT (справа). Диспетчер pack старается разместить данный элемент управления как можно ближе к указанной стороне фрейма в рамках ограничений, налагаемых наличием других элементов управления. Если, к примеру, задать side=TOP для всех элементов управления в фрейме, они выстроятся в цепочку в вертикальном направлении, будучи выровнены по центру в горизонтальном; первым будет элемент управления, для которого вызов pack иерархически самый старший.

Еще более сложные конструкции можно создать, используя фреймы внутри фреймов. Директивы для диспетчера pack, относящиеся к элементу управления, влияют на его расположение внутри родительского фрейма, у которого, в свою очередь, есть свои опции диспетчера pack для размещения его внутри другого фрейма.

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

from Tkinter import *

def report():
    print "A button has been pressed"

outer_frame = Frame(None)
outer_frame.pack()

first_line = Frame(outer_frame)
first_line.pack(side=TOP)
button1 = Button(first_line, text='Button 1', command=report)
button1.pack(side=LEFT)
button2 = Button(first_line, text='Button 2', command=report)
button2.pack(side=LEFT)

second_line = Frame(outer_frame)
second_line.pack(side=TOP)
button3 = Button(second_line, text='Button 3', command=report)
button3.pack(side=LEFT)
button4 = Button(second_line, text='Button 4', command=report)
button4.pack(side=LEFT)

outer_frame.mainloop()

Обратите внимание, что в роли родителя первого фрейма выступает None (никакой), т.е. корневое окно. Внутренние фреймы, представляющие прямые линии, в качестве родителя имеют внешний фрейм (outer frame), а кнопки присоединяются к внутренним фреймам.

Создание собственных элементов управления

В принципе любой GUI можно запрограммировать как в последнем примере, конструируя один объект за другим и задавая опции диспетчера pack. Однако, такая программа быстро становится громоздкой и нечитаемой, а ее модификация - почти невозможной. Объекты графического интерфейса пользователя должны быть тем или иным способом структурированы. Как правило этого добиваются за счет введения специально предназначенных для данного конкретного приложения элементов управления более высокого уровня. Скажем, в вышеприведенном примере окажется полезным элемент управления "ряд кнопок".

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

С элементом управления "ряд кнопок" пример с четырьмя кнопками может быть изменен следующим образом:

from Tkinter import *

def report():
    print "A button has been pressed"

class ButtonLine(Frame):

    def __init__(self, master, button_data):
	Frame.__init__(self, master)
	for text, command in button_data:
	    Button(self, text=text, command=command).pack(side=LEFT)

outer_frame = Frame(None)
outer_frame.pack()

first_line = ButtonLine(outer_frame, [('Button 1', report),
				      ('Button 2', report)])
first_line.pack(side=TOP)

second_line = ButtonLine(outer_frame, [('Button 3', report),
				       ('Button 4', report)])
second_line.pack(side=TOP)

outer_frame.mainloop()

Элементы управления Tk: обзор

Здесь нет возможности приводить подробное описание всех элементов управления, но их краткий обзор с указанием обычной сферы применения будет хорошим вступлением перед дальнейшим знакомством с руководством по Tk.

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

Более подробную информацию об этих элементах управления можно найти в документации Tkinter и в справке Tk в формате man.

Не столь простой пример

Следующая программа выводит на экран небольшое окно, в котором пользователь может напечатать имя файла и потом вывести текст либо в окно, либо на печать. Вместо того чтобы набирать имя файла, можно использовать средство для просмотра файлов (файловый браузер). На рисунке показано главное окно, текстовое окно и файловый браузер:

[screen snapshot]

В программе используются три элемента управления, которые относятся не к стандартным элементам Tk, а к стандартной библиотеке Python: "диалоговое окно" для вывода сообщений об ошибке, файловый браузер и текстовый элемент управления с линейками прокрутки. Также в программе определяются два собственных элемента управления: элемент управления "ввод имени файла", содержащий метку, поле ввода для имени файла и кнопку запуска файлового браузера, а также "панель кнопок", содержащая группу кнопок, выравненных по левому краю, и группу кнопок, выравненных по правому краю. Также с помощью классов определяются два полнооконных элемента управления, один для текстового окна и другой для главного.

В программе применяется опция упаковщика, о которой еще не упоминалось. Как правило, все элементы управления имеют фиксированный размер, который зависит от их содержимого. Для некоторых элементов управления, тем не менее, не существует явного корректного размера, например в случае полей ввода. Поля ввода могут принимать произвольный размер, устанавливаемый упаковщиком. Опция fill=X предписывает упаковщику увеличить данный элемент управления в направлении x, так чтобы он заполнил все пространство, доступное в родительском элементе управления.

Обратите внимание, что действия, связанные с кнопками, определяются как методы элементов управления. Это нужно, для того чтобы можно было оперировать атрибутами элемента управления, что часто требуется на практике. В реальных Tk-программах почти все действия определяются через методы элементов управления.

from Tkinter import *
from Dialog import Dialog
from FileDialog import LoadFileDialog
from ScrolledText import ScrolledText

class FilenameEntry(Frame):

    def __init__(self, master, text):
	Frame.__init__(self, master)
	Label(self, text=text).pack(side=LEFT)
	self.filename = StringVar()
	Entry(self, textvariable=self.filename).pack(side=LEFT, fill=X)
	Button(self, text="Browse...", command=self.browse).pack(side=RIGHT)

    def browse(self):
	file = LoadFileDialog(self).go(pattern='*')
	if file:
	    self.filename.set(file)

    def get(self):
	return self.filename.get()


class ButtonBar(Frame):

    def __init__(self, master, left_button_list, right_button_list):
	Frame.__init__(self, master, bd=2, relief=SUNKEN)
	for button, action in left_button_list:
	    Button(self, text=button, command=action).pack(side=LEFT)
	for button, action in right_button_list:
	    Button(self, text=button, command=action).pack(side=RIGHT)


class FileNotFoundMessage(Dialog):

    def __init__(self, master, filename):
	Dialog.__init__(self, master, title = 'File not found',
			text = 'File ' + filename + ' does not exist',
			bitmap = 'warning', default = 0,
			strings = ('Cancel',))


class TextWindow(Frame):

    def __init__(self, master, text):
	Frame.__init__(self, master)
	text_field = ScrolledText(self)
	text_field.insert(At(0,0), text)
	text_field.pack(side=TOP)
	text_field.config(state=DISABLED)
	ButtonBar(self, [],
		  [('Close', self.master.destroy)]).pack(side=BOTTOM, fill=X)


class MainWindow(Frame):

    def __init__(self, master):
	Frame.__init__(self, master)
	Label(self, text="Enter a filename and " +
                         "select an action").pack(side=TOP)
	self.filename_field = FilenameEntry(self, "Filename: ")
	self.filename_field.pack(side=TOP, fill=X)
	ButtonBar(self, [('Show', self.show), ('Print', self.lpr)],
		  [('Quit', self.quit)]).pack(side=BOTTOM, fill=X)

    def show(self):
	filename = self.filename_field.get()
	try:
	    text = open(filename).read()
	except IOError:
	    FileNotFoundMessage(self, filename)
	else:
	    new_window = Toplevel()
	    new_window.title(filename)
	    TextWindow(new_window, text).pack()

    def lpr(self):
	filename = self.filename_field.get()
	import os
	os.system('lpr ' + filename)


mw = MainWindow(None)
mw.pack()
mw.mainloop()

Дальнейшее изучение Tk

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


Упражнения

Решение упражнений, часть 1

Вывод на экран центра масс:

from Scientific.Geometry import Vector

a = Vector(0, 1, 0)
b = Vector(4.3, -2.4, 0.005)
c = Vector(-3.2, 5.1, -3.)

print (a+b+c)/3
Правильным, но крайне неудачным решением будет сложение всех трех компонент вектора по отдельности.

Вывод на экран углов:

from Scientific.Geometry import Vector

a = Vector(0, 1, 0)
b = Vector(4.3, -2.4, 0.005)
c = Vector(-3.2, 5.1, -3.)

print "Angle a-b-c:", (a-b).angle(c-b)
print "Angle c-a-b:", (c-a).angle(b-a)
print "Angle b-c-a:", (b-c).angle(a-c)
Проверить работу программы очень просто: для этого достаточно просуммировать все углы; сумма должна равняться числу пи.

Решение упражнений, часть 2

Подсчет строк и слов в файле

from Scientific.IO.TextFile import TextFile
from string import split

lines = 0
words = 0
for line in TextFile('some_file'):
    lines = lines + 1
    words = words + len(split(line))

print lines, " lines, ", words, " words."
Эту программу можно проверить, запустив для сравнения утилиту Unix wc.

Подсчет атомов углерода в PDB-файле

from Scientific.IO.TextFile import TextFile
from Scientific.IO.FortranFormat import FortranFormat, FortranLine
from string import strip

generic_format = FortranFormat('A6')

atom_format = FortranFormat('A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,' +
                            '3X,3F8.3,2F6.2,7X,A4,2A2')

carbons = 0
for line in TextFile('protein.pdb'):
    record_type = FortranLine(line, generic_format)[0]
    if record_type == 'ATOM  ' or record_type == 'HETATM':
        data = FortranLine(line, atom_format)
        atom_name = strip(data[2])
	if atom_name[0] == 'C':
	    carbons = carbons + 1

print carbons, " carbon atoms"
Обратите внимание, как используется strip для удаления всех пробелов рядом с названием атома, присутствующих в правильно отформатированном PDB-файле (но не в PDB-файлах, генерируемых большинством программ).

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


Преобразование из формата PDB в формат XYZ

from Scientific.IO.TextFile import TextFile
from Scientific.IO.FortranFormat import FortranFormat, FortranLine
from string import join, strip

generic_format = FortranFormat('A6')

atom_format = FortranFormat('A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,' +
                            '3X,3F8.3,2F6.2,7X,A4,2A2')

atom_lines = []

for line in TextFile('protein.pdb'):
    record_type = FortranLine(line, generic_format)[0]
    if record_type == 'ATOM  ' or record_type == 'HETATM':
        data = FortranLine(line, atom_format)
        atom_name = strip(data[2])
	element = atom_name[0]
	atom_lines.append(join([element, str(data[8]),
				str(data[9]), str(data[10])]) + '\n')

output_file = TextFile('protein.xyz', 'w')
output_file.write(`len(atom_lines)` + '\n\n')
output_file.write(join(atom_lines, ''))
output_file.close()
Есть несколько способов избежать повторного считывания файла для нахождения количества атомов, но общая идея такова: собрать информацию в список и на следующем шаге ее записать. В примере, приведенном выше, список содержит текстовые строки выходных файлов, но его также можно использовать для хранения вводимых строк или какого-либо промежуточного представления даных.

Организация выходных строк в список имеет то преимущество, что отпадает необходимость во втором цикле для генерации выходного файла; достаточно одной операции записи и вызова join. Фактически эта комбинация встречается так часто, что имеет смысл ввести соответствующее сокращение, что даже более эффективно: вместо output_file.write(join(atom_lines, '')) записать output_file.writelines(atom_lines).


Решение упражнений, часть 3

Нахождение общей массы пептидной цепочки

Сначала модуль, содержащийся в файле под названием PeptideChain.py:
import string

# Таблица перевода из однобуквенного кода в трехбуквенный.
residue_name = {'a': 'ala', 'r': 'arg', 'n': 'asn', 'd': 'asp', 'c': 'cys',
		'q': 'gln', 'e': 'glu', 'g': 'gly', 'h': 'his', 'i': 'ile',
		'l': 'leu', 'k': 'lys', 'm': 'met', 'f': 'phe', 'p': 'pro',
		's': 'ser', 't': 'thr', 'w': 'trp', 'y': 'tyr', 'v': 'val'}

# Таблица масс остатков.
residue_mass = {'ala': 71.0790184654,
		'arg': 157.196106437,
		'asn': 114.104059219,
		'asp': 114.080688839,
		'cys': 103.143406585,
		'gln': 128.131048075,
		'glu': 128.107677695,
		'gly': 57.0520296093,
		'his': 137.141527428,
		'ile': 113.159985034,
		'leu': 113.159985034,
		'lys': 129.182660137,
		'met': 131.197384297,
		'phe': 147.17714379,
		'pro': 97.1170442246,
		'ser': 87.0783231891,
		'thr': 101.105312045,
		'trp': 186.213916723,
		'tyr': 163.176448514,
		'val': 99.1329961777}

# Дополнительная ("концевая") масса пептидной цепочки.
terminus_mass = 18.0152566767

# Функция, выполняющая расчет
def totalMass(peptide_chain):
    mass = terminus_mass
    for residue in peptide_chain:
	residue = string.lower(residue)
	if len(residue) == 1:
	    residue = residue_name[residue]
	mass = mass + residue_mass[residue]
    return mass
Заметьте, что для хранения информации об остатках применяются два словаря.

Теперь простое приложение, использующее этот модуль:

from PeptideChain import totalMass

print totalMass('AEG')
print totalMass(['ala', 'arg', 'gly', 'his'])

Чтение и запись матриц

Модуль MatrixIO.py:

from Scientific.IO.TextFile import TextFile
import string
import Numeric

def readMatrix(filename):
    rows = []
    for line in TextFile(filename):
        columns = []
        for number in string.split(line):
            columns.append(string.atof(number))
	rows.append(columns)
    return Numeric.array(rows)

def writeMatrix(a, filename):
    file = TextFile(filename, 'w')
    for line in a:
	for number in line:
	    file.write(`number` + ' ')
	file.write('\n')
    file.close()

Функция string.atof в первой функции применяется ко всем элементам списка. Эта операция встречается настолько часто, что имеет смысл ввести специальное сокращение: функция map(function,sequence) использует произвольную функцию по отношению к каждому элементу последовательности и возвращает список результатов. Благодаря этому можно заменить строки

        columns = []
        for number in string.split(line):
            columns.append(string.atof(number))
одной единственной строкой columns = map(string.atof, string.split(line)).

Обработка данных

Сначала генерация тестовых данных:

import MatrixIO
import Numeric

for i in range(1, 51):
    time = Numeric.arrayrange(10.)
    values = 4.*Numeric.ones(10)
    parameters = (i, -3.5)
    data = Numeric.zeros((11, 2), Numeric.Float)
    data[0] = parameters
    data[1:, 0] = time
    data[1:, 1] = values
    MatrixIO.writeMatrix(data, "data"+`i`)
Заметьте, как перед выводом все было помещено в один большой массив; этот трюк дает возможность использовать матричную функцию вывода, написанную в последнем упражнении. Результатом программы анализа должен быть постоянный временной ряд, значение которого равно двойной сумме целых чисел от 1 до 50, т.е. 50*51=2550.

А сейчас программа, которая будет "анализировать" файлы данных.

import MatrixIO
import Numeric
import Gnuplot

total = 0.

for i in range(1, 51):
    data = MatrixIO.readMatrix("data"+`i`)
    time = data[1:, 0]
    values = data[1:, 1]
    parameters = data[0]
    total = total + parameters[0]*Numeric.sqrt(values)

plot_data = Numeric.transpose(Numeric.array([time, total]))
Gnuplot.plot(plot_data)

Решение упражнений, часть 4

Матричные функции

import Numeric; N = Numeric
import LinearAlgebra; LA = LinearAlgebra

# Вычислить произвольную матричную функцию 
def matrixFunction(function, matrix):
    eigenvalues, eigenvectors = LA.eigenvectors(matrix)
    eigenvalues = function(eigenvalues)
    return N.dot(eigenvalues*N.transpose(eigenvectors), eigenvectors)

# Вычислить экспоненту по разложению в ряд Тейлора
def matrixExponential(matrix):
    sum = N.identity(matrix.shape[0])
    product = sum
    for i in range(1, 50):
	product = N.dot(product, matrix)/i
	sum = sum + product
    return sum

# Создать симметричную матрицу
m = N.reshape(N.arrayrange(9.), (3, 3))**2/10.
m = (m + N.transpose(m))/2.

# Сравнить результаты
print matrixExponential(m)
print matrixFunction(N.exp, m)

Визуализация отличий в структуре

from Scientific.IO.TextFile import TextFile
from Scientific.IO.FortranFormat import FortranFormat, FortranLine
from Scientific.Geometry import Vector
import string
from Scientific.Visualization import VRML

generic_format = FortranFormat('A6')
atom_format = FortranFormat('A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,' +
                            '3X,3F8.3,2F6.2,7X,A4,2A2')

# Считать PDB-файл и составить список всех положений C-альфа
def readCAlphaPositions(filename):
    c_alpha_positions = []
    for line in TextFile(filename):
	record_type = FortranLine(line, generic_format)[0]
	if record_type == 'ATOM  ' or record_type == 'HETATM':
	    data = FortranLine(line, atom_format)
	    atom_name = string.strip(data[2])
	    position = Vector(data[8:11])
	    if atom_name == 'CA':
		c_alpha_positions.append(position)
    return c_alpha_positions

# Считать обе структуры
conf1 = readCAlphaPositions('conf1.pdb')
conf2 = readCAlphaPositions('conf2.pdb')

# Создать сцену VRML и материал для стрелок
scene = VRML.Scene([])
arrow_material = VRML.DiffuseMaterial('black')
arrow_radius = 0.2

# Создать сферы и расположить их на сцене
for i in range(len(conf1)):
    scene.addObject(VRML.Arrow(conf1[i], conf2[i], arrow_radius,
			       material=arrow_material))

# Просмотр сцены
scene.view()

Решение упражнений, часть 5

Матричные классы

import Numeric;  N = Numeric
import LinearAlgebra; LA = LinearAlgebra

class Matrix:

    def __init__(self, array):
	if type(array) != type(N.array([0])) or len(array.shape) != 2:
	    raise ValueError, "argument not a rank-2 array"
	self.array = array

    def __str__(self):
	return str(self.array)

    def __repr__(self):
	return 'Matrix:\n' + str(self.array)

    def __add__(self, other_matrix):
	return Matrix(self.array+other_matrix.array)

    def __sub__(self, other_matrix):
	return Matrix(self.array-other_matrix.array)

    def __mul__(self, other_matrix):
	return Matrix(N.dot(self.array, other_matrix.array))

    def generalized_inverse(self):
	return Matrix(LA.generalized_inverse(self.array))


class SquareMatrix(Matrix):

    def __init__(self, array):
	Matrix.__init__(self, array)
	if self.array.shape[0] != self.array.shape[1]:
	    raise ValueError, "argument not a square matrix"

    def inverse(self):
	return SquareMatrix(LA.inverse(self.array))

    def eigenvalues(self):
	return LA.eigenvalues(self.array)

Расчет массы пептида

import string

class Residue:

    residue_name = {'a': 'ala', 'r': 'arg', 'n': 'asn', 'd': 'asp', 'c': 'cys',
		    'q': 'gln', 'e': 'glu', 'g': 'gly', 'h': 'his', 'i': 'ile',
		    'l': 'leu', 'k': 'lys', 'm': 'met', 'f': 'phe', 'p': 'pro',
		    's': 'ser', 't': 'thr', 'w': 'trp', 'y': 'tyr', 'v': 'val'}

    residue_mass = {'ala': 71.0790184654, 'arg': 157.196106437,
		    'asn': 114.104059219, 'asp': 114.080688839,
		    'cys': 103.143406585, 'gln': 128.131048075,
		    'glu': 128.107677695, 'gly': 57.0520296093,
		    'his': 137.141527428, 'ile': 113.159985034,
		    'leu': 113.159985034, 'lys': 129.182660137,
		    'met': 131.197384297, 'phe': 147.17714379,
		    'pro': 97.1170442246, 'ser': 87.0783231891,
		    'thr': 101.105312045, 'trp': 186.213916723,
		    'tyr': 163.176448514, 'val': 99.1329961777}

    def __init__(self, name):
	name = string.lower(name)
	if len(name) == 1:
	    name = self.residue_name[name]
	self.name = name

    def mass(self):
	return self.residue_mass[self.name]


class PeptideChain:

    terminus_mass = 18.0152566767

    def __init__(self, sequence):
	self.sequence = map(Residue, sequence)

    def mass(self):
	total_mass = self.terminus_mass
	for residue in self.sequence:
	    total_mass = total_mass + residue.mass()
	return total_mass


# Пример:
print PeptideChain('AEG').mass()
print PeptideChain(['ala', 'arg', 'gly', 'his']).mass()