Метка: расчеты на python

Коронавирус: предсказание на 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 👈 

Decimal числа. Отличия от float

После рассказа про float меня просили рассказать про Decimal. Узнаем же, что это за зверь, как он устроен внутри и как с ним работать. Итак, Decimal – это класс из стандартного модуля decimal. Он представляет собой число с плавающей точкой, как и float. Да, именно с плавающей, потому что некоторые, я слышал, думают, что это число с фиксированной точкой.

Однако, Decimal имеет ряд существенных отличий от float.

Цель

Тип Decimal создан, чтобы операции над рациональными числами в компьютере выполнялись также, как они выполняются людьми, как их преподают в школе. Иными словами, чтобы все-таки 0.1 + 0.1 + 0.1 == 0.3. Из-за ошибок представления, float приводит к утере точности, и такие простые на первый взгляд равенства не выполняются. А это может быть критично в высокоточных научных вычислениях, и главное в сфере бизнеса и финансов!

Внутреннее устройство

float – реализован по стандарту IEEE-754 как число с плавающей запятой двойной точности (64 бита) на основании 2. Реализация таких чисел заложена прямо в железо почти любого современного процессора. Поэтому float в Python работает примерно также, как и double в С, С++, Java и прочих языках. И имеет такие же ограничения и «странности». Так как поддержка float имеет аппаратный характер, то его быстродействие сравнительно велико.

Decimal – число с плавающей точкой с основанием экспоненты – 10, отсюда и название (decima лат. – десятая часть, десятина).

Он реализован по стандарту IBM: General Decimal Arithmetic Specification Version 1.70 – 7 Apr 2009, который в свою очередь основан на стандартах IEEE. По поводу реализации, в исходниках CPython я нашел два варианта: на чистом Python и с помощью Си-библиотеки libmpdec. Обе реализации есть в кодовой базе, хотя в последних версиях Python 3 используется именно Си-версия, очевидно, она в разы быстрее! Видите букву Си?

Python 3.7.5 (default, Nov 13 2019, 14:05:23)
[Clang 11.0.0 (clang-1100.0.33.12)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import decimal
>>> help(decimal)

Help on module decimal:

NAME
    decimal - C decimal arithmetic module
...

Поэтому первый важный вывод:

Хоть Decimal и написан на Си, он в разы медленнее, чем float, так как реализован программно, а float – аппаратно.

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

Теперь самое главное – основание 10. Оно позволяет записывать десятичные дроби точно, без ошибок представления.

Decimal_Число = ±мантисса * 10 экcпонента

Мантисса и экспоненты – целые числа.

Помните, что мы не могли представить 0.1 в float с основанием 2? С основанием 10 – это элементарно:

0.1 = 1 * 10-1, и таким образом, 0.3 = 3 * 10-1 = (1 + 1 + 1) * 10-1 = 0.1 + 0.1 + 0.1

Как мы в школе учили десятичные дроби и знаем, как оперировать ими, так и здесь. Все точно и привычно. Ну почти все. Если мы разделим единицу на тройку, то получим бесконечную периодическую дробь 0.33333333…, либо по другому ее пишут 0.(3) – три в периоде. Естественно, что бесконечных чисел, записанных цифрами в памяти компьютера быть не может, иначе бы потребовалась бесконечная память. Поэтому количество троек в записи может быть большим, но обязано быть конечным.

Decimal оперирует с числами с произвольной (задаваемой пользователем), но конечной точностью.

По умолчанию точность – 28 десятичных знаков.

Еще одно следствие того, что Decimal реализовано программно – то, что его можно на ходу настраивать, как угодно пользователю. Для этого есть контекст – объект, содержащий настройки для выполнения операций и флаги. Операции выполняемые в этом контексте следуют правилам, заданным в нем. В отличии от float, где все правила фиксированы на аппаратном или низшим программным уровнях. Настроить можно:

  • Точность выполнения операций в количестве десятичных знаках
  • Режимы округления (их целых 8 штук)
  • Пределы по экспоненте
  • Режимы обработки исключительных ситуаций – настройки сигналов (например, деление на ноль, переполнение и прочее).

Флаги в контексте устанавливаются со стороны модуля decimal, если при последнем вычислении случился какой-то из сигналов. (Это отдельная тема, о ней потом.)

Сам же объект Decimal содержит знак, мантиссу (коэффициент перед экспонентой) и саму экспоненту (степень). Лишние нули в мантиссе на обрезаются, чтобы сохранять значимость числа (1.20 * 2.40 = 2.8800).

Decimal – иммутабельный (неизменяемый) тип. Операции над ним приводят к созданию новых объектов, а старые не меняются.

Поработаем с Decimal

Начинаем с импорта и посмотрим, каков контекст по умолчанию:

>>> from decimal import *
>>> getcontext()
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999, Emax=999999, capitals=1, clamp=0, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])

Мы видим здесь, что точность 28 знаков, округление к ближайшему четному, пределы по экспоненте ±999999, capitals – это про заглавную Е при печати, clamp не будем трогать пока что, флаги все сброшены, а включенные ловушки – неправильная операция, деление на ноль, переполнение. Если ловушка включена, это значит, что при возникновении соответствующего сигнала будет брошено исключение. Если нет ловушки, то при сигнале будет только втихую установлен флаг. Я оставлю тему ловушек на следующую статью.

Создание Decimal

Создать Decimal можно из обычного целого числа, из float, из строки или кортежа. С обычным числом все просто – int представлены и так точно:

>>> Decimal(1)
Decimal('1')
>>> Decimal(-1)
Decimal('-1')
>>> Decimal(10002332)
Decimal('10002332')

Из float – надо быть очень аккуратным. Потому что, float округляется внутри до ближайшего возможного, а Decimal не знает о ваших первоначальных намерениях, поэтому копирует содержимое float. К примеру, числа 0.1 в представлении float просто не существует. Python считывает 0.1 из кода как строку, потому ищет наиболее близкий к нему возможный float, а из него уже копируется содержимое в Decimal, как есть – уже с ошибкой:

>>> Decimal(0.1)
Decimal('0.1000000000000000055511151231257827021181583404541015625')

Не рекомендуется создавать Decimal из float. В Decimal попадет уже неправильно округленное число. Создавайте Decimal из целых чисел, либо из строк!

Логически правильно создавать Decimal сразу из строки, избегая фазу с превращением его в float! Что есть в строке – попадет в Decimal. Может показаться, что это немного криво – хранить числах в строках, но теперь вы знаете о представлении двоичного float, и строки обретают реальный смысл.

>>> Decimal('0.1')
Decimal('0.1')
>>> Decimal('3.14')
Decimal('3.14')
>>> Decimal('1.2e+10')
Decimal('1.2E+10')
>>> Decimal('10_000_000_000')  # c версии Python 3.6 можно подчеркивания
Decimal('10000000000')

Можно строкой еще задавать бесконечности и NaN (не число). Примеры:

>>> Decimal('Inf')
Decimal('Infinity')
>>> Decimal('-Inf')
Decimal('-Infinity')
>>> Decimal('nan')
Decimal('NaN')

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

  1. Знак, как число: 0 – это плюс, 1 – это минус.
  2. Кортеж из значащих цифр мантиссы
  3. Число – показатель экспоненты

Вообще кортеж для Decimal использует редко. Но вот вам пример:

>>> Decimal((0, (1, 2, 3, 4, 5), -1))
Decimal('1234.5')
>>> Decimal((1, (7, 7, 7), 3))
Decimal('-7.77E+5')

Если число слишком большое, то будет сигнал – неправильная операция. А так как на этом сигнале ловушка по умолчание – то будет исключение:

>>> Decimal("1e9999999999999999999")
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
decimal.InvalidOperation: [<class 'decimal.InvalidOperation'>]

Точность представление Decimal задается исключительно длиной задающего числа (или длиной строки). Настройки точности и режимов округления из контекста в ступают в игру только во время совершения математических операций.

>>> с = Context(prec=3)  # точность 3
>>> Decimal('5.643434231', c)  # но число целиком сохраняется
Decimal('5.643434231')

>>> Decimal('5.643434231', c) * 2  # после операции уже применяется округление до нужной точности
Decimal('11.287')

>>> +Decimal('5.643434231', c)  # трюк: унарный плюс применит контекст
Decimal('5.6434')

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

Работайте с Decimal как с обычными числами: складывайте, вычитайте, умножайте, делите и прочее. Можете, миксовать их с целыми числами. Но не рекомендуется миксовать их с float.

>>> x = Decimal('1.2')
>>> y = Decimal('2.3')
>>> x + y
Decimal('3.5')
>>> x - y
Decimal('-1.1')
>>> x * y
Decimal('2.76')
>>> x / y
Decimal('0.52174')

>>> y // x  # деление нацело
Decimal('1')
>>> y % x  # остаток
Decimal('1.1')

>>> Decimal('2.2') * 2
Decimal('4.4')
>>> Decimal('2.2') - 1
Decimal('1.2')

Дополнительно еще доступны некоторые математические функции:

>>> getcontext().prec = 10  # просто точность задали

>>> Decimal(2).sqrt()  # корень квадратный
Decimal('1.414213562')

>>> Decimal(2).ln()  # логарифм натуральный
Decimal('0.6931471806')

>>> Decimal(100).log10()  # логарифм десятичный
Decimal('2')

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

Кстати, Decimal можно передавать как аргументы функций, ожидающих float. Тогда они будут преобразованы во float:

>>> math.sin(Decimal(1))
0.8414709848078965

Метод quantize округляет число до фиксированной экспоненты, полезно для финансовых операций, когда нужно округлить копейки (центы). Первый аргумент – Decimal – что-то вроде шаблона округления. Смотрите примеры:

>>> Decimal('10.4266').quantize(Decimal('.01'), rounding=ROUND_DOWN)
Decimal('10.42')
>>> Decimal('10.4266').quantize(Decimal('.01'), rounding=ROUND_UP)
Decimal('10.43')
>>> Decimal('10.4266').quantize(Decimal('1.'), rounding=ROUND_UP)
Decimal('11')

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

>>> x = Decimal('0.1')
>>> x + x + x == Decimal('0.3')
True

Можно сортировать списки Decimal, искать минимум и максимум. А также преобразовывать Decimal обратно в обычные типы int, float, str. Пример из документации:

>>> data = list(map(Decimal, '1.34 1.87 3.45 2.35 1.00 0.03 9.25'.split()))
>>> max(data)
Decimal('9.25')
>>> min(data)
Decimal('0.03')
>>> sorted(data)
[Decimal('0.03'), Decimal('1.00'), Decimal('1.34'), Decimal('1.87'),
 Decimal('2.35'), Decimal('3.45'), Decimal('9.25')]
>>> sum(data)
Decimal('19.29')
>>> a, b, c = data[:3]
>>> str(a)
'1.34'
>>> float(a)
1.34
>>> round(a, 1)
Decimal('1.3')
>>> int(a)
1

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

Не все операции над Decimal абсолютно точные, если результат неточен, то возникает сигнал Inexact.

>>> c = getcontext()
>>> c.clear_flags()

>>> Decimal(1) / Decimal(3)
Decimal('0.3333333333')

>>> c.flags[Inexact]
True

Выводы

Выбор между Decimal и float – это поиск компромисса с учетом условий конкретной задачи.

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

Сравнительная таблица

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

Еще

В этой статье я не осветил полностью вопросы:

  • Сигналы, флаги и ловушки
  • Обзор режимов округления
  • Управление контекстами
  • Контексты и многопоточность

Если сообществу будет интересно, то я продолжу тему. Голосование будет на канале!

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

LRU-кэш в одну строчку

Кэш нужен, чтобы запоминать результаты каких-то тяжелых операций: вычислений, доступа к диску или запросов в сеть. В Python есть отличный декоратор, чтобы элегантно снабдить вашу функцию кэшированием: @functools.lru_cache(maxsize=128, typed=False)