Метка: математика

3D-бродилка с трассировкой лучей в 184 строки на Python

184 строки – это еще с комментариями! И без использования сторонних движков! На самом деле все будет очень просто, в качестве средства отрисовки мы будем использовать любимую ASCII графику в консоли через библиотеку curses. Изображение будет строиться по принципам игры Wolfenstein 3D.

Готовая бродилка - иллюстрация
Что получим в итоге!

Библиотека curses – это интерфейс к Ncurses, консольной библиотеки для отрисовки всяческих квадратно-текстовых окошек, кнопок и прочих элементов управления в стиле Turbo Pascal или QBasic, кто помнит… Из нее нам понадобится только способность получать размеры терминала и рисовать в нужном месте символ нужного цвета.

Пользователи Windows, к сожалению, в вашей версии Python скорее всего не встроен модуль curses, поэтому вам придется установить пакет. Я использовал такой вариант:

pip install windows-curses

К сожалению, у мне не удалось добиться поддержки цвета на Windows, но я не слишком старался. Это не беда, потому что код был изначально заточен под имитацию оттенков серого через выбор символов разной плотности закраски. На Linux и Macos все должно работать на голом Python и сразу в цвете.

Код, который я вам представлю, является моим портом проекта 3D-Walk, что в свою очередь является портом проекта CommandLineFPS (по ссылке – видео) от javidx9 aka OneLoneCoder. Я внес в код небольшие модификации, исправления и поддержку цвета.

Итак, побежали по коду. Все что, нам надо импортировать:

import curses
import locale
from math import pi, cos, sin

Затем идут некоторые константы:

POS_X, POS_Y, POS_A = 2, 2, 0  # Положение и поворот игрока на карте (начальные)
ROTATION_SPEED = 0.1  # скорость поворота игрока в радианах
SPEED = 0.3  # Скорость игрока вперед назад за одно нажатие

FOV = pi / 2  # Ширина угла обзор в радинах
RESOLUTION = 0.1  # разрешение шага луча
DEPTH = 16  # Максимальная глубина прорисовки

# Наша карта строчкой
MAP = """
################
#..............#
#..............#
#...########...#
#.......#..#...#
#..........#...#
#..............#
#...########...#
#.......#..#...#
#.......####...#
#..............#
#....##..##....#
#...#...#..#...#
#....###...#...#
#..............#
################
"""

Карта представляет собой двумерную схему лабиринта, вид сверху. Карта записана многострочной строкой, где символом «решетка» обозначены непроходимые стены, а «точками» – пустое пространство. Позиция игрока – число блоков (вероятно, дробное) от угла карты. Угол поворота игрока измеряется в радианах, а широта взгляда по горизонтали FOV – это 90 градусов от левого края экрана до правого. В этим числом будет весело поиграть. Графически выглядит примерно так. Правда, здесь ошибка в точке отсчета, думаю, что внимательный читатель заметит.

Координаты и поворот игрока относительно карты
Положение игрока на карте, вид сверху

Главная функция получает параметром объект экрана (рабочего окна) и запускается через curses.wrapper:

def main_3dwalk(screen):
   ...

curses.wrapper(main_3dwalk)

В начале ее работы мы настраиваем curses:

# для корректного отображение юникода
locale.setlocale(locale.LC_ALL, '')

curses.noecho()  # нажатые клавиши не печатаются на экране
curses.curs_set(0)  # курсор убран
curses.start_color()  # цветной режим
curses.use_default_colors()  # стандартная палитра
# инициализация всех цветов!
for i in range(0, curses.COLORS):
    curses.init_pair(i, i, -1)

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

def make_map(string_map):
    rows = string_map.strip().split('\n')
    h = len(rows)
    w = len(rows[0])
    return string_map.replace('\n', ''), w, h

# форматируем карту и получаем ее размеры
level_map, map_width, map_height = make_map(MAP)

На выходе будет:

('#################..............##..............##...########...##.......#..#...##..........#...##..............##...########...##.......#..#...##.......####...##..............##....##..##....##...#...#..#...##....###...#...##..............#################', 16, 16)

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

def get_block(x, y):
    x, y = int(x), int(y)
    if 0 <= x < map_width and 0 <= y < map_height:
        return level_map[y * map_width + x]
    else:
        return '#'

Если координаты попадут по какой-либо причине за пределы карты, что считается, что там всегда глухая стена непроходимого вещества, иными словами – символ '#'. Больше от карты ничего не требуется, только знание, есть ли стена в этой точке пространства или нет!

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

# текущие положение и угол
pos_x, pos_y, pos_a = POS_X, POS_Y, POS_A

exit_flag = False  # флаг выхода
while not exit_flag:
   # получаем размер экрана в каждом кадре, чтобы не глючить, если юзер изменил размер терминала
   screen_height, screen_width = screen.getmaxyx()
   ...

Изображение строится по столбцам. Пока номеру столбца мы находим угол отклонения луча от прямого взгляда. Левая колонка соответствует минимальному углу, то есть pos_a - FOV / 2, а самая правая колонка – pos_a + FOV / 2. Таким образом, вычисляя синус и косинус того угла, мы получаем вектор направления взгляда.

for col in range(screen_width):
    # 1. определим направление луча
    # угол сканирует от pos_a - FOV / 2 до pos_a + FOV / 2
    ray_angle = (pos_a - FOV / 2) + (col / screen_width) * FOV

    # вектор, куда смотрит луч на карте
    eye_x, eye_y = sin(ray_angle), cos(ray_angle)

Вдоль этого направления испускается луч, стартуя с позиции игрока pos_x, pos_y по направлению eye_x, eye_y. Небольшими шажками (размер шага – константа RESOLUTION), мы продвигаемся вдоль луча и проверяем, пользуясь картой, нет ли в этой точке стены. Как только луч натыкается на стену, в этот момент фиксируется дистанция, цикл прерывается, и алгоритм переходит к следующей колонке. Правильнее было бы называть его не трассировкой лучшей, а чем-то сродни ray casting. Вот код для определения дистанции:

# 2. Ищем ближайшую стену и дистанцию до нее
distance = 0.0
# пока не достигли стены и дистанция менее предельной
while distance < DEPTH:
    # луч делает шаг вперед
    distance += RESOLUTION
    
    # "текущее" положение на луче
    test_x = int(pos_x + eye_x * distance)
    test_y = int(pos_y + eye_y * distance)
    
    # смотрим карту, есть ли там стена или край
    if get_block(test_x, test_y) == '#':
        break  # стена. расстояние в distance

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

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

ceiling = int(screen_height / 2 - screen_height / distance)  # высота потолка
floor = int(screen_height - ceiling)  # высота пола
Сканирование колонки

Теперь остается только заполнить колонку сверху вниз. Потолок мы заливаем символом одиночной кавычки красного цвета. А ниже потолка идет серая стена. Цвет и символ заливки стены зависит от дистанции. Сюда прекрасно подойдут квадраты разной плотности, благо они есть в Unicode. Под стеной начинается пол зеленого цвета. Для выразительности символ, который образует пол, тоже зависит от дистанции от нижней кромки экрана, что косвенно отвечает дистанции этого кусочка пола от игровой камеры.

# рисуем вертикальную линию
for row in range(screen_height):
    if row <= ceiling:  # Ряд выше или равен границе потолка
        shade = '`'
        color = curses.COLOR_RED
    elif floor >= row > ceiling:  # Кусок стены
        if distance <= DEPTH / 4:  # совсем близко
            shade = "█"
        elif distance <= DEPTH / 3:  # ближе
            shade = "▓"
        elif distance <= DEPTH / 2:  # дальше
            shade = "▒"
        elif distance <= DEPTH:  # еще дальше
            shade = "░"
        else:
            shade = " "  # совсем далеко
        # оттенок цвета, нормированный на предельную дистанцию
        color = color_by_distance(1 - (distance / DEPTH))
    else:
        # Оттенок пола, чем ближе к низу экрана, тем гуще заливка
        b = 1 - (row - screen_height / 2) / (screen_height / 2)
        if b < 0.25:
            shade = '#'
        elif b < 0.5:
            shade = "x"
        elif b < 0.75:
            shade = "."
        else:
            shade = ' '
        color = curses.COLOR_GREEN
    # заменяем символ в row/col на shade с цветом color
    screen.insstr(row, col, shade, curses.color_pair(color))

После циклов остается только отрисовать все изменения на экране:

screen.refresh()  # отрисуем все на экране

Дальше мы ожидаем нажатия игроком клавиш управления. Для перемещения используются клавиши WASD, W/S – вперед и назад, A/D – повороты влево и вправо. Esc – выход.

Не забудьте переключиться на английскую раскладку. На русской раскладке игра не будет реагировать!

Преобразовав код клавиши в символ, мы поворачиваем камеру либо смещаем игрока назад или вперед вдоль направления взгляда. Если новое положение игрока оказалось «внутри» стены, то такое движение отменяется. Так просто происходит обработка столкновений со стенами. Вот код, решающий задачу движения игрока:

key_code = screen.getch()  # ждем клавишу и обрабатываем
key = chr(key_code) if 0 < key_code < 256 else 0
if key in ('w', 's'):
    # шаг вперед или назад
    dx, dy = sin(pos_a) * SPEED, cos(pos_a) * SPEED
    if key == 's':  # назад - обратим вектор
        dx, dy = -dx, -dy

    # сдвинем игрока в направлении
    pos_x += dx
    pos_y += dy
    if get_block(pos_x, pos_y) == '#':  # упс, мы в стене
        # отменим движение
        pos_x -= dx
        pos_y -= dy
elif key == 'a':  # поворот налево
    pos_a -= ROTATION_SPEED
elif key == 'd':  # поворот направо
    pos_a += ROTATION_SPEED
elif key_code == 27:  # esc
    break  # выход из игры

В конце не забудем завершить работу curses корректно, восстановив все настройки терминала:

curses.endwin()

Вот и все! Наша бродилка готова! Теперь мы знаем немного больше про 3D графику!

Код программы я залил на gist.github.com. Там два файла: цветная версия для macOS и Linux, и черно-белая для Windows. Наслаждайтесь. Возможно, кто-то из читателей модифицирует этот код, добавив больше цветов, текстур, возможно, противников 🙂

Специально для канала @pyway. Подписывайтесь на мой канал в Телеграм @pyway 👈 

Комплексные числа в Python. Бонус: фрактал

В Python есть встроенный тип данных complex, который моделируют комплексные числа. По-моему, теория комплексных чисел – настоящий прорыв в математике, оказавший колоссальное влияние на современную физику. Неудивительно, что комплексные числа оказались в стандартной библиотеке такого языка, как Python.

Фрактал

Ко́мпле́ксные чи́сла — числа вида a + bi, где a, b — вещественные числа, i — мнимая единица, то есть число, для которого выполняется равенство: i2 = -1.

Графическое представление комплексного числа
Графическое представление комплексного числа

Цель этой статьи – рассказать про комплексные числа с точки зрения программирования на Python, если вам хочется узнать математическую теорию или просто освежить воспоминания о комплексных числах, то смело переходите по ссылке: Комплексные числа (откроется в новом окне), а потом возвращайтесь сюда.

В Python комплексное число состоит из пары чисел с плавающей запятой, которые отвечают за реальную и мнимые части. В исходнике на Си – это структура из пары чисел типа double. (Вы же помните, что float в Python это числа двойной точности?):

typedef struct {
    double real;
    double imag;
} Py_complex;

✅ Давайте посмотрим, какими способами мы можем задать комплексное число. Во-первых, с помощью встроенной функции complex(real[, imag]):

>>> complex()   # нуль! (комплексный)
0j
>>> complex(1)   # из обычного числа int или float (мнимая часть будет 0)
(1+0j)
>>> complex(2, 3)   # из пары чисел (реальная и мнимая части)
(2+3j)
>>> complex('2+3j')  # из строки (без пробелов)!
(2+3j)

❌ Когда создаете комплексное число из строки, то в ней не должно быть пробелов, иначе будет ошибка ValueError!

✅ Во-вторых, что очень круто, комплексное число можно задать особым синтаксисом в форме a+bj. Примеры:

>>> 8+5j
(8+5j)
>>> -1-1j
(-1-1j)
>>> 0j
0j
>>> -4.4e+5 + 1.5e+6j
(-440000+1500000j)

❌ Но! Нельзя задавать их так:

>>> 1+j
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'j' is not defined
>>> a = 5
>>> 2 + aj
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'aj' is not defined

В первом случае надо писать 1j, вместо j, потому что j – это может быть название какой-то иной переменной, ведь имя j не зарезервировано под мнимую единицу. Чтобы задать мнимое число, обязательно требуется перед j ставить цифры.

Во втором случае, мы желаем задать мнимую часть через переменную, но Python думает, что у нас должна быть переменная с именем aj, которой нет. Правильно писать в таком случае a*1j.

✅ Правильные варианты:

>>> 1+1j
(1+1j)
>>> a = 5
>>> 2 + a * 1j
(2+5j)

Кстати, буква j может быть и заглавной J.

Какие свойства есть у типа complex?

Можно извлекать из комплексного числа его мнимую (imag) и реальную (real) части – это обычные float:

>>> z = -3+7j
>>> z.real, z.imag
(-3.0, 7.0)
>>> z == z.real + z.imag * 1j
True

Комплексно-сопряженное число – то же самое, только с другим знаком у мнимой части:

>>> z.conjugate()
(-3-7j)
>>> (4-1j).conjugate()
(4+1j)

Модуль комплексного числа – фактически длина вектора на комплексной плоскости – вычисляется обычной функцией abs:

>>> abs(4+3j)
5.0

Операции

✅ К комплексным числам применимы обычные арифметический операторы, такие как +, -, *, /.

Как вы догадались из формы синтаксиса 1+2j, можно без проблем складывать комплексные числа с обычными вещественными float или целыми int. Т.е. полное комплексное число в этой форме представляется суммой действительного числа и чисто мнимого. Аналогично, умножение комплексного числа на вещественное просто масштабирует пропорционально его компоненты на это число.

Пару комплексных чисел можно складывать и вычитать. Это просто, они работают подобно двумерным векторам на плоскости. Комплексные числа можно умножать и делить. Тут математика несколько сложнее, не буду повторять ее, читайте вики. А примеры кода вот:

>>> z1 = 3 + 4j
>>> z2 = 5 - 2j
>>> z1 + z2
(8+2j)
>>> z1 - z2
(-2+6j)
>>> z1 * z2
(23+14j)
>>> z1 / z2
(0.24137931034482757+0.896551724137931j)

Не забывайте ставить скобки:

>>> 1+2j * 5-4j
(1+6j)
>>> (1+2j) * (5-4j)
(13+6j)

❌ Целочисленные деления, взятие остатка и подобное не применимы к комплексным числам.

✅ Комплексные числа можно сравнивать на равенство и неравенство.

>>> 1+2j == 2j+1
True
>>> 3 + 4j != 2 - 2j
True

❌ Но нельзя к ним применять знаки больше, меньше.

>>> 1j > 2j
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: '>' not supported between instances of 'complex' and 'complex'

✅ Можно возводить одно комплексное число в степень другого следующими способами*:

>>> z1 ** z2
(3046.7438186304985+19732.04597993193j)
>>> pow(z1, z2)
(3046.7438186304985+19732.04597993193j)

Комплексный 0 в степени комплексного 0 даст… барабанная дробь… единицу:

>>> 0j ** 0j
(1+0j)

*) О многозначности замолвлено будет ниже, математики приберегите свои помидоры, пока не дочитаете до конца.

cmath

Модуль, который отвечает за стандартные функции над комплексными числами называется cmath (документация на английском). Обычный math не умеет извлекать корень из минус единицы, а cmath – запросто:

>>> import math, cmath
>>> math.sqrt(-1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: math domain error
>>> cmath.sqrt(-1)
1j

Рассмотрим основные функции из cmath:

cmath.phase(x) – фаза (или у нас она более известна, как аргумент Arg z) – эквивалент math.atan2(x.imag, x.real) – иными словами, это угол поворота φ вектора на комплексной плоскости. Результат лежит в промежутке [-ππ]. При этом разрез на комплексной области (т. е. луч, через который результат функции разрывается и перепрыгивает) выбирается вдоль отрицательной части реальной оси).

Представление комплексного числа в полярных координатах
Представление комплексного числа в полярных координатах
>>> phase(complex(-1.0, 0.0))
3.141592653589793
>>> phase(complex(-1.0, -0.0))  # минус ноль! подходим с другой стороны к точке разреза
-3.141592653589793

За модуль отвечает обычная abs (без cmath).

cmath.polar(x) – перевод комплексного числа в полярные координаты (кортеж из двух float): (abs(x), cmath.phase(x)).

cmath.rect(rphi) – наоборот выдает комплексное число по его полярным координатам – эквивалент: r * (math.cos(phi) + math.sin(phi)*1j).

cmath.sqrt(x) – корень квадратный из комплексного числа.

>>> cmath.sqrt(1j)
(0.7071067811865476+0.7071067811865475j)

Те, кто разбирается в математике сразу заметят, что многие функции от комплексной переменной – многозначны. (Видео о многозначных функциях) В частности у корня квадратного из комплексного числа всегда ровно два ответа. Но! Python всегда возвращает ровно одно значение. Он берет всегда «главную ветвь», а всякие вращения на 2πk/n остаются на совести программиста. Вот теория, как брать корень из комплексного числа. Я тщетно пытался нагуглить функцию, возвращающую множество значений корня n-й степени, и в итоге написал свою:

import cmath

def roots(z: complex, n):
    assert isinstance(n, int) and n > 1
    r, phi = cmath.polar(z)
    r **= 1 / n
    for k in range(n):
        yield cmath.rect(r, (phi + 2 * cmath.pi * k) / n)

z1, z2 = roots(1j, 2)
print(z1, z2)
print(cmath.isclose(z1 * z1, 1j))
print(cmath.isclose(z2 * z2, 1j))

# (0.7071067811865476+0.7071067811865475j) (-0.7071067811865477-0.7071067811865475j)
# True
# True
Разрез комплексной плоскости
Разрез комплексной плоскости

Функция cmath.isclose(ab*rel_tol=1e-09abs_tol=0.0) проверяет близки ли два комплексных числа между собой с некоторой точностью (относительной или абсолютной). Так как complex строится на базе чисел с плавающей точкой, то из-за ошибок округления редко можно получить точное равенство, приходится проверять близость. Код эквивалентен: abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol). Числа NaN не близки ни к кому, включая самих себя. А inf и -inf близки только каждое само к себе.

cmath.exp(x),cmath.log(x[, base]), cmath.log10(x), cmath.acos(x), cmath.asin(x), cmath.atan(x), cmath.cos(x), cmath.sin(x), cmath.tan(x), cmath.acosh(x), cmath.asinh(x), cmath.atanh(x), cmath.cosh(x), cmath.sinh(x), cmath.tanh(x) –   обычный набор функций, только для комплексных чисел. Каждая из функций возвращает один результат. Думаю, если вам действительно нужны в работе комплексные функции, значит вы и так знаете, как достать все значения функций.

Множество Мандельброта

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

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

Реккурентное соотношение
Реккурентное соотношение

Иными словами мы берем каждую точку с = x + yj на выбранном участке комплексной плоскости и возводим в квадрат, потом снова прибавляем c, опять возводим в квадрат и так несколько раз. Смотрим, улетает ли результат в бесконечность. Если не улетает, значит точка принадлежит множеству, красим ее в черный, а если улетает, то красим в белый.

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

 pip install pillow tqdm 

Пишем код:

from PIL import Image
from tqdm import tqdm

W, H = 1024, 768  # размеры картинки
ITER = 100  # максимальное число итераций, чтобы убедиться расходится или нет формула в данной точке
LIMIT = 2.0  # предельное значение, выше которого уже наверняка расходится

img = Image.new('RGB', (W, H))

for px in tqdm(range(W)):
    for py in range(H):
        # преобразование координат
        x = px / W * 3 - 2  # x = -2..1
        y = py / H * 2 - 1  # y = -1..1

        color = (0, 0, 0)  # черный
        c = x + 1j * y  # смещение из координат
        z = 0j  # начальная точка
        for n in range(ITER):
            z = z ** 2 + c
            if abs(z) > LIMIT:  # разошлось
                color = (255, 255, 255)  # белый цвет
                break
        img.putpixel((px, py), color)

img.save('mand0.png')  # сохраним
img.show()  # покажем

Результат:

Фрактал в ЧБ
Фрактал в ЧБ

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

from PIL import Image
from tqdm import tqdm
import math

W, H = 1024, 768  # размеры картинки
ITER = 1000  # максимальное число итераций, чтобы убедиться расходися или нет формула в данной точке
LIMIT = 2.0  # предельное значение, выше которого уже расходится

img = Image.new('RGB', (W, H))

# создадим палитру от числа итераций
palette = [
    (
        int(255 * math.sin(i / 50.0 + 1.0) ** 2),
        int(255 * math.sin(i / 50.0 + 0.5) ** 2),
        int(255 * math.sin(i / 50.0 + 1.7) ** 2)
    ) for i in range(ITER - 1)
]
palette.append((0, 0, 0))  # последняя итерация - значит мы внутри - черный

for px in tqdm(range(W)):
    for py in range(H):
        # преобразование координат
        x = px / W * 3 - 2  # x = -2..1
        y = py / H * 2 - 1  # y = -1..1

        c = x + 1j * y  # смещение из координат
        z = 0j  # начальная точка
        for n in range(ITER):
            z = z ** 2 + c
            if abs(z) > LIMIT:  # разошлось
                break
        img.putpixel((px, py), palette[n])


img.save('mand1.png')  # сохраним
img.show()  # покажем

Вот итог:

Цветной фрактал
Цветной фрактал

Вот такая красота получается из простейшей формулы!

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

                if (z * z.conjugate()).real > 4.0:  # вместо abs(z) > 2.0!
                    break

От кадра к кадру будем менять область отображения. Генератор фрактала для заданной области описан функцией:

def mandelbrot(w, h, palette, x1=-2, y1=-1, x2=1, y2=1):
    img = Image.new('RGB', (w, h))

    dx = x2 - x1
    dy = y2 - y1
    iters = len(palette)

    n = 0
    for px in range(w):
        for py in range(h):
            # преобразование координат
            x = px / w * dx + x1
            y = py / h * dy + y1

            c = x + 1j * y  # смещение из координат
            z = 0j  # начальная точка
            for n in range(iters):
                z = z ** 2 + c
                if (z * z.conjugate()).real > 4.0:  # разошлось
                    break
            img.putpixel((px, py), palette[n])
    return img

Остаетя задать закон движения камеры и сохранить пачку кадров в GIF.

x, y, r = 0, 0, 5  # от
x_tar, y_tar, r_tar = -0.74529, 0.113075, 1.5e-6  # до

frames = []
for _ in tqdm(range(FRAMES)):
    frames.append(mandelbrot(W, H, palette, x - r, y - r, x + r, y + r))
    x += (x_tar - x) * 0.1
    y += (y_tar - y) * 0.1
    r += (r_tar - r) * 0.1

frames[0].save('mandel.gif', save_all=True, append_images=frames[1:], duration=60, loop=0)

Полный код генератора анимации здесь.

Анимация погружения в фрактал
Анимация погружения в фрактал

Специально для канала @pyway. Подписывайтесь на мой канал в Телеграм @pyway 👈 

NumPy-бродкастинг

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

Бродкастинг (broadcasting) – автоматическое расширение размерности (ndim) и размеров (shape) массивов, при совершении операций (сложение, умножение и подобные) над массивами с разными размерами или размерностями, при условии, что они совместимы с правилами бродкастинга.

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

>>> y = np.array([2] * 3)
>>> y
array([2, 2, 2])
>>> x * y
array([2, 4, 6])

>>> x = np.array([1, 2, 3])
>>> x * 2
array([2, 4, 6])

NumPy не заставляет нас вручную превращать скаляр (2) в вектор [2, 2, 2], он сам добавляет размерность и клонирует содержимое нужно число раз, а потом уже производит поэлементное умножение.

Бродкастинг скалара в вектор

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

Бродкастинг работает и в более сложных ситуациях по четким правилам.

Правила бродкастинга

Буду сразу объяснять на примере. Есть два массива, которые мы желаем сложить:

>>> a = np.ones((8, 1, 6, 1))
>>> b = np.ones((7, 1, 5))
>>> a.shape
(8, 1, 6, 1)
>>> b.shape
(7, 1, 5)

# (a + b) = ?

Сначала размеры (shape) массивов выстраивается друг над другом, выравнивая по правому краю. Напомню, что справа у нас самая «глубокая» размерность.

A         (4d массив):  8 x 1 x 6 x 1
B         (3d массив):      7 x 1 x 5

Затем NumPy идет справа налево, поэлементно сравнивая каждый размер операндов. Два размера считаются совместимыми, если они равны или один из них равен единице (1). Если два размера несовместимы, бродкастинг не пройдет, возникнет ошибка.
ValueError: operands could not be broadcast together with shapes

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

A         (4d массив):  8 x 1 x 6 x 1

B         (3d массив):      7 x 1 x 5
B'        (4d массив):  1 x 7 x 1 x 5
B' = [ B ] 

Мы видим, что в примере два массива полностью совместимы для бродкастинга – (8 над 1, 1 над 7, 6 над 1, 1 над 5): в каждом из столбиков есть единичка.

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

A         (4d массив):  8 x 1 x 6 x 1
B         (3d массив):      7 x 1 x 5
Результат (4d массив):  8 x 7 x 6 x 5

Т. е. у A на глубоком уровне по одному числу, а у B – по 5 штук. Примерно так:

A = [ [ [ [123], ... ] ], ... ]
B = [ [ [456, 456, 456, 456, 456] ], ... ]

A' = [ [ [ [123, 123, 123, 123, 123], ... ] ], ... ] 

Значит внутренний подмассив [123] у A тоже раскопируется в 5 значений [123, 123, 123, 123] и, таким образом, станет совместим с внутренним подмассивом B, где уже было 5 чисел.

Как только все размерности выровнены путем «копирования», то можно делать любую операцию поэлементно. Форма результата будет равна форме операндов. В итоге:

>>> (a + b).shape
(8, 7, 6, 5)

Еще примеры

Примеры привожу по следам документации. Такой код:

import numpy as np

a = np.array([[0, 0, 0],
              [10, 10, 10],
              [20, 20, 20],
              [30, 30, 30]])
b = np.array([1, 2, 3])

print(a + b)
# [[ 1  2  3]
#  [11 12 13]
#  [21 22 23]
#  [31 32 33]]

Работает по следующей схеме:

Сложение матрицы и вектора

Еще больше примеров, как получается финальный размер после операции:

Картинка  (3d массив)	256 x	256 x	3
Масштаб   (1d массив)	 	 	3
Результат (3d массив)	256 x	256 x	3

A      (2d array):  5 x 4
B      (1d array):      1
Result (2d array):  5 x 4

A      (2d array):  5 x 4
B      (1d array):      4
Result (2d array):  5 x 4

A      (3d array):  15 x 3 x 5
B      (3d array):  15 x 1 x 5
Result (3d array):  15 x 3 x 5

A      (3d array):  15 x 3 x 5
B      (2d array):       3 x 5
Result (3d array):  15 x 3 x 5

A      (3d array):  15 x 3 x 5
B      (2d array):       3 x 1
Result (3d array):  15 x 3 x 5

Примеры несовместимости

А вот примеры несовместимости:

A      (1d array):  3
B      (1d array):  4 # не совпадают 3 и 4 (и ни одна из них не 1)

A      (2d array):      2 x 1
B      (3d array):  8 x 4 x 3 # второй столбик справа не совпадает (2 и 4) 

Такой код на практике даст ошибку:

import numpy as np

a = np.array([[0, 0, 0],
              [10, 10, 10],
              [20, 20, 20],
              [30, 30, 30]])
b = np.array([0, 1, 2, 3])

print(a + b)
# ValueError: operands could not be broadcast together with shapes (4,3) (4,)

Потому что тут бродкастинг не работает, так как нарушены его правила:

Тут бродкастинг не работает

Опасность бродкастинга

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

Я пока не нашел опции запретить бродкастинг в NumPy, а ответы со Stackoverflow вроде [1], [2] оказались НЕРАБОЧИМИ. Как всегда, будьте осторожны!

Специально для канала @pyway. Подписывайтесь на мой канал в Телеграм @pyway 👈 

Коронавирус: предсказание на Python

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

КДПВ: вирус + Китай + график

Нам понадобятся библиотеки:

  • pandas – для загрузки данных
  • matplotlib – для построения графика
  • scipy – для построения предсказания
  • numpy – для работы с массивами

Установите их, если они еще не установлены, набрав в терминале:

pip install pandas matplotlib scipy numpy

Подключим модули:

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from datetime import timedelta

Я собрал данные по заболевшим в последние дни в CSV файл и закачал его на свой сайт по адресу: https://tirinox.ru/wp-content/uploads/2020/01/corona.csv

Pandas умеет читать данные прямо из интернета по URL:

df = pd.read_csv('https://tirinox.ru/wp-content/uploads/2020/01/corona.csv', parse_dates=['day'])

print(df.head())

""" Вывод:
         day  infected  dead
0 2020-01-20       278     4
1 2020-01-21       326     6
2 2020-01-22       547     8
3 2020-01-23       639    14
4 2020-01-24       916    25
"""

Мы будем экстраполировать (то есть предсказывать) количество заболевших на будущие дни с учетом текущей динамики. Но какой функцией нам воспользоваться, чтобы она максимально точно вписалась в данные о заболевших? Построим ка график:

df.plot(kind='bar', x='day', y=['infected', 'dead'])
plt.show()
График роста вируса по дням
infected = зараженные люди, dead = погибшие

Очевидно, что прирост с каждым днем все больше и больше. Явно нелинейная зависимость! Это и понятно. Представим, что 1 человек в среднем заражает двоих. Те двое, каждый заразив еще 2, добавят 4 новых человека. А те 4 заразят 8. Затем 16, 32, 64. К этому нужно не забыть добавить уже больных.

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

def fit_func_exp(x, a, b, c):
    return a * np.exp(c * x) + b

Теперь вычислим оптимальные параметры, обеспечивающие наименьшую ошибку:

# зависимая переменная - ее будем предсказывать
infected = df['infected']

# дни - преобразуем их в целые числа от 0 до максимального
days = range(len(df['day']))

# у нас 3 параметра в функции: a, b, c – начальное приближение - единицы
p0 = np.ones(3)

# впишем кривую в наши данные
(a, b, c), _ = curve_fit(fit_func_exp, days, infected, p0=p0)

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

# предскажем динамику вируса на 60 дней (начало соответствует 20 января)
MAX_DAY = 60

x_days = np.linspace(0, MAX_DAY - 1, MAX_DAY)
y_infect = fit_func_exp(x_days, a, b, c)

Построим график и убедимся, что он хорошо описывает экспериментальные данные:

plt.xlabel('Дни')
plt.ylabel('Больные')

# график в log шкале
plt.yscale('log')

# это данные текущей статистики - нарисуем их синими точками
plt.scatter(days, infected, marker='D', label='Реальные')

# это красная линия – предсказание (первые 22 дня)
plt.plot(x_days[:22], y_infect[:22], 'r', label='Предсказание')
plt.legend()

plt.show()

Вот, что у нас получится:

График предсказания.

Ого! Рост поражает: сотни, тысячи, десятки тысяч! Узнаем, на какой день число зараженных людей достигнет населения всей Земли:

# население Земли
EARTH_POPULATION = 7_530_000_000

# найдем номер дня, когда количество зараженных достигнет популяции Земли
doom_index = np.argmax(y_infect >= EARTH_POPULATION)
doom_day = x_days[doom_index]

# вычислим дату
day0 = df['day'][0]
doom_date = day0 + timedelta(days=int(doom_day))

# дата конца!
print(f'Doom date: {doom_date:%d.%m.%Y}')

Doom date: 13.03.2020

13 марта…

Полный код загрузил в gist.

UPD: добавил табличку:

+------------+------------------+----------------+
|    Дата    | Число заболевших | Число погибших |
+------------+------------------+----------------+
| 20.01.2020 |        54        |       7        |
| 21.01.2020 |       239        |       9        |
| 22.01.2020 |       492        |       11       |
| 23.01.2020 |       839        |       15       |
| 24.01.2020 |       1314       |       21       |
| 25.01.2020 |       1963       |       30       |
| 26.01.2020 |       2853       |       45       |
| 27.01.2020 |       4072       |       68       |
| 28.01.2020 |       5740       |      104       |
| 29.01.2020 |       8024       |      161       |
| 30.01.2020 |      11150       |      252       |
| 31.01.2020 |      15431       |      394       |
| 01.02.2020 |      21293       |      618       |
| 02.02.2020 |      29317       |      971       |
| 03.02.2020 |      40304       |      1528      |
| 04.02.2020 |      55347       |      2405      |
| 05.02.2020 |      75942       |      3787      |
| 06.02.2020 |      104139      |      5965      |
| 07.02.2020 |      142745      |      9398      |
| 08.02.2020 |      195601      |     14807      |
| 09.02.2020 |      267968      |     23331      |
| 10.02.2020 |      367047      |     36763      |
| 11.02.2020 |      502698      |     57930      |
| 12.02.2020 |      688423      |     91288      |
| 13.02.2020 |      942703      |     143854     |
| 14.02.2020 |     1290844      |     226690     |
| 15.02.2020 |     1767494      |     357228     |
| 16.02.2020 |     2420088      |     562938     |
| 17.02.2020 |     3313572      |     887108     |
| 18.02.2020 |     4536865      |    1397953     |
| 19.02.2020 |     6211706      |    2202971     |
| 20.02.2020 |     8504777      |    3471566     |
| 21.02.2020 |     11644280     |    5470689     |
| 22.02.2020 |     15942657     |    8621023     |
| 23.02.2020 |     21827679     |    13585498    |
| 24.02.2020 |     29885019     |    21408803    |
| 25.02.2020 |     40916535     |    33737214    |
| 26.02.2020 |     56020077     |    53165030    |
| 27.02.2020 |     76698735     |    83780495    |
| 28.02.2020 |    105010434     |   132026095    |
| 29.02.2020 |    143772730     |   208054274    |
| 01.03.2020 |    196843216     |   327863828    |
| 02.03.2020 |    269503423     |   516666580    |
| 03.03.2020 |    368984436     |   814192761    |
| 04.03.2020 |    505186524     |   1283051543   |
| 05.03.2020 |    691664408     |   2021906043   |
| 06.03.2020 |    946976216     |   3186235247   |
| 07.03.2020 |    1296530371    |   5021051835   |
| 08.03.2020 |    1775114218    |   7912460814   |
| 09.03.2020 |    2430356032    |  12468908548   |
| 10.03.2020 |    3327464945    |  19649219634   |
| 11.03.2020 |    4555720506    |  30964364746   |
| 12.03.2020 |    6237357709    |  48795417935   |
| 13.03.2020 |    8539731720    |  76894611953   |
| 14.03.2020 |   11691972927    |  121174929895  |
| 15.03.2020 |   16007789810    |  190954388899  |
| 16.03.2020 |   21916688952    |  300916853606  |
| 17.03.2020 |   30006719188    |  474201998218  |
| 18.03.2020 |   41082993743    |  747274645537  |
| 19.03.2020 |   56247814448    | 1177598150075  |
+------------+------------------+----------------+

Предостережение

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

Узнавать оперативно о статистике вируса можно через телеграм бот @NovelCoronaVirusBot.

Желаю вам здоровья и не подхватить даже обычной простуды! Меньше бывайте в людных местах, часто мойте руки с мылом, носите маски, при любых симптомах обращайтесь к врачу!

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

Специально для канала @pyway. Подписывайтесь на мой канал в Телеграм @pyway 👈 

​​Анимация Jupyter Notebook

Сегодня мы будем анимировать график прямо внутри Jupyter Notebook. Сперва сделаем плавную отрисовку графика. Переключим режим отображения графиков в notebook:

%matplotlib notebook

Импортируем все, что нужно:

import matplotlib.pyplot as plt
from matplotlib import animation
import numpy as np

Сгенерируем наши данные:

# время (200 точек)
t = np.linspace(0, 2 * np.pi, 200)
x = np.sin(t)  # синусоида

Создадим пустой график:

fig, ax = plt.subplots()
# пределы отображения
ax.axis([0, 2 * np.pi, -2, 2])
l, = ax.plot([], [])

Функция animate будет вызываться при отрисовка каждого кадра, аргумент i – номер кадра:

def animate(i):
    # рисуем данные только от 0 до i
    # на первом кадре будет 0 точек, 
    # а на последнем - все
    l.set_data(t[:i], x[:i])

Запускаем анимацию:

fps = 30  # карды в сек
# frames - число кадров анимации
ani = animation.FuncAnimation(fig, animate, frames=len(t), interval=1000.0 / fps)

Если мы хотим анимировать сами данные, например, заставить синусоиду «плясать», то на каждом шаге перегенерируем данные заново, используя переменную i:

def animate(i):
    x = np.sin(t - i / len(t) * np.pi * 2) * np.sin(t * 15)
    l.set_data(t, x)

Можно сохранить в GIF:

ani.save('myAnimation.gif', writer='imagemagick', fps=30)

Сам ноутбук я загрузил на GitHub, но поиграться онлайн с ним не получится, надо скачать себе и запустить локально. Анимированные графики отрисовываются в реальном времени, поэтому требуют достаточно много ресурсов. Пример 3D анимации:

Специально для канала @pyway. Подписывайтесь на мой канал в Телеграм @pyway 👈