Деплой

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

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

• Разделение ответственности в ряде случаев ответственность за деплой оказывается высокой.

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

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

Как это происходит у нас

В Учи.ру был создан собственный инструмент для выкатки сервисов и настройки стейджей под названием Шаман по сути кастомный интерфейс к HashiCorp Nomad. Несколько лет назад он помог уйти от исполнения скриптов на локальных машинах, автоматизировать процесс деплоя и позволил разработчикам создавать тестовые конфигурации максимально быстро. Об этом мои коллеги из команды инфраструктуры написали отдельную статью.

Немного истории

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

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

Да придет спаситель

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

А что делать, если тебе нужно поднять копию среды для тестирования второго релиза? Не создавать же всю связку целиком заново? А если ее незначительно видоизменить? Для подобного в Шамане есть возможность создать шаблон конфигурации, в котором с помощью YAML описана вся россыпь машин, сервисов и эндпоинтов. В результате, если ваша команда поддерживает свой шаблон, поднятие нового инстанса развесистой связки сервисов становится тривиальной задачей.

И, кажется, коллеги со мной согласны:

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

Схема связи Шамана с кодом

На данный момент почти все наши GitHub-репозитории собираются при помощи GitHub Actions, в результате чего докер-образ пушится в общую приватную докер-репу. Этот workflow можно триггерить как автоматически (пуш в релизную ветку), так и руками, если хочется потестить какой-то сиюминутный фикс. Дальше Шаман подцепляет свежий образ и по кнопке раскатывает его на стейдж. Ну не красота, а?

И так как этот процесс находится в ведении разработчиков чуть более чем полностью, у нас есть возможность его упрощать. Меня, например, все-таки расстраивает необходимость сначала идти в GitHub, чтобы посмотреть статус билда, а потом в Шаман чтобы нажать на заветную зеленую кнопку для выкатки образа. После незначительной встряски коллег из инфраструктуры выяснилось, что последний предоставляет API, ручку которого можно дернуть из Github Actions с адресом для деплоя и идентификатором образа для деплоя. А это значит, что деплоить код можно полностью автоматически!

Когда что то идёт не так

Понятное дело, что после выкатки кода на прод надо убедится, что он ведёт себя так, как задумано. И тут у нас чуть хуже с автоматизацией, но тем не менее всё под контролем. Для мониторинга мы используем сервис Sentry, который нотифицирует как о новых индивидуальных ошибках, так и об изменениях общего фона. Ну и продуктовые дашборды с метриками никто не отменял. А Шаман, в свою очередь, быстро и непринужденно позволяет откатить бажную версию фронтового приложения в случае серьезных проблем.

Так нужно ли деплоить разработчику?

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

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

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

• Сервисы приходят и уходят, соответственно, шаблон связки сервисов нужно держать актуальным, иначе поднятие новых тестовых стендов превратится в боль.

• Стоит держать в голове, что мы работаем с in-house-решением со всеми вытекающими: оно может сломаться. Поэтому здесь нужно быть достаточно открытым и помогать коллегам чинить проблему, а не бранить решение.

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

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

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

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

Несколько лет назад я работал над системой, которая позволила нам ускорить процесс деплоя конфигов на 1000+ серверов с минуты до нескольких секунд.

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

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

Случай из жизни

За прошедшие годы в Badoo было написано более ста различных сервисов, и их количество продолжает расти. При этом у каждого сервиса может быть от двух-трёх до нескольких сотен инстансов, а значит, в нашем коде должна быть возможность оперативно убирать запросы с определённого инстанса. Ведь как гласит закон Мёрфи, если что-то может пойти не так, оно пойдёт не так.

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

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

Для отключения сервисов мы используем свою систему с говорящим названием Выключалка хостов (disable hosts). Принцип её работы довольно прост:

• выбираем в веб-интерфейсе сервисы, которые нужно отключить (или, наоборот, включить);

• нажимаем на кнопку Deploy;

• изменения сохраняются в базе данных и затем доставляются на все машины, на которых выполняется PHP-код.

В коде при подключении к сервису стоит проверка вида:

if (\DownChecker\Host::isDisabled($host)) {   $this->errcode = self::ERROR_CONNECT_FAILED;   return false;}

В первой версии Выключалки хостов для доставки изменений на серверы мы решили использовать нашу основную систему деплоя конфигов под названием mcode. Она по SSH выполняет определённый набор команд с дополнительной обвязкой этого процесса.

Процесс деплоя разбит на несколько шагов:

• упаковываем конфиги в tar-архив;

• копируем архив на сервер через rsync или scp;

• распаковываем архив в отдельную директорию;

• переключаем симлинк на новую директорию.

Особенность mcode в том, что она ждёт завершения выполнения текущего шага на каждом сервере, прежде чем перейти к следующему. С одной стороны, это даёт больше контроля на каждом этапе и гарантирует, что на 99% серверов изменения доедут примерно в один момент времени (псевдоатомарность). С другой проблемы с одним из серверов приводят к тому, что процесс деплоя зависает, ожидая завершения выполнения текущей команды по тайм-ауту. Если на каком-то сервере несколько раз подряд не удалось выполнить команду, то он временно исключается из раскладки. Это позволяет уменьшить влияние проблемных серверов на общую продолжительность деплоя.

Таким образом, mcode не даёт никаких гарантий высокой скорости доставки изменений. Можно лишь спрогнозировать её максимальную продолжительность, как сумму таймаутов всех шагов.

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

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

В поисках альтернативного транспорта

Тут может возникнуть резонный вопрос: зачем изобретать велосипед, если можно использовать классическую схему с базой данных (БД) и кешированием?

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

• использование БД и кеша это дополнительная завязка на внешний сервис, а значит, ещё одна потенциальная точка отказа;

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

• ещё одним плюсом чтения из файла является то, что весь конфиг по сути является статическим массивом, а значит, после первого исполнения он осядет в OPCache, что даст небольшой прирост производительности.

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

Но в этой схеме нам не понравилась идея с запросами в цикле. У нас около 2000 серверов, на которых может выполняться PHP-код. Если мы будем делать запрос в базу/кеш хотя бы один раз в секунду, то это будет создавать довольно большой фон запросов (2k rps), а сами данные при этом обновляются не так часто. На этот случай есть решение событийная модель, например шаблон проектирования Publisher-Subscriber (PubSub). Из популярных решений можно было использовать Redis, но у нас он не прижился (для кеширования мы используем Memcache, а для очередей у нас есть своя отдельная система).

Зато прижился Consul, у которого есть механизм отслеживания изменений (watches) на базе блокирующих запросов. Это в целом похоже на PubSub и вписывается в нашу схему с обновлением конфига по событию. Мы решили сделать прототип нового транспорта на базе Consul, который со временем эволюционировал в отдельную систему под названием AutoConfig.

Как работает AutoConfig

Общая схема работы выглядит следующим образом:

• разработчик вызывает специальный метод, который обновляет значение ключа в базе данных и добавляет ключ в очередь на деплой;

• в облаке работает скрипт, который исключает пачку ключей из очереди и записывает их в Consul; помимо самих ключей, в Consul хранится специальная индексная карта со списком всех ключей (о ней я расскажу чуть позже);

• на каждом сервере запущен Consul watch, который отслеживает изменения индексной карты и вызывает специальный скрипт-обработчик, когда она обновляется;

• обработчик записывает изменения в файл.

Со стороны кода работа с системой выглядит довольно просто:

Обновление ключа

$record = \AutoConfig\AutoConfigRecord::initByKeyData('myKey', 'Hello, Habr!', 'Eugene Tupikov');$storage = new \AutoConfig\AutoConfigStorage();$storage->deployRecord($record);

Чтение ключа

$reader = new \AutoConfig\AutoConfigReader();$config = $reader->getFromCurrentSpace('myKey');

Удаление ключа

$storage = new \AutoConfig\AutoConfigStorage();$storage->removeKey('example');

Подробнее про Consul watch

Consul watch реализован с использованием блокирующих запросов в HTTP API, работу которых лучше продемонстрировать на примере отслеживания изменений ключа.

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

curl -X PUT --data 'hello, harb!' http://127.0.0.1:8500/v1/kv/habr-key

Затем отправляем запрос на чтение нашего ключа

curl -v http://127.0.0.1:8500/v1/kv/habr-key

API обрабатывает запрос и возвращает ответ, содержащий HTTP-заголовок X-Consul-Index с уникальным идентификатором, который соответствует текущему состоянию нашего ключа.

......< X-Consul-Index: 266834870< X-Consul-Knownleader: true......<[  {    "LockIndex": 0,    "Key": "habr-key",    "Flags": 0,    "Value": "dXBkYXRlZCAy",    "CreateIndex": 266833109,    "ModifyIndex": 266834870  }]

Мы отправляем новый запрос на чтение и дополнительно передаем значение из заголовка X-Consul-Index в параметре запроса index

curl http://127.0.0.1:8500/v1/kv/habr-key?index=266834870Ждем изменений ключа...

При виде параметра index API начинает ждать изменений нашего ключа и ждёт до тех пор, пока время ожидания не превысит заданного тайм-аута.

Затем открываем новую вкладку терминала и отправляем запрос на обновление ключа

curl -X PUT --data 'updated value' http://127.0.0.1:8500/v1/kv/habr-key

Возвращаемся на первую вкладку и видим, что запрос на чтение вернул обновленное значение (изменились ключи Value и ModifyIndex):

[  {    "LockIndex": 0,    "Key": "habr-key",    "Flags": 0,    "Value": "dXBkYXRlZA==",    "CreateIndex": 266833109,    "ModifyIndex": 266835734  }]

При вызове команды

consul watch -type=key -key=habr_key <handler>

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

Зачем нужна индексная карта

Помимо отслеживания изменений конкретного ключа, можно подписаться на изменения ключей с определённым префиксом. Предположим, что все наши ключи имеют общий префикс auto_config. Если мы запустим команду

consul watch -type=keyprefix -prefix=auto_config/ <handler>

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

По этому поводу на GitHub уже довольно давно открыт Issue и, судя по комментариям, лёд тронулся. Разработчики Consul начали работу над улучшением подсистемы блокирующих запросов, что должно решить описанную выше проблему.

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

Она имеет следующий формат:

return [    'value' => [        'version' => 437036,        'keys' => [            'my/awesome/key' => '80003ff43027c2cc5862385fdf608a45',            ...            ...        ],        'created_at' => 1612687434    ]]

В случае если карта обновилась, обработчик:

• считывает текущее состояние карты с диска;

• находит изменившиеся ключи (для этого и нужен хеш значения);

• вычитывает через HTTP API актуальные значения изменившихся ключей и обновляет нужные файлы на диске;

• сохраняет новую индексную карту на диск.

И ещё немного про Consul

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

Ограничение размера значения ключа (и не только)

Consul это распределённая система, и для достижения согласованности используется протокол Raft. Для стабильной работы протокола в Consul установлен максимальный размер значения одного ключа 512 Кб. Его можно изменить, воспользовавшись специальной опцией, но делать это крайне не рекомендуется, так как изменение может привести к непредсказуемому поведению всей системы.

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

Чтобы обойти эти ограничения, мы сделали следующее:

• разбили индексную карту на несколько частей (по 1000 ключей в каждой) и обновляем только части, содержащие изменённые ключи;

• ограничили максимальный размер одного ключа AutoConfig 450 Кб, чтобы оставить место для шардов индексной карты (значение выбрано опытным путём);

• доработали скрипт, обрабатывающий очередь на деплой таким образом, что он

  • сначала вычитывает N ключей из очереди и проверяет их суммарный размер;

  • если размер не превышает заданный предел, то деплоит все ключи разом, а в противном случае деплоит ключи по одному.

Отсутствие встроенной репликации

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

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

В качестве заключения

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

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

За последние пару лет система стала популярна у наших разработчиков и сегодня фактически является стандартом при работе с конфигами, для которых не требуется атомарная раскладка. Например, через AutoConfig у нас деплоятся параметры A/B-тестов, настройки промокампаний, на его базе реализован функционал Service Discovery и многое другое.

Общее количество ключей на данный момент порядка 16 000, а их суммарный размер примерно 120 Мб.

Спасибо за внимание!

Расскажите в комментариях, как вы деплоите конфиги в своих проектах.

Ссылка на статью в моем блоге

Тропические леса и фикусы-душители

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

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

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

Рефакторинг сервиса приложения доставки продуктов

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

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

Допустим имеются следующие действия, которые у нас хранятся в одной таблице:

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

• Можно оформитьвозврат товара. Если вам не понравился кефир - вы оформляете возврат и вам возвращают его цену.

• Можносписать бонусысо счета. В таком случае часть стоимости оплачивается этими бонусами.

• Начисляются бонусы. Каким-либо алгоритмом нам не важно каким конкретно.

• Также заказ может бытьзарегистрирован в некотором приложении-партнере(ExternalOrder)

Все перечисленная информация по заказам и пользователям хранится в таблице (пусть она будет называтьсяOrderHistory):

Логика такой организации данных вполне справедлива на раннем этапе разработки системы. Ведь наверняка пользователь может посмотреть историю своих действий. Где он одним списком видит что он заказывал, как начислялись и списывались бонусы. В таком случае мы просто выводим записи, относящиеся к нему за указанный диапазон. Организовать в виде одной таблицы банальная экономия на создании дополнительных таблиц, их поддержании. Однако, по мере роста бизнес-логики и добавления новых типов операций число столбцов с null значениями начало расти. Записей в таблице сотни миллионов. Причем распределены они очень неравномерно. В основном это операции открытия и закрытия заказов. Но вот операции начисления бонусов составляют 0.1% от общего числа, однако эти записи используются при расчете новых бонусов, что происходит регулярно.В итоге логика расчета бонусов работает медленнее, чем если бы эти записи хранились в отдельной таблице. Ну и расширять таблицу новыми столбцами не хотелось бы в дальнейшем. Кроме того заказы в закрытом статусе с датой создания более 2 месяцев для бизнес-логики интереса не представляют. Они нужны только для отчетов не более.

И вот возникает идея.Разделить таблицу на две, три или даже больше.

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

Изменение структуры хранения в три этапа

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

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

Оба экземпляра работают с одной базой данных. Реализуя паттернShared Database.

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

Отдельная новая база данных вполне может появиться. Однако не всегда. Ввиду сложностей обеспечения транзакционности между двумя БД. Все в конечном счете зависит от реализации и от ограничений бизнес-логики.

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

BonusOperations:

Отдельную таблицу для данных из внешних систем -ExternalOrders:

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

Для оставшихся типов записей -OrderHistoryArchive(старше 2х недель). Где теперь также можно удалить несколько лишних столбцов.

Выделение таких архивных данных часто бывает удобным. Если оперативная часть очень требовательна к производительности она вполне может себе размещается на быстрых SSD дисках. В то время как архивные данные могут использоваться один раз в месяц для отчета. И их больше в разы. Поэтому размещая их на дешевых дисках мы экономим иногда вполне приличную сумму.

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

Монолит версии 1 и монолит версии 2 вполне могут работать совместно. Лишь для тех запросов, которые обрабатывались монолитом версии 1 в новой базе данных будут пробелы. Эти пробелы, а также недостающие данные можно будет в дальнейшем скопировать отдельным скриптом или утилитой.

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

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

Итого. Внешне система никогда не менялась. Однако внутренняя организация радикально преобразилась. Возможно под капотом теперь работает новая система. Которая лишена недостатков предыдущей. Не напоминает фикусов-душителей? Что-то похожее есть. Поэтому именно такое название паттерн и получил Strangler.

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

Выводы

• ПаттернStranglerпозволяет совершенствовать системы с высокими требованиями к SLA.

• Для обновления системы без простоя нужно сделать как минимум 3 последовательных развертования на продакшен. Это одна из причин, почему системы требовательные к показателям общего простоя заметно дороже.

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

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

Всем привет! Меня зовут Борис, я выпускник программы Науки о данных ФКН ВШЭ, работаю ML Инженером и преподаю в OTUS на курсах ML Professional, DL Basic, Computer Vision.

В первых числах января 2021 я узнал про британский штамм коронавируса,прогнозы о новой волне в США. Я подумал: аналитик данных я или кто? Мне захотелось забить гвоздик своим микроскопом и узнать, вызовет ли британский штамм волну заражений в Москвеи стоит ли покупать авиабилеты на лето.

Выглядело как приключение на две недели, но превратилось в исследование на три месяца. В процессе я выяснил, что хороших материалов по созданию эпидемиологических моделей практически нет. Банально авторы статей по моделированию COVID-19 в топовых журналах даже не делают train-test split.

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

Ноутбук к туториалу.

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

Препринт статьи: TBA

Какой британский штамм?

COVID-19 в представлении не нуждается, к нему мы уже привыкли. Однако новые мутации (штаммы) вируса способны изменить ход эпидемии. Сейчас ВОЗвыделяеттри таких штамма (variants of concern). Среди них штамм B.1.1.7, так называемый британский штамм.

Почему меня заинтересовал именно он? Во-первых,исследования показывают, что этот штамм распространяется на 40% - 90% быстрее, чем обычный коронавирус (для продвинутых: R0 на 40% - 90% больше). Во-вторых, случай заражения этим штаммом уже былобнаруженв России. В-третьих, непохоже, чтобы кто-то в России готовился к его появлению.

По результатам исследования пришел к выводу, что B.1.1.7 действительно способен привести к новой волне заражений COVID-19 и может унести 30 тысяч жизней.

Исходные данные

Нам понадобится статистика по заражениям, смертям и выздоровлениям в Москве. Её можно скачать со страницы.
На момент исследования был доступен промежуток дат: 2020.03.10 - 2021.03.23.

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

df.columns = ['date', 'region',             'total_infected', 'total_recovered', 'total_dead',            'deaths_per_day', 'infected_per_day', 'recovered_per_day']df = df[df.region == 'Москва'].reset_index()df['date'] = pd.to_datetime(df['date'], format='%d.%m.%Y')df['infected'] = df['total_infected'] - df['total_recovered'] - df['total_dead']df = df.drop(columns=['index', 'region'])df = df.sort_values(by='date')df.index = df.datedf['date'] = df.date.asfreq('D')

Добавим сглаженные скользящим окном в 7 дней версии всех колонок. Они пригодятся нам позже.

df_smoothed = df.rolling(7).mean().round(5)df_smoothed.columns = [col + '_ma7' for col in df_smoothed.columns]full_df = pd.concat([df, df_smoothed], axis=1)for column in full_df.columns:    if column.endswith('_ma7'):        original_column = column.strip('_ma7')        full_df[column] = full_df[column].fillna(full_df[original_column])

В итоге работаем со следующими колонками, а так же их сглаженными версиями:

• infected_per_day новых заражений в данный день.

• recovered_per_day новых выздоровлений в данный день.

• deaths_per_day погибших от инфекции в данный день.

• total_infected всего заражений к этому моменту, кумулятивная суммаinfected_per_day.

• total_dead всего погибших к этому моменту, кумулятивная сумма отdeaths_per_day.

• total_recovered всего выздоровлений к этому моменту, кумулятивная сумма отrecovered_per_day.

Основа модели: SEIRD

SIRэто очень простая модель, которая симулирует развитие эпидемии во времени. Популяция разделяется на три группы: Susceptible, Infected, Recovered.
Идея на пальцах:

• Есть замкнутая популяция и изначальное количество зараженных (Infected).

• В течение дня каждый зараженный с некоторой вероятностью инфицирует кого-то из группы Susceptible. Новые зараженные переходят в группу Infected.

• Находящиеся в группе Infected выздоравливают, когда проходит срок длительности болезни и переходят в группу Recovered.

• Запускаем этот процесс на много дней вперед и смотрим, как меняется численность групп.

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

• Susceptible могут быть заражены.

• Exposed заражены вирусом, но не распространяют его, то есть проходят инкубационный период.

• Infectious распространяют вирус.

• Recovered переболели и не могут быть заражены.

• Dead погибли.

Изменение численности групп за день определяют дифференциальные уравнения:

Например, третье уравнение описывает изменение группы Infected за день. С вероятностьюу человека из группыEзаканчивается инкубационный период и он переходит в группуI. Так же с вероятностьючеловек из группыIвыздоравливает или погибает, переходя в группуRилиD. Куда он попадет определяется параметром, который соответствует доле смертности (Infection Fatality Rate). Таким образом, каждый день группу пополняют новые зараженные, а переболевшие уходят из группы.

Параметры модели:
case fatality rate.
количество человек, которых один переносчик заражает за день.
1 делить на среднюю длину инкубационного периода.
y 1 делить на среднюю длительность болезни.
R0 = /y базовое репродуктивное число, то есть количество человек, которых один переносчик заражает за время своей болезни.

Есть три причины использовать именно SEIRD. Во-первых, параметры SEIRD это просто характеристики болезни, а значит мы сможем прикинуть их значения на основании медицинских исследований COVID-19 как болезни. Во-вторых, эту модель очень легко модифицировать, что мы и будем делать дальше, добавляя меры карантина, неполноту статистики и второй штамм. В-третьих, эта модель позволяет легко покрутить разные сценарии. Например, можно смоделировать ситуацию: Что будет, если в этой же популяции появится второй штамм и параметрR0у него на 90% больше? Сделать подобное с помощью, например, градиентного бустинга, гораздо сложнее.

Главный минус SEIRD: это очень простая и полностью детерминированная модель. Она чувствительна к изначальным значениям и параметрам. Предсказания такой модели могут на порядки отличаться от реальных значений. Для нас это не страшно: мы ведь не пытаемся сделать самый точный прогноз по количеству зараженных, мы просто пытаемся понятьможет ли штамм B.1.1.7 вызвать новую волну. Интересно, что, несмотря на это, лучшая на текущий момент модель по прогнозированию COVID-19 этоочень навороченная SEIR модель.

Реализуем классический SEIRD

Привожу минимальную реализацию SEIRD.

class BarebonesSEIR:    def init(self, params=None):        self.params = params<span class="hljs-function" style="box-sizing: border-box;"><span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">def</span> <span class="hljs-title" style="box-sizing: border-box; color: rgb(121, 93, 163);">get_fit_params</span>(<span class="hljs-params" style="box-sizing: border-box;">self</span>):</span>    params = lmfit.Parameters()    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"population"</span>, value=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">12_000_000</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">False</span>)    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"epidemic_started_days_ago"</span>, value=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">10</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">False</span>)    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"r0"</span>, value=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">4</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">min</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">3</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">max</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">5</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">True</span>)    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"alpha"</span>, value=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0.0064</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">min</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0.005</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">max</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0.0078</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">True</span>)  <span class="hljs-comment" style="box-sizing: border-box; color: rgb(150, 152, 150);"># CFR</span>    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"delta"</span>, value=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>/<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">3</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">min</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>/<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">14</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">max</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>/<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">2</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">True</span>)  <span class="hljs-comment" style="box-sizing: border-box; color: rgb(150, 152, 150);"># E -&gt; I rate</span>    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"gamma"</span>, value=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>/<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">9</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">min</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>/<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">14</span>, <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">max</span>=<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>/<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">7</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">True</span>)  <span class="hljs-comment" style="box-sizing: border-box; color: rgb(150, 152, 150);"># I -&gt; R rate</span>    params.add(<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">"rho"</span>, expr=<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'gamma'</span>, vary=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">False</span>)  <span class="hljs-comment" style="box-sizing: border-box; color: rgb(150, 152, 150);"># I -&gt; D rate</span>    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> params<span class="hljs-function" style="box-sizing: border-box;"><span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">def</span> <span class="hljs-title" style="box-sizing: border-box; color: rgb(121, 93, 163);">get_initial_conditions</span>(<span class="hljs-params" style="box-sizing: border-box;">self, data</span>):</span>    <span class="hljs-comment" style="box-sizing: border-box; color: rgb(150, 152, 150);"># Simulate such initial params as to obtain as many deaths as in data</span>    population = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'population'</span>]    epidemic_started_days_ago = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'epidemic_started_days_ago'</span>]    t = np.arange(epidemic_started_days_ago)    (S, E, I, R, D) = self.predict(t, (population - <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>, <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0</span>, <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>, <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0</span>, <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0</span>))    I0 = I[-<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>]    E0 = E[-<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>]    Rec0 = R[-<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>]    D0 = D[-<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>]    S0 = S[-<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>]    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> (S0, E0, I0, Rec0, D0)<span class="hljs-function" style="box-sizing: border-box;"><span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">def</span> <span class="hljs-title" style="box-sizing: border-box; color: rgb(121, 93, 163);">step</span>(<span class="hljs-params" style="box-sizing: border-box;">self, initial_conditions, t</span>):</span>    population = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'population'</span>]    delta = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'delta'</span>]    gamma = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'gamma'</span>]    alpha = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'alpha'</span>]    rho = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'rho'</span>]        rt = self.params[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'r0'</span>].value    beta = rt * gamma    S, E, I, R, D = initial_conditions    new_exposed = beta * I * (S / population)    new_infected = delta * E    new_dead = alpha * rho * I    new_recovered = gamma * (<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span> - alpha) * I    dSdt = -new_exposed    dEdt = new_exposed - new_infected    dIdt = new_infected - new_recovered - new_dead    dRdt = new_recovered    dDdt = new_dead    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">assert</span> S + E + I + R + D - population &lt;= <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1e10</span>    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">assert</span> dSdt + dIdt + dEdt + dRdt + dDdt &lt;= <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1e10</span>    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> dSdt, dEdt, dIdt, dRdt, dDdt<span class="hljs-function" style="box-sizing: border-box;"><span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">def</span> <span class="hljs-title" style="box-sizing: border-box; color: rgb(121, 93, 163);">predict</span>(<span class="hljs-params" style="box-sizing: border-box;">self, t_range, initial_conditions</span>):</span>    ret = odeint(self.step, initial_conditions, t_range)    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> ret.T

Давайте разберемся.

Методget_fit_paramsпросто создает хранилище парамтеров. Мы используем библиотекуlmfit, но пока что не оптимизируем параметры, а просто используем структуруParametersкак продвинутый словарь. Для задания диапазонов параметров я использовалмета-анализы характеристик COVID-19.

Параметрepidemic_started_days_agoпозволяет задать дату появления первого зараженного. В моём предположении это2 марта 2020.

Методstepсодержит дифференциальные уравнения модели. Он получает на вход изначальные численности всех групп в виде кортежаinitial_conditionsи номер дняt, возвращает изменения всех групп за день.

Методpredictполучает на вход начальные условия и массивt_range. Он просто запускаетscipy.integrate.odeintна методеstep, чтобы проинтегрировать дифференциальные уравнения и получить численность групп в каждый день из массиваt_range.

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

Запустим модель и сравним её предсказания с реальностью:

model = BarebonesSEIR()model.params = model.get_fit_params()train_initial_conditions = model.get_initial_conditions(train_subset)train_t = np.arange(len(train_subset))(S, E, I, R, D) = model.predict(train_t, train_initial_conditions)plt.figure(figsize=(10, 7))plt.plot(train_subset.date, train_subset['total_dead'], label='ground truth')plt.plot(train_subset.date, D, label='predicted', color='black', linestyle='dashed' )plt.legend()plt.title('Total deaths')plt.show()

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

Взглянем на динамику всех групп SEIRD:

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

Добавляем сдерживающие меры

Одна из причин появления нескольких волн это сдерживающие меры. Предположим, что карантинные меры в деньt, снижаютR0болезни на процентq(t). Тогда для каждого дня можно расчитать репродуктивное число с учетом карантина:Rt = R0 - R0 * q(t). Отсюда мы можем вычислить(t) = Rt * yи использовать её в уравнениях.

Осталось только задать функцииq(t). В своём исследовании я использовал простую ступенчатую функцию. Делим весь промежуток на отрезки по 60 дней, для каждого отрезка добавляем новый параметр: уровень карантина в данный период. Переходы между отрезками стоит сделать плавными с помощью, иначе функция не будет гладкой и подбирать параметры будет сложно.

Реализуем небольшой пример для наглядности.

def sigmoid(x, xmin, xmax, a, b, c, r):    x_scaled = (x - xmin) / (xmax - xmin)    out = (a * np.exp(c * r) + b * np.exp(r * x_scaled)) / (np.exp(c * r) + np.exp(x_scaled * r))    return outdef stepwise_soft(t, coefficients, r=20, c=0.5):    t_arr = np.array(list(coefficients.keys()))min_index = np.<span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">min</span>(t_arr)max_index = np.<span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">max</span>(t_arr)<span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">if</span> t &lt;= min_index:    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> coefficients[min_index]<span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">elif</span> t &gt;= max_index:    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> coefficients[max_index]<span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">else</span>:    index = np.<span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">min</span>(t_arr[t_arr &gt;= t])<span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">if</span> <span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">len</span>(t_arr[t_arr &lt; index]) == <span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">0</span>:    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> coefficients[index]prev_index = np.<span class="hljs-built_in" style="box-sizing: border-box; color: rgb(0, 92, 197);">max</span>(t_arr[t_arr &lt; index])<span class="hljs-comment" style="box-sizing: border-box; color: rgb(150, 152, 150);"># sigmoid smoothing</span>q0, q1 = coefficients[prev_index], coefficients[index]out = sigmoid(t, prev_index, index, q0, q1, c, r)<span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> outt_range = np.arange(100)coefficients = {    0: 0,    30: 0.5,    60: 1,    100: 0.4,}plt.title('Пример функции уровня карантина')plt.scatter(coefficients.keys(), coefficients.values(), label='Моменты изменения уровня карантина')plt.plot(t_range, [stepwise_soft(t, coefficients, r=20, c=0.5) for t in t_range], label='Ступенчатая функция с плавными переходами')plt.xlabel('t')plt.ylabel('Уровень канартина')plt.legend(loc='lower center', bbox_to_anchor=(0.5, -0.6),)plt.show()

Чтобы добавить это в нашу SEIRD модель нужно сделать несколько шагов:

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

def get_fit_params(self, data):    ...    params.add(f"t0_q", value=0, min=0, max=0.99, brute_step=0.1, vary=False)    piece_size = self.stepwise_size    for t in range(piece_size, len(data), piece_size):        params.add(f"t{t}_q", value=0.5, min=0, max=0.99, brute_step=0.1, vary=True)    return params

2. Добавляемновый методget_step_rt_beta, который вычисляет(t)иRtиспользуя ступенчатую функцию.

3. В методеstepиспользуемget_step_rt_beta,(t)с учетом карантина в деньt.

Учитываем неполноту статистики

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

Изменим группы Infected, Recovered и Dead следующим образом:

• Iv(t) распространяют вирус, видны в статистике.

• I(t) распространяют вирус, не видны в статистике.

• Rv(t) переболели, видны в статистике.

• R(t) переболели, не видны в статистике.

• Dv(t) погибли, видны в статистике.

• D(t) погибли, не видны в статистике.

Добавим два новых параметра:

• pi вероятность попадания зараженного в статистику, то есть в группу Iv.

• pd вероятность регистрации смерти зараженного, невидимого в статистике, в группу Dv. Зараженные из группы Iv всегда попадают в Dv, но этот параметр позволяет учесть, что часть смертей незамеченных зараженных так же попадает в статистику по смертности от COVID-19.

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

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

Оптимизируем параметры

Модель готова, осталось только подобрать оптимальные значения параметров так, чтобы выдаваемые моделью Iv(t), Rv(t) и Dv(t) были как можно ближе к историческим данным. Это можно сделать с помощью метода наименьших квадратов.

Конкретнее, используем для этогоlmfit.minimize. Эта функция получает на вход callable, возвращающий список отклонений предсказаний от истины (остатки, residuals), и подбирает такие параметры, чтобы сумма отклонений была минимальна. Используем алгоритм Levenberg-Marquardt (method=leastsq), который двигает каждый параметр на очень маленькое значение и проверяет, уменьшается ли ошибка. Важно знать, что по этой причине алгоритм не может оптимизировать дискретные параметры, поэтому, например, не сможет подобратьstepwise_size.

Я написалнебольшую оберткувокруг этой функции под названиемBaseFitter, которая сохраняет историю параметров и делает другие служебные вещи. Но в конце концов всё сводится к вызовуminimize.

Рассмотрим получение остатков подробнее:

def smape_resid_transform(true, pred, eps=1e-5):    return (true - pred) / (np.abs(true) + np.abs(pred) + eps)class HiddenCurveFitter(BaseFitter):...    def residual(self, params, t_vals, data, model):        model.params = params    initial_conditions = model.get_initial_conditions(data)    (S, E, I, Iv, R, Rv, D, Dv), history = model.predict(t_vals, initial_conditions, history=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">False</span>)    (new_exposed,     new_infected_invisible, new_infected_visible,     new_recovered_invisible,     new_recovered_visible,     new_dead_invisible, new_dead_visible) = model.compute_daily_values(S, E, I, Iv, R, Rv, D, Dv)    new_infected_visible = new_infected_visible    new_dead_visible = new_dead_visible    new_recovered_visible = new_recovered_visible    true_daily_cases = data[self.new_cases_col].values[<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>:]    true_daily_deaths = data[self.new_deaths_col].values[<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>:]    true_daily_recoveries = data[self.new_recoveries_col].values[<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>:]    resid_I_new = smape_resid_transform(true_daily_cases, new_infected_visible)    resid_D_new = smape_resid_transform(true_daily_deaths, new_dead_visible)    resid_R_new = smape_resid_transform(true_daily_recoveries, new_recovered_visible)    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">if</span> self.weights:        residuals = np.concatenate([            self.weights[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'I'</span>] * resid_I_new,            self.weights[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'D'</span>] * resid_D_new,            self.weights[<span class="hljs-string" style="box-sizing: border-box; color: rgb(223, 80, 0);">'R'</span>] * resid_R_new,        ]).flatten()    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">else</span>:        residuals = np.concatenate([            resid_I_new,            resid_D_new,            resid_R_new,        ]).flatten()    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> residuals

Основная идея простая. Для промежуткаt_valsполучаем предсказания заражений в день, смертей в день и выздоровлений в день. Вычисляем отклонения предсказанных значений от реальных и кладем всё это в один массив.

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

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

Во-вторых, можно предположить, что статистика по смертям достовернее статистики по заражениям и выздоровлениям. Чтобы внести это предположение в оптимизацию вводятся веса (self.weights) для остатков. Хорошо сработали такие веса:0.5для остатков по смертям и0.25для остальных.

Тренируем модель

Наконец, натренируем SEIRD со скрытыми состояниями. Статистика вносится с некоторым лагом по времени. Чтобы это меньше влияло на нашу модель будем использовать сглаженные версии колонок для расчета остатков.

from sir_models.fitters import HiddenCurveFitterfrom sir_models.models import SEIRHiddenstepwize_size = 60weights = {    'I': 0.25,    'R': 0.25,    'D': 0.5,}model = SEIRHidden(stepwise_size=stepwize_size)fitter = HiddenCurveFitter(     new_deaths_col='deaths_per_day_ma7',     new_cases_col='infected_per_day_ma7',     new_recoveries_col='recovered_per_day_ma7',     weights=weights,     max_iters=1000,)fitter.fit(model, train_subset)result = fitter.resultresult

Полученные параметры:

Проверим фит модели на тренировочных данных:

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

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

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

Верифицируем модель c помощью кросс-валидации

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

В статьях по моделированию COVID-19 стандартом является проверять модель на предсказаниях кумулятивного числа погибших, поэтому сделаем так же.

Процесс верификации:

• Выбираем несколько дат из тренировочных данных, например каждый 20-ый день.

• Для каждой даты:

  • Обучаем модель на всех данных до этой даты.

  • Делаем прогноз суммарного количества погибших на 30 дней вперед.

  • Считаем ошибку предсказания.

• Вычисляем среднюю ошибку модели.

• Сравниваем ошибку с ошибкой бейзлайна: предсказанием предыдущего дня.

from sir_models.utils import eval_on_select_dates_and_k_days_aheadfrom sir_models.utils import smapefrom sklearn.metrics import mean_absolute_errorK = 30last_day = train_subset.date.iloc[-1] - pd.to_timedelta(K, unit='D')eval_dates = pd.date_range(start='2020-06-01', end=last_day)[::20]def eval_hidden_moscow(train_df, t, train_t, eval_t):    weights = {        'I': 0.25,        'R': 0.25,        'D': 0.5,    }    model = SEIRHidden()    fitter = HiddenCurveFitter(        new_deaths_col='deaths_per_day_ma7',        new_cases_col='infected_per_day_ma7',        new_recoveries_col='recovered_per_day_ma7',        weights=weights,        max_iters=1000,        save_params_every=500)    fitter.fit(model, train_df)train_initial_conditions = model.get_initial_conditions(train_df)train_states, history = model.predict(train_t, train_initial_conditions, history=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">False</span>)test_initial_conds = [compartment[-<span class="hljs-number" style="box-sizing: border-box; color: rgb(0, 134, 179);">1</span>] <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">for</span> compartment <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">in</span> train_states]test_states, history = model.predict(eval_t, test_initial_conds, history=<span class="hljs-literal" style="box-sizing: border-box; color: rgb(0, 134, 179);">False</span>)    <span class="hljs-keyword" style="box-sizing: border-box; color: rgb(167, 29, 93);">return</span> model, fitter, test_states(models, fitters, model_predictions, train_dfs, test_dfs) = eval_on_select_dates_and_k_days_ahead(train_subset,                                        eval_func=eval_hidden_moscow,                                        eval_dates=eval_dates,                                        k=K)model_pred_D = [pred[7] for pred in model_predictions]true_D = [tdf.total_dead.values for tdf in test_dfs]baseline_pred_D = [[tdf.iloc[-1].total_dead]*K for tdf in train_dfs]overall_errors_model = [mean_absolute_error(true, pred) for true, pred in zip(true_D, model_pred_D)]overall_errors_baseline = [mean_absolute_error(true, pred) for true, pred in zip(true_D, baseline_pred_D)]print('Mean overall error baseline', np.mean(overall_errors_baseline).round(3))print('Mean overall error model', np.mean(overall_errors_model).round(3))overall_smape_model = [smape(true, pred) for true, pred in zip(true_D, model_pred_D)]np.median(overall_smape_model)

Результат:

• Mean Absolute Arror бейзлайна: 714.

• Mean Absolute Arror модели: 550.

• Symmetric Mean Absolute Percentage Error: 4.6%

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

Расширение модели для двух штаммов

Ситуацию с двумя штаммами можно смоделировать как две параллельные эпидемии, которые заражают людей из одной популяции. При этом зараженные одним штаммом не могут заразиться другим. В контексте SEIR моделей это означает, например, что вместо группыI(t)для всех зараженных у нас теперь будет две группы для разных штаммов:I1(t)иI2(t).

Модель эпидемии обычного SARS-CoV-2 мы только что обучили. Осталось получить модель для британского штамма. СогласноисследованиямB.1.1.7 отличается только параметромR0: он в 1.4 - 1.9 раз больше. Значит мы можем взять обученную модель для обычного SARS-CoV-2, увеличитьR0и получить модель для B.1.1.7.

Полная схема модели с двумя штаммами выглядит так:

Главное отличиереализации этой модели в кодеот предыдущих в том, что эта модель не обучается, а создается на основе обученной модели для одного штамма. Добавляется параметрbeta2_mult, который задает, на сколько надо умножить1(t)первого штамма, чтобы получить2(t)второго.

class SEIRHiddenTwoStrains(SEIRHidden):    ...    @classmethod    def from_strain_one_model(cls, model):        strain1_params = model.params        strain1_params.add("beta2_mult", value=1.5, min=1, max=2, vary=False)        return cls(params=deepcopy(strain1_params))

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

Сценарии развития эпидемии при появлении B.1.1.7

Наконец, мы получили модель для двух штаммов. Неизвестно какое именноR0у B.1.1.7, а так же сколько сейчас в Москве носителей этого штамма. Рассмотрим несколько сценариев.

ЗначениеR0влияет на дату и высоту пика заражений новым штаммом. В лучшем случае B.1.1.7 просто продлевает эпидемию в текущем состоянии на полгода-год. В худшем случае случается большая волна заражений.

Изначальное количество носителей британского штамма влияет только на дату пика. Во всех случаях новая волна случается в течение года после появления нового штамма.

Ослабление карантина приводит к значительному количеству заражений.

Автор статьи: Борис Цейтлин

Узнайте подробнее о курсе "Machine Learning. Professional".

Также приглашаем всех желающих на бесплатный интенсив по теме: Деплой ML модели: от грязного кода в ноутбуке к рабочему сервису.

Перевод подготовлен в рамках курса "Machine Learning. Professional".

Также приглашаем всех желающих принять участие в двухдневном онлайн-интенсиве Деплой ML модели: от грязного кода в ноутбуке к рабочему сервису.

Ошибки наблюдения и различия в подгруппах вызывают статистические парадоксы

Ошибки наблюдения и различия в подгруппах могут легко привести к статистическим парадоксам в любом прикладном решении data science. Игнорирование этих элементов может полностью дискредитировать заключения нашего анализа.

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

1. Парадокс Берксона

Первым ярчайшим примером является обратная корреляция между степенью тяжести заболевания COVID-19 и курением сигарет (см., например, обзор Европейской комиссии Wenzel 2020). Курение сигарет широко известный фактор риска респираторных заболеваний, так как же объяснить это противоречие?

Работа Griffith 2020, недавно опубликованная в Nature, предполагает, что это может быть случай ошибки коллайдера (Collider Bias), также называемой парадоксом Берксона. Чтобы понять этот парадокс, давайте рассмотрим следующую графическую модель, в которую мы включили третью случайную переменную: госпитализация.

Парадокс Берксона: госпитализация это переменная-коллайдер для курения сигарет, и для тяжести течения COVID-19. (Изображение автора)

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

Давайте рассмотрим следующий пример набора данных. На левом рисунке у нас есть данные по всему населению, а на правом рисунке мы рассматриваем только подмножество госпитализированных людей (то есть мы используем переменную-коллайдер).

Парадокс Берксона: если мы добавим условие в соответствии с коллайдером госпитализация, мы увидим обратную связь между курением и COVID-19! (Изображение автора)

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

Но на правом рисунке где мы рассматриваем только пациентов больниц мы видим противоположную тенденцию! Чтобы понять это, обратите внимание на следующие моменты.

1.Тяжелая форма COVID-19 увеличивает шансы на госпитализацию. То есть, если степень тяжести заболевания выше 1, то требуется госпитализация.

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

3. Таким образом, если у пациента легкая форма COVID-19, он имеет больше шансов оказаться курильщиком! Более того, в отличие от COVID-19 причиной для госпитализации станет наличие у пациента какого-либо заболевания, которое может быть вызвано курением (например, сердечно-сосудистые заболевания, рак, диабет).

Этот пример очень похож на оригинальную работу Berkson 1946, где автор заметил отрицательную корреляцию между холециститом и диабетом у пациентов больниц, несмотря на то, что диабет является фактором риска холецистита.

2. Скрытые (латентные) переменные

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

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

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

Парадокс скрытой переменной: степень тяжести пожара это скрытая переменная для n задействованных пожарных и для n пострадавших. (Изображение автора)

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

Давайте рассмотрим следующий пример с набором данных. На левом рисунке у нас отражены общие данные по всем видам пожаров, а на правом рисунке мы рассматриваем только сведения, соответствующие трем фиксированным степеням тяжести пожара (т.е. мы обусловливаем наши данные наблюдений скрытой переменной).

Скрытые переменные: если мы примем за условие скрытую переменную степень тяжести пожара, мы увидим обратную корреляцию между количеством задействованных пожарных и количеством раненых! (Изображение автора)

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

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

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

3. Парадокс Симпсона

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

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

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

• С другой стороны, общий процент приема среди абитуриентов женского пола был ниже, чем у абитуриентов мужского пола.

Чтобы понять, как как такое может быть, давайте рассмотрим следующий набор данных с двумя факультетами: Факультет A и Факультет B.

• Из 100 абитуриентов мужского пола: 80 подали заявки на Факультет A, из которых 68 были приняты (85%), а 20 подали заявки на Факультет В, из которых приняты были 12 человек (60%).

• Из 100 абитуриентов женского пола: 30 подали заявки на Факультет А, из которых 28 были приняты (93%), в то время как 70 подали заявки на Факультет B, из которых были приняты 46 (66%).

Парадокс Симпсона: женщины-абитуриенты с большей вероятностью будут приняты в каждом факультете, но общий процент приема женщин в сравнении с мужчинами ниже! (Изображение автора)

Парадокс выражается следующими неравенствами.

Парадокс Симпсона: неравенство, лежащее в основе очевидного противоречия. (Изображение автора)

Теперь мы можем понять происхождение наших, казалось бы, противоречивых наблюдений. Дело в том, что существует ощутимый классовый гендерный дисбаланс среди абитуриентов на каждом из двух факультетов (Факультет A: 8030, Факультет B: 2070). Действительно, большинство студентов женского пола подали заявку на более конкурентный Факультет B (который имеет низкие показатели приема), в то время как большинство студентов мужского пола подали документы на менее конкурентный Факультет А (который имеет более высокие показатели приема). Это обусловливает противоречивые данные, которые мы получили.

Заключение

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

Узнать подробнее о курсе "Machine Learning. Professional"

Участвовать в онлайн-интенсиве Деплой ML модели: от грязного кода в ноутбуке к рабочему сервису