Рассказы о математике с примерами на языках Python и C - читать бесплатно онлайн полную версию книги автора Дмитрий Дмитриевич Елисеев (ч. 2)

20. Нейронные сети - распознавание текста

Рассмотрим практический пример - распознавание выше написанной нейронной сетью рукописного текста. Для тренировки сети воспользуемся имеющейся в открытом доступе базой MNIST, содержащей 60000 черно-белых рукописных изображений цифр размером 28x28пкс.

Выглядят изображения примерно так:

Массив изображений каждой цифры хранится в виде чисел float в диапазоне от 0 до 1, длина массива равна 28х28 = 784 значения. Кстати, при выводе массива в консоль форма цифры вполне узнаваема:

(картинка соответствует числу “5”)

Всего таких чисел, как было сказано, 60000, что вполне хватит для тренировки нашей сети.

Воспользуемся уже готовым классом MLP, изменим лишь некоторые параметры:

- Число входов будет равно 28х28 = 784, т.к. на вход нейросети будет подаваться массив целиком.

- На выходе нейросеть должна давать значения от 0 до 9, однако мы не можем получить значение > 1, так что воспользуемся другим путем: нейросеть будет иметь 10 выходов, где “1” будет соответствовать указанной цифре. Т.е. для числа “0” мы получим [1,0,0,0,0,0,0,0,0,0], для “1” [0,1,0,0,0,0,0,0,0,0], и так далее.

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

В итоге, мы имеем такую структуру сети:

Рассмотрим подробнее программу на языке Python.

Тренировка сети

Создаем нейронную сеть и подготавливаем 2 массива тестовых данных (исходные и целевых значения). Массив “digits” загружается из зараннее подготовленного текстового файла базы MNIST. Для каждой цифры содержится исходный массив в 784 цифры, и правильное значение числа, которое и будет использоваться для тренировки.

mlp = MLP(n_in=28 * 28, n_hidden=48, n_out=10)

print "1. Загрузка данных"

with open("digitsMnist.txt", 'rb') as fp:

digits = cPickle.load(fp)

print "2. Подготовка данных"

batch = 50

inputs = []

targets = []

for p in range(batch):

data = digits[p]["data"]

resulVal = digits[p]["result"]

# Конвертация целевого числа в массив: “3” => [0,0,0,1,0,0,0,0,0,0]

result_flat = [0.0] * 10

result_flat[resulVal] = 1.0

inputs.append(data)

targets.append(result_flat)

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

print "3. Обучение"

mlp.load()

try:

N = 10000

for p in range(N):

errSum = 0.0

for d in range(batch):

mlp.forwardPass(inputs[d])

mlp.backPass(inputs[d], targets[d])

errSum += mlp.calcError(inputs[d], targets[d])

print "Step", p, "error", errSum

except KeyboardInterrupt:

pass

mlp.save()

Всего было проведено около 50000 итераций на блоке размером 50 цифр, весь массив в 60000 цифр не использовался. Обучение сети заняло порядка 2х часов на компьютере с Core i7. Разумеется, этот проект учебный, и программа никак не оптимизирована по быстродействию. С использованием специальных фреймворков типа Keras тренировка такой сети занимает не более 5 минут.

Тестирование сети

И наконец, тестирование. В качестве теста, “от руки” была нарисована цифра “3”, которая получилась в виде массива примерно такой:

t = [[ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 ] +

[ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 ]]

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

for digit in t:

printDigit(digit)

mlp.forwardPass(digit)

mlp.printOutput()

Результат весьма интересен:

[0.019443, 0.0102247, 0.0091067, 0.467039, 0.0000103, 0.0150096, 0.0008709, 0.0004573, 0.655264, 0.0867915]

На тренировочных значениях сеть должна была выдавать “1” в том разряде, который соответствует “загаданному” числу. И хотя нигде явным образом это не планировалось и специально не программировалось, сеть выдала максимальные значения в тех разрядах массива, где “по ее мнению”, степень схожести максимальна. Как можно видеть, эти разряды соответствуют числам “3” и “8”, что визуально действительно похоже на нарисованное число. Строго говоря, результат не совсем верный, картинка больше похожа на “3” чем на ”8”, но и тренировка не была закончена полностью.

Для следующего теста число нейронов было увеличено до 128, а количество итераций обучения было увеличено до 75000.

Результаты теста на том же массиве:

[0.0044085, 0.219943, 0.0118179, 0.816858, 0.0002372, 0.06267359, 0.0008674, 0.00020319, 0.55439027, 0.00061541]

Максимальное значение теперь действительно соответствует числу “3”, можно сказать что сеть работает как задумывалось.

Для тех кто захочет поэкспериментировать самостоятельно, исходный текст программы и база чисел MNIST находятся в архиве, ссылка на который есть в начале книги.

21. Использование библиотеки TensorFlow

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

Для установки TensorFlow для Windows необходимо:

- Установить 64-битную версию Python (дистрибутив можно скачать по ссылке https://www.python.org/downloads/release/python-365/)

- Запустить из консоли команду pip3 install tensorflow

Рассмотрим использование Tensorflow на простом примере: умножение двух массивов. Исходный код приведен ниже.

import tensorflow as tf

# Initialize two constants

x1 = tf.constant([1,2,3,4])

x2 = tf.constant([5,6,7,8])

# Multiply

task = tf.multiply(x1, x2)

# Intialize the Session

session = tf.Session()

# Print the result

print(session.run(task))

# Close the session

session.close()

Сохраним программу в файле tf1.py и запустим, набрав в консоли python3 tf1.py (обратим внимание на цифру “3” в версии Python).

Разберем код программы подробнее. Функция tf.constant позволяет создать константу, в нашем случае это массив со значениями x1 = [1,2,3,4]. Аналогично создается массив x2. Наконец, команда task = tf.multiply(x1,x2) умножает два массива. Казалось бы все пока просто - но не совсем. Если мы напишем print(task), то никакого результата не увидим, на экран отобразится лишь что-то типа Tensor(“Mul:0”, dtype=int32). На самом деле, в памяти TensorFlow формируется так называемый “граф вычислений” (computation graph). Как уже упоминалось, Python - это интерпретатор, и работает довольно-таки медленно, поэтому библиотека TensorFlow формирует в памяти собственную модель данных для вычислений. Для того, чтобы получить реальный результат, мы должны создать объект Session, и запустить в нем вычисления командой session.run(task). Когда вычисления произведены, сессию можно закрыть командой session.close().

Разумеется, можно вычислять не только константы. Рассмотрим пример использования переменных.

import tensorflow as tf

k = tf.constant(2)

a = tf.placeholder(tf.int32)

b = tf.placeholder(tf.int32)

m = tf.multiply(k, a)

task = tf.add(m, b)

session = tf.Session()

print(session.run(task, feed_dict={a: 2, b: 3}))

print(session.run(task, feed_dict={a: [10,11], b: [3,4]}))

session.close()

Здесь мы задаем граф вычислений по формуле a*k + b, где a и b - это переменные, для указания этого используется специальный тип placeholder. Теперь, при вычислении командой session.run мы можем передать различные параметры, например a=2 и b=3:

print(session.run(task, feed_dict={a: 2, b: 3}))

>> 7

Можно даже передать массивы в качестве параметра:

print(session.run(task, feed_dict={a: [10,11], b: [3,4]}))

> [23, 26]

Библиотека TensorFlow поддерживает множество различных операций, например, вычисления с матрицами, ну и разумеется, большое количество функций предназначено специально для обработки данных и машинного обучения. Вернемся к предыдущей главе, и перепишем нашу нейронную сеть, распознающую текст. За основу исходного кода был взят готовый пример из https://github.com/aymericdamien/TensorFlow-Examples/tree/master/examples желающие могут изучить приведенные там примеры более подробно.

Код целиком приведен ниже:

import tensorflow as tf

import _pickle as cPickle

# Prepare training data

training_epochs = 500

learning_rate = 0.001

inputs = []

targets = []

with open("digitsMnist.txt", 'rb') as fp:

digits = cPickle.load(fp)

print("Digits loaded:", len(digits))

inputs = list(map(lambda x: x["data"], digits))

def arrayFromDigit(data):

result = data["result"]

# Convert value to array: "3" => [0,0,0,1,0,0,0,0,0,0]

result_flat = [0.0]*10

result_flat[result] = 1.0

return result_flat

targets = list(map(arrayFromDigit, digits))

# Network Parameters

n_hidden_1 = 48 # 1st layer number of neurons

n_input = 784 # MNIST data input (img shape: 28*28)

n_classes = 10 # MNIST total classes (0-9 digits)

model_path = "model_mlp.ckpt"

# tf Graph input

X = tf.placeholder("float", [None, n_input])

Y = tf.placeholder("float", [None, n_classes])

# Store layers weight & bias

weights = {

'h1': tf.Variable(tf.random_normal([n_input, n_hidden_1])),

'out': tf.Variable(tf.random_normal([n_hidden_1, n_classes]))

}

biases = {

'b1': tf.Variable(tf.random_normal([n_hidden_1])),

'out': tf.Variable(tf.random_normal([n_classes]))

}

# Create model

def multilayer_perceptron(x):

# Hidden fully connected layer with 256 neurons

layer_1 = tf.add(tf.matmul(x, weights['h1']), biases['b1'])

# Output fully connected layer with a neuron for each class

out_layer = tf.matmul(layer_1, weights['out']) + biases['out']

return out_layer

# Construct model

logits = multilayer_perceptron(X)

# Define loss and optimizer

loss_op = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits_v2(logits=logits, labels=Y))

optimizer = tf.train.AdamOptimizer(learning_rate=learning_rate)

train_op = optimizer.minimize(loss_op)

# Initializing the variables

init = tf.global_variables_initializer()

# Save/restore model

saver = tf.train.Saver()

session = tf.Session()

session.run(init)

# Training cycle

for p in range(training_epochs):

_, c = session.run([train_op,loss_op], feed_dict={X:inputs, Y:targets})

if p % 50 == 0:

print("Step {} of {}".format(p, training_epochs))

print("Training Finished!")

# Save model (optional)

# save_path = saver.save(sess, './' + model_path)

# print("Model saved in file: %s" % save_path)

# Restore model weights from previously saved model (optional)

# saver.restore(sess, model_path)

# print("Model restored from file: %s" % save_path)

# Use model

pred = tf.nn.softmax(logits) # Apply softmax to logits

print("")

t = [ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,

0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,

0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,

0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,

0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 ]

res = session.run(pred, feed_dict={X: [t]})

print("Result:", res[0])

session.close()

Как можно видеть, суть практически не отличается от описанного в главе 19. Мы создаем модель в функции multilayer_perceptron, состоящую из входных данных x, скрытого слоя layer_1 и выходного слоя out_layer. Мы также имеем массивы weights и biases. Далее, весь громоздкий код тренировки сети, который был написан в главе 19, здесь не нужен - он уже реализован в TensorFlow, нам достаточно вызвать функцию optimizer.minimize в цикле тренировки.

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

Результаты работы программы показаны на рисунке:

Обучение сети заняло … 20 секунд вместо 2х часов - что показывает кардинальный прирост в быстродействии по сравнению с неоптимизированным кодом на Python, выходные данные соответствуют “нарисованной” цифре “3”.

22. Фракталы

Фракталы - интересные математические объекты, которые долгое время не замечались математиками. В 1872 году, объекты такого типа впервые описал немецкий математик Карл Вейерштрасс, “придумавший” функцию следующего вида:

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

В Википедии можно найти ее график:

Во времена Вейерштрасса математики назвали такую функцию “оскорблением здравого смысла”, практической применимости не мог найти никто. Лишь через 100 лет, в историческом плане совсем недавно, в 1975 году, математик Мандельброт ввел термин “фрактал”. Фракталы стали широко известны после выхода в 1977 году его книги “Фрактальная геометрия природы” - как оказалось, такие структуры распространены в природе гораздо чаще, чем это можно представить. От ветвей деревьев или рисунка береговой линии, до сосудов человека или колебаний курса акций - все эти процессы имеют природу, схожую с фрактальной.

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

Другим примером служит так называемый “парадокс береговой линии”. При измерении длины государственной границы между Испанией и Португалией, выяснилось, что по измерениям из Португалии, длина границы составила 987км, а по измерениям с Испанской стороны - 1214км. Т.к. противоречие было налицо, им заинтересовался английский математик Льюис Ричардсон, который выяснил, что все дело в масштабе - чем меньше масштаб, тем больше получается длина. Поразительным оказалось то, что при теоретически бесконечном масштабе, длина линии также будет бесконечной!

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

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

Если замкнуть 3 такие линии, получится так называемая “снежинка Коха”. Нарисовать ее несложно на языке Python, благодаря встроенной библиотеке рисования Turtle. Она позволяет управляя виртуальной “черепашкой”, рисовать примитивы простыми командами, типа “повернуть на угол” или “сдвинуться на расстояние”.

Для рисования линии Коха, нам потребуются всего 3 команды: сдвинуться вперед, повернуть на 60 градусов, и повернуть на -60 градусов. Повторив процедуру 3 раза, мы получим “снежинку Коха”.

Код на языке Python выглядит так:

import turtle

def koch(size, order):

if order > 0:

koch(size/3, order - 1)

turtle.left(60)

koch(size/3, order - 1)

turtle.right(120)

koch(size/3, order - 1)

turtle.left(60)

koch(size/3, order - 1)

else:

turtle.forward(size)

turtle.ht()

turtle.speed(20)

for i in range(3):

koch(size=400, order=4)

turtle.right(120)

turtle.exitonclick()

В результате работы программы мы получаем следующую картинку:

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

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

S = 0.9 + 0.09 + 0.009 + …

Как нетрудно догадаться, сам по себе ряд бесконечный, но его сумма стремится ко вполне определенному числу 0.999999… = 1.

“Снежинка Коха” интересна своей простотой, но есть и более сложные варианты фракталов. Рассмотрим множество Мандельброта, которое строится по следующему алгоритму: для каждой точки c=x + j*y на комплексной плоскости, эта точка попадает в множество при условии что ряд Z_n+1 = Z_n² + c не расходится при z₀ = 0.

Рассмотрим пример.

При c = (0,0) = 0 + 0j, ряд очевидно, всегда равен нулю, т.е. не расходится, значит точка (0,0) попадает в множество.

При с = (0.5, 0.5) = 0.5 + 0.5j, получаем следующие значения:

Z₀ = 0

Z₁ = Z₀ + c = 0.5+0.5j

Z₂ = Z₁ + c = 0.5+1j

Z₃ = Z₂ + c = -0.25+1.5j

Z₄ = Z₃ + c = -1.6875-0.25j

Z₅ = Z₄ + c = 3.28515625+1.34375j

Z₆ = Z₅ + c = 9.48658752+9.32885j

Легко видеть, что ряд расходящийся, числа увеличиваются, значит (0.5, 0.5) не попадает в множество.

При с = (0.25, 0.25) = 0.25 + 0.25j, имеем следующий ряд:

Z₀ = 0

Z₁ = Z₀ + c = 0.25+0.25j

Z₂ = Z₁ + c = 0.25+0.375j

Z₃ = Z₂ + c = 0.171875+0.4375j

Z₄ = Z₃ + c = 0.088134+0.400390625j

Z₅ = Z₄ + c = 0.09745508+0.32057666j

Z₆ = Z₅ + c = 0.15672809+0.31248365j

Z₇ = Z₆ + c = 0.17691766+0.34794993j

Z₈ = Z₇ + c = 0.160230+0.37311697j

Z₉ = Z₈ + c = 0.1364575+0.36956959j

Как можно видеть, ряд не расходится, т.е. точка (0.25, 0.25) попадает в множество.

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

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

Сама функция имеет следующий вид:

def countMandelbrotIterationsForPt(pt):

z = 0

c = pt

threshold = 64

for iteration in xrange(threshold):

z = z*z + c

if abs(z) > 4:

break

return iteration

Здесь вычисляется 64 итерации, которых вполне достаточно чтобы определить, является ли ряд расходящимся или нет.

Чтобы сохранить весь фрактал, нужно перебрать все точки в диапазоне [-2, 2], программа, сохраняющая фрактал в файл, приведена ниже:

from PIL import Image

points = 1000

img = Image.new('RGB', (2*points,2*points), "black")

pixels = img.load()

# Range.X: -2..2,

# Range.Y: -2..2

for ix in xrange(-points, points, 1):

for iy in xrange(-points, points, 1):

pt = complex(2.0*ix/points, 2.0*iy/points)

i = countMandelbrotIterationsForPt(pt)

if i > 16:

colR, colG, colB = 4*i, 4*i, 0

if colR >= 255: colR = 255

if colB >= 255: colB = 255

img_x = points + ix

img_y = points + iy

pixels[img_x, img_y] = (colR, colG, colB)

if ix % 10 is 0:

print "Done: {}%".format(100.0*(ix + points)/(2*points))

img.save("fractal.png")

Как можно видеть, программа содержит 2 вложенных цикла ix, iy, затем значение преобразуются в комплексное число, для которого и выполняется описанная выше функция проверки.

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

Если в несколько раз увеличить количество точек, можно увидеть структуру верхней части более детально:

Немного изменив формулу, можно получить другие виды фракталов. Например, множество Жюлиа, описывается также, но Z₀ = pt, и c = const. Казалось бы, небольшое отличие, приводит к совершенно другому рисунку фрактала:

Кстати, этот рисунок напоминает сорт цветной капусты “Романеско”:

Наконец, если в формуле фрактала Мандельброта Z_n+1 = Z_n² + c, вычислять абсолютное значение Z*Z в виде (|a| + j|b|)*(|a| + j|b|), мы получим фрактал с названием “горящий корабль”:

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

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

В youtube можно найти интересные трехмерные визуализации фракталов, например, набрав в поиске “mandelbrot 3d”:

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

23. Горн Гавриила

В математике есть некоторые парадоксы, связанные с бесконечностью.

Так называемый “Горн Гавриила” (Gabriel's Horn) - это фигура, образованная функцией y = 1/x, если создать из нее 3х-мерную фигуру:

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

Действительно, формулы площади (V) и объема (A) “горна” легко найти в Википедии:

Из них видно, что при a = ∞, V становится равным ℼ, а A действительно стремится к бесконечности.

24. Построение графиков функций

Очень многие закономерности гораздо легче исследовать графически. Очень просто и эффективно строить графики с помощью языка Python и библиотеки matplotlib.

Выведем графики функций y = x² и у = 0.5x².

import matplotlib.pyplot as plt

plt.title('Graph')

plt.xlabel('X')

plt.ylabel('Y')

x_values = range(0,10)

y_values1 = map(lambda x: 0.5*x*x, x_values)

y_values2 = map(lambda x: x*x, x_values)

plt.axis([0, 10, 0, 20]) # xmin, xmax, ymin, ymax

plt.plot(x_values, y_values1, 'ro', label='0.5*x*x')

plt.plot(x_values, y_values2, marker='o', linestyle='--', label='x*x')

plt.legend()

plt.show()

Результат работы программы выглядит так:

Как можно видеть, мы имеем массив входных данных x_values = range(0,20), и 2 массива выходных данных y_values1 и y_values2, вычисляемых по соответствующим формулам. Все остальное за нас делает библиотека matplotlib. Команда plt.plot выводит график на экран, причем можно задать разные варианты оформления. Команда выводит подсказку (“легенду”), и является опциональной. Команды plt.title, plt.xlabel и plt.ylabel задают названия графика и осей. Команда plt.show выводит окно с графиком на экран.

Рассмотрим более сложный пример.

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

Программа вывода графика выглядит так:

import matplotlib.pyplot as plt

import math

def is_prime(n):

if n % 2 == 0 and n > 2:

return False

for i in xrange(3, int(math.sqrt(n)) + 1, 2):

if n % i == 0:

return False

return True

def get_primes(n):

cnt = 0

for i in xrange(1, n):

if is_prime(i):

cnt += 1

return cnt

N = 400

x_values = range(2, N)

y_values = map(lambda n: get_primes(n), x_values)

plt.axis([0, N, 0, N/4])

plt.plot(x_values, y_values, 'ro', markersize=1)

plt.show()

Кстати, согласно Википедии, еще в 18м веке Гауссом и Лежандром было высказано предположение что функция распределения простых чисел выглядит так:

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

plt.plot(x_values,map(lambda n:n/math.log(n),x_values), marker='o', linestyle='--')

Окончательный график выглядит так:

С помощью библиотеки numpy строить графики еще проще.

График функции sin(t) в диапазоне 0..2𝜋:

import matplotlib.pyplot as plt

import numpy as np

t = np.arange(0.0, 2.0, 0.01)

s = np.sin(2*np.pi*t)

plt.plot(t, s)

plt.show()

Важно заметить, что t здесь - не одно число, а целый массив - в данном примере используется возможность numpy работать непосредственно с массивами.

И наконец, с помощью matplotlib можно строить даже 3х-мерные графики функции двух переменных, для этого используется модуль mplot3d:

from mpl_toolkits.mplot3d import axes3d

import matplotlib.pyplot as plt

import math

fig = plt.figure()

ax = fig.gca(projection='3d')

# Axis

ax.set_xlabel('X')

ax.set_xlim(0, 20)

ax.set_ylabel('Y')

ax.set_ylim(0, 20)

ax.set_zlabel('Z')

ax.set_zlim(0, 20)

# Graph

for x in range(0,20):

for y in range(0,20):

z = 10*math.sin(0.1*x + 0.1*y)

ax.scatter(x, y, z, c='r', s=1)

plt.show()

Результат показан на рисунке:

Интересно отметить кросплатформенность кода на Python - приведенная выше программа работает и на Windows и на OSX без каких-либо изменений.

25. Точность компьютерных вычислений

Запустим интерпретатор Python и введем простую программу:

a1 = 0.5 - 0.4

a2 = 1.5 - 1.4

print a1

print a2

Получим результаты 0.1 и еще раз 0.1. Пока все логично.

Введем новую команду:

print a1 == a2

В результате получаем … False. 0.1 != 0.1? Как такое может быть?

Чтобы убедиться, что значения не равны, введем:

print a1 - a2

Получаем вовсе не 0, как хотелось бы, а -1.1102230^-16. Почему так? Разберемся.

Полученный результат - вовсе не ошибка, а вполне известная особенность компьютерных вычислений с вещественными числами. Любое число в языке программирования имеет свой определенный тип. Обычно это или целое число (1,2,3…), или так называемое число с плавающей точкой (1.1, 2.5, 3.8,...).

С целыми числами все просто - в памяти компьютера они хранятся в двоичном виде как суммы степеней 2. Например для числа “5”:

5 = 0*128 + 0*64 + 0*32 + 0*16 + 0*8 + 1*4 + 0*2 + 1*1 = 00000101b.

Любое десятичное число можно перевести в двоичное “туда” и “обратно”, соответствие тут однозначно, никакой потери точности нет. Ограничения лишь в количестве используемых байт, например для 1 байта (8 бит) максимальное число будет 128+64+32+16+8+4+2+1 = 255. Для 2х бит максимально хранимое значение равно 65535, для 4х бит 2147483647, и так далее.

C вещественными числами все сложнее. Для их хранения используется так называемый стандарт IEEE-754. Согласно этому стандарту, значения хранятся в так называемом формате “чисел с плавающей точкой”, состоящем из трех величин - знака (+ или -), мантиссы и экспоненты.

value = (1 + b₀/2 + b₁/4 + b₂/8 + b₃/16 + … ) * 2^e-127

(b - биты мантиссы, e - экспонента)

К примеру, значение 3.14 в двоичном виде будет храниться так: 01000000010010001111010111000011. Первый 0 - это знак (+), 1000000 = 128 - это экспонента, а 10010001111010111000011 - это мантисса.

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

(1 + 1/2 + 1/16 + ...) * 2^128-127 = 3.14

Ключевая особенность таких чисел - это их ограниченная точность, ведь на мантиссу и экспоненту отводится вполне определенное число бит. Ошибка весьма невелика, но тем не менее она есть. Например, число 0.5 будет храниться как 00111111000000000000000000000000, что действительно дает 0.5. А вот 0.6 будет храниться как 00111111000110011001100110011010, что в результате дает 0.60000002384185791015625.

Разумеется, записывая в языке программирования строку a=0.6, пользователь не задумывается о битах. Но как в поговорке про суслика, мы его не видим, а он есть.

Введем команды в языке Python:

a = 0.8

print a

Получим 0.8 - команда print достаточно “умна” чтобы округлить результат. Но сделаем вывод с большим числом знаком после запятой:

print("{0:.20f}".format(a))

Получим 0.8000000000000004441, что как говорится, и требовалось доказать.

Что с этим делать? Ну в общем-то ничего. Реальная ошибка в точности достаточно мала, чтобы это было проблемой. Но при желании можно использовать типы данных большей точности, правда ценой увеличения объема памяти и замедления скорости вычислений. Например, в С++ можно выбрать между типами float и double, второй имеет большую точность, но занимает больше памяти.

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

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

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

a = 0.5

b = 0.4

if a-b == 0.1:

print "0.1!"

Вместо этого достаточно убедиться, что результат мал:

if abs(a-b-0.1) < 1e-8:

print "0.1!"

3) Для денежных расчетов, где ошибки округления весьма нежелательны и баланс должен сойтись, существуют специальные библиотеки и типы данных, например библиотека python-money для Python.

4) Если все-таки желательно повысить точность, в Python можно использовать тип Decimal. Он хранит, по возможности, все значимые биты числа.

from decimal import Decimal

a = Decimal("0.5")

b = Decimal("0.4")

print a-b == Decimal("0.1")

Данная программа действительно вернет True. Но достигается это за счет значительного (до 10 раз или даже выше) замедления скорости расчета, да и не для всех вычислений это работает, например для уже упомянутого числа 1/3 точное значение сохранить невозможно. Но для таких случаев можно хранить числа в виде дробей, сохраняя и числитель и знаменатель. В общем, единого решения тут нет, каждый случай нужно рассматривать отдельно.

Но опять же, стоит повторить, что в большинстве случаев можно не беспокоиться - точности стандартных типов данных вполне достаточно. К примеру, если сохранить в переменной расстояние от Москвы до Владивостока, равное 9136.9км, то ошибка округления составит 10^-11м, что лишь немногим больше размера одного атома. Скорее всего, такой ошибкой можно пренебречь.

26. Математика и календарь

Еще одна большая тема, где соединяются математика, астрономия и реальная жизнь - это наш обычный календарь. Например, часто можно услышать что високосный год считается “тяжелым”. Может ли такое быть?

Для начала, вспомним что Земля вращается по эллиптической орбите вокруг Солнца (что обеспечивает смену времен года), и в то же время, вокруг своей оси (что обеспечивает чередование дней). Второй период всем известен, и составляет 24 часа, а вот первый период известен уже не всем, и составляет 365.242187 суток. А в году, как известно, 365 дней.

На картинке это выглядит примерно так:

Первую неточность в календаре знали уже в древнем Риме, еще во времена Юлия Цезаря. Каждые 4 года добавлялся один день, чтобы скомпенсировать те самые недостающие 1/4 суток. В какой именно год добавлять - несложно догадаться, что принципиальной разницы тут нет, главное чтобы число лет было кратно 4. Можно даже добавлять по 2 дня каждые 8 лет, принципиально это ни на что не повлияет. Может ли високосный год быть “тяжелее” другого, соседнего? Очевидно что нет, т.к. все равно, какой именно год считать високосным - это лишь общепринятая условность, и на орбиту Земли выбор года никак не влияет.

Но если с орбитой все более-менее понятно, то с календарем не все оказалось так просто. Внимательный читатель наверное заметил, что разница в 0.242187 суток вовсе не кратна ¼. И совершенно справедливо. В юлианском календаре год оказывался длиннее астрономического на 11.2 минуты. Казалось бы немного, но за 128 лет накапливался целый день разницы. Впрочем, в древние времена ритм жизни был не столь быстр как сейчас, и ошибку заметили … только спустя 1.5 тысячи лет, когда “набежало” уже более 10 дней разницы. Это стало заметно не только астрономам - такой сдвиг календаря влиял и на даты посевов, и на празднование дней равноденствия. Тогда, в 1582 году папа Григорий XIII ввел новый календарь, названный теперь “григорианским” - в нем убрано часть високосных годов. Современное правило звучит так: год високосный, если он делится на четыре без остатка, но если он делится на 100 без остатка, это не високосный год, но если он делится без остатка на 400, это високосный год. Таким образом, 1996, 2000, 2004 - високосные годы, но 1900й год високосным не был. Правило выглядит запутанным, но зная, что день разницы накапливается за 128 лет, нетрудно понять почему именно так. Каждый 100й год 1 день “убирается” чтобы скомпенсировать разницу, набегающую за 100 лет, но каждый 400й год делается еще одно исключение, чтобы скомпенсировать оставшиеся 28 лет.

Кстати, в России увы, исчисление европейской точностью не отличалось. Юлианский календарь использовался в России аж до 1918 года, суммарная “набежавшая” разница составила 13 суток, что и составляет разницу между “старым” и “новым” стилем (лишь Турция и Египет приняли новый календарь еще позже, в 1927 году).

Так что теперь мы знаем причину, почему все отмечают “старый новый год” именно 13 января - это тот самый промежуток, который накопился за почти 2000 лет из “всего лишь” 11 минутной разницы между реальным и календарным годами.

Является ли григорианский календарь абсолютно точным? Нет, но дополнительных високосных дней никто уже не вводит, это было бы слишком запутанно. Вместо этого, по необходимости, просто чуть корректируют текущее время, например 1 июня 2015 года, по международному соглашению, день длился на 1 секунду дольше - секунда была добавлена для компенсации разницы календарного и астрономического времени. С 1972 года такие поправки вводились 25 раз.

27. Можно ли выиграть в азартные игры?

Одним из “классических” примеров использования теории вероятности являются азартные игры. Собственно, многие из статистических расчетов изначально и возникли из-за желания понять такие процессы, как например, бросание кубика. Посмотрим, какова вероятность выиграть в ту или иную азартную игру.

Игральные кости

Самый наверное, интуитивно простой и понятный вариант - есть кубик с метками от 1 до 6, и вероятность выпадения того или иного числа равна 1/6.

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

Каждый ход игры состоит из бросания двух кубиков, набранные очки суммируются. Как написано в статье http://www.vokrugsveta.ru/article/215452/ “Правила игры незамысловаты: игрок кидает две кости, и, если сумма очков на них равна 7 или 11, он выигрывает, если 2, 3 или 12 — проигрывает. Когда на кубиках выпадает другая сумма, шутер бросает их до выигрышной или проигрышной комбинаций”.

Посмотрим, сколько можно выиграть таким способом. Для этого необязательно идти в казино, для симуляции игры воспользуемся Python. Напишем функцию для одного броска:

import random

def shoot():

return random.randint(1,6) + random.randint(1,6)

Напишем функцию симуляции одного хода по вышеописанным правилам.

def move():

while True:

val = shoot()

print("Dice roll:", val)

if val == 7 or val == 11:

return True

if val == 2 or val == 3 or val == 12:

return False

Зададим нашему виртуальному игроку начальную сумму в 100 у.е., и запустим процесс игры. Пусть наш игрок сделает 100 ставок, каждый раз по 1у.е.

money_total = 100

win = 0

loose = 0

for p in range(100):

bet = 1

step = move()

if step is True:

money_total += bet

win += 1

else:

money_total -= bet

loose += 1

print("Win", win, "Loose", loose, "Money", money_total)

Запускаем симуляцию 100 игр и удивляемся результату, игрок выиграл, причем с заметным отрывом побед от поражений: Win 63, Loose 37, Money 126.

Увеличиваем число игр до 1000 и запускаем еще раз, игрок опять выиграл: Win 680, Loose 320, Money 460.

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

Интуитивно кажется, что при бросании кубика вероятность выпадения любой грани равновероятна. И это действительно так, но в случае одного кубика. Если кубиков два, то все становится гораздо сложнее. К примеру, число 7 может выпасть как 3+4, 2+5, 1+6, а вот число 12 может выпасть только в виде комбинации 6+6.

Построим в Jupyter notebook график выпадения сумм от 1 до 12 для 100 бросков:

from matplotlib import pyplot as plt

%matplotlib inline

y = [ shoot() for v in range(100) ]

plt.hist(y)

Предположение подтвердилось, и суммы из центра действительно выпадают чаще. Таким образом, числа 7 и 11 действительно выпадают чаще чем 2,3 или 12. И вероятность получить выигрышную комбинацию “7 или 11” действительно выше.

Как такое может быть? Увы, ответ прост - автор процитированной выше статьи просто не разобрался досконально в правилах игры. Текст “правила игры незамысловаты: игрок кидает две кости, и, если сумма очков на них равна 7 или 11, он выигрывает, если 2, 3 или 12 — проигрывает” весьма далек от правды, и правила совсем не так незамысловаты как кажутся.

Реальные правила для ставки на pass line оказались несколько сложнее (есть и другие виды ставок, желающие могут разобраться самостоятельно).

Ход-1: Делается бросок. Если выпадает 7 или 11, игрок выигрывает, если 2, 3 или 12, игрок проигрывает. Если выпадает другое число, оно запоминается под названием point.

Ход-2: Делается бросок. Если выпадает 7, то игрок проиграл. Если выпадает point то игрок выиграл. Если выпадают другие числа, ход повторяется (в это время игроки могут также делать ставки на другие числа).

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

def move():

point = 0

while True:

val = shoot()

if point == 0:

# First move

if val == 7 or val == 11:

return True

if val == 2 or val == 3 or val == 12:

return False

point = val

else:

# 2..N move

if val == 7:

return False

if val == point:

return True

Результат теперь более похож на правду: за 100 игр игрок выиграл 43 раза, проиграл 57 раз, баланс в конце игры составил 86у.е. от изначальных 100. Интересно и то, что число выигрышей оказалось довольно-таки велико, и составляет лишь немногим менее 50%. Это грамотная стратегия с точки зрения казино - она позволяет поддерживать интерес участника к игре (все время проигрывать было бы неинтересно), но в то же время баланс прибыли казино остается положительным, а баланс денег игрока - соответственно, отрицательным.

Посмотрим подробнее, что получается для симуляции в 100 игр.

- Шанс выиграть на первом шаге выпал примерно в 20 случаях.

- Шанс проиграть сразу на первом шаге выпал в 15 случаях

- В остальных 65 случаях игра продолжается, и тут все хитро: выбор происходит из двух чисел, 7 и point, но как было видно из графика выше, вероятность выпадения “проигрышной” цифры 7 максимальна, что в общем-то и требовалось доказать.

Интересно заметить, что шанс выигрыша в 45% - довольно-таки высок. Так можно ли выиграть? В краткосрочном периоде, да, например, в другой симуляции игроку “повезло”, и он за 100 игр увеличил свой виртуальный капитал со 100 до 112уе.

Но уже следующая симуляция показала отрицательный баланс: игрок уменьшил свое состояние со 100 до 88уе, потеряв, кстати, те же самые 12уе, “выигранные” в предыдущий раз.

Если увеличить число итераций до 1000, то можно увидеть как может выглядеть денежный баланс игрока в долгосрочной перспективе:

Понятно, что при шансе выигрыша каждой игры менее 50%, результирующая сумма денег на счету игрока будет постепенно уменьшаться, а сумма прибыли казино постепенно увеличиваться. На графике кстати, видны всплески и падения, и может возникнуть резонный вопрос - можно ли их предсказать? Увы нет, т.к. бросания кубика - это независимые друг от друга события, и предыдущие результаты никак не влияют на следующие. Можно еще раз повторить главную мысль - можно выиграть один или даже несколько раз, но в долгосрочном периоде остаться в плюсе у казино невозможно, правила игры составлены так что баланс будет не в пользу игрока.

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

Американская рулетка

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

Игровое поле рулетки делится на 38 ячеек: 36 зон с цифрами + 2 зоны “zero” и “double zero”. Брошенный на рулетку шарик очевидно, остановится в одной из зон. Игрок может делать разнообразные ставки, видов которых более 10, рассмотрим некоторые из них.

Черное-белое (или чет-нечет)

Игрок выигрывает, если названная им ставка совпала. Очевидно, что вероятность черного или белого была бы 50/50, если бы не два поля zero - при попадании на них ставка проигрывает. Как и в случае с крэпсом, это делает вероятность выигрыша лишь чуть менее 50% - но этого “чуть” достаточно, чтобы быть в минусе.

Напишем функцию симуляции хода с помощью случайных чисел от 1 до 38, последние 2 цифры будем считать за “зеро”.

def move_roulette1():

val = random.randint(1, 38)

if val == 37 or val == 38:

return False

return val % 2 != 0

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

money_total = 100

win = 0

loose = 0

for p in range(100):

bet = 1

step = move_roulette1()

if step is True:

money_total += bet

win += 1

else:

money_total -= bet

loose += 1

print("Win", win, "Loose", loose, "Money", money_total)

Результат: игрок выиграл 46 раз и проиграл 54 раза. На графике видно, что у игрока были и “взлеты” и “падения”, но итоговый баланс все равно негативный.

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

Ставка на конкретный номер

Игрок также может поставить на определенный номер, ставка при выигрыше составляет 35:1. Это кажется большим, но нетрудно догадаться, что шанс выпадения определенного номера рулетки 1:38, т.е. опять же, чуть меньше.

Допишем функцию ставки на конкретный номер:

def move_roulette2(num):

val = random.randint(1,38)

return val == num

Симуляция, будем считать что игрок ставит на число 10:

money_total = 100

win = 0

loose = 0

for p in range(100):

bet = 1

step = move_roulette2(10)

if step is True:

money_total += 35*bet

win += 1

else:

money_total -= bet

loose += 1

print("Win", win, "Loose", loose, "Money", money_total)

В результате, игрок выиграл 2 раза и проиграл 98 раз, итоговый баланс -28у.е.

Ставка на два номера

Можно поставить на два номера - шанс выигрыша выше, но зато ставка меньше и составляет 17:1.

Напишем функцию:

def move_roulette3(num1, num2):

val = random.randint(1,38)

return val == num1 or val == num2

За 100 попыток нашей симуляции игрок выиграл 3 раза и проиграл 97 раз, баланс составил -46у.е.

Существуют и другие виды ставок, например, на 4 номера с коэффициентом 1:8, желающие могут поэкспериментировать самостоятельно. Как нетрудно догадаться, все коэффициенты рассчитаны так, чтобы игрок оказался в минусе. Кажется заманчивым поставить 1уе на номер, чтобы выиграть целых 35уе. Но сумма выигрыша увеличивается в 35 раз, а шанс выигрыша уменьшается в 38 раз - итоговый баланс все равно будет в пользу казино.

Лото 6 из 45

Следующее, что интересно проверить, это лото. Принцип игры довольно прост - в барабане находятся 45 шаров, выпадают случайным образом 6 из них. Цена билета составляет 100р, а выигрыш зависит от количества угаданных шаров. Примерный порядок сумм выигрыша таков: 2 угаданных шара дают выигрыш в 100р, 3 угаданных шара дают 300р, 4 шара - 3000р, 5 шаров - 300.000р и 6 шаров - суперприз порядка 10.000.000р.

Для начала напишем программу выбрасывания шаров:

def lottery(values):

balls = range(1, 45+1)

b1 = balls.pop(random.randint(0, len(balls)-1))

b2 = balls.pop(random.randint(0, len(balls)-1))

b3 = balls.pop(random.randint(0, len(balls)-1))

b4 = balls.pop(random.randint(0, len(balls)-1))

b5 = balls.pop(random.randint(0, len(balls)-1))

b6 = balls.pop(random.randint(0, len(balls)-1))

s = [b1,b2,b3,b4,b5,b6]

# print s

res = list(set(s) & set(values))

return len(res)

Из массива balls 6 раз “достается” случайный элемент, затем определяется число элементов пересечений двух множеств. Затем построим график суммарного выигрыша от количества купленных билетов. Для простоты будем считать что игрок ставит на одни и те же числа (тем более, что это ни на что не влияет).

money = []

money_total = 0

for p in xrange(N):

val = lottery([3,7,12,18,33,28])

if val == 2:

money_total += 100

if val == 3:

money_total += 300

if val == 4:

money_total += 3000

if val == 5:

money_total += 300000

if val == 6:

money_total += 3000000

money.append(money_total)

x = range(0, N)

price = map(lambda x: ticket_price*x, x)

from matplotlib import pyplot as plt

%matplotlib inline

plt.plot(price, money)

Интересно проанализировать работу программы. Если купить 100 билетов (суммарная потраченная сумма будет 10.000р), то это дает около 14 угаданных “двойных” шаров и один угаданный “тройной”. Суммарный выигрыш составляет около 2.000р при потраченных 10.000р.

График выигрыша от потраченной суммы получился практически линейным:

Получается, что если купить билетов на миллион, выигрыш составит 250тыс. “Суперприз” в симуляции так не разу и не выпал, хотя теоретически, он конечно возможен. Кстати, как написано в правилах, призовой фонд составляет 50% с проданных билетов, но “суперприз” выпадает далеко не всегда, так что как и в случае с казино, организаторы разумеется, в выигрыше.

В завершении этой главы, интересно сравнить принципиально разный психологический подход к играм. В лотерее выигрыши потенциально велики, но весьма редки. В казино подход обратный — ставки настроены так, что человек будет выигрывать максимально часто. Условно, сделав 10 игр в казино, человек выиграет 4 раза и проиграет 6 раз. Это позволяет игроку не терять интерес к игре, но в любом случае общий баланс остается негативным — человек будет много раз выигрывать, но и чуть-чуть больше проигрывать.

В итоге, может ли существовать «выигрышная стратегия» в таких случайных играх? Очевидно нет, т.к. ни кость, ни шарик, ни лотерейные билеты, не имеют памяти, и их поведение не зависит от предыдущих итераций. Кстати, этот момент важно помнить - интуитивно, проиграв несколько раз, человек может решить, что вот сейчас-то он “точно” выиграет. Увы нет - рулетка или кубик не имеют памяти, и “не знают” о количестве предыдущих попыток, каждая игра по сути начинается с чистого листа. Можно ли выиграть в азартные игры? Как показывает, симуляция, в принципе можно, теория вероятности допускает. Но недолго — стоит начать играть 2й, 3й,… Nй раз, как баланс пойдет вниз. В долгосрочной перспективе выиграть у казино невозможно.

28. Теорема Пифагора

Теорему Пифагора, можно надеяться, знает каждый школьник. Она была известна еще до нашей эры, и говорит о том, что для прямоугольного треугольника со сторонами a,b,c выполняется равенство a² + b² = c².

Ничего сложного, в этом разумеется, нет. Интерес тут представляет другое - так называемые “пифагоровы тройки” - наборы из целых чисел, удовлетворяющие теореме Пифагора. Если например, мы возьмем треугольник со сторонами a=1 и b=2, то с будет равна 2.236, что не является целым числом. А вот треугольник со сторонами a=3 и b=4 удовлетворяет данному условию (с=5).

Таких наборов достаточно много, например (3,4,5), (5,12,13), (8,15,17) и пр. Интересно посмотреть на них графически - для чего отметим соответствующие точки (3,4, 5,12 и пр) на плоскости. Для вывода напишем простую программу на языке Python - просто переберем все пары X,Y, и отметим точками те, для которых выполняется целочисленное равенство.

import math

import matplotlib.pyplot as plt

r_max = 20

pt_x, pt_y = [], []

cnt = 0

for x in range(1, r_max):

for y in range(1, r_max):

d = math.sqrt(x*x + y*y)

if d.is_integer():

print x, y, d

pt_x.append(x)

pt_y.append(y)

cnt += 1

print("CNT:", cnt)

plt.plot(pt_x, pt_y, 'ro', markersize=1)

plt.axhline(0, color='black')

plt.axvline(0, color='black')

plt.show()

Полученная картинка для чисел от 1 до 20 показана на рисунке. Ей соответствуют треугольники со сторонами 3,4,5, 6,8,10, 5,12,13, 9,12,15, 8,15,17 и 12,16,20. Некоторые треугольники, например 6,8,10 являются умножением на 2 предыдущей тройки 3,4,5. Те тройки, которые не имеют таких пар, называются примитивными.

Если увеличить количество точек по оси, то мы увидим вполне красивый периодический узор:

Разумеется, изображение является симметричным, его можно продлить и в сторону отрицательных чисел.

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

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

Кстати, для более высоких степеней N > 2 выражение aⁿ+ bⁿ = cⁿ уже не имеет решений для целых чисел. Это так называемая “Великая теорема Ферма”, которая занимала умы математиков не одну сотню лет. Теорема была полностью доказана лишь в 1995г Эндрю Уайлсом, а доказательство занимает 130 страниц.

29. ABC-гипотеза

Еще один интересный пример из теории чисел - так называемая ABC-гипотеза. Она была сформулирована в 1985г, и до сих пор является одной из нерешенных задач.

Гипотеза формируется довольно-таки просто: для любого числа s>1 существует ограниченное множество троек взаимно простых чисел A+B = C, таких что (rad(A)*rad(B)*rad(C))^s < C.

Операция rad(N) называется радикалом, и определяется как произведение всех простых чисел, из которых можно составить число N:

rad(7) = 7 (7 = 1*7)

rad(8) = 2 (8 = 2*2*2)

rad(10) = 10 (10 = 5*2)

rad(17) = 17 (17 = 17*1)

rad(30) = 30 (30 = 2*3*5)

Числа называются взаимно простыми, если они не имеют одинаковых простых множителей, например 3 и 20 взаимно-простые (3=3, 20=2*2*5), а 12 и 15 нет (12=2*2*3, 15=3*5).

Так вот, ABC-гипотеза говорит о том, что если взять некое число S>1, например 1.2, то количество равенств A,B,C, удовлетворяющих описанному условию, будет ограниченным. Пример такой тройки чисел: 1+8=9.

Действительно, rad(1)=1, rad(8)=2, rad(9)=3, и 3 > (1*2)^1.2.

Интерес гипотезы, во-первых, в том, что таких равенств мало, для большинства чисел неравенство не выполняется. Например, если взять наугад, равенство 10+6=16, то (rad(10)*rad(6)*rad(16))^1.2 = (10*6*2)^1.2 = 312.6, > 16.

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

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

Разумеется, мы здесь не претендуем на доказательство гипотезы, но можем проверить некоторые значения на языке Python.

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

def prime_factors(n):

factors = []

while n % 2 == 0:

factors.append(int(2))

n = n / 2

for i in range(3, int(math.sqrt(n)) + 1, 2):

# while i divides n , print i ad divide n

while n % i == 0:

factors.append(int(i))

n = n / i

if n > 2:

factors.append(int(n))

return factors

Как можно видеть, функция последовательно делит числа, и если остаток от деления равен 0, результат заносится в массив. Отдельно проверяются два частных случая - если число четное (делится на 2) и если число простое (ни на что не делится кроме себя). В результате функция возвращает массив чисел, напримем для 8 это будет [2,2,2] (2*2*2 = 8).

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

def rad(n):

result = 1

for num in set(prime_factors(n)):

result *= num

return result

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

def not_mutual_primes(a,b,c):

fa, fb, fc = set(prime_factors(a)), set(prime_factors(b)), set(prime_factors(c))

return len(fa.intersection(fb)) == 0 and len(fa.intersection(fc)) == 0 and

len(fb.intersection(fc)) == 0

Здесь мы пользуемся функцией нахождения пересечения множеств языка Python.

Написанных 3х функций достаточно, чтобы написать функцию вывода троек чисел A,B,C:

def calculate(N):

S = 1.2

cnt = 0

for a in range(1, N):

for b in range(1, N):

if b < a: continue

c = a+b

if not_mutual_primes(a, b, c):

if c > (rad(a*b*c))**S:

print("{} + {} = {}".format(a, b, c))

cnt += 1

print("N: {}, CNT: {}".format(N, cnt))

return cnt

Как написано в условии задачи, количество таких троек должно быть ограниченно. Мы можем вывести график количества троек от максимального значения N, тогда будет примерно видно, как быстро он растет.

import matplotlib.pyplot as plt

x_values = range(1,2000,25)

y_values = list(map(calculate, x_values))

plt.plot(x_values, y_values, 'ro', label='count(N), S=1.2')

plt.legend()

plt.show()

Как можно видеть, количество троек реально мало - в диапазоне до 2000 их всего 10:

1 + 8 = 9

1 + 80 = 81

1 + 242 = 243

1 + 288 = 289

1 + 512 = 513

3 + 125 = 128

5 + 1024 = 1029

13 + 243 = 256

49 + 576 = 625

81 + 1250 = 1331

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

Приложение 1 - Вычисления с помощью видеокарты

Еще 20 лет назад, во времена процессоров 80386, пользователям приходилось покупать математический сопроцессор, позволяющий быстрее выполнять вычисления с плавающей точкой. Сейчас такой сопроцессор покупать уже не надо - благодаря прогрессу в игровой индустрии, даже встроенная видеокарта компьютера имеет весьма неплохую вычислительную мощность. Например, даже бюджетный видеочип Intel Graphics 4600 имеет 20 вычислительных блоков, что превышает количество ядер “основного” процессора. Разумеется, каждое ядро GPU по отдельности слабее CPU, но здесь как раз тот случай, когда количество дает преимущество над качеством. Вычисления с помощью GPU сейчас очень популярны - от майнинга биткоинов до научных расчетов, диапазон ценовых решений также различен, от “бесплатной” встроенной видеокарты до NVIDIA Tesla ценой более 100тыс рублей. Поэтому интересно посмотреть, как же это работает.

Есть две основные библиотеки для GPU-расчетов - NVidia CUDA и OpenCL. Первая обладает большими возможностями, однако работает только с картами NVIDIA. Библиотека OpenCL работает с гораздо большим числом графических карт, поэтому мы рассмотрим именно ее.

Основной принцип GPU-расчетов - параллельность вычислений. Данные, хранящиеся в “глобальной памяти” (global & constant memory) устройства, обрабатываются модулями (каждый модуль называется “ядром”, или “kernel”), каждый из которых работает параллельно с другими. Модуль имеет и свою собственную память для промежуточных данных (private memory). Так это выглядит в виде блок-схемы:

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

Для запуска OpenCL-программы из Python, необходимо поставить библиотеку pyopencl, скачать ее дистрибутив можно со страницы https://wiki.tiker.net/PyOpenCL/Installation/Windows. Для установки библиотеки под Windows необходимо скачать файл и ввести команду pip install pyopencl-2018.1.1+cl12-cp27-cp27m-win_amd64.whl для 64-разрядной версии Python, или pip install pyopencl‑2018.1.1+cl12‑cp27‑cp27m‑win32.whl для 32-разрядной.

Для тестирования pyopencl можно запустить программу, выводящую информацию о системе:

import pyopencl as cl

for plat in cl.get_platforms():

print 'Platform: {}'.format(plat.name)

print 'Version: ' + plat.version

devices = plat.get_devices(cl.device_type.ALL)

print 'Devices:'

for dev in devices:

print('\t{} ({})'.format(dev.name, dev.vendor))

flags = [('Version', dev.version),

('Type', cl.device_type.to_string(dev.type)),

('Memory (global), MB', str(dev.global_mem_size/(1024*1024))),

('Memory (local), KB', str(dev.local_mem_size/1024)),

('Max work item dims', str(dev.max_work_item_dimensions)),

('Max work group size', str(dev.max_work_group_size)),

('Max compute units', str(dev.max_compute_units)),

('Driver version', dev.driver_version),

('Device available', str(bool(dev.available))),

('Compiler available', str(bool(dev.compiler_available)))]

for name, flag in flags:

print '\t\t{0:<25}{1:<10}'.format(name + ':', flag)

Если pyopencl и драйвер видеокарты установлен корректно, мы увидим на экране примерно такой вывод (программа запускалась на Macbook Pro):

Platform: Apple

Version: OpenCL 1.2 (Sep 12 2017 16:28:17)

Devices:

Iris Pro (Intel)

Version: OpenCL 1.2

Type: GPU

Memory (global), MB: 1536

Memory (local), KB: 64

Max work item dims: 3

Max work group size: 512

Max compute units: 40

Driver version: 1.2 (Oct 4 2017 01:28:36)

Device available: True

Compiler available: True

Программа удобна своей кросс-платформенностью, один и тот же код может работать и на OSX и на Windows без каких-либо изменений.

Рассмотрим пример: сформировать массив простых чисел от 1 до N. Для решения такой задачи можно использовать алгоритм “решето Эратосфена”, но мы в учебных целях будем решать задачу напрямую, проверяя каждое число отдельно.

Код программы:

import pyopencl as cl

import pyopencl.array as cl_array

import numpy, time

type = cl.device_type.GPU # cl.device_type.CPU

platform = cl.get_platforms()[0]

devices = platform.get_devices(device_type=type)

ctx = cl.Context(devices)

queue = cl.CommandQueue(ctx)

prg = cl.Program(ctx, """

__kernel void primes(__global int *results)

{

int gid = get_global_id(0);

if (gid < 2) return;

int num = gid + 1;

for(int p=2; p<=num/2 + 1; p++) {

if ((num % p) == 0) {

return;

}

int val = atomic_add(&results[0], 1);

results[val+1] = num;

}

""").build()

N = 500000

results = numpy.zeros(N, dtype=numpy.int32)

dest_dev = cl_array.to_device(queue, results)

size = (results.shape[0]-1,)

prg.primes(queue, size, None, dest_dev.data)

res = dest_dev.get()

count = res[0]

primes = numpy.sort(res[1:count+1])

print "Found prime numbers:", count

# for p in primes:

# print p

Разберем текст программы подробнее.

1. Мы создаем контекст устройства , передав ему в качестве параметра тип устройства device_type.GPU. Переменная типа CommandQueue хранит очередь команд, которые будут посланы на устройство.

2. Класс cl.Program получает в качестве параметра программу ядра (kernel) и компилирует ее с помощью вызова функции build. Функция primes написана на языке Си, и будет выполняться параллельно на всех ядрах видеокарты.

3. С помощью функции numpy.zeros мы создаем заполненный нулями массив размера N, затем копируем его на видеокарту с помощью вызова dest_dev = cl_array.to_device. Важно понимать, что существуют две разных копии массива - один в памяти компьютера, второй на видеокарте, с которым и будут выполняться вычисления.

4. Вызовом функции prg.primes выполняется параллельный расчет на видеокарте. При этом в качестве параметра функция получает массив, общий для всех экземпляров функции. Библиотека OpenCL сама распараллеливает вызовы функций, общее количество вызовов равно числу элементов массива, а для получения конкретного идентификатора мы используем функцию get_global_id. В нашем случае идентификатор соответствует числу, которое мы хотим проверить.

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

6. Когда выполнение программы завершено, мы загружаем данные обратно с помощью функции dest_dev.get(). Т.к. ядра запускались параллельно, то простые числа лежат в массиве вперемешку, чтобы их отсортировать, мы вызываем функцию numpy.sort.

Как можно видеть, процесс довольно-таки громоздкий, но оно того стоит. Однопоточная программа, написанная на Си с таким же кодом, выполнялась 14 секунд. Вышеприведенная программа, запускаемая на видеокарте, выполнила расчет за 1.2 секунды.

Разумеется, еще раз стоит повторить, что “игра стоит свеч” лишь в том случае, если задача хорошо распараллеливается на небольшие блоки, в таком случае выигрыш будет заметен.

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

import pycuda.driver as cuda

import pycuda.autoinit

from pycuda.compiler import SourceModule

import numpy, time

# Define kernel

mod = SourceModule('''

__global__ void primes(int *results, int maxVal) {

int gid = blockDim.x*blockIdx.x + threadIdx.x;

if (gid < 2 || gid >= maxVal) return;

int num = gid + 1;

for(int p=2; p<=num/2 + 1; p++) {

if ((num % p) == 0) {

return;

}

int val = atomicAdd(&results[0], 1);

results[val+1] = num;

}''')

N = 500000

n, grid = 1 + N/512, 512

# Host array

h_array = numpy.zeros(n*grid, dtype=numpy.int32)

# Copy host array to device array

d_array = cuda.mem_alloc(h_array.nbytes)

cuda.memcpy_htod(d_array, h_array)

# Invoke kernel

func = mod.get_function("primes")

func(d_array, numpy.int32(N), block=(n, 1, 1), grid=(grid, 1))

# Copy device array to host array

cuda.memcpy_dtoh(h_array, d_array)

count = h_array[0]

primes = numpy.sort(h_array[1:count+1])

print "Found prime numbers:", count

# for p in primes:

# print p

Как можно видеть, принцип тот же. Единственная разница состоит в том, что “вручную” пришлось указать размер блока в 512 значений, и весь массив разбить на блоки, используя параметр grid - библиотека OpenCL делала этот процесс автоматически. И разумеется, синтаксис и названия функций разных библиотек тоже разные, хотя суть остается примерно та же.

Приложение 2 - Тестирование скорости языка Python

Язык Python очень удобен своей краткостью и лаконичностью, возможностью использования большого количества сторонних библиотек. Однако, один из его минусов, который может быть ключевым для математических расчетов - это быстродействие. Python это интерпретатор, он не создает exe-файл, что разумеется, сказывается на скорости выполнения программы.

Рассмотрим простой пример: рассчитаем сумму квадратов чисел от 1 до 1000000. Также выведем время выполнения программы.

Программа на языке Python выглядит так:

import time

start_time = time.time()

s = 0

for x in xrange(1,1000001):

s += x*x

print("Sum={}, T={}s".format(s, time.time() - start_time))

Результаты работы:

Sum = 333333833333500000, T = 0.47s

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

import time

start_time = time.time()

l = xrange(1000001)

s = sum(x*x for x in l)

print("Sum = {}, T = {}s".format(s, time.time() - start_time))

Результаты работы:

Sum = 333333833333500000, T = 0.32s

Быстрее, но лишь чуть-чуть.

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

import numpy as np

start_time = time.time()

a = np.arange(1,1000001, dtype=np.uint64)

s = np.sum(a ** 2)

print("Sum = {}, T = {}s".format(s, time.time() - start_time))

Результат 0.008с - в 40 раз быстрее “обычного” Python-кода.

И наконец, призываем “тяжелую артиллерию”: перепишем программу на языке Cи.

Код выглядит так:

#include

int main() {

clock_t start = clock();

unsigned long long int sum = 0, i;

for (i = 1; i < 1000001; i + +) {

sum += i * i;

}

clock_t end = clock();

printf("Sum = %llu, T = %fs", sum, (float)(end - start)/CLOCKS_PER_SEC);

return 0;

}

Как можно видеть, он ненамного сложнее python-версии. Перед запуском программы, ее надо скомпилировать, выполнив команду C:\GCC\bin\gcc.exe "Appendix-2 - speedTest.c" -o"Appendix-2 - speedTest". Результат очевиден: T = 0.007 секунд. И еще чуть-чуть: добавляем флаг оптимизации по скорости, выполнив команду C:\GCC\bin\gcc.exe "Appendix-2 - speedTest.c" -o"Appendix-2 - speedTest" -O3. Результат: 0.0035 секунд, разница в быстродействии более 100 раз!

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

a = xrange(1000001)

s = 0

for x in a:

s += x

print(s)

Можно использовать встроенную функцию sum:

a = range(1000001)

s = sum(a)

print(s)

Данный код выполняется за 0.02 секунды, т.е. в 5 раз быстрее первого варианта.

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

Для тех, кто захочет максимально широко использовать Python в математических расчетах, крайне рекомендуется ознакомиться с библиотекой numpy, она является стандартом де-факто для всех серьезных расчетов.

Теперь посмотрим, как еще можно ускорить Python-программу с помощью многопоточности и других хитростей.

Приложение 3 - Пример ускорения Python-программы

Рассмотрим практический пример сравнения скорости вычислений на Python. Допустим, есть массив целых чисел, для каждого числа надо найти количество его делителей.

Функцию нахождения делителей числа будем использовать “как есть”, без какой-либо оптимизации: перебираем все числа, если число делится без остатка, увеличиваем счетчик.

def get_dividers(n):

divs = 0

for i in xrange(1, n+1):

if n % i == 0:

divs += 1

return divs

Рассмотрим способы решения задачи. Исходные данные: массив целых чисел:

values = [41212317, 4672313, 4342311, 46512319, 51212317, 5672313, 5342311, 56512319]

Способ-1

Решение “в лоб” (заодно показан способ измерения времени выполнения).

import timeit

start_time = timeit.default_timer()

res = []

for v in values:

res.append(get_dividers(v))

print(res)

print("T =", timeit.default_timer() - start_time)

Время выполнения кода на компьютере с процессором Core i7: 10c.

Способ-2

Используем функцию map, позволяющую применить функцию сразу к массиву.

res = map(get_dividers, values)

print(list(res))

Время выполнения: те же 10с, пока мы выиграли разве что в краткости записи.

Способ-3

Используем многопоточность, класс multiprocessing.Pool, позволяющий разбить вычисления над массивом на несколько потоков (способ работает только в Python 3).

if __name__ == '__main__':

p = multiprocessing.Pool(processes=4)

res = p.map(get_dividers, values)

print(list(res))

Время выполнения: 3.7с, что уже лучше. Интересно, что несмотря на разбивку на 4 процесса, реальное время выполнения уменьшилось лишь вдвое. Многопоточные вычисления - достаточно дорогостоящая в плане накладных расходов операция, множество ресурсов процессора тратится на синхронизацию и передачу данных между процессами.

Способ-4

Используем библиотеку параллельных вычислений ipyparallel. Она имеет большие возможности, даже позволяет производить вычисления параллельно на нескольких компьютерах.

Перед тестом запустим сервер вычислений на том же компьютере, командой

ipcluster start -n 4.

Код запуска вычислений не намного больше:

import ipyparallel as ipp

rc = ipp.Client()

dv = rc[:]

res = dv.map_sync(get_dividers, values)

print(list(res))

Как можно видеть, синтаксис практически тот же, библиотека сама делает всю работу по распараллеливанию вычислений. Впрочем, время выполнения практически то же - 3.4с, ведь мы выполняем вычисления на одном компьютере. При желании читатели могут повторить эксперимент с несколькими компьютерами самостоятельно, описание процесса есть на странице http://ipyparallel.readthedocs.io/en/latest/process.html.

Способ-5

Как говорилось еще в начале книги, Python - это интерпретатор, и код выполняется достаточно медленно. Кардинально повысить скорость можно, используя компиляцию программы перед выполнением, что было реализовано в библиотеке Numba. Ее можно установить, выполнив команду pip install numba.

Метод очень простой и эффективный, и практически не требует изменения кода: всего лишь добавляем перед описанием функции ключевое слово @numba.jit.

Данный код скопирован в отдельный файл Appendix3_test2.py:

import numba

@numba.jit

def get_dividers(n):

divs = 0

for i in xrange(1, n+1):

if n % i == 0:

divs += 1

return divs

Запуск ничем не отличается от предыдущего варианта:

res = map(test2.get_dividers, values)

print(list(res))

Время выполнения в обычном, однопоточном режиме - 1.5с. Это уже заметная разница!

И наконец, повторяем вышеописанные процедуры, но функцию заменяем на новый вариант.

Многопоточность:

p = multiprocessing.Pool(processes=8)

res = p.map(test2.get_dividers, values)

Время выполнения: 1.04c.

И последний вариант, библиотека ipyparallel:

res = dv.map_sync(test2.get_dividers, values)

print(list(res))

Время выполнения: 0.57c.

Таким образом, по сравнению с изначальным вариантом, мы увеличили скорость расчетов почти в 20 раз, при этом сам код изменять или переписывать практически не пришлось.

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

1) По возможности, стараться использовать встроенные функции, имеющиеся в библиотеках math, numpy и др. Они уже написаны на C или C++, и работают максимально эффективно.

2) Использовать библиотеку numba, указав в начале функции префикс @numba.jit. Это самый простой и “безболезненный” способ, он чуть увеличит время запуска программы, но скорость расчетов может вырасти в 10 раз. Правда, этот способ не поможет, если код уже состоит только из вызовов встроенных функций из п1. Если же код содержит циклы, условия, другие проверки, то скорость выполнения вырастет значительно.

3) Если скорость расчета нужно еще увеличить, переходим к более тяжелой артиллерии - многопоточности. Используя библиотеку multiprocessing, можно разбить вычисления на несколько ядер процессора - это может увеличить скорость в 2-4 раза. Но важно иметь в виду, что межпроцессное взаимодействие - сама по себе довольно-таки долгая операция, так что делать ее имеет смысл для действительно долго выполняющихся фрагментов кода. Небольшая функция, вызывающаяся много раз в разных потоках, будет работать в итоге даже медленнее за счет дополнительных “накладных расходов”.

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

5) Если необходимо обработать большой объем данных, целесообразно рассмотреть вычисления на видеокарте (GPU). За счет параллельности обработки, можно запустить программу одновременно на нескольких сотнях ядер видеокарты, что может дать огромный прирост в скорости. Такие библиотеки (OpenCL и CUDA) есть практически под все видеокарты, имеющиеся сегодня в продаже.

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

Приложение 4 - Тестирование программы на языке Fortran-90

Язык Фортран всегда считался специализированным для математических вычислений. Его первая версия появилась в IBM еще в 1957 году, более-менее современными версиями являются Fortran 90 и Fortran 2003.

Поскольку и сейчас существует версия gfortran под Linux и OSX, стало интересно сравнить быстродействие программ. Отдельного компьютера с Linux у меня не было, так что программа была запущена на Raspberry Pi. Один и тот же код без проблем компилировался и под Linux и под OSX.

В качестве теста использовалась несложная программа вычисления числа Пи по формуле Валлиса.

Программа на С выглядит так:

Appendix4_test1.c

#include

double wallis_formula(unsigned long int count)

{

double pi = 1.0, num;

unsigned long int i;

for (i=1; i<=count; i++)

{

num = 4.0*i*i;

pi *= num/(num - 1);

}

return pi*2;

}

int main() {

unsigned long int input_count = 100000000;

double result = wallis_formula(input_count);

printf("Calc done: %.20f\n", result);

return 0;

}

Программа, переписанная на Fortran 90 выглядит немного непривычно, если не сказать, архаично:

Appendix4_test1.f90:

function wallisSeries() result(pi)

real(8) :: pi

real(8) :: num

integer :: ITERATIONS = 100000000

pi = 1.0

do i = 1,ITERATIONS

num = 4.0*i*i

pi = pi*num/(num - 1)

end do

pi = 2*pi

end function

program hello

real(8) :: wallisSeries

print *, "Calc done:", wallisSeries()

end program hello

Впрочем, несмотря на другой синтаксис, основная структура кода та же.

Ключи компиляции также похожи:

gcc -o Appendix4_test1_c Appendix4_test1.c -Ofast

gfortran -o Appendix4_test1_f Appendix4_test1.f90 -Ofast

Результаты запуска обеих программ оказались весьма интересными.

Raspberry Pi

$ gfortran -o Appendix4_test1_f Appendix4_test1.f90 -Ofast

$ time ./Appendix4_test1_f

Calc done: 3.1415926364809641

real 2.222s

$ gcc -o Appendix4_test1_c Appendix4_test1.c -Ofast

$ time ./Appendix4_test1_c

Calc done: 3.14159264306626218044

real 3.441s

Как можно видеть, программа на Си оказалась где-то на 30% медленнее. Все-таки более простой код Fortran позволяет легче получить более эффективный код. Но на OSX с gcc последней версии все оказалось по-другому.

OSX

$ time ./Appendix4_test1_f

Calc done: 3.1415926364809641

real 0.788s

$ time ./Appendix4_test1_c

Calc done: 3.14159264306626218044

real 0.434s

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

Из всего этого следует несколько выводов. Фортран все-таки устаревший язык, и в общем случае, переписывать программу с Си на Фортран, надеясь ее ускорить, смысла нет - прироста производительности не будет. По читаемости и удобству написания кода Фортран тоже проигрывает - все-таки прогресс в создании компиляторов за последние 20 лет весьма значителен. Но с другой стороны, использование Фортран может быть целесообразно на каких-то специализированных устройствах, если современные компиляторы недоступны. Или в тех случаях, когда уже есть большое количество готовых библиотек, которые в противном случае, пришлось бы переписывать.

Заключение

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

Продолжение следует.

Обо найденных неточностях или дополнениях просьба писать на электронную почту dmitryelj@gmail.com. Наличие новой версии можно проверить на странице https://cloud.mail.ru/public/4SGE/oE2EGEWnp.