{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "RVU1Q900Dl2r"
},
"source": [
"##### Copyright 2020 The TensorFlow Authors.\n",
"##### Copyright 2020 Sen Pei (Columbia University).\n",
"\n",
"Licensed under the Apache License, Version 2.0 (the \"License\");"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "QKUGkj9ZDmfg"
},
"outputs": [],
"source": [
"#@title Licensed under the Apache License, Version 2.0 (the \"License\"); { display-mode: \"form\" }\n",
"# you may not use this file except in compliance with the License.\n",
"# You may obtain a copy of the License at\n",
"#\n",
"# https://www.apache.org/licenses/LICENSE-2.0\n",
"#\n",
"# Unless required by applicable law or agreed to in writing, software\n",
"# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
"# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
"# See the License for the specific language governing permissions and\n",
"# limitations under the License."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "KCZNiAlqEAr1"
},
"source": [
"# Substantial Undocumented Infection Facilitates the Rapid Dissemination of Novel Coronavirus (SARS-CoV2)\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "iBdidAV1Rw-G"
},
"source": [
"This is a TensorFlow Probability port of the eponymous 16 March 2020 paper by Li et al. We faithfully reproduce the original authors' methods and results on the TensorFlow Probability platform, showcasing some of TFP's capabilities in the setting of modern epidemiology modeling. Porting to TensorFlow gives us a ~10x speedup relative to the original Matlab code, and, since TensorFlow Probability pervasively supports vectorized batch computation, also favorably scales to hundreds of independent replications.\n",
"\n",
"### Original paper\n",
"\n",
"Ruiyun Li, Sen Pei, Bin Chen, Yimeng Song, Tao Zhang, Wan Yang, and Jeffrey Shaman. Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV2). (2020), doi:\n",
"https://doi.org/10.1126/science.abb3221 .\n",
"\n",
"*Abstract:* \"Estimation of the prevalence and contagiousness of undocumented novel coronavirus (SARS-CoV2) infections is critical for understanding the overall prevalence and pandemic potential of this disease. Here we use observations of reported infection within China, in conjunction with mobility data, a networked dynamic metapopulation model and Bayesian inference, to infer critical epidemiological characteristics associated with SARS-CoV2, including the fraction of undocumented infections and their contagiousness. We estimate 86% of all infections were undocumented (95% CI: [82%–90%]) prior to 23 January 2020 travel restrictions. Per person, the transmission rate of undocumented infections was 55% of documented infections ([46%–62%]), yet, due to their greater numbers, undocumented infections were the infection source for 79% of documented cases. These findings explain the rapid geographic spread of SARS-CoV2 and indicate containment of this virus will be particularly challenging.\"\n",
"\n",
"[Github link](https://github.com/SenPei-CU/COVID-19) to the code and data."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "0yETSzjhoVBN"
},
"source": [
"## Overview"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "HtZGdFItx4Z3"
},
"source": [
"The model is a [compartmental disease model](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology), with compartments for \"susceptible\", \"exposed\" (infected but not yet infectious), \"never-documented infectious\", and \"eventually-documented infectious\". There are two noteworthy features: separate compartments for each of 375 Chinese cities, with an assumption about how people travel from one city to another; and delays in reporting infection, so that a case that becomes \"eventually-documented infectious\" on day $t$ doesn't show up in the observed case counts until a stochastic later day.\n",
"\n",
"The model assumes that the never-documented cases end up undocumented by being milder, and thus infect others at a lower rate. The main parameter of interest in the original paper is the proportion of cases that go undocumented, to estimate both the extent of existing infection, and the impact of undocumented transmission on the spread of the disease.\n",
"\n",
"This colab is structured as a code walkthrough in bottom-up style. In order, we will\n",
"- Ingest and briefly examine the data,\n",
"- Define the state space and dynamics of the model,\n",
"- Build up a suite of functions for doing inference in the model following Li et al, and\n",
"- Invoke them and examine the results. Spoiler: They come out the same as the paper.\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "Brvp_GH-ItsG"
},
"source": [
"## Installation and Python Imports"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "iUEZ07T3yuBj"
},
"outputs": [],
"source": [
"!pip3 install -q tf-nightly tfp-nightly"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "-Za3hnj2y5M9"
},
"outputs": [],
"source": [
"import collections\n",
"import io\n",
"import requests\n",
"import time\n",
"import zipfile\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"import tensorflow.compat.v2 as tf\n",
"import tensorflow_probability as tfp\n",
"from tensorflow_probability.python.internal import samplers\n",
"\n",
"tfd = tfp.distributions\n",
"tfes = tfp.experimental.sequential"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "KYfPl6uqzIs1"
},
"source": [
"## Data Import\n",
"\n",
"Let's import the data from github and inspect some of it."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "1q22N4rPzXIg"
},
"outputs": [],
"source": [
"r = requests.get('https://raw.githubusercontent.com/SenPei-CU/COVID-19/master/Data.zip')\n",
"z = zipfile.ZipFile(io.BytesIO(r.content))\n",
"z.extractall('/tmp/')\n",
"raw_incidence = pd.read_csv('/tmp/data/Incidence.csv')\n",
"raw_mobility = pd.read_csv('/tmp/data/Mobility.csv')\n",
"raw_population = pd.read_csv('/tmp/data/pop.csv')"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "IC3Nt8lJSllh"
},
"source": [
"Below we can see the raw incidence count per day. We are most interested in the first 14 days (January 10th to January 23rd), as the travel restrictions were put in place on the 23rd. The paper deals with this by modeling Jan 10-23 and Jan 23+ separately, with different parameters; we will just restrict our reproduction to the earlier period."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {
"height": 1000
},
"colab_type": "code",
"id": "yFP-Tah4zfne",
"outputId": "f02e7cf1-e4ae-4a56-ad1c-8bbfd61fde4a"
},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
" City Population\n",
"0 Beijing 21540000.0\n",
"1 Tianjin 11760000.0\n",
"2 Shijiazhuang 10151200.0\n",
"3 Tangshan 7843600.0\n",
"4 Qinhuangdao 3094600.0\n",
"5 Handan 9492800.0\n",
"6 Xintai 7319900.0\n",
"7 Baoding 10425300.0\n",
"8 Zhangjiakou 4425100.0\n",
"9 Chengde 3531800.0\n",
"10 Cangzhou 7505500.0\n",
"11 Langfang 4615000.0\n",
"12 Hengshui 4453100.0\n",
"13 Taiyuan 4421400.0\n",
"14 Datong 3456000.0\n",
"15 Yangquan 1414400.0\n",
"16 Changzhi 3468200.0\n",
"17 Jincheng 2343100.0\n",
"18 Shuozhou 1781200.0\n",
"19 Jinzhong 3381600.0\n",
"20 Yuncheng 5359600.0\n",
"21 Xinzhou 3171200.0\n",
"22 Linfen 4500300.0\n",
"23 Lvliang 3885600.0\n",
"24 Huhehaote 2437900.0\n",
"25 Baotou 2297400.0\n",
"26 Wuhai 434900.0\n",
"27 Chifeng 4518000.0\n",
"28 Tongliao 3083500.0\n",
"29 Eerduosi 2078400.0\n",
".. ... ...\n",
"345 Urumuqi 2330000.0\n",
"346 Kelamayi 270000.0\n",
"347 Tulufan 590000.0\n",
"348 Hami 550000.0\n",
"349 Changji 1590000.0\n",
"350 Bozhou (Xinjiang) 470000.0\n",
"351 Bazhou (Xinjiang) 1170000.0\n",
"352 Akesu 2260000.0\n",
"353 Kezhou 470000.0\n",
"354 Kashi 3650000.0\n",
"355 Hetian (Xinjiang) 1800000.0\n",
"356 Yili 2600000.0\n",
"357 Tacheng 990000.0\n",
"358 Aletai 640000.0\n",
"359 Shihezi 660000.0\n",
"360 Ala'er 326800.0\n",
"361 Tumushuke 255600.0\n",
"362 Wujiaqu (Xinjiang) 125600.0\n",
"363 Beitun 80000.0\n",
"364 Wujiaqushi (Xinjiang) 125600.0\n",
"365 Tiemenguan 50000.0\n",
"366 Shuanghe 53800.0\n",
"367 Kekedala 75000.0\n",
"368 Taiwan 23580000.0\n",
"369 Hongkong 7450000.0\n",
"370 Macao 630000.0\n",
"371 Hetianxian (Xinjiang) 360000.0\n",
"372 Moyu 646200.0\n",
"373 Yutian 289500.0\n",
"374 Gejiu 472000.0\n",
"\n",
"[375 rows x 2 columns]"
]
},
"execution_count": 0,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"raw_population"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "Xf6xLARfsBwe"
},
"source": [
"Let's also check and record which entry is Wuhan."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {
"height": 35
},
"colab_type": "code",
"id": "phQQFI401Pr8",
"outputId": "b0d8334e-c714-4a8e-813f-de31aa6947e0"
},
"outputs": [
{
"data": {
"text/plain": [
"'Wuhan'"
]
},
"execution_count": 0,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"raw_population['City'][169]"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "uS5vWXkBsQLG"
},
"outputs": [],
"source": [
"WUHAN_IDX = 169"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "3RWGcEfL65sS"
},
"source": [
"And here we see the mobility matrix between different cities. This is a proxy for the number of people moving between different cities on the first 14 days. It's dervied from GPS records provided by Tencent for the 2018 Lunar New Year season. Li et al model mobility during the 2020 season as some unknown (subject to inference) constant factor $\\theta$ times this."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {
"height": 1000
},
"colab_type": "code",
"id": "BvPonRBJ1AIQ",
"outputId": "95eec3c7-cf53-4f30-f4b9-fe6742fbb5fe"
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Day
\n",
"
Origin
\n",
"
Destination
\n",
"
Mobility Index
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
1
\n",
"
Beijing
\n",
"
Tianjin
\n",
"
17585
\n",
"
\n",
"
\n",
"
1
\n",
"
1
\n",
"
Beijing
\n",
"
Shijiazhuang
\n",
"
21289
\n",
"
\n",
"
\n",
"
2
\n",
"
1
\n",
"
Beijing
\n",
"
Handan
\n",
"
27676
\n",
"
\n",
"
\n",
"
3
\n",
"
1
\n",
"
Beijing
\n",
"
Baoding
\n",
"
19007
\n",
"
\n",
"
\n",
"
4
\n",
"
1
\n",
"
Beijing
\n",
"
Zhangjiakou
\n",
"
20327
\n",
"
\n",
"
\n",
"
5
\n",
"
1
\n",
"
Beijing
\n",
"
Langfang
\n",
"
31180
\n",
"
\n",
"
\n",
"
6
\n",
"
1
\n",
"
Beijing
\n",
"
Ha'erbin
\n",
"
24390
\n",
"
\n",
"
\n",
"
7
\n",
"
1
\n",
"
Beijing
\n",
"
Shanghai
\n",
"
66799
\n",
"
\n",
"
\n",
"
8
\n",
"
1
\n",
"
Beijing
\n",
"
Nanjing
\n",
"
17719
\n",
"
\n",
"
\n",
"
9
\n",
"
1
\n",
"
Beijing
\n",
"
Hangzhou
\n",
"
28536
\n",
"
\n",
"
\n",
"
10
\n",
"
1
\n",
"
Beijing
\n",
"
Wenzhou
\n",
"
16767
\n",
"
\n",
"
\n",
"
11
\n",
"
1
\n",
"
Beijing
\n",
"
Xinyang
\n",
"
17522
\n",
"
\n",
"
\n",
"
12
\n",
"
1
\n",
"
Beijing
\n",
"
Wuhan
\n",
"
32517
\n",
"
\n",
"
\n",
"
13
\n",
"
1
\n",
"
Beijing
\n",
"
Changsha
\n",
"
35671
\n",
"
\n",
"
\n",
"
14
\n",
"
1
\n",
"
Beijing
\n",
"
Chongqing
\n",
"
66639
\n",
"
\n",
"
\n",
"
15
\n",
"
1
\n",
"
Beijing
\n",
"
Xi'an
\n",
"
17088
\n",
"
\n",
"
\n",
"
16
\n",
"
1
\n",
"
Tianjin
\n",
"
Beijing
\n",
"
13942
\n",
"
\n",
"
\n",
"
17
\n",
"
1
\n",
"
Tianjin
\n",
"
Tangshan
\n",
"
6631
\n",
"
\n",
"
\n",
"
18
\n",
"
1
\n",
"
Tianjin
\n",
"
Handan
\n",
"
5140
\n",
"
\n",
"
\n",
"
19
\n",
"
1
\n",
"
Tianjin
\n",
"
Baoding
\n",
"
6398
\n",
"
\n",
"
\n",
"
20
\n",
"
1
\n",
"
Tianjin
\n",
"
Cangzhou
\n",
"
7030
\n",
"
\n",
"
\n",
"
21
\n",
"
1
\n",
"
Tianjin
\n",
"
Langfang
\n",
"
7363
\n",
"
\n",
"
\n",
"
22
\n",
"
1
\n",
"
Tianjin
\n",
"
Changchun
\n",
"
4258
\n",
"
\n",
"
\n",
"
23
\n",
"
1
\n",
"
Tianjin
\n",
"
Ha'erbin
\n",
"
12675
\n",
"
\n",
"
\n",
"
24
\n",
"
1
\n",
"
Tianjin
\n",
"
Shanghai
\n",
"
8698
\n",
"
\n",
"
\n",
"
25
\n",
"
1
\n",
"
Tianjin
\n",
"
Hangzhou
\n",
"
4364
\n",
"
\n",
"
\n",
"
26
\n",
"
1
\n",
"
Tianjin
\n",
"
Fuzhou (Fujian)
\n",
"
14632
\n",
"
\n",
"
\n",
"
27
\n",
"
1
\n",
"
Tianjin
\n",
"
Dezhou
\n",
"
6857
\n",
"
\n",
"
\n",
"
28
\n",
"
1
\n",
"
Tianjin
\n",
"
Wuhan
\n",
"
4669
\n",
"
\n",
"
\n",
"
29
\n",
"
1
\n",
"
Tianjin
\n",
"
Changsha
\n",
"
6579
\n",
"
\n",
"
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
\n",
"
\n",
"
92218
\n",
"
14
\n",
"
Gejiu
\n",
"
Changsha
\n",
"
52
\n",
"
\n",
"
\n",
"
92219
\n",
"
14
\n",
"
Gejiu
\n",
"
Hengyang
\n",
"
42
\n",
"
\n",
"
\n",
"
92220
\n",
"
14
\n",
"
Gejiu
\n",
"
Changde
\n",
"
53
\n",
"
\n",
"
\n",
"
92221
\n",
"
14
\n",
"
Gejiu
\n",
"
Chenzhou
\n",
"
36
\n",
"
\n",
"
\n",
"
92222
\n",
"
14
\n",
"
Gejiu
\n",
"
Yongzhou
\n",
"
42
\n",
"
\n",
"
\n",
"
92223
\n",
"
14
\n",
"
Gejiu
\n",
"
Guangzhou
\n",
"
48
\n",
"
\n",
"
\n",
"
92224
\n",
"
14
\n",
"
Gejiu
\n",
"
Shenzhen
\n",
"
77
\n",
"
\n",
"
\n",
"
92225
\n",
"
14
\n",
"
Gejiu
\n",
"
Huizhou
\n",
"
43
\n",
"
\n",
"
\n",
"
92226
\n",
"
14
\n",
"
Gejiu
\n",
"
Nannig
\n",
"
30
\n",
"
\n",
"
\n",
"
92227
\n",
"
14
\n",
"
Gejiu
\n",
"
Beihai
\n",
"
30
\n",
"
\n",
"
\n",
"
92228
\n",
"
14
\n",
"
Gejiu
\n",
"
Baise
\n",
"
37
\n",
"
\n",
"
\n",
"
92229
\n",
"
14
\n",
"
Gejiu
\n",
"
Chongqing
\n",
"
132
\n",
"
\n",
"
\n",
"
92230
\n",
"
14
\n",
"
Gejiu
\n",
"
Neijiang
\n",
"
47
\n",
"
\n",
"
\n",
"
92231
\n",
"
14
\n",
"
Gejiu
\n",
"
Liupanshui
\n",
"
41
\n",
"
\n",
"
\n",
"
92232
\n",
"
14
\n",
"
Gejiu
\n",
"
Bijie
\n",
"
32
\n",
"
\n",
"
\n",
"
92233
\n",
"
14
\n",
"
Gejiu
\n",
"
Kunming
\n",
"
452
\n",
"
\n",
"
\n",
"
92234
\n",
"
14
\n",
"
Gejiu
\n",
"
Qujing
\n",
"
359
\n",
"
\n",
"
\n",
"
92235
\n",
"
14
\n",
"
Gejiu
\n",
"
Yuxi
\n",
"
170
\n",
"
\n",
"
\n",
"
92236
\n",
"
14
\n",
"
Gejiu
\n",
"
Shaotong
\n",
"
94
\n",
"
\n",
"
\n",
"
92237
\n",
"
14
\n",
"
Gejiu
\n",
"
Pu'er
\n",
"
50
\n",
"
\n",
"
\n",
"
92238
\n",
"
14
\n",
"
Gejiu
\n",
"
Chuxiong
\n",
"
62
\n",
"
\n",
"
\n",
"
92239
\n",
"
14
\n",
"
Gejiu
\n",
"
Wenshan
\n",
"
138
\n",
"
\n",
"
\n",
"
92240
\n",
"
14
\n",
"
Gejiu
\n",
"
Dali
\n",
"
50
\n",
"
\n",
"
\n",
"
92241
\n",
"
14
\n",
"
Gejiu
\n",
"
Xi'an
\n",
"
55
\n",
"
\n",
"
\n",
"
92242
\n",
"
14
\n",
"
Gejiu
\n",
"
Xianyang
\n",
"
38
\n",
"
\n",
"
\n",
"
92243
\n",
"
14
\n",
"
Gejiu
\n",
"
Lanzhou
\n",
"
60
\n",
"
\n",
"
\n",
"
92244
\n",
"
14
\n",
"
Gejiu
\n",
"
Pingliang
\n",
"
73
\n",
"
\n",
"
\n",
"
92245
\n",
"
14
\n",
"
Gejiu
\n",
"
Dongxi
\n",
"
126
\n",
"
\n",
"
\n",
"
92246
\n",
"
14
\n",
"
Gejiu
\n",
"
Xining
\n",
"
31
\n",
"
\n",
"
\n",
"
92247
\n",
"
14
\n",
"
Gejiu
\n",
"
Guyuan
\n",
"
33
\n",
"
\n",
" \n",
"
\n",
"
92248 rows × 4 columns
\n",
"
"
],
"text/plain": [
" Day Origin Destination Mobility Index\n",
"0 1 Beijing Tianjin 17585\n",
"1 1 Beijing Shijiazhuang 21289\n",
"2 1 Beijing Handan 27676\n",
"3 1 Beijing Baoding 19007\n",
"4 1 Beijing Zhangjiakou 20327\n",
"5 1 Beijing Langfang 31180\n",
"6 1 Beijing Ha'erbin 24390\n",
"7 1 Beijing Shanghai 66799\n",
"8 1 Beijing Nanjing 17719\n",
"9 1 Beijing Hangzhou 28536\n",
"10 1 Beijing Wenzhou 16767\n",
"11 1 Beijing Xinyang 17522\n",
"12 1 Beijing Wuhan 32517\n",
"13 1 Beijing Changsha 35671\n",
"14 1 Beijing Chongqing 66639\n",
"15 1 Beijing Xi'an 17088\n",
"16 1 Tianjin Beijing 13942\n",
"17 1 Tianjin Tangshan 6631\n",
"18 1 Tianjin Handan 5140\n",
"19 1 Tianjin Baoding 6398\n",
"20 1 Tianjin Cangzhou 7030\n",
"21 1 Tianjin Langfang 7363\n",
"22 1 Tianjin Changchun 4258\n",
"23 1 Tianjin Ha'erbin 12675\n",
"24 1 Tianjin Shanghai 8698\n",
"25 1 Tianjin Hangzhou 4364\n",
"26 1 Tianjin Fuzhou (Fujian) 14632\n",
"27 1 Tianjin Dezhou 6857\n",
"28 1 Tianjin Wuhan 4669\n",
"29 1 Tianjin Changsha 6579\n",
"... ... ... ... ...\n",
"92218 14 Gejiu Changsha 52\n",
"92219 14 Gejiu Hengyang 42\n",
"92220 14 Gejiu Changde 53\n",
"92221 14 Gejiu Chenzhou 36\n",
"92222 14 Gejiu Yongzhou 42\n",
"92223 14 Gejiu Guangzhou 48\n",
"92224 14 Gejiu Shenzhen 77\n",
"92225 14 Gejiu Huizhou 43\n",
"92226 14 Gejiu Nannig 30\n",
"92227 14 Gejiu Beihai 30\n",
"92228 14 Gejiu Baise 37\n",
"92229 14 Gejiu Chongqing 132\n",
"92230 14 Gejiu Neijiang 47\n",
"92231 14 Gejiu Liupanshui 41\n",
"92232 14 Gejiu Bijie 32\n",
"92233 14 Gejiu Kunming 452\n",
"92234 14 Gejiu Qujing 359\n",
"92235 14 Gejiu Yuxi 170\n",
"92236 14 Gejiu Shaotong 94\n",
"92237 14 Gejiu Pu'er 50\n",
"92238 14 Gejiu Chuxiong 62\n",
"92239 14 Gejiu Wenshan 138\n",
"92240 14 Gejiu Dali 50\n",
"92241 14 Gejiu Xi'an 55\n",
"92242 14 Gejiu Xianyang 38\n",
"92243 14 Gejiu Lanzhou 60\n",
"92244 14 Gejiu Pingliang 73\n",
"92245 14 Gejiu Dongxi 126\n",
"92246 14 Gejiu Xining 31\n",
"92247 14 Gejiu Guyuan 33\n",
"\n",
"[92248 rows x 4 columns]"
]
},
"execution_count": 0,
"metadata": {
"tags": []
},
"output_type": "execute_result"
}
],
"source": [
"raw_mobility"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "vSY98S14Oaqs"
},
"source": [
"Finally, let's preprocess all this into numpy arrays that we can consume."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "9txbadTC14YP"
},
"outputs": [],
"source": [
"# The given populations are only \"initial\" because of intercity mobility during\n",
"# the holiday season.\n",
"initial_population = raw_population['Population'].to_numpy().astype(np.float32)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "w-pKDK461zFX"
},
"source": [
"Convert the mobility data into an [L, L, T]-shaped Tensor, where L is the number of locations, and T is the number of timesteps."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "f-2zASNxw_wi"
},
"outputs": [],
"source": [
"daily_mobility_matrices = []\n",
"for i in range(1, 15):\n",
" day_mobility = raw_mobility[raw_mobility['Day'] == i]\n",
" \n",
" # Make a matrix of daily mobilities.\n",
" z = pd.crosstab(\n",
" day_mobility.Origin, \n",
" day_mobility.Destination, \n",
" values=day_mobility['Mobility Index'], aggfunc='sum', dropna=False)\n",
" \n",
" # Include every city, even if there are no rows for some in the raw data on\n",
" # some day. This uses the sort order of `raw_population`.\n",
" z = z.reindex(index=raw_population['City'], columns=raw_population['City'], \n",
" fill_value=0)\n",
" # Finally, fill any missing entries with 0. This means no mobility.\n",
" z = z.fillna(0)\n",
" daily_mobility_matrices.append(z.to_numpy())\n",
"\n",
"mobility_matrix_over_time = np.stack(daily_mobility_matrices, axis=-1).astype(\n",
" np.float32)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "GaJYQFT91-15"
},
"source": [
"Finally take the observed infections and make an [L, T] table."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "qs4TKTnv176D"
},
"outputs": [],
"source": [
"# Remove the date parameter and take the first 14 days.\n",
"observed_daily_infectious_count = raw_incidence.to_numpy()[:14, 1:]\n",
"observed_daily_infectious_count = np.transpose(\n",
" observed_daily_infectious_count).astype(np.float32)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "B9WeGRug2HLQ"
},
"source": [
"And double-check that we got the shapes the way we wanted. As a reminder, we're working with 375 cities and 14 days."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {
"height": 69
},
"colab_type": "code",
"id": "dymKgF8p1hgi",
"outputId": "040b68e0-2d29-462a-a208-dd82c75c02b9"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Mobility Matrix over time should have shape (375, 375, 14): (375, 375, 14)\n",
"Observed Infectious should have shape (375, 14): (375, 14)\n",
"Initial population should have shape (375): (375,)\n"
]
}
],
"source": [
"print('Mobility Matrix over time should have shape (375, 375, 14): {}'.format(\n",
" mobility_matrix_over_time.shape))\n",
"print('Observed Infectious should have shape (375, 14): {}'.format(\n",
" observed_daily_infectious_count.shape))\n",
"print('Initial population should have shape (375): {}'.format(\n",
" initial_population.shape))"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "LGTcVy2K2mxv"
},
"source": [
"## Defining State and Parameters\n",
"\n",
"\n",
"Let's start defining our model. The model we are reproducing is a variant of an [SEIR model](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SEIR_model). In this case we have the following time-varying states:\n",
"* $S$: Number of people susceptible to the disease in each city.\n",
"* $E$: Number of people in each city exposed to the disease but not infectious yet. Biologically, this corresponds to contracting the disease, in that all exposed people eventually become infectious.\n",
"* $I^u$: Number of people in each city who are infectious but undocumented. In the model, this actually means \"will never be documented\".\n",
"* $I^r$: Number of people in each city who are infectious and documented as such. Li et al model reporting delays, so $I^r$ actually corresponds to something like \"case is severe enough to be documented at some point in the future\".\n",
"\n",
"As we will see below, we will be inferring these states by running an Ensemble-adjusted Kalman Filter (EAKF) forward in time. The state vector of the EAKF is one city-indexed vector for each of these quantities.\n",
"\n",
"The model has the following inferrable global, time-invariant parameters:\n",
"\n",
"* $\\beta$: The transmission rate due to documented-infectious individuals.\n",
"* $\\mu$: The relative transmission rate due to undocumented-infectious\n",
" individuals. This will act through the product $\\mu \\beta$.\n",
"* $\\theta$: The intercity mobility factor. This is a factor greater than\n",
" 1 correcting for underreporting of mobility data (and for population growth\n",
" from 2018 to 2020).\n",
"* $Z$: The average incubation period (i.e., time in the \"exposed\" state).\n",
"* $\\alpha$: This is the fraction of infections severe enough to be (eventually) documented.\n",
"* $D$: The average duration of infections (i.e., time in either \"infectious\" state).\n",
"\n",
"We will be inferring point estimates for these parameters with an Iterative-Filtering loop around the EAKF for the states.\n",
"\n",
"The model also depends on un-inferred constants:\n",
"* $M$: The intercity mobility matrix. This is time-varying and presumed given. Recall that it's scaled by the inferred parameter $\\theta$ to give the actual population movements between cities.\n",
"* $N$: The total number of people in each city. The initial populations are taken as given, and the time-variation of population is computed from the mobility numbers $\\theta M$.\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "--4jYF-Q9maG"
},
"source": [
"First, we give ourselves some data structures for holding our states and parameters.\n"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "VacihFdC25l4"
},
"outputs": [],
"source": [
"SEIRComponents = collections.namedtuple(\n",
" typename='SEIRComponents',\n",
" field_names=[\n",
" 'susceptible', # S\n",
" 'exposed', # E\n",
" 'documented_infectious', # I^r\n",
" 'undocumented_infectious', # I^u\n",
" # This is the count of new cases in the \"documented infectious\" compartment.\n",
" # We need this because we will introduce a reporting delay, between a person\n",
" # entering I^r and showing up in the observable case count data.\n",
" # This can't be computed from the cumulative `documented_infectious` count,\n",
" # because some portion of that population will move to the 'recovered'\n",
" # state, which we aren't tracking explicitly.\n",
" 'daily_new_documented_infectious'])\n",
"\n",
"ModelParams = collections.namedtuple(\n",
" typename='ModelParams',\n",
" field_names=[\n",
" 'documented_infectious_tx_rate', # Beta\n",
" 'undocumented_infectious_tx_relative_rate', # Mu\n",
" 'intercity_underreporting_factor', # Theta\n",
" 'average_latency_period', # Z\n",
" 'fraction_of_documented_infections', # Alpha\n",
" 'average_infection_duration' # D\n",
" ]\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "kBgUMzTL_2gc"
},
"source": [
"We also code Li et al's bounds for the values of the parameters."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "H_HWdcXZ_1Cw"
},
"outputs": [],
"source": [
"PARAMETER_LOWER_BOUNDS = ModelParams(\n",
" documented_infectious_tx_rate=0.8,\n",
" undocumented_infectious_tx_relative_rate=0.2,\n",
" intercity_underreporting_factor=1.,\n",
" average_latency_period=2.,\n",
" fraction_of_documented_infections=0.02,\n",
" average_infection_duration=2.\n",
")\n",
"\n",
"PARAMETER_UPPER_BOUNDS = ModelParams(\n",
" documented_infectious_tx_rate=1.5,\n",
" undocumented_infectious_tx_relative_rate=1.,\n",
" intercity_underreporting_factor=1.75,\n",
" average_latency_period=5.,\n",
" fraction_of_documented_infections=1.,\n",
" average_infection_duration=5.\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "-yizvLfw-CiT"
},
"source": [
"## SEIR Dynamics\n",
"\n",
"Here we define the relationship between the parameters and state.\n",
"\n",
"The time-dynamics equations from Li et al (supplemental material, eqns 1-5) are as follows:\n",
"\n",
"$\\frac{dS_i}{dt} = -\\beta \\frac{S_i I_i^r}{N_i} - \\mu \\beta \\frac{S_i I_i^u}{N_i} + \\theta \\sum_k \\frac{M_{ij} S_j}{N_j - I_j^r} - + \\theta \\sum_k \\frac{M_{ji} S_j}{N_i - I_i^r}$\n",
"\n",
"$\\frac{dE_i}{dt} = \\beta \\frac{S_i I_i^r}{N_i} + \\mu \\beta \\frac{S_i I_i^u}{N_i} -\\frac{E_i}{Z} + \\theta \\sum_k \\frac{M_{ij} E_j}{N_j - I_j^r} - + \\theta \\sum_k \\frac{M_{ji} E_j}{N_i - I_i^r}$\n",
"\n",
"$\\frac{dI^r_i}{dt} = \\alpha \\frac{E_i}{Z} - \\frac{I_i^r}{D}$\n",
"\n",
"$\\frac{dI^u_i}{dt} = (1 - \\alpha) \\frac{E_i}{Z} - \\frac{I_i^u}{D} + \\theta \\sum_k \\frac{M_{ij} I_j^u}{N_j - I_j^r} - + \\theta \\sum_k \\frac{M_{ji} I^u_j}{N_i - I_i^r}$\n",
"\n",
"$N_i = N_i + \\theta \\sum_j M_{ij} - \\theta \\sum_j M_{ji}$\n",
"\n",
"As a reminder, the $i$ and $j$ subscripts index cities. These equations model the time-evolution of the disease through\n",
"- Contact with infectious individuals leading to more infection;\n",
"- Disease progression from \"exposed\" to one of the \"infectious\" states;\n",
"- Disease progression from \"infectious\" states to recovery, which we model by removal from the modeled population;\n",
"- Inter-city mobility, including exposed or undocumented-infectious persons; and\n",
"- Time-variation of daily city populations through inter-city mobility.\n",
"\n",
"Following Li et al, we assume that people with cases severe enough to eventually be reported do not travel between cities.\n",
"\n",
"Also following Li et al, we treat these dynamics as subject to term-wise Poisson noise, i.e., each term is actually the rate of a Poisson, a sample from which gives the true change. The Poisson noise is term-wise because subtracting (as opposed to adding) Poisson samples does not yield a Poisson-distributed result.\n",
"\n",
"We will evolve these dynamics forward in time with the classic fourth-order Runge-Kutta integrator, but first let's define the function that computes them (including sampling the Poisson noise)."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "_Cpj5Lyo-AGw"
},
"outputs": [],
"source": [
"def sample_state_deltas(\n",
" state, population, mobility_matrix, params, seed, is_deterministic=False):\n",
" \"\"\"Computes one-step change in state, including Poisson sampling.\n",
" \n",
" Note that this is coded to support vectorized evaluation on arbitrary-shape\n",
" batches of states. This is useful, for example, for running multiple\n",
" independent replicas of this model to compute credible intervals for the\n",
" parameters. We refer to the arbitrary batch shape with the conventional\n",
" `B` in the parameter documentation below. This function also, of course,\n",
" supports broadcasting over the batch shape.\n",
"\n",
" Args:\n",
" state: A `SEIRComponents` tuple with fields Tensors of shape\n",
" B + [num_locations] giving the current disease state.\n",
" population: A Tensor of shape B + [num_locations] giving the current city\n",
" populations.\n",
" mobility_matrix: A Tensor of shape B + [num_locations, num_locations] giving\n",
" the current baseline inter-city mobility.\n",
" params: A `ModelParams` tuple with fields Tensors of shape B giving the\n",
" global parameters for the current EAKF run.\n",
" seed: Initial entropy for pseudo-random number generation. The Poisson\n",
" sampling is repeatable by supplying the same seed.\n",
" is_deterministic: A `bool` flag to turn off Poisson sampling if desired.\n",
"\n",
" Returns:\n",
" delta: A `SEIRComponents` tuple with fields Tensors of shape\n",
" B + [num_locations] giving the one-day changes in the state, according\n",
" to equations 1-4 above (including Poisson noise per Li et al).\n",
" \"\"\"\n",
" undocumented_infectious_fraction = state.undocumented_infectious / population\n",
" documented_infectious_fraction = state.documented_infectious / population\n",
"\n",
" # Anyone not documented as infectious is considered mobile\n",
" mobile_population = (population - state.documented_infectious)\n",
" def compute_outflow(compartment_population):\n",
" raw_mobility = tf.linalg.matvec(\n",
" mobility_matrix, compartment_population / mobile_population)\n",
" return params.intercity_underreporting_factor * raw_mobility\n",
" def compute_inflow(compartment_population):\n",
" raw_mobility = tf.linalg.matmul(\n",
" mobility_matrix,\n",
" (compartment_population / mobile_population)[..., tf.newaxis],\n",
" transpose_a=True)\n",
" return params.intercity_underreporting_factor * tf.squeeze(\n",
" raw_mobility, axis=-1)\n",
"\n",
" # Helper for sampling the Poisson-variate terms.\n",
" seeds = samplers.split_seed(seed, n=11)\n",
" if is_deterministic:\n",
" def sample_poisson(rate):\n",
" return rate\n",
" else:\n",
" def sample_poisson(rate):\n",
" return tfd.Poisson(rate=rate).sample(seed=seeds.pop())\n",
"\n",
" # Below are the various terms called U1-U12 in the paper. We combined the\n",
" # first two, which should be fine; both are poisson so their sum is too, and\n",
" # there's no risk (as there could be in other terms) of going negative.\n",
" susceptible_becoming_exposed = sample_poisson(\n",
" state.susceptible *\n",
" (params.documented_infectious_tx_rate *\n",
" documented_infectious_fraction +\n",
" (params.undocumented_infectious_tx_relative_rate *\n",
" params.documented_infectious_tx_rate) *\n",
" undocumented_infectious_fraction)) # U1 + U2\n",
"\n",
" susceptible_population_inflow = sample_poisson(\n",
" compute_inflow(state.susceptible)) # U3\n",
" susceptible_population_outflow = sample_poisson(\n",
" compute_outflow(state.susceptible)) # U4\n",
"\n",
" exposed_becoming_documented_infectious = sample_poisson(\n",
" params.fraction_of_documented_infections *\n",
" state.exposed / params.average_latency_period) # U5\n",
" exposed_becoming_undocumented_infectious = sample_poisson(\n",
" (1 - params.fraction_of_documented_infections) *\n",
" state.exposed / params.average_latency_period) # U6\n",
"\n",
" exposed_population_inflow = sample_poisson(\n",
" compute_inflow(state.exposed)) # U7\n",
" exposed_population_outflow = sample_poisson(\n",
" compute_outflow(state.exposed)) # U8\n",
"\n",
" documented_infectious_becoming_recovered = sample_poisson(\n",
" state.documented_infectious /\n",
" params.average_infection_duration) # U9\n",
" undocumented_infectious_becoming_recovered = sample_poisson(\n",
" state.undocumented_infectious /\n",
" params.average_infection_duration) # U10\n",
"\n",
" undocumented_infectious_population_inflow = sample_poisson(\n",
" compute_inflow(state.undocumented_infectious)) # U11\n",
" undocumented_infectious_population_outflow = sample_poisson(\n",
" compute_outflow(state.undocumented_infectious)) # U12\n",
"\n",
" # The final state_deltas\n",
" return SEIRComponents(\n",
" # Equation [1]\n",
" susceptible=(-susceptible_becoming_exposed +\n",
" susceptible_population_inflow +\n",
" -susceptible_population_outflow),\n",
" # Equation [2]\n",
" exposed=(susceptible_becoming_exposed +\n",
" -exposed_becoming_documented_infectious +\n",
" -exposed_becoming_undocumented_infectious +\n",
" exposed_population_inflow +\n",
" -exposed_population_outflow),\n",
" # Equation [3]\n",
" documented_infectious=(\n",
" exposed_becoming_documented_infectious +\n",
" -documented_infectious_becoming_recovered),\n",
" # Equation [4]\n",
" undocumented_infectious=(\n",
" exposed_becoming_undocumented_infectious +\n",
" -undocumented_infectious_becoming_recovered +\n",
" undocumented_infectious_population_inflow +\n",
" -undocumented_infectious_population_outflow),\n",
" # New to-be-documented infectious cases, subject to the delayed\n",
" # observation model.\n",
" daily_new_documented_infectious=exposed_becoming_documented_infectious)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "K7T_QQTRKqOE"
},
"source": [
"Here's the integrator. This is completely standard, except for passing the PRNG seed through to the `sample_state_deltas` function to get independent Poisson noise at each of the partial steps that the Runge-Kutta method calls for."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "2qx4tX8sCAZv"
},
"outputs": [],
"source": [
"@tf.function(autograph=False)\n",
"def rk4_one_step(state, population, mobility_matrix, params, seed):\n",
" \"\"\"Implement one step of RK4, wrapped around a call to sample_state_deltas.\"\"\"\n",
" # One seed for each RK sub-step\n",
" seeds = samplers.split_seed(seed, n=4)\n",
"\n",
" deltas = tf.nest.map_structure(tf.zeros_like, state)\n",
" combined_deltas = tf.nest.map_structure(tf.zeros_like, state)\n",
"\n",
" for a, b in zip([1., 2, 2, 1.], [6., 3., 3., 6.]):\n",
" next_input = tf.nest.map_structure(\n",
" lambda x, delta, a=a: x + delta / a, state, deltas)\n",
" deltas = sample_state_deltas(\n",
" next_input,\n",
" population,\n",
" mobility_matrix,\n",
" params,\n",
" seed=seeds.pop(), is_deterministic=False)\n",
" combined_deltas = tf.nest.map_structure(\n",
" lambda x, delta, b=b: x + delta / b, combined_deltas, deltas)\n",
"\n",
" return tf.nest.map_structure(\n",
" lambda s, delta: s + tf.round(delta),\n",
" state, combined_deltas)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "3gVKuTqi32Bn"
},
"source": [
"## Initialization\n",
" \n",
"Here we implement the initialization scheme from the paper.\n",
"\n",
"Following Li et al, our inference scheme will be an ensemble adjustment Kalman filter inner loop, surrounded by an iterated filtering outer loop (IF-EAKF). Computationally, that means we need three kinds of initialization:\n",
"- Initial state for the inner EAKF\n",
"- Initial parameters for the outer IF, which are also the initial\n",
" parameters for the first EAKF\n",
"- Updating parameters from one IF iteration to the next, which serve\n",
" as the initial parameters for each EAKF other than the first.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "rdRf_UmdvPLo"
},
"outputs": [],
"source": [
"def initialize_state(num_particles, num_batches, seed):\n",
" \"\"\"Initialize the state for a batch of EAKF runs.\n",
" \n",
" Args:\n",
" num_particles: `int` giving the number of particles for the EAKF.\n",
" num_batches: `int` giving the number of independent EAKF runs to\n",
" initialize in a vectorized batch.\n",
" seed: PRNG entropy.\n",
" \n",
" Returns:\n",
" state: A `SEIRComponents` tuple with Tensors of shape [num_particles,\n",
" num_batches, num_cities] giving the initial conditions in each\n",
" city, in each filter particle, in each batch member.\n",
" \"\"\"\n",
" num_cities = mobility_matrix_over_time.shape[-2]\n",
" state_shape = [num_particles, num_batches, num_cities]\n",
" susceptible = initial_population * np.ones(state_shape, dtype=np.float32)\n",
" documented_infectious = np.zeros(state_shape, dtype=np.float32)\n",
" daily_new_documented_infectious = np.zeros(state_shape, dtype=np.float32)\n",
"\n",
" # Following Li et al, initialize Wuhan with up to 2000 people exposed\n",
" # and another up to 2000 undocumented infectious.\n",
" rng = np.random.RandomState(seed[0] % (2**31 - 1))\n",
" wuhan_exposed = rng.randint(\n",
" 0, 2001, [num_particles, num_batches]).astype(np.float32)\n",
" wuhan_undocumented_infectious = rng.randint(\n",
" 0, 2001, [num_particles, num_batches]).astype(np.float32)\n",
" \n",
" # Also following Li et al, initialize cities adjacent to Wuhan with three\n",
" # days' worth of additional exposed and undocumented-infectious cases,\n",
" # as they may have traveled there before the beginning of the modeling\n",
" # period.\n",
" exposed = 3 * mobility_matrix_over_time[\n",
" WUHAN_IDX, :, 0] * wuhan_exposed[\n",
" ..., np.newaxis] / initial_population[WUHAN_IDX]\n",
" undocumented_infectious = 3 * mobility_matrix_over_time[\n",
" WUHAN_IDX, :, 0] * wuhan_undocumented_infectious[\n",
" ..., np.newaxis] / initial_population[WUHAN_IDX]\n",
"\n",
" exposed[..., WUHAN_IDX] = wuhan_exposed\n",
" undocumented_infectious[..., WUHAN_IDX] = wuhan_undocumented_infectious\n",
"\n",
" # Following Li et al, we do not remove the initial exposed and infectious\n",
" # persons from the susceptible population.\n",
" return SEIRComponents(\n",
" susceptible=tf.constant(susceptible),\n",
" exposed=tf.constant(exposed),\n",
" documented_infectious=tf.constant(documented_infectious),\n",
" undocumented_infectious=tf.constant(undocumented_infectious),\n",
" daily_new_documented_infectious=tf.constant(daily_new_documented_infectious))\n",
" \n",
"def initialize_params(num_particles, num_batches, seed):\n",
" \"\"\"Initialize the global parameters for the entire inference run.\n",
"\n",
" Args:\n",
" num_particles: `int` giving the number of particles for the EAKF.\n",
" num_batches: `int` giving the number of independent EAKF runs to\n",
" initialize in a vectorized batch.\n",
" seed: PRNG entropy.\n",
" \n",
" Returns:\n",
" params: A `ModelParams` tuple with fields Tensors of shape\n",
" [num_particles, num_batches] giving the global parameters\n",
" to use for the first batch of EAKF runs.\n",
" \"\"\"\n",
" # We have 6 parameters. We'll initialize with a Sobol sequence,\n",
" # covering the hyper-rectangle defined by our parameter limits.\n",
" halton_sequence = tfp.mcmc.sample_halton_sequence(\n",
" dim=6, num_results=num_particles * num_batches, seed=seed)\n",
" halton_sequence = tf.reshape(\n",
" halton_sequence, [num_particles, num_batches, 6])\n",
" halton_sequences = tf.nest.pack_sequence_as(\n",
" PARAMETER_LOWER_BOUNDS, tf.split(\n",
" halton_sequence, num_or_size_splits=6, axis=-1))\n",
" def interpolate(minval, maxval, h):\n",
" return (maxval - minval) * h + minval\n",
" return tf.nest.map_structure(\n",
" interpolate,\n",
" PARAMETER_LOWER_BOUNDS, PARAMETER_UPPER_BOUNDS, halton_sequences)\n",
"\n",
"def update_params(num_particles, num_batches,\n",
" prev_params, parameter_variance, seed):\n",
" \"\"\"Update the global parameters between EAKF runs.\n",
"\n",
" Args:\n",
" num_particles: `int` giving the number of particles for the EAKF.\n",
" num_batches: `int` giving the number of independent EAKF runs to\n",
" initialize in a vectorized batch.\n",
" prev_params: A `ModelParams` tuple of the parameters used for the previous\n",
" EAKF run.\n",
" parameter_variance: A `ModelParams` tuple specifying how much to drift\n",
" each parameter.\n",
" seed: PRNG entropy.\n",
" \n",
" Returns:\n",
" params: A `ModelParams` tuple with fields Tensors of shape\n",
" [num_particles, num_batches] giving the global parameters\n",
" to use for the next batch of EAKF runs.\n",
" \"\"\"\n",
" # Initialize near the previous set of parameters. This is the first step\n",
" # in Iterated Filtering.\n",
" seeds = tf.nest.pack_sequence_as(\n",
" prev_params, samplers.split_seed(seed, n=len(prev_params)))\n",
" return tf.nest.map_structure(\n",
" lambda x, v, seed: x + tf.math.sqrt(v) * tf.random.stateless_normal([\n",
" num_particles, num_batches, 1], seed=seed),\n",
" prev_params, parameter_variance, seeds)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "jKWm64gF0JkZ"
},
"source": [
"## Delays\n",
"\n",
"One of the important features of this model is taking explicit account of the fact that infections are reported later than they begin. That is, we expect that a person who moves from the $E$ compartment to the $I^r$ compartment on day $t$ may not show up in the observable reported case counts until some later day.\n",
"\n",
"We assume the delay is gamma-distributed. Following Li et al, we use 1.85 for the shape, and parameterize the rate to produce an average reporting delay of 9 days."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "rFNy6L491BFD"
},
"outputs": [],
"source": [
"def raw_reporting_delay_distribution(gamma_shape=1.85, reporting_delay=9.):\n",
" return tfp.distributions.Gamma(\n",
" concentration=gamma_shape, rate=gamma_shape / reporting_delay)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "yQTwUeKr8dVP"
},
"source": [
"Our observations are discrete, so we will round the raw (continuous) delays up to the nearest day. We also have a finite data horizon, so the delay distribution for a single person is a categorical over the remaining days. We can therefore compute the per-city predicted observations more efficiently than sampling $O(I^r)$ gammas, by pre-computing multinomial delay probabilities instead."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "0YRv0FjX9DfB"
},
"outputs": [],
"source": [
"def reporting_delay_probs(num_timesteps, gamma_shape=1.85, reporting_delay=9.):\n",
" gamma_dist = raw_reporting_delay_distribution(gamma_shape, reporting_delay)\n",
" multinomial_probs = [gamma_dist.cdf(1.)]\n",
" for k in range(2, num_timesteps + 1):\n",
" multinomial_probs.append(gamma_dist.cdf(k) - gamma_dist.cdf(k - 1))\n",
" # For samples that are larger than T.\n",
" multinomial_probs.append(gamma_dist.survival_function(num_timesteps))\n",
" multinomial_probs = tf.stack(multinomial_probs)\n",
" return multinomial_probs"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "7WAHi68lBFqz"
},
"source": [
"Here's the code for actually applying these delays to the new daily documented infectious counts:"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "vrPHDvUEBPsS"
},
"outputs": [],
"source": [
"def delay_reporting(\n",
" daily_new_documented_infectious, num_timesteps, t, multinomial_probs, seed):\n",
" # This is the distribution of observed infectious counts from the current\n",
" # timestep.\n",
"\n",
" raw_delays = tfd.Multinomial(\n",
" total_count=daily_new_documented_infectious,\n",
" probs=multinomial_probs).sample(seed=seed)\n",
"\n",
" # The last bucket is used for samples that are out of range of T + 1. Thus\n",
" # they are not going to be observable in this model.\n",
" clipped_delays = raw_delays[..., :-1]\n",
"\n",
" # We can also remove counts that are such that t + i >= T.\n",
" clipped_delays = clipped_delays[..., :num_timesteps - t]\n",
" # We finally shift everything by t. That means prepending with zeros.\n",
" return tf.concat([\n",
" tf.zeros(\n",
" tf.concat([\n",
" tf.shape(clipped_delays)[:-1], [t]], axis=0),\n",
" dtype=clipped_delays.dtype),\n",
" clipped_delays], axis=-1)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "DRsd3wsW9WYQ"
},
"source": [
"## Inference"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "2lO5MVT6M2lc"
},
"source": [
"First we'll define some data structures for inference.\n",
"\n",
"In particular, we'll be wanting to do Iterated Filtering, which packages the\n",
"state and parameters together while doing inference. So we'll define\n",
"a `ParameterStatePair` object.\n",
"\n",
"We also want to package any side information to the model."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "6JrCjAhXM6TL"
},
"outputs": [],
"source": [
"ParameterStatePair = collections.namedtuple(\n",
" 'ParameterStatePair', ['state', 'params'])\n",
"\n",
"# Info that is tracked and mutated but should not have inference performed over.\n",
"SideInfo = collections.namedtuple(\n",
" 'SideInfo', [\n",
" # Observations at every time step.\n",
" 'observations_over_time',\n",
" 'initial_population',\n",
" 'mobility_matrix_over_time',\n",
" 'population',\n",
" # Used for variance of measured observations.\n",
" 'actual_reported_cases',\n",
" # Pre-computed buckets for the multinomial distribution.\n",
" 'multinomial_probs',\n",
" 'seed',\n",
" ])\n",
"\n",
"# Cities can not fall below this fraction of people\n",
"MINIMUM_CITY_FRACTION = 0.6\n",
"\n",
"# How much to inflate the covariance by.\n",
"INFLATION_FACTOR = 1.1\n",
"\n",
"INFLATE_FN = tfes.inflate_by_scaled_identity_fn(INFLATION_FACTOR)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "eoZPQ9YDM9mj"
},
"source": [
"Here is the complete observation model, packaged for the Ensemble Kalman Filter.\n",
"\n",
"The interesting feature is the reporting delays (computed as previously). The upstream model emits the `daily_new_documented_infectious` for each city at each time step."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "lBrt1BkO9Xje"
},
"outputs": [],
"source": [
"# We observe the observed infections.\n",
"def observation_fn(t, state_params, extra):\n",
" \"\"\"Generate reported cases.\n",
" \n",
" Args:\n",
" state_params: A `ParameterStatePair` giving the current parameters\n",
" and state.\n",
" t: Integer giving the current time.\n",
" extra: A `SideInfo` carrying auxiliary information.\n",
"\n",
" Returns:\n",
" observations: A Tensor of predicted observables, namely new cases\n",
" per city at time `t`.\n",
" extra: Update `SideInfo`.\n",
" \"\"\"\n",
" # Undo padding introduced in `inference`.\n",
" daily_new_documented_infectious = state_params.state.daily_new_documented_infectious[..., 0]\n",
" # Number of people that we have already committed to become\n",
" # observed infectious over time.\n",
" # shape: batch + [num_particles, num_cities, time]\n",
" observations_over_time = extra.observations_over_time\n",
" num_timesteps = observations_over_time.shape[-1]\n",
"\n",
" seed, new_seed = samplers.split_seed(extra.seed, salt='reporting delay')\n",
" \n",
" daily_delayed_counts = delay_reporting(\n",
" daily_new_documented_infectious, num_timesteps, t,\n",
" extra.multinomial_probs, seed)\n",
" observations_over_time = observations_over_time + daily_delayed_counts\n",
"\n",
" extra = extra._replace(\n",
" observations_over_time=observations_over_time,\n",
" seed=new_seed)\n",
"\n",
" # Actual predicted new cases, re-padded.\n",
" adjusted_observations = observations_over_time[..., t][..., tf.newaxis]\n",
" # Finally observations have variance that is a function of the true observations:\n",
" return tfd.MultivariateNormalDiag(\n",
" loc=adjusted_observations,\n",
" scale_diag=tf.math.maximum(\n",
" 2., extra.actual_reported_cases[..., t][..., tf.newaxis] / 2.)), extra"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "ugNrSlDAyUc6"
},
"source": [
"Here we define the transition dynamics. We've done the semantic work already; here we just package it for the EAKF framework, and, following Li et al, clip city populations to prevent them from getting too small."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "95S7LAj4NB_5"
},
"outputs": [],
"source": [
"def transition_fn(t, state_params, extra):\n",
" \"\"\"SEIR dynamics.\n",
"\n",
" Args:\n",
" state_params: A `ParameterStatePair` giving the current parameters\n",
" and state.\n",
" t: Integer giving the current time.\n",
" extra: A `SideInfo` carrying auxiliary information.\n",
"\n",
" Returns:\n",
" state_params: A `ParameterStatePair` predicted for the next time step.\n",
" extra: Updated `SideInfo`.\n",
" \"\"\"\n",
" mobility_t = extra.mobility_matrix_over_time[..., t]\n",
" new_seed, rk4_seed = samplers.split_seed(extra.seed, salt='Transition')\n",
" new_state = rk4_one_step(\n",
" state_params.state,\n",
" extra.population,\n",
" mobility_t,\n",
" state_params.params,\n",
" seed=rk4_seed)\n",
"\n",
" # Make sure population doesn't go below MINIMUM_CITY_FRACTION.\n",
" new_population = (\n",
" extra.population + state_params.params.intercity_underreporting_factor * (\n",
" # Inflow\n",
" tf.reduce_sum(mobility_t, axis=-2) - \n",
" # Outflow\n",
" tf.reduce_sum(mobility_t, axis=-1)))\n",
" new_population = tf.where(\n",
" new_population < MINIMUM_CITY_FRACTION * extra.initial_population,\n",
" extra.initial_population * MINIMUM_CITY_FRACTION,\n",
" new_population)\n",
"\n",
" extra = extra._replace(population=new_population, seed=new_seed)\n",
"\n",
" # The Ensemble Kalman Filter code expects the transition function to return a distribution.\n",
" # As the dynamics and noise are encapsulated above, we construct a `JointDistribution` that when\n",
" # sampled, returns the values above.\n",
"\n",
" new_state = tfd.JointDistributionNamed(\n",
" model=tf.nest.map_structure(lambda x: tfd.VectorDeterministic(x), new_state))\n",
" params = tfd.JointDistributionNamed(\n",
" model=tf.nest.map_structure(lambda x: tfd.VectorDeterministic(x), state_params.params))\n",
" \n",
" state_params = tfd.JointDistributionNamed(\n",
" model=ParameterStatePair(state=new_state, params=params))\n",
"\n",
" return state_params, extra"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "nsdNsvzuyaN-"
},
"source": [
"Finally we define the inference method. This is two loops, the outer loop\n",
"being Iterated Filtering while the inner loop is Ensemble Adjustment Kalman Filtering."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "EUWinMkUNFGJ"
},
"outputs": [],
"source": [
"# Use tf.function to speed up EAKF prediction and updates.\n",
"ensemble_kalman_filter_predict = tf.function(\n",
" tfes.ensemble_kalman_filter_predict, autograph=False)\n",
"ensemble_adjustment_kalman_filter_update = tf.function(\n",
" tfes.ensemble_adjustment_kalman_filter_update, autograph=False)\n",
"\n",
"def inference(\n",
" num_ensembles,\n",
" num_batches,\n",
" num_iterations,\n",
" actual_reported_cases,\n",
" mobility_matrix_over_time,\n",
" seed=None,\n",
" # This is how much to reduce the variance by in every iterative\n",
" # filtering step.\n",
" variance_shrinkage_factor=0.9,\n",
" # Days before infection is reported.\n",
" reporting_delay=9.,\n",
" # Shape parameter of Gamma distribution.\n",
" gamma_shape_parameter=1.85):\n",
" \"\"\"Inference for the Shaman, et al. model.\n",
"\n",
" Args:\n",
" num_ensembles: Number of particles to use for EAKF.\n",
" num_batches: Number of batches of IF-EAKF to run.\n",
" num_iterations: Number of iterations to run iterative filtering.\n",
" actual_reported_cases: `Tensor` of shape `[L, T]` where `L` is the number\n",
" of cities, and `T` is the timesteps.\n",
" mobility_matrix_over_time: `Tensor` of shape `[L, L, T]` which specifies the\n",
" mobility between locations over time.\n",
" variance_shrinkage_factor: Python `float`. How much to reduce the\n",
" variance each iteration of iterated filtering.\n",
" reporting_delay: Python `float`. How many days before the infection\n",
" is reported.\n",
" gamma_shape_parameter: Python `float`. Shape parameter of Gamma distribution\n",
" of reporting delays.\n",
"\n",
" Returns:\n",
" result: A `ModelParams` with fields Tensors of shape [num_batches],\n",
" containing the inferred parameters at the final iteration.\n",
" \"\"\"\n",
" print('Starting inference.')\n",
" num_timesteps = actual_reported_cases.shape[-1]\n",
" params_per_iter = []\n",
"\n",
" multinomial_probs = reporting_delay_probs(\n",
" num_timesteps, gamma_shape_parameter, reporting_delay)\n",
"\n",
" seed = samplers.sanitize_seed(seed, salt='Inference')\n",
"\n",
" for i in range(num_iterations):\n",
" start_if_time = time.time()\n",
" seeds = samplers.split_seed(seed, n=4, salt='Initialize')\n",
" if params_per_iter:\n",
" parameter_variance = tf.nest.map_structure(\n",
" lambda minval, maxval: variance_shrinkage_factor ** (\n",
" 2 * i) * (maxval - minval) ** 2 / 4.,\n",
" PARAMETER_LOWER_BOUNDS, PARAMETER_UPPER_BOUNDS)\n",
" params_t = update_params(\n",
" num_ensembles,\n",
" num_batches,\n",
" prev_params=params_per_iter[-1],\n",
" parameter_variance=parameter_variance,\n",
" seed=seeds.pop())\n",
" else:\n",
" params_t = initialize_params(num_ensembles, num_batches, seed=seeds.pop())\n",
"\n",
" state_t = initialize_state(num_ensembles, num_batches, seed=seeds.pop())\n",
" population_t = sum(x for x in state_t)\n",
" observations_over_time = tf.zeros(\n",
" [num_ensembles,\n",
" num_batches,\n",
" actual_reported_cases.shape[0], num_timesteps])\n",
"\n",
" extra = SideInfo(\n",
" observations_over_time=observations_over_time,\n",
" initial_population=tf.identity(population_t),\n",
" mobility_matrix_over_time=mobility_matrix_over_time,\n",
" population=population_t,\n",
" multinomial_probs=multinomial_probs,\n",
" actual_reported_cases=actual_reported_cases,\n",
" seed=seeds.pop())\n",
"\n",
" # Clip states\n",
" state_t = clip_state(state_t, population_t)\n",
" params_t = clip_params(params_t, seed=seeds.pop())\n",
"\n",
" # Accrue the parameter over time. We'll be averaging that\n",
" # and using that as our MLE estimate.\n",
" params_over_time = tf.nest.map_structure(\n",
" lambda x: tf.identity(x), params_t)\n",
"\n",
" state_params = ParameterStatePair(state=state_t, params=params_t)\n",
"\n",
" eakf_state = tfes.EnsembleKalmanFilterState(\n",
" step=tf.constant(0), particles=state_params, extra=extra)\n",
"\n",
" for j in range(num_timesteps):\n",
" seeds = samplers.split_seed(eakf_state.extra.seed, n=3)\n",
" \n",
" extra = extra._replace(seed=seeds.pop())\n",
" \n",
" # Predict step.\n",
"\n",
" # Inflate and clip.\n",
" new_particles = INFLATE_FN(eakf_state.particles)\n",
" state_t = clip_state(new_particles.state, eakf_state.extra.population)\n",
" params_t = clip_params(new_particles.params, seed=seeds.pop())\n",
" eakf_state = eakf_state._replace(\n",
" particles=ParameterStatePair(params=params_t, state=state_t))\n",
"\n",
" eakf_predict_state = ensemble_kalman_filter_predict(eakf_state, transition_fn)\n",
"\n",
" # Clip the state and particles.\n",
" state_params = eakf_predict_state.particles\n",
" state_t = clip_state(\n",
" state_params.state, eakf_predict_state.extra.population)\n",
" state_params = ParameterStatePair(state=state_t, params=state_params.params)\n",
"\n",
" # We preprocess the state and parameters by affixing a 1 dimension. This is because for\n",
" # inference, we treat each city as independent. We could also introduce localization by\n",
" # considering cities that are adjacent.\n",
" state_params = tf.nest.map_structure(lambda x: x[..., tf.newaxis], state_params)\n",
" eakf_predict_state = eakf_predict_state._replace(particles=state_params)\n",
"\n",
" # Update step.\n",
" \n",
" eakf_update_state = ensemble_adjustment_kalman_filter_update(\n",
" eakf_predict_state,\n",
" actual_reported_cases[..., j][..., tf.newaxis],\n",
" observation_fn)\n",
" \n",
" state_params = tf.nest.map_structure(\n",
" lambda x: x[..., 0], eakf_update_state.particles)\n",
"\n",
" # Clip to ensure parameters / state are well constrained.\n",
" state_t = clip_state(\n",
" state_params.state, eakf_update_state.extra.population)\n",
" \n",
" # Finally for the parameters, we should reduce over all updates. We get\n",
" # an extra dimension back so let's do that.\n",
" params_t = tf.nest.map_structure(\n",
" lambda x, y: x + tf.reduce_sum(y[..., tf.newaxis] - x, axis=-2, keepdims=True),\n",
" eakf_predict_state.particles.params, state_params.params)\n",
" params_t = clip_params(params_t, seed=seeds.pop())\n",
" params_t = tf.nest.map_structure(lambda x: x[..., 0], params_t)\n",
"\n",
" state_params = ParameterStatePair(state=state_t, params=params_t)\n",
" eakf_state = eakf_update_state\n",
" eakf_state = eakf_state._replace(particles=state_params)\n",
"\n",
" # Flatten and collect the inferred parameter at time step t.\n",
" params_over_time = tf.nest.map_structure(\n",
" lambda s, x: tf.concat([s, x], axis=-1), params_over_time, params_t)\n",
"\n",
" est_params = tf.nest.map_structure(\n",
" # Take the average over the Ensemble and over time.\n",
" lambda x: tf.math.reduce_mean(x, axis=[0, -1])[..., tf.newaxis],\n",
" params_over_time)\n",
" params_per_iter.append(est_params)\n",
" print('Iterated Filtering {} / {} Ran in: {:.2f} seconds'.format(\n",
" i, num_iterations, time.time() - start_if_time))\n",
"\n",
" return tf.nest.map_structure(\n",
" lambda x: tf.squeeze(x, axis=-1), params_per_iter[-1])"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "UvALK7e1fgKU"
},
"source": [
"Final detail: clipping the parameters and state consists of making sure they are within range, and non-negative."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "aQj7o0nU6-ru"
},
"outputs": [],
"source": [
"def clip_state(state, population):\n",
" \"\"\"Clip state to sensible values.\"\"\"\n",
" state = tf.nest.map_structure(\n",
" lambda x: tf.where(x < 0, 0., x), state)\n",
"\n",
" # If S > population, then adjust as well.\n",
" susceptible = tf.where(state.susceptible > population, population, state.susceptible)\n",
" return SEIRComponents(\n",
" susceptible=susceptible,\n",
" exposed=state.exposed,\n",
" documented_infectious=state.documented_infectious,\n",
" undocumented_infectious=state.undocumented_infectious,\n",
" daily_new_documented_infectious=state.daily_new_documented_infectious)\n",
"\n",
"def clip_params(params, seed):\n",
" \"\"\"Clip parameters to bounds.\"\"\"\n",
" def _clip(p, minval, maxval):\n",
" return tf.where(\n",
" p < minval,\n",
" minval * (1. + 0.1 * tf.random.stateless_uniform(p.shape, seed=seed)),\n",
" tf.where(p > maxval,\n",
" maxval * (1. - 0.1 * tf.random.stateless_uniform(\n",
" p.shape, seed=seed)), p))\n",
" params = tf.nest.map_structure(\n",
" _clip, params, PARAMETER_LOWER_BOUNDS, PARAMETER_UPPER_BOUNDS)\n",
"\n",
" return params"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "c7ACJedFAyHQ"
},
"source": [
"## Running it all together"
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {
"height": 207
},
"colab_type": "code",
"id": "-T603NTSA0f3",
"outputId": "690136df-c409-4b1b-fa06-97e46ced19ba"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Starting inference.\n",
"Iterated Filtering 0 / 10 Ran in: 26.65 seconds\n",
"Iterated Filtering 1 / 10 Ran in: 28.69 seconds\n",
"Iterated Filtering 2 / 10 Ran in: 28.06 seconds\n",
"Iterated Filtering 3 / 10 Ran in: 28.48 seconds\n",
"Iterated Filtering 4 / 10 Ran in: 28.57 seconds\n",
"Iterated Filtering 5 / 10 Ran in: 28.35 seconds\n",
"Iterated Filtering 6 / 10 Ran in: 28.35 seconds\n",
"Iterated Filtering 7 / 10 Ran in: 28.19 seconds\n",
"Iterated Filtering 8 / 10 Ran in: 28.58 seconds\n",
"Iterated Filtering 9 / 10 Ran in: 28.23 seconds\n"
]
}
],
"source": [
"# Let's sample the parameters.\n",
"#\n",
"# NOTE: Li et al. run inference 1000 times, which would take a few hours.\n",
"# Here we run inference 30 times (in a single, vectorized batch).\n",
"best_parameters = inference(\n",
" num_ensembles=300,\n",
" num_batches=30,\n",
" num_iterations=10,\n",
" actual_reported_cases=observed_daily_infectious_count,\n",
" mobility_matrix_over_time=mobility_matrix_over_time)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "ZWjY1KpJ3JuT"
},
"source": [
"The results of our inferences. We plot the maximum-likelihood values for all the global paramters to show their variation across our `num_batches` independent runs of inference. This corresponds to Table S1 in the supplemental materials."
]
},
{
"cell_type": "code",
"execution_count": 0,
"metadata": {
"colab": {
"height": 298
},
"colab_type": "code",
"id": "TZqHM3nIQ8Zj",
"outputId": "b2d7e3f7-9ce6-401a-99b1-6d2a895c5f6d"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAEZCAYAAAAt5touAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X1QVOfdPvDryKrRoMgC2wTXp0io\noihsdFEkhpqkvG0IddVxwKbKMGaDL2mdZBImMzWptJlqmpmk0SR0GxJHayQzpkqjwGAqmGgTcEde\nJERLHFBB5cVINETDsvD7w3F/rizsgruce3evzwzz7HnD7+a5ey7Ofe5zH6m/v78fREREAhojdwFE\nRESDYUgREZGwGFJERCQshhQREQmLIUVERMJiSBERkbAYUkREJCyGFBERCYsh5WJNTU1ITU1FYGAg\npk6dig8//FDukojIi1RXV+ORRx7BxIkTsWDBApw/f17uktyKIeViK1asQGJiIjo7O/GPf/wDf/7z\nn+UuiYi8REtLC3Q6HXJzc3HlyhWEh4d7/TmGIeVCdXV1uHLlCp5//nn4+fkBAEJCQmSuikT12muv\nYd26ddblq1evYuzYsbh586aMVZHIXnjhBTzzzDNIT0/HhAkTkJGRgRMnTshdllsp5C7Amxw/fhyL\nFy9GX18fqqur8fzzz+Pll1+WuywS1KlTp/DYY49Zl2tqajBz5kzcd999MlZForp27RqKiorwv//9\nz7qur6/P69sLr6RcqKamBlqtFo899hi0Wi0mTpyIZcuWyV0WCerUqVPQaDTW5ZqaGsTExMhYEYns\nP//5D8xmM6KjozFlyhRMmTIFv/nNb/Dzn/9c7tLciiHlQjU1NYiNjUV5eTm+/fZbKJVKvPTSS3KX\nRQLq6enB2bNnMXfuXOu62tpam9AiulNzczPS09PR1dVl/XnssceQkpIid2luxZByEYvFgm+++QYP\nP/wwxowZg4ceegiPPPKI3GWRoBoaGjB16lRMnDgRANDf34+KigpeSdGgfvrpJ2t7AW6NJDaZTEhP\nT5exKvdjSLnImTNn8OOPP6KkpAQWiwU1NTUoKCjAmjVr5C6NBHTq1Cm0t7fj7NmzuHHjBjZv3oxz\n584hLCxM7tJIULGxsTh69CguXryICxcuYNWqVXjttdegVCrlLs2tGFIuUl1djdmzZ+OFF17AlClT\nkJWVhbfffhtxcXFyl0YCOnXqFJKTk5GamoqIiAj87Gc/Q3h4OF577TW5SyNBPf7443jqqacwY8YM\nLF68GL/97W/xzDPPyF2W20l8M69rvPjii1AqlRzNR05JTU3F2rVrsXz5crlLIRIar6RcpLq6GrNm\nzZK7DPIQp06dYnshcgJDykVqa2sRGRkpdxnkAa5evYr29nb84he/kLsUIuGxu4+IiITFKykiIhIW\nQ4qIiIQl5Nx9wcHBfF5EZs3Nzejs7JS7DKewvYjBU9oM24sYnG0vQoZUWFgYTCaT3GX4NK1WK3cJ\nTmN7EYOntBm2FzE4217Y3UdERMJiSBERkbAYUi4iSdKQP76ktLQUM2fOREREBLZu3Tpge0VFBQIC\nAqDRaKDRaJCXl2fdFhYWhrlz50Kj0XhM99FIOGovvtZmyDFfbS9C3pPyRHc+biZJEnz18TOLxYIN\nGzbg8OHDUKvViI2NRXp6OmbPnm2z36OPPoqDBw/a/R3l5eUIDg4ejXJlc3f78OU2Q87x1XMMr6TI\npaqqqhAREYHw8HCMGzcOGRkZKCoqkrssIvJQDClyqdbWVkybNs26rFar0draOmC/L7/8EjExMUhN\nTcXXX39tXS9JEpKSkjB//nwYjcZRqZmIxMXuPnIpe10Qd/eXz5s3D+fOnYO/vz+Ki4uxdOlSNDY2\nAgCOHz+O0NBQtLe3IzExEZGRkUhISBjwO41GozXEOjo63PBNiEgEvJIil1Kr1bhw4YJ1uaWlBaGh\noTb7TJ48Gf7+/gAAnU4Hs9lsfajv9r4qlQp6vR5VVVV2/x2DwQCTyQSTyYSQkBB3fBUaJdnZ2VCp\nVJgzZ86Q+504cQJ+fn7Yt2+fdZ2jQTrk+RhS5FKxsbFobGxEU1MTenp6UFhYOOD11pcvX7ZecVVV\nVaGvrw9BQUHo7u7G9evXAQDd3d0oKytzeOIiz5eVlYXS0tIh97FYLMjNzUVycrLNug0bNqCkpAQN\nDQ3Yu3cvGhoa3F0ujTJ295FLKRQK7NixA8nJybBYLMjOzkZUVBTy8/MBADk5Odi3bx/ee+89KBQK\nTJgwAYWFhZAkCW1tbdDr9QCA3t5erFq1CikpKXJ+HRoFCQkJaG5uHnKf7du3Y/ny5Thx4oR13Z2D\ndABYB+ncPZKUPBtDilxOp9NBp9PZrMvJybF+3rhxIzZu3DjguPDwcNTW1rq9PvIsra2t2L9/P44c\nOWITUvYG6VRWVspRIrmRw+4+R/3Fp0+fxqJFizB+/Hi88cYbNtvYX0xE92rTpk3Ytm0b/Pz8bNY7\nM0jnNqPRCK1WC61Wy4E2HsbhlVRWVhY2btyI1atX292uVCrx9ttv48CBAzbrnX2ok4hoKCaTCRkZ\nGQCAzs5OFBcXQ6FQODVI5zaDwQCDwQDAcybCpVschpSj/mKVSgWVSoVDhw7ZrGd/MRG5QlNTk/Vz\nVlYW0tLSsHTpUvT29loH6UydOhWFhYX46KOPZKyU3MFt96TYX0xEzsjMzERFRQU6OzuhVquxZcsW\nmM1mALb3Mu822CAd8i5uC6nh9BcDfDiTyFft3bvX6X137txps2xvkA55F7c9JzWc/mLA8x7OVCqV\nQ85GPNRsxUqlUubqiYg8g9tCypmHOj3Z1atX0d/fP6Kfq1evyl0+EZFHcNjd56i/+PLly9Bqtbh2\n7RrGjBmDt956Cw0NDZg8eTL7i4mI6J44DClH/cUPPPAAWlpa7G5jfzEREd0Lzt1HRETCYkgREZGw\nGFJERCQshhQREQmLIUVERMJiSBERkbAYUkREJCyGFBERCYshRUREwmJIEREJiJNY3+K2V3UQkS2l\nUjnk5MKDvcomMDAQ3333nbvKIkHdnsR6JIZ6LZKnYUgRjZKRnnS86YRDNFzs7iMiImExpIiISFgM\nKXK50tJSzJw5ExEREdi6deuA7RUVFQgICIBGo4FGo0FeXp7TxxKRb2FIkUtZLBZs2LABJSUlaGho\nwN69e9HQ0DBgv0cffRQ1NTWoqanBK6+8MqxjybtkZ2dDpVJhzpw5drcXFRUhOjoaGo0GWq0Wx44d\ns24LCwvD3LlzrdvI+zCkyKWqqqoQERGB8PBwjBs3DhkZGSgqKnL7seS5srKyUFpaOuj2J554ArW1\ntaipqcEHH3yAtWvX2mwvLy9HTU0NTCaTu0slGXB03wj1vzoZ+GPAyI/1Uq2trZg2bZp1Wa1Wo7Ky\ncsB+X375JWJiYhAaGoo33ngDUVFRTh9L3iUhIQHNzc2Dbvf397d+7u7u5mhHH8OQGiFpy7V7eoah\n/4+urUcU9v6b3H1SmTdvHs6dOwd/f38UFxdj6dKlaGxsdOrY24xGI4xGIwCgo6PDBZWTyPbv34+X\nX34Z7e3tOHTokHW9JElISkqCJEl49tlnYTAY7B7P9uK52N1HLqVWq3HhwgXrcktLC0JDQ232mTx5\nsvWvY51OB7PZjM7OTqeOvc1gMMBkMsFkMiEkJMQN34REotfrcfr0aRw4cACbN2+2rj9+/DhOnjyJ\nkpISvPPOO/j888/tHs/24rkYUuRSsbGxaGxsRFNTE3p6elBYWIj09HSbfS5fvmy9aqqqqkJfXx+C\ngoKcOpZ8W0JCAs6ePYvOzk4AsP4Ro1KpoNfrUVVVJWd55Abs7iOXUigU2LFjB5KTk2GxWJCdnY2o\nqCjk5+cDAHJycrBv3z689957UCgUmDBhAgoLCyFJ0qDHkm/79ttv8dBDD0GSJJw8eRI9PT0ICgpC\nd3c3+vr6MGnSJHR3d6OsrMw6UpS8B0OKXE6n00Gn09msy8nJsX7euHEjNm7c6PSx5N0yMzNRUVFh\n7fLdsmULzGYzgFvt5pNPPsGuXbswduxYTJgwAR9//DEkSUJbWxv0ej0AoLe3F6tWrUJKSoqcX4Xc\ngCFFRLLau3fvkNtzc3ORm5s7YH14eDhqa2vdVRYJgvekiIhIWAwpIiISFkOKiIiExZAiIiJhMaSI\niEhYDCkiIhIWQ4qIiITF56SIiATENy3cwpAiIhIQ37RwC7v7iIhIWAwpIiISlsPuvuzsbBw8eBAq\nlQr19fUDtvf39+P3v/89iouLMXHiROzcuRPz5s0DAISFhWHSpEnw8/ODQqHg653Jp430HoM33V8g\nGi6HIZWVlYWNGzdi9erVdreXlJSgsbERjY2NqKysxLp162xe+V1eXo7g4GDXVUzkoUZ6j8Gb7i8Q\nDZfD7r6EhAQolcpBtxcVFWH16tWQJAlxcXHo6urCpUuXXFokERH5pnu+J9Xa2opp06ZZl9VqNVpb\nWwHc+gswKSkJ8+fPh9FovNd/ioiIfMw9D0G3130hSRIA4Pjx4wgNDUV7ezsSExMRGRmJhIQEu7/H\naDRag6yjo+NeyyIiIi9wz1dSarUaFy5csC63tLQgNDQUAKz/V6VSQa/Xo6qqatDfYzAYYDKZYDKZ\nEBIScq9lERGRF7jnkEpPT8euXbvQ39+Pr776CgEBAXjwwQfR3d2N69evAwC6u7tRVlaGOXPm3HPB\nRETkOxyGVGZmJhYtWoQzZ85ArVajoKAA+fn5yM/PBwDodDqEh4cjIiICzzzzDN59910AQFtbGxYv\nXoyYmBgsWLAATz75JFJSUtz7bYjI42RnZ0OlUg36R2xRURGio6Oh0Wig1Wpx7Ngx67bS0lLMnDkT\nERER2Lp162iVTKPI4T2pvXv3DrldkiS88847A9aHh4ejtrZ25JURkU9w9JjLE088gfT0dEiShLq6\nOqxcuRKnT5+GxWLBhg0bcPjwYajVasTGxiI9PR2zZ88e5W9A7sQZJ4hIVo4ec/H397cOxuru7rZ+\nrqqqQkREBMLDwzFu3DhkZGSgqKhoVGqm0cOQIiLh7d+/H5GRkXjyySfxwQcfABj68Ze7GY1GaLVa\naLVajh72MAwpcjln7xOcOHECfn5+2Ldvn3VdWFgY5s6da73/QAQAer0ep0+fxoEDB7B582YAQz/+\ncjeOHvZcfFUHuZSz9wksFgtyc3ORnJw84Hd481Rag51EhxIYGOiGSjxTQkICzp49i87OziEffyHv\nwSspciln7xNs374dy5cvh0qlkqFKefT39w/6M9T27777TubK5fXtt99a/xudPHkSPT09CAoKQmxs\nLBobG9HU1ISenh4UFhYiPT1d5mrJ1XgldQ9G8lcx4N1/Gdu7T3DnhMO399m/fz+OHDmCEydO2Gy7\nPZWWJEl49tlnYTAYRqVukk9mZiYqKiqsV0dbtmyB2WwGAOTk5OCTTz7Brl27MHbsWEyYMAEff/wx\nJEmCQqHAjh07kJycDIvFguzsbERFRcn8bcjVGFIjNNRs1pIkjfiNmp7OmfsEmzZtwrZt2+Dn5zdg\nX2en0uI0Wt7D0WMuubm5yM3NtbtNp9NBp9O5oywSBEOKXMqZ+wQmkwkZGRkAgM7OThQXF0OhUGDp\n0qV2p9KyF1IGg8F6lcUBFuSt2FvDe1LkYs7cJ2hqakJzczOam5uxYsUKvPvuu1i6dCmn0iK6w0jv\nYXrbfUxeSZFLDXaf4PY0Wjk5OYMe29bWBr1eDwDo7e3FqlWrOJUWkY9jSJHL2btPMFg47dy50/qZ\nU2kR0d3Y3UdERMJiSBERkbAYUkREJCyGFBERCYshRUREwmJIERGRsBhSREQkLIYUEREJiyFFRETC\nYkgREZGwGFJERCQshhQREQmLIUVERMLiLOgucvfLye5e9tU39ZJ99l5mxzZDQ/HVcwyvpFxkqBeQ\neWvjoZFz1F58qc1kZ2dDpVIN+oLLPXv2IDo6GtHR0YiPj7d5nUtYWBjmzp0LjUbj9W9o9tX2wpAi\nIlllZWWhtLR00O3Tp0/H0aNHUVdXh82bN8NgMNhsLy8vR01NDUwmk7tLJRmwu4+IZJWQkIDm5uZB\nt8fHx1s/x8XFoaWlZRSqIlHwSoqIPEZBQQFSU1Oty5IkISkpCfPnz4fRaBz0OKPRCK1WC61Wi46O\njtEolVxE6hewMzM4OBhhYWFylzFiHR0dCAkJkbuMe9Lc3IzOzk65y3CKp7cXgG2mubkZaWlpqK+v\nH3Sf8vJyrF+/HseOHUNQUBAA4OLFiwgNDUV7ezsSExOxfft2JCQkDPlvsb2Iwdn2ImR3n6ecHAej\n1WrZPz6KPL29AGwzjtTV1WHt2rUoKSmxBhQAhIaGAgBUKhX0ej2qqqochhTbi2dhdx8RCe38+fNY\ntmwZdu/ejRkzZljXd3d34/r169bPZWVlg44QJM8l5JUUEfmOzMxMVFRUoLOzE2q1Glu2bIHZbAYA\n5OTkIC8vD1euXMH69esBAAqFAiaTCW1tbdDr9QCA3t5erFq1CikpKbJ9D3IPIe9JeTqj0ThgmCzR\nUNhmaDh8qb0wpIiISFi8J0VERMJiSLmQo+ldiO7GNkPD4YvthSHlQo6mdyG6G9sMDYcvtheGlAsl\nJCRAqVTKXQZ5ELYZGg5fbC8MKSIiEhZDioiIhMWQIiIiYTGkiIhIWAwpF8rMzMSiRYtw5swZqNVq\nFBQUyF0SCY5thobDF9sLZ5wgIiJh8UqKiIiExZAiIiJhMaSIiEhYDCkiIhIWQ4qIiITFkCIiImEx\npIiISFgMKSIiEhZDioiIhMWQIiIiYTGkiIhIWAwpIiISFkPKRfLy8uDv72/zc99990GSJHz88cdy\nl0dEHu7q1auQJMl6fvm///s/rFixAvX19XKX5lYMKRd55ZVX8MMPP1h/Ll++DI1GA51Oh2XLlsld\nHhF5uJqaGiiVSus5prq6GjExMVi4cCFOnz4td3luw5Bygxs3biAtLQ33338/PvnkE4wdO1bukkhA\nvb29+NOf/oSwsDAEBQXho48+wuuvv47XXntN7tJIQDU1NdBoNNbloKAgbN68GfPmzfPq90op5C7A\n2/T09GDZsmXo6elBWVkZ7rvvPrlLIkH94Q9/gMlkQm1tLT7//HO89NJLkCQJlZWVcpdGAqqurrYJ\nqdsiIyPR2toqQ0Wjg1dSLtTb24uMjAy0tbWhpKQE/v7+cpdEgrp27RreeustGI1GBAQEWLtsnn76\naUyaNEnu8khAd19J3fb9998jJCREhopGB0PKRfr6+pCVlYUzZ86grKwMAQEBcpdEAjty5AhmzJiB\n8PBwALeuwAMCAvDcc8/JXBmJ6KeffsI333yDmJgYm/UWiwX//e9/8ctf/lKmytyP3X0usm7dOnz1\n1Vf44osvEBwcLHc5JLiLFy8iNDTUumw0GjF16lReRZFd9fX1GDNmDGbNmmWzPj8/H+PGjcNTTz0l\nU2Xux5Bygeeffx4lJSX44osv8OCDD8pdDnkAtVqNmpoaXLp0CefPn8fu3bvxww8/oKenB+PGjZO7\nPBJMdXU1oqKirIOwLly4gL///e/Iz89HcXGxVw/OYkjdo/r6erz55psYO3YsoqKibLbdf//9uHTp\nEsaMYa8q2UpJSUFSUhJmzZoFpVKJf/3rX3jppZfw+OOP49ixY3KXR4KpqalBXV0dJk2aBIVCAZVK\nhV/96lcwmUwICwuTuzy3kvr7+/vlLoKIiMge/olPRETCYkgREZGwGFJERCQshhQREQmLIUVERMIS\ncgh6cHCw1w+rFF1zczM6OzvlLsMpbC9i8JQ2w/YiBmfbi5AhFRYWBpPJJHcZPk2r1cpdgtPYXsTg\nKW2G7UUMzrYXdvcREZGwGFJERCQsIbv7PJEkSUNu58QedCdH7QVgmyFbvnqO4ZWUi/T391t/7l72\n1sYzmNLSUsycORMRERHYunXrgO179uxBdHQ0oqOjER8fj9raWqeP9Rb22ocvtxlyzFfbC0OKXMpi\nsWDDhg0oKSlBQ0MD9u7di4aGBpt9pk+fjqNHj6Kurg6bN2+GwWBw+lgi8i0MKXKpqqoqREREIDw8\nHOPGjUNGRgaKiops9omPj0dgYCAAIC4uDi0tLU4fS0S+hSFFLtXa2opp06ZZl9VqNVpbWwfdv6Cg\nAKmpqcM+1mg0QqvVQqvVoqOjw0XVE5FoOHCCXMpe3/hgN3zLy8tRUFBgfX/ScI41GAzWbkJPeT6H\niIaPIUUupVarceHCBetyS0uLzWvSb6urq8PatWtRUlKCoKCgYR1LRL6D3X3kUrGxsWhsbERTUxN6\nenpQWFiI9PR0m33Onz+PZcuWYffu3ZgxY8awjiUi3+J0SFksFjz88MNIS0sbsK2/vx+/+93vEBER\ngejoaJw8edK6zVeGFNMtCoUCO3bsQHJyMmbNmoWVK1ciKioK+fn5yM/PBwDk5eXhypUrWL9+PTQa\njbW7brBjich3Od3d97e//Q2zZs3CtWvXBmwrKSlBY2MjGhsbUVlZiXXr1qGystI6pPjw4cNQq9WI\njY1Feno6Zs+e7dIvQWLR6XTQ6XQ263Jycqyf33//fbz//vtOH0tEvsupK6mWlhYcOnQIa9eutbu9\nqKgIq1evhiRJiIuLQ1dXFy5dusQhxUREdE+cCqlNmzbh9ddfx5gx9ncfbOgwhxQTEdG9cBhSBw8e\nhEqlwvz58wfdZ7Chw8MdUmwymWAymRASEuKoLCIi8gEO70kdP34c//73v1FcXIybN2/i2rVrePrp\np/HPf/7Tus9gQ4d7eno4pJiIiEbM4ZXUX/7yF7S0tKC5uRmFhYV4/PHHbQIKANLT07Fr1y709/fj\nq6++QkBAAB588EEOKSYionsy4od5bw8nzsnJgU6nQ3FxMSIiIjBx4kR8+OGHt375HUOKLRYLsrOz\nOaSYiIicNqyQWrJkCZYsWQLAdkixJEl455137B7DIcVERDRSnHGCiIiExZAiGiVKpRKSJNn9ATDo\nNqVSKXPlRPLhBLNEo+Tq1asjeoOqM6+aJ/JWvJIiIiJhMaSIiEhYDCkiIhIWQ4qIiITFkCIiId28\neRMLFixATEwMoqKi8Oqrrw6674kTJ+Dn54d9+/aNYoU0Gji6j4iENH78eBw5cgT+/v4wm81YvHgx\nUlNTERcXZ7OfxWJBbm4ukpOTZaqU3IlXUkQkJEmS4O/vDwAwm80wm812h+Nv374dy5cvh0qlGu0S\naRQwpIhIWBaLBRqNBiqVComJiVi4cKHN9tbWVuzfv99mmjZ7+L46z8WQIiJh+fn5oaamBi0tLaiq\nqkJ9fb3N9k2bNmHbtm3w8/Mb8vfwfXWei/ekRkipVOLq1auDbh9qloDAwEB899137iiLyCtNmTIF\nS5YsQWlpKebMmWNdbzKZkJGRAQDo7OxEcXExFAoFli5dKlep5GIMqREa6RQ3AKe5IXJGR0cHxo4d\niylTpuDGjRv47LPPkJuba7NPU1OT9XNWVhbS0tIYUF6GIUVEQrp06RLWrFkDi8WCvr4+rFy5Emlp\naTbvsiPvx5AiIiFFR0ejurp6wPrBwmnnzp1urmh08ZbCLQwpIiIB8ZbCLQ5D6ubNm0hISMBPP/2E\n3t5erFixAlu2bLHZ569//Sv27NkDAOjt7cU333yDjo4OKJVKhIWFYdKkSfDz84NCoYDJZHLPNyEi\nIq/jMKSceer7xRdfxIsvvggA+PTTT/Hmm2/avKitvLwcwcHBbiifiIi8mcPnpJx96vu2vXv3IjMz\n03UVEhGRz3LqYV5HT33f9uOPP6K0tBTLly+3rpMkCUlJSZg/fz6MRqNrqiYiIp/g1MCJ2099d3V1\nQa/Xo76+3uaButs+/fRTPPLIIzZdfcePH0doaCja29uRmJiIyMhIJCQkDDjWaDRaQ4zTlhARETDM\naZHufOrbnsLCwgFdfaGhoQAAlUoFvV6Pqqoqu8dy2hIiIrqbw5Dq6OhAV1cXAFif+o6MjByw3/ff\nf4+jR4/i17/+tXVdd3c3rl+/bv1cVlZm9wqMiIjIHofdfc4+9b1//34kJSXh/vvvtx7b1tYGvV4P\n4NbQ9FWrViElJcUd32PU9b86GfhjwMiPJSIihxyGlLNPfWdlZSErK8tmXXh4OGpra++tQkFJW67d\n04N2/X90bT1ERN6Ir+ogIiJhMaSIiEhYDClyudLSUsycORMRERHYunXrgO2nT5/GokWLMH78eLzx\nxhs228LCwjB37lxoNBpotdrRKpmIBMUJZsmlLBYLNmzYgMOHD0OtViM2Nhbp6emYPXu2dR+lUom3\n334bBw4csPs7OI0WEd3GkCKXqqqqQkREBMLDwwEAGRkZKCoqsgkplUoFlUqFQ4cOyVUmkfA4gvgW\nhhS5VGtrK6ZNm2ZdVqvVqKysdPr429NoSZKEZ599FgaDwe5+nKGEvB1HEN/CkCKXsvc/quG828bZ\nabQMBoM1wHjvish7MaTIpdRqNS5cuGBdbmlpsU6N5Qx702jZCylPNNLuG2/quiEaLoYUuVRsbCwa\nGxvR1NSEqVOnorCwEB999JFTx3Z3d6Ovrw+TJk2yTqP1yiuvuLni0TPS7htv6rohGi6GFLmUQqHA\njh07kJycDIvFguzsbERFRdlMo3X58mVotVpcu3YNY8aMwVtvvYWGhgZ0dnZ67TRaRDQyDClyOZ1O\nB51OZ7Puzmm0HnjgAbS0tAw4bvLkyV47jRYRjQwf5iUiImHxSuoeDGfU2p0CAwNdXAkRkXdiSI3Q\nUDfAJUka8fMNRET0/7G7j4iIhMUrKaJRNJIuYnYPky9jSBGNEnYREw0fu/uIiEhYDkPq5s2bWLBg\nAWJiYhAVFYVXX311wD4VFRUICAiARqOBRqNBXl6edZujdwsREdnjzLlnz549iI6ORnR0NOLj4/mc\nnRdy2N03fvx4HDlyBP7+/jCbzVi8eDFSU1MRFxdns9+jjz6KgwcP2qxz5t1CRET2OHPumT59Oo4e\nPYrAwECUlJTAYDAMa9Z90fExFydCSpIk+Pv7AwDMZjPMZrPT/+GcebcQEZE9zpx74uPjrZ/j4uLs\nzmTiqXgP8xan7klZLBZoNBqoVCokJiZi4cKFA/b58ssvERMTg9TUVHz99dcA7L9bqLW11e6/YTQa\nodVqodVq+X4gIgLg3LnntoLxQbz2AAAE2ElEQVSCAqSmptrdxvOL53IqpPz8/FBTU4OWlhZUVVWh\nvr7eZvu8efNw7tw51NbW4rnnnsPSpUsBDO/dQgaDASaTCSaTCSEhIcP9HkTkhRyde24rLy9HQUEB\ntm3bZnc7zy+ea1ij+6ZMmYIlS5agtLTUZv3kyZOtl+U6nQ5msxmdnZ33/G4hIiJg8HMPANTV1WHt\n2rUoKipCUFCQDNWROzkMqY6ODnR1dQEAbty4gc8++wyRkZE2+1y+fNl61VRVVYW+vj4EBQXZvFuo\np6cHhYWFSE9Pd8PXICJv48y55/z581i2bBl2796NGTNmyFEmuZnDgROXLl3CmjVrYLFY0NfXh5Ur\nVyItLc3m/UD79u3De++9B4VCgQkTJqCwsBCSJA36biEiIkecOffk5eXhypUrWL9+PYBb7zMzmUxy\nlk0uJvULOEREq9V6dEPzhpE3nvT/A0+qdTBsM6PHU+ocii+1F844QUREwmJIERGRsBhSREQkLIYU\nEREJiyFFRETCYkgREZGwGFJERCQshhQREQmLIUVERMJiSBERkbAYUkREJCyGFBERCYshRUREwmJI\nERGRsBhSREQkLIYUEREJy+Gbeck5kiQNuezpLygj17q7fdhbxzZDd/LVc4zDK6mbN29iwYIFiImJ\nQVRUFF599dUB++zZswfR0dGIjo5GfHw8amtrrdvCwsIwd+5caDQaaLVa11YvkP7+/iF/iO7kqL2w\nzdDdfLW9OLySGj9+PI4cOQJ/f3+YzWYsXrwYqampiIuLs+4zffp0HD16FIGBgSgpKYHBYEBlZaV1\ne3l5OYKDg93zDYiIyGs5DClJkuDv7w8AMJvNMJvNAy4z4+PjrZ/j4uLQ0tLi4jKJiMgXOTVwwmKx\nQKPRQKVSITExEQsXLhx034KCAqSmplqXJUlCUlIS5s+fD6PROOhxRqMRWq0WWq0WHR0dw/gKRETk\nraT+YXRmdnV1Qa/XY/v27ZgzZ86A7eXl5Vi/fj2OHTuGoKAgAMDFixcRGhqK9vZ2JCYmYvv27UhI\nSBjy3wkODkZYWNjwvolAOjo6EBISIncZ96S5uRmdnZ1yl+EUT28vANvMaGJ7EYOz7WVYo/umTJmC\nJUuWoLS0dEBI1dXVYe3atSgpKbEGFACEhoYCAFQqFfR6PaqqqhyGlCc09KFotVqYTCa5y/AZnt5e\nALaZ0cT24lkcdvd1dHSgq6sLAHDjxg189tlniIyMtNnn/PnzWLZsGXbv3o0ZM2ZY13d3d+P69evW\nz2VlZXavwIiIiOxxeCV16dIlrFmzBhaLBX19fVi5ciXS0tKQn58PAMjJyUFeXh6uXLmC9evX3/ql\nCgVMJhPa2tqg1+sBAL29vVi1ahVSUlLc+HWIiMibDOueFDnHaDTCYDDIXQZ5ELYZGg5fai8MKSIi\nEhbn7iMiImExpFwoOzsbKpWKg0PIaWwzNBy+2F4YUi6UlZWF0tJSucsgD8I2Q8Phi+2FIeVCCQkJ\nUCqVcpdBHoRthobDF9sLQ4qIiITFkCIiImExpIiISFgMKSIiEhZDyoUyMzOxaNEinDlzBmq1GgUF\nBXKXRIJjm6Hh8MX2whkniIhIWLySIiIiYTGkiIhIWAwpIiISFkOKiIiExZAiIiJhMaSIiEhYDCki\nIhIWQ4qIiIT1/wBaRBLoa7Z0JwAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"tags": []
},
"output_type": "display_data"
}
],
"source": [
"fig, axs = plt.subplots(2, 3)\n",
"axs[0, 0].boxplot(best_parameters.documented_infectious_tx_rate,\n",
" whis=(2.5,97.5), sym='')\n",
"axs[0, 0].set_title(r'$\\beta$')\n",
"\n",
"axs[0, 1].boxplot(best_parameters.undocumented_infectious_tx_relative_rate,\n",
" whis=(2.5,97.5), sym='')\n",
"axs[0, 1].set_title(r'$\\mu$')\n",
"\n",
"axs[0, 2].boxplot(best_parameters.intercity_underreporting_factor,\n",
" whis=(2.5,97.5), sym='')\n",
"axs[0, 2].set_title(r'$\\theta$')\n",
"\n",
"axs[1, 0].boxplot(best_parameters.average_latency_period,\n",
" whis=(2.5,97.5), sym='')\n",
"axs[1, 0].set_title(r'$Z$')\n",
"\n",
"axs[1, 1].boxplot(best_parameters.fraction_of_documented_infections,\n",
" whis=(2.5,97.5), sym='')\n",
"axs[1, 1].set_title(r'$\\alpha$')\n",
"\n",
"axs[1, 2].boxplot(best_parameters.average_infection_duration,\n",
" whis=(2.5,97.5), sym='')\n",
"axs[1, 2].set_title(r'$D$')\n",
"plt.tight_layout()"
]
}
],
"metadata": {
"colab": {
"collapsed_sections": [],
"machine_shape": "hm",
"name": "Substantial Undocumented Infection Facilitates the Rapid Dissemination of Novel Coronavirus (SARS-CoV2)",
"toc_visible": true
},
"kernelspec": {
"display_name": "Python 3",
"name": "python3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
![](\"https://www.tensorflow.org/images/tf_logo_32px.png\"){kind=logo}
View on TensorFlow.org\n",
" ![](\"https://www.tensorflow.org/images/colab_logo_32px.png\"){kind=logo}
Run in Google Colab\n",
" ![](\"https://www.tensorflow.org/images/GitHub-Mark-32px.png\")
View source on GitHub\n",
" ![](\"https://www.tensorflow.org/images/download_logo_32px.png\"){kind=logo}
Download notebook\n",
"