dli_gpr_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GP regression with heteroscedasic noise\n",
    "## Pyro implementation demo notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import pyro\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# import GP regression code\n",
    "from importlib import reload\n",
    "import dli_gpr\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading input data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 233,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time</th>\n",
       "      <th>proportion(y)</th>\n",
       "      <th>numcellsincluster(x)</th>\n",
       "      <th>totalcellsinsample(n)</th>\n",
       "      <th>patientID</th>\n",
       "      <th>sampleID</th>\n",
       "      <th>tumorburden</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-30</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>1</td>\n",
       "      <td>329</td>\n",
       "      <td>DFCI-5310</td>\n",
       "      <td>1</td>\n",
       "      <td>94.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>180</td>\n",
       "      <td>0.0931</td>\n",
       "      <td>21</td>\n",
       "      <td>1310</td>\n",
       "      <td>DFCI-5310</td>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-42</td>\n",
       "      <td>0.1725</td>\n",
       "      <td>127</td>\n",
       "      <td>2834</td>\n",
       "      <td>DFCI-5309</td>\n",
       "      <td>5</td>\n",
       "      <td>11.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>90</td>\n",
       "      <td>0.0808</td>\n",
       "      <td>88</td>\n",
       "      <td>2030</td>\n",
       "      <td>DFCI-5309</td>\n",
       "      <td>6</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>0.3389</td>\n",
       "      <td>7</td>\n",
       "      <td>3184</td>\n",
       "      <td>DFCI-5311</td>\n",
       "      <td>9</td>\n",
       "      <td>44.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>90</td>\n",
       "      <td>0.2036</td>\n",
       "      <td>32</td>\n",
       "      <td>3762</td>\n",
       "      <td>DFCI-5311</td>\n",
       "      <td>11</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>270</td>\n",
       "      <td>0.0487</td>\n",
       "      <td>24</td>\n",
       "      <td>8662</td>\n",
       "      <td>DFCI-5311</td>\n",
       "      <td>12</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0</td>\n",
       "      <td>0.4367</td>\n",
       "      <td>116</td>\n",
       "      <td>1358</td>\n",
       "      <td>DFCI-5312</td>\n",
       "      <td>21</td>\n",
       "      <td>38.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>270</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>1</td>\n",
       "      <td>13</td>\n",
       "      <td>DFCI-5312</td>\n",
       "      <td>22</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>-5</td>\n",
       "      <td>0.1644</td>\n",
       "      <td>10</td>\n",
       "      <td>1174</td>\n",
       "      <td>DFCI-5317</td>\n",
       "      <td>23</td>\n",
       "      <td>80.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>345</td>\n",
       "      <td>0.0228</td>\n",
       "      <td>5</td>\n",
       "      <td>4390</td>\n",
       "      <td>DFCI-5317</td>\n",
       "      <td>24</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>-30</td>\n",
       "      <td>0.2759</td>\n",
       "      <td>5</td>\n",
       "      <td>1566</td>\n",
       "      <td>DFCI-5314</td>\n",
       "      <td>25</td>\n",
       "      <td>94.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1095</td>\n",
       "      <td>0.1840</td>\n",
       "      <td>38</td>\n",
       "      <td>1853</td>\n",
       "      <td>DFCI-5314</td>\n",
       "      <td>26</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>-416</td>\n",
       "      <td>0.4462</td>\n",
       "      <td>4</td>\n",
       "      <td>995</td>\n",
       "      <td>DFCI-5311</td>\n",
       "      <td>29</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>-90</td>\n",
       "      <td>0.6127</td>\n",
       "      <td>52</td>\n",
       "      <td>1544</td>\n",
       "      <td>DFCI-5311</td>\n",
       "      <td>30</td>\n",
       "      <td>16.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>-152</td>\n",
       "      <td>0.3670</td>\n",
       "      <td>65</td>\n",
       "      <td>1447</td>\n",
       "      <td>DFCI-5314</td>\n",
       "      <td>31</td>\n",
       "      <td>78.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>29</td>\n",
       "      <td>0.2993</td>\n",
       "      <td>4</td>\n",
       "      <td>1216</td>\n",
       "      <td>DFCI-5314</td>\n",
       "      <td>32</td>\n",
       "      <td>98.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>-416</td>\n",
       "      <td>0.0164</td>\n",
       "      <td>4</td>\n",
       "      <td>731</td>\n",
       "      <td>DFCI-5317</td>\n",
       "      <td>38</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>-50</td>\n",
       "      <td>0.1313</td>\n",
       "      <td>4</td>\n",
       "      <td>1051</td>\n",
       "      <td>DFCI-5317</td>\n",
       "      <td>39</td>\n",
       "      <td>80.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>-416</td>\n",
       "      <td>0.0427</td>\n",
       "      <td>68</td>\n",
       "      <td>1430</td>\n",
       "      <td>DFCI-5309</td>\n",
       "      <td>40</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    time  proportion(y)  numcellsincluster(x)  totalcellsinsample(n)  \\\n",
       "0    -30         0.0000                     1                    329   \n",
       "1    180         0.0931                    21                   1310   \n",
       "2    -42         0.1725                   127                   2834   \n",
       "3     90         0.0808                    88                   2030   \n",
       "4      0         0.3389                     7                   3184   \n",
       "6     90         0.2036                    32                   3762   \n",
       "7    270         0.0487                    24                   8662   \n",
       "8      0         0.4367                   116                   1358   \n",
       "9    270         0.0000                     1                     13   \n",
       "10    -5         0.1644                    10                   1174   \n",
       "11   345         0.0228                     5                   4390   \n",
       "12   -30         0.2759                     5                   1566   \n",
       "13  1095         0.1840                    38                   1853   \n",
       "14  -416         0.4462                     4                    995   \n",
       "15   -90         0.6127                    52                   1544   \n",
       "16  -152         0.3670                    65                   1447   \n",
       "17    29         0.2993                     4                   1216   \n",
       "18  -416         0.0164                     4                    731   \n",
       "19   -50         0.1313                     4                   1051   \n",
       "20  -416         0.0427                    68                   1430   \n",
       "\n",
       "    patientID  sampleID  tumorburden  \n",
       "0   DFCI-5310         1         94.0  \n",
       "1   DFCI-5310         2          0.0  \n",
       "2   DFCI-5309         5         11.0  \n",
       "3   DFCI-5309         6          0.0  \n",
       "4   DFCI-5311         9         44.0  \n",
       "6   DFCI-5311        11          0.0  \n",
       "7   DFCI-5311        12          0.0  \n",
       "8   DFCI-5312        21         38.0  \n",
       "9   DFCI-5312        22          0.0  \n",
       "10  DFCI-5317        23         80.0  \n",
       "11  DFCI-5317        24          0.0  \n",
       "12  DFCI-5314        25         94.0  \n",
       "13  DFCI-5314        26          0.0  \n",
       "14  DFCI-5311        29          0.0  \n",
       "15  DFCI-5311        30         16.0  \n",
       "16  DFCI-5314        31         78.0  \n",
       "17  DFCI-5314        32         98.0  \n",
       "18  DFCI-5317        38          0.0  \n",
       "19  DFCI-5317        39         80.0  \n",
       "20  DFCI-5309        40          0.0  "
      ]
     },
     "execution_count": 233,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import os\n",
    "\n",
    "# read input data\n",
    "test_fname_r = os.path.join(\"data\", \"v3_MetaclusterTE_R_gprinput.txt\")\n",
    "test_fname_nr = os.path.join(\"data\", \"v3_MetaclusterTE_NR_gprinput.txt\")\n",
    "test_df_r = pd.read_csv(test_fname_r, sep=\"\\t\", header=0)\n",
    "test_df_nr = pd.read_csv(test_fname_nr, sep=\"\\t\", header=0)\n",
    "\n",
    "# get rid of rows with NaN\n",
    "test_df_r.dropna(inplace=True)\n",
    "\n",
    "# get rid of rows with NaN\n",
    "test_df_nr.dropna(inplace=True)\n",
    "\n",
    "test_df_r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time</th>\n",
       "      <th>proportion(y)</th>\n",
       "      <th>numcellsincluster(x)</th>\n",
       "      <th>totalcellsinsample(n)</th>\n",
       "      <th>patientID</th>\n",
       "      <th>sampleID</th>\n",
       "      <th>tumorburden</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-30</td>\n",
       "      <td>0.0506</td>\n",
       "      <td>1</td>\n",
       "      <td>791</td>\n",
       "      <td>DFCI-5322</td>\n",
       "      <td>3</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>165</td>\n",
       "      <td>0.0935</td>\n",
       "      <td>15</td>\n",
       "      <td>759</td>\n",
       "      <td>DFCI-5322</td>\n",
       "      <td>4</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>0.1418</td>\n",
       "      <td>241</td>\n",
       "      <td>1692</td>\n",
       "      <td>DFCI-5324</td>\n",
       "      <td>7</td>\n",
       "      <td>93.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-14</td>\n",
       "      <td>0.0771</td>\n",
       "      <td>1</td>\n",
       "      <td>1804</td>\n",
       "      <td>DFCI-5325</td>\n",
       "      <td>17</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>730</td>\n",
       "      <td>0.1688</td>\n",
       "      <td>736</td>\n",
       "      <td>3756</td>\n",
       "      <td>DFCI-5325</td>\n",
       "      <td>18</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>-4</td>\n",
       "      <td>0.1195</td>\n",
       "      <td>39</td>\n",
       "      <td>1515</td>\n",
       "      <td>DFCI-5326</td>\n",
       "      <td>19</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>120</td>\n",
       "      <td>0.0466</td>\n",
       "      <td>36</td>\n",
       "      <td>5409</td>\n",
       "      <td>DFCI-5326</td>\n",
       "      <td>20</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0</td>\n",
       "      <td>0.3891</td>\n",
       "      <td>49</td>\n",
       "      <td>1627</td>\n",
       "      <td>DFCI-5318</td>\n",
       "      <td>27</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>182</td>\n",
       "      <td>0.3800</td>\n",
       "      <td>4</td>\n",
       "      <td>771</td>\n",
       "      <td>DFCI-5318</td>\n",
       "      <td>28</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>-543</td>\n",
       "      <td>0.3027</td>\n",
       "      <td>3</td>\n",
       "      <td>2709</td>\n",
       "      <td>DFCI-5322</td>\n",
       "      <td>33</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>-540</td>\n",
       "      <td>0.1500</td>\n",
       "      <td>30</td>\n",
       "      <td>3599</td>\n",
       "      <td>DFCI-5322</td>\n",
       "      <td>34</td>\n",
       "      <td>34.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>-428</td>\n",
       "      <td>0.3142</td>\n",
       "      <td>48</td>\n",
       "      <td>1200</td>\n",
       "      <td>DFCI-5322</td>\n",
       "      <td>35</td>\n",
       "      <td>97.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>-274</td>\n",
       "      <td>0.7735</td>\n",
       "      <td>87</td>\n",
       "      <td>2945</td>\n",
       "      <td>DFCI-5324</td>\n",
       "      <td>36</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>285</td>\n",
       "      <td>0.0216</td>\n",
       "      <td>1</td>\n",
       "      <td>878</td>\n",
       "      <td>DFCI-5318</td>\n",
       "      <td>41</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>-543</td>\n",
       "      <td>0.0094</td>\n",
       "      <td>12</td>\n",
       "      <td>1493</td>\n",
       "      <td>DFCI-5321</td>\n",
       "      <td>42</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>0</td>\n",
       "      <td>0.1209</td>\n",
       "      <td>1</td>\n",
       "      <td>430</td>\n",
       "      <td>DFCI-5321</td>\n",
       "      <td>43</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    time  proportion(y)  numcellsincluster(x)  totalcellsinsample(n)  \\\n",
       "0    -30         0.0506                     1                    791   \n",
       "1    165         0.0935                    15                    759   \n",
       "2      0         0.1418                   241                   1692   \n",
       "4    -14         0.0771                     1                   1804   \n",
       "5    730         0.1688                   736                   3756   \n",
       "6     -4         0.1195                    39                   1515   \n",
       "7    120         0.0466                    36                   5409   \n",
       "8      0         0.3891                    49                   1627   \n",
       "9    182         0.3800                     4                    771   \n",
       "10  -543         0.3027                     3                   2709   \n",
       "11  -540         0.1500                    30                   3599   \n",
       "12  -428         0.3142                    48                   1200   \n",
       "13  -274         0.7735                    87                   2945   \n",
       "14   285         0.0216                     1                    878   \n",
       "15  -543         0.0094                    12                   1493   \n",
       "16     0         0.1209                     1                    430   \n",
       "\n",
       "    patientID  sampleID  tumorburden  \n",
       "0   DFCI-5322         3        100.0  \n",
       "1   DFCI-5322         4        100.0  \n",
       "2   DFCI-5324         7         93.0  \n",
       "4   DFCI-5325        17        100.0  \n",
       "5   DFCI-5325        18        100.0  \n",
       "6   DFCI-5326        19        100.0  \n",
       "7   DFCI-5326        20        100.0  \n",
       "8   DFCI-5318        27        100.0  \n",
       "9   DFCI-5318        28        100.0  \n",
       "10  DFCI-5322        33          0.0  \n",
       "11  DFCI-5322        34         34.0  \n",
       "12  DFCI-5322        35         97.0  \n",
       "13  DFCI-5324        36          3.0  \n",
       "14  DFCI-5318        41        100.0  \n",
       "15  DFCI-5321        42          0.0  \n",
       "16  DFCI-5321        43        100.0  "
      ]
     },
     "execution_count": 234,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_df_nr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show some descriptive statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 241,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean tumor burden, responders:  31.65\n",
      "Mean tumor burden, non-responders:  76.6875\n",
      "Variance of tumor burden, responders:  1494.1274999999998\n",
      "Variance of tumor burden, non-responders:  1570.46484375\n",
      "Overall variance:  2028.8888888888885\n",
      "Standard error:  56.35802469135801\n"
     ]
    }
   ],
   "source": [
    "tb_r = test_df_r[\"tumorburden\"].values\n",
    "tb_nr = test_df_nr[\"tumorburden\"].values\n",
    "\n",
    "print(\"Mean tumor burden, responders: \", np.mean(tb_r))\n",
    "print(\"Mean tumor burden, non-responders: \", np.mean(tb_nr))\n",
    "\n",
    "print(\"Variance of tumor burden, responders: \", np.var(tb_r))\n",
    "print(\"Variance of tumor burden, non-responders: \", np.var(tb_nr))\n",
    "\n",
    "print(\"Overall variance: \", np.var(np.concatenate([tb_r, tb_nr])))\n",
    "print(\"Standard error: \", np.var(np.concatenate([tb_r, tb_nr]))/(len(tb_r) + len(tb_nr)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 243,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3828.891966759003\n"
     ]
    }
   ],
   "source": [
    "diffs = tb_r[1:] - tb_r[:-1]\n",
    "print(np.var(diffs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "78.63828947368421"
      ]
     },
     "execution_count": 248,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.var(tb_r[np.argsort(test_df_r[\"time\"].values)])/len(diffs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing model input\n",
    "\n",
    "Below we prepare two models: one to describe tumor burden and one to describe the proportion of cells found in that metacluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 358,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare torch tensor for time points\n",
    "y = torch.FloatTensor(test_df_r.time.values)\n",
    "\n",
    "# reponse variables (one for metclusters, one for tumor burden)\n",
    "metacluster = torch.FloatTensor(test_df_r[\"proportion(y)\"].values)\n",
    "tumorburden = torch.FloatTensor(tb_r / 100.)\n",
    "\n",
    "# number of cells size in each time point\n",
    "cluster_sizes = torch.FloatTensor(test_df_r[\"totalcellsinsample(n)\"].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model initialization and optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 373,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variances\n",
      "Tumor burden:  0.007863828912377357\n",
      "Metacluster:  0.0015130525252896516\n"
     ]
    }
   ],
   "source": [
    "reload(dli_gpr)\n",
    "\n",
    "# tumor burden model\n",
    "tumor_model = dli_gpr.gpr(y, tumorburden)\n",
    "\n",
    "# metacluster model\n",
    "metacluster_model = dli_gpr.dli_gpr(y, metacluster, cluster_sizes)\n",
    "\n",
    "# initialize kernel and center the response variable\n",
    "tumor_model.initialize_variables()\n",
    "metacluster_model.initialize_variables()\n",
    "\n",
    "print(\"Variances\")\n",
    "print(\"Tumor burden: \", tumor_model.sigma.item())\n",
    "print(\"Metacluster: \", 1./metacluster_model.gam.item())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimize each model\n",
    "We will run 2000 iterations of SVI for each model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 376,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 5000/5000 [00:05<00:00, 998.12it/s] \n",
      "100%|██████████| 2000/2000 [00:03<00:00, 543.23it/s]\n"
     ]
    }
   ],
   "source": [
    "pyro.clear_param_store()\n",
    "losses1 = tumor_model.optimize(n_steps=5000)\n",
    "\n",
    "pyro.clear_param_store()\n",
    "losses2 = metacluster_model.optimize(n_steps=2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 368,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x13428ad90>]"
      ]
     },
     "execution_count": 368,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZhU1Zn48e/LKi6IKCCCikai4i5ESUxmjIpgTMQkY4Y4KpnJDPk5JqOPyRhI4qiJRrMZl0QS44YrIqKiCIogi4i0jew0Lc3aDU1309D73vX+/qjb3dXVte9V9/08Tz1VdereW+di+95T7zn3HFFVjDHGuEOvdFfAGGNM6ljQN8YYF7Ggb4wxLmJB3xhjXMSCvjHGuEifdFcgnBNOOEFHjRqV7moYY0xWWbt27UFVHeJfnvFBf9SoUeTn56e7GsYYk1VEZE+gckvvGGOMi1jQN8YYF7Ggb4wxLmJB3xhjXMSCvjHGuIgFfWOMcREL+sYY4yIW9OPw1vp91Da1prsaxhgTMQv6Mdqyv5rbZ69n+uub0l0VY4yJmAX9GDW2tANwoKYpzTUxxpjIWdA3xhgXsaBvjDEuYkHfGGNcxIK+Mca4iAV9Y4xxEQv6xhjjIhb0jTHGRSzoG2OMi1jQN8YYF7Ggb4wxLmJB3xhjXMSCvjHGuEhEQV9EBonIXBHZJiIFIvJlERksIotFZLvzfJzP9jNEpEhECkVkok/5WBHZ5Hz2mIhIMk7KGGNMYJG29B8FFqnqWcAFQAEwHViiqqOBJc57RGQMMAU4B5gEPCEivZ3jzASmAaOdx6QEnYcxxpgIhA36IjIQ+CfgaQBVbVHVKmAyMMvZbBZwvfN6MjBbVZtVdRdQBFwiIsOBgaq6WlUVeN5nH2OMMSkQSUv/dKACeFZE1onIUyJyFDBMVUsBnOehzvYjgGKf/UucshHOa//yHkRkmojki0h+RUVFVCdkjDEmuEiCfh/gYmCmql4E1OOkcoIIlKfXEOU9C1WfVNVxqjpuyJAhEVTRGGNMJCIJ+iVAiaqucd7PxXsRKHNSNjjP5T7bn+yz/0hgv1M+MkC5McaYFAkb9FX1AFAsImc6RVcCW4H5wFSnbCrwlvN6PjBFRPqLyGl4O2zznBRQrYiMd0bt3OKzjzHGmBToE+F2PwFeEpF+wE7g3/FeMOaIyA+BvcANAKq6RUTm4L0wtAG3qWq7c5xbgeeAAcBC52GMMSZFIgr6qroeGBfgoyuDbP8A8ECA8nzg3GgqaIwxJnHsjtwYBeyBNsaYDGdB3xhjXMQVQV9VeSVvL9WNrQk7ps0fYYzJRq4I+pv2VTNj3iZ+PndjuqtijDFp5Yqg39TqAaCyvjnNNTHGmPRyRdA3xhjjZUHfGGNcxIK+Mca4iAV9Y4xxEQv6xhjjIq4P+is+r6C+uS3d1TDGmJRwddDfW9nALc/kcZeN3zfGuISrg36d08LfUVEX9b42944xJhu5Iuh7l+Q1xhjjiqCfDDb3jjEmG7ki6HsX6jLGGOOKoG/pHWOM8XJF0DfGGOPliqBv6R1jjPFyRdA3xhjj5YqgHyynbz8AjDFuE1HQF5HdIrJJRNaLSL5TNlhEFovIduf5OJ/tZ4hIkYgUishEn/KxznGKROQxSXPexfp3jTFuE01L/+uqeqGqjnPeTweWqOpoYInzHhEZA0wBzgEmAU+ISG9nn5nANGC085gU/ymEZzl9Y4zxiie9MxmY5byeBVzvUz5bVZtVdRdQBFwiIsOBgaq6Wr35lud99kmqWIdstnuUptb2BNfGGGPSJ9Kgr8D7IrJWRKY5ZcNUtRTAeR7qlI8Ain32LXHKRjiv/ct7EJFpIpIvIvkVFRURVjE88buPNtwPgDteXc9Zdy8K+Jllhowx2ahPhNtdpqr7RWQosFhEtoXYNlAo1RDlPQtVnwSeBBg3blzC4qtGGarf3rA/UV9tjDEZIaKWvqrud57LgTeAS4AyJ2WD81zubF4CnOyz+0hgv1M+MkB50gXL6cfTkWu9BMaYbBQ26IvIUSJyTMdr4GpgMzAfmOpsNhV4y3k9H5giIv1F5DS8HbZ5TgqoVkTGO6N2bvHZJ6lsGgZjjPGKJL0zDHjDaS33AV5W1UUi8ikwR0R+COwFbgBQ1S0iMgfYCrQBt6lqR2/orcBzwABgofMwxhiTImGDvqruBC4IUF4JXBlknweABwKU5wPnRl/N+ARL79hITmOM20TakZvVEpneaWxpZ9uBmoQdzxhjUskV0zB08B+yGcu14H/nbuDbT3xMeW1zgmpljDGp46qgH+2QzUA276sGoN5ZX3d9cRXf+/tqmtvsJi5jTOZzVdD3l4icfrtHydt1iB3l9fEfzBhjksxVQd8/vWOMMW7jqqCfiPRO17GMMSb7uCLoh5tlc9uB2uiPZVHfGJOFXBH07Y5cY4zxckXQ75DQnL7foRKZOjLGmGRxVdBPaGC2GG+MyUKuCvrGGON2rgr6yUzvGGNMNnBV0E8muwfAGJMNXBX0k5nTt45cY0w2cFXQTwRrzxtjspmrgr6lYIwxbueqoG8pGGOM27kq6CdCq8cDQHO7J801McaY6Lkq6CcivVN8qBGAp1fujPtYxhiTaq4K+onU0NJ90RSb3scYkw1cG/SbWttpbos9RWOLqhtjslHEQV9EeovIOhF5x3k/WEQWi8h25/k4n21niEiRiBSKyESf8rEissn57DEJN+dxEl38m8Vc/9dVne/boszR+7fsozmTz/YetuUVjTFpEU1L/3agwOf9dGCJqo4GljjvEZExwBTgHGAS8ISI9Hb2mQlMA0Y7j0lx1T4O/umZPYcaoto/1mzOzoo6vvPEx/z67a0xHsEYY2IXUdAXkZHAtcBTPsWTgVnO61nA9T7ls1W1WVV3AUXAJSIyHBioqqvVO8H98z77ZIRZH+/mwYUF4TeMw+GGVgC2ltYk9XuMMSaQSFv6jwB3Ab45kGGqWgrgPA91ykcAxT7blThlI5zX/uU9iMg0EckXkfyKiooIqxi/e+Zv4e/LYxuVE6ojt6XN02Mhl5LDjTF9jzHGxCNs0BeRbwLlqro2wmMGym5riPKehapPquo4VR03ZMiQCL82cRK50lZNUytf/NVCHl9a1K28orY5Yd9hjDGRiqSlfxlwnYjsBmYDV4jIi0CZk7LBeS53ti8BTvbZfySw3ykfGaA8aeZ9VsL/vbU56v1eyy8Jv1GEDte3ADB3beKOaYwxsQob9FV1hqqOVNVReDtol6rqTcB8YKqz2VTgLef1fGCKiPQXkdPwdtjmOSmgWhEZ74zaucVnn6S4c84Gnl+9p+tcQnS/+v4MOVDTlMRaGWNM+vSJY9+HgDki8kNgL3ADgKpuEZE5wFagDbhNVTuGytwKPAcMABY6j5TJjAnXui48Ta3tHNG3d4htjTEmsaIK+qq6DFjmvK4Ergyy3QPAAwHK84Fzo61kJkpE2v/dTaV85+KR4Tc0xpgEcdUduZHOshlZQI8u6nd9dyb82jDGuJWrgn5iRRa8MyOlZIwxXvHk9LOOILy35UDniJpUCPXrwiZpM8akmquCvqL86IXIbjeob24LezRjjMk2rgr6ofjO/fbCJ7v58wefJ+a4PdI7drEwxqSPq3L6bZ7IAu7ButSlf4wxJpVcFfQPpTCXb4wxmchVQT+Rou2EtU5bY0wmsKCfcjaE0xiTPq4I+slYoCvaQ9ryisaYTOCKoN8hlXHX41H+6fcfMn/DPiBwescyPsaYVHNV0E9kkA2Xo29p97D3UAN/fN9/6KeFemNM+rgq6CeSx3pmjTFZyIK+Ix0pd0vzG2NSzYJ+jBLROWy/FYwxqeaqoL+nsiFhx4o95Fv73hiTPq4K+sYY43YW9GMUbWqmq9/XkjrGmPSxoJ9GaiOAjDEpZkE/AsWH4u8LSMUdua3tHkZNX8ALn+xJ/pcZY7KSK4J+JAE31DZf+/2HPcqibaXXN7dR1dB9ls9ETg+xds9hisrrAPjT+4UJO64xJreEDfoicoSI5InIBhHZIiL3OeWDRWSxiGx3no/z2WeGiBSJSKGITPQpHysim5zPHpNkTIoTQKKyKO0+8/H7H9L/O/zfH25o5cJfL05MRQL47syPuebRlUk7vjEmN0TS0m8GrlDVC4ALgUkiMh6YDixR1dHAEuc9IjIGmAKcA0wCnhCR3s6xZgLTgNHOY1ICzyXpXssvTujxov21oKp85cElvL62JOR2VQ2tPLtqVzxVM8bkqLBBX73qnLd9nYcCk4FZTvks4Hrn9WRgtqo2q+ouoAi4RESGAwNVdbV6o93zPvskVaJ+T9T5rJsb7pDJ+A3T5lH2Vzfx89c3ht32vre3xv19re0emlrb4z6OMSZzRJTTF5HeIrIeKAcWq+oaYJiqlgI4z0OdzUcAvk3iEqdshPPavzzQ900TkXwRya+oqIjmfAK64W+rw27Tcy3b0FI17ubzslpm5+1N0bd1N/HPKzjr7kVp+W5jTHJEFPRVtV1VLwRG4m21nxti80DRU0OUB/q+J1V1nKqOGzJkSCRVjJtm6Pj5q/+8gunzNnUrS1VNdx6sT9E3GWNSJarRO6paBSzDm4svc1I2OM/lzmYlwMk+u40E9jvlIwOUZ41Q/c4i8OG2cr7y4BKaWtuTsjxiuiZweGBB/KkiY0xmiGT0zhARGeS8HgBcBWwD5gNTnc2mAm85r+cDU0Skv4ichrfDNs9JAdWKyHhn1M4tPvukXSTpnVAdr6pw39tb2F/dxPz1ybmWJbuFX3K4gVHTF/Dmun3dyv+x0jqFjckVkbT0hwMfishG4FO8Of13gIeACSKyHZjgvEdVtwBzgK3AIuA2Ve3oDbwVeApv5+4OYGECzyWlQrXk73p9Y1Yuj/h5WS0Ab63fF2bLzFB4oNbuajYmSn3CbaCqG4GLApRXAlcG2ecB4IEA5flAqP6AjJai2wqCf39avz2zLCss5wfPfsof/uV8bhh3cvgdjDGAS+7IjUQ64nmsbVRr3dJ593FBaW2aa2JMdrGgn4U8Cm3tnnRXwxiThSzoJ4j/kM9kNMZ9D3nGL7O2OyShMnWorTGZyoJ+GmVSjn5PZT37qprSXY2Ipbt/xZhsFbYj1wTmn1evbmzt9j5v96Hwx0hojeLzz39Ylu4qGGNSwFr6jo6OwWCKDzVQXhO8JXzz03nd3k99Ji/IloHNyS9m6/6aqPYxxphoWUvfsbigLOTngebUj0WwXP9dc72TqO1+6Nqo9020TPoFYoxJLGvpO9w4CrK+uY1lheXhN8xgbvzvZkw8LOg7Glvawm+UAN36H9McsO6au5EfPPspeyq7T6yWDV2k2VBHYzJRzgb9aG9gejPK+XICHb2pNfzY+dc/y5wpDnZUePsx6pvDz5mfaTeEZVZtjMkeORv0V24/mNwvCBB1DoTo6O3w8hqfufGtuWqMSbGcDfr1zalJ18Tjrrkb8Xgib7MGuhHp4x0Hw448ilY2tKLtemlMbGz0TqwSFHVqm9o49si+Me9/4z/WJKYiYaimZ34iY0xi5WxLP9dkWk49lGc+2sWo6QuobWoNv3GM7AJkTGxyNujPXVsSfqN4JDAGN7akbvHxgtIaWsNM1hbxupZBPL96NwAH61qi2Cs6WXQNNCaj5GzQX7ItO8afz9+4n8MN4YOjKjS3xTez5s6KOq55dCW/X7QtruOEY/HYmMyVs0E/W9z95uaItlPg12/Ht1ZtRW0zAOuLq/yO7TdDaFzf0tUKT2YGxtI7xsTGgn6M0tGa3R7nKJ2OOkeyHnCPfWPIpwQKzNvLalOazjLGdGdBP0ap7liN9fvydh2ivLbJOYZT6BeMY7kIxKKxpZ0Jf17B/8xel7BjZlMHtzGZwIJ+jvve31dz7WMfdSvrCPHBOlrj7cgNtrBJi9MnsWZnZRRHC8yyO8bExoJ+jFK9iIdC1ENW/rZ8B9CVy/cNxjVNrRysaw7+XXHoyukn79/I2vfGxCZs0BeRk0XkQxEpEJEtInK7Uz5YRBaLyHbn+TiffWaISJGIFIrIRJ/ysSKyyfnsMcni5Y/qknTHb1F5Lb95Z2uPtEUsWYyHFvqN0nGOsWbXIc6/9/0Ya2iMyWaRtPTbgJ+q6tnAeOA2ERkDTAeWqOpoYInzHuezKcA5wCTgCRHp7RxrJjANGO08JiXwXHLCD579lKc/2kXJ4UaaWv06POO8RsbTOs601HnWthaMSbOwQV9VS1X1M+d1LVAAjAAmA7OczWYB1zuvJwOzVbVZVXcBRcAlIjIcGKiqq9XbjH3eZx/j6EyNSPeWuqLsjHf0TqoWYfE5h6R/V/K/IqHufHU9T63cme5qGBeLKqcvIqOAi4A1wDBVLQXvhQEY6mw2Aij22a3EKRvhvPYvD/Q900QkX0TyKyoqoqli1vNN63SMuvGWQ20cKaWSww0Rbxu4Izezwmu2ZAZHTV/Aj17I73w/b90+7l9QkMYaGbeLOOiLyNHA68AdqhpqMddggz8iHhSiqk+q6jhVHTdkyJBIq5i1OjpafSU6qD23anfEgTtZ4T2RF45sGqr53pbQS3Eak0oRBX0R6Ys34L+kqvOc4jInZYPz3DHvQQlwss/uI4H9TvnIAOWuN/mvq3qUJTqoKZmXlzfGpF4ko3cEeBooUNWHfT6aD0x1Xk8F3vIpnyIi/UXkNLwdtnlOCqhWRMY7x7zFZx+TZE9/tItbnsmLef9gF4zdB+u7paG823o39v+xEssQzpY2T+f4/m7HSmJ658kVOyJKh725bh8b/Ka0MCbTRdLSvwy4GbhCRNY7j28ADwETRGQ7MMF5j6puAeYAW4FFwG2q2jEM5VbgKbyduzuAhYk8mVziH9QytZV++R+XcckDS7qVdd74m4DAPO7+xZx7z3txHydSpdWN/PbdbfzHc5+G3faOV9cH/JVmTCYLu4iKqn5E8BFyVwbZ5wHggQDl+cC50VTQbYLF9mR2pPofO1Ft6MP1LRw7oC9H9+8T8HsiUdPU1nms447q1+PzRF8M252VzCJZN9iYbGR35GYo/8BbcrgxZd+dqDj6zcc/4vI/fBjRtm+u2xfyhreLfrM4QbWKTDZ1FBsTDQv6Gaa02pkcza/86j+vSH1lYuQbLw/WtfBK3l5Om7GgsxXtb0NxFXe8up5fvrEpAd+tfLT9oAVtY4KwoJ+hGprbeHfTgZR8V3VD92UNlxX2vDdClZ53CEfotwsKUIWGIFMqd7TwS6uaAn4eSLDugtfyS7jp6TXM+2xf1PXsfvzsuA/AmGhZ0M9QE1LYsr/xqZ6Lq++s6H737+b91Zx19yLe3xL+QhRt7v7lNXsByNt9KKr9Ail2Rt2ES4epKj+fu5G1ew4H/dyYXGRB3wR0xZ+Wd3u/bq83OC7/PPF3SAeb7TMSsXZC17e082p+Mbc83f2CF0sL/2BdM7sO1ke9nzHpEHb0jjEQ3Xw60TaSeyUhlZLKaSO+/OASWtuV3Q9dm7LvNCZW1tI3UUlG1sM/5t/64lp+9tqGyPbtsQxY6nPxre2WCjLZw4K+Sbh4Q+DCzQeYu7aE+ggmmAvWoo/34rS/OvJOZWOyiQV9E5Guu2wTf+xg6R1PiMgdrBqxVm/F5xX8Y4VNeWxyn+X0TdJ1hO45+cUBP0/GhSTahn7HvETfOH944itjTAaxoJ/FPi+rTdl3dTS6X/xkL8OPHcCJA4/gzfWBx8IHa6A/vrQI6JpaIRlseL0xoVnQz2KpvEvXN3f+h/cKw27tK1QcbmnzsHL7wdgrFrQKsSX1I71mTP7LRzEd35h0s5y+iUg8HaOhdq1tag3xaQTH9jt4LNM3x2JDSXXQzyrrmtlfFdtcSbVNrbS295xK2phEsaBvEq4tyBw78Xp7g8+aO2HyOOkcRDn2/g/4ykNLY9r3vHvfZ9rz+eE3NCZGFvRNwlX5zeUTKjy/vzXypQR/8so69laGXtwk2px+fZD5gNLpwwBzHxmTKBb0TdKFanXPmNd9Zs33wszt09wWWZD2T/u8s3E/r366N6J9s53HYzONmuAs6JuM8qMX1sa1f7CG/o9fXsfPX49/6uZYLdxUyvjfLgm/YQI8s2oXNz29hsVR/Ioy7mFB30RkeYakHCJtuz65MvE3Wr3wyR6ueXQleyqjn1zt1pc+40BNau7y3eOkwHy/r7qxlVHTF/BBHBeCFZ9XUNXQEnf9THpZ0DcRKThQE3ab1nZPxOmXSEQS4P236cjpB1pMHaDkcAOb9wUfeROqT+DuNzdTUFrDxhAjdzLVdueejpnLd0S13+odleytbKCmqZVbnsnjv2LoZK5tau0xVbdJHxunbyITQQSe8PBydofpaE2UWAdmfvV33uUbN983MeDnuZQGD3Qu0eb5v/+PTwBY+6urANhREf2vnBv+tpptB2ptFtIMEbalLyLPiEi5iGz2KRssIotFZLvzfJzPZzNEpEhECkVkok/5WBHZ5Hz2mNjSRFmlNoLJz4IF/EBr31Y3hh+f/6s3NvcoCxezIv2z2hKitR9Osq4Lh+u7UieNcYwqCvRPkM7/27YdSN2d4ya8SNI7zwGT/MqmA0tUdTSwxHmPiIwBpgDnOPs8ISK9nX1mAtOA0c7D/5jGRb73t9Vht5nvOy4/wfYcCnyBiic4zpi3MaLzCubmZ7oWdKlviX2qilAXxhz6IWNiFDboq+oKwH8du8nALOf1LOB6n/LZqtqsqruAIuASERkODFTV1er9ffm8zz4mB4yaviCq7QtjnDeoYzqIWNfr9TlQVHwXdV8RZPWwV/KK41rysaA0sS3i7hcw+2FtvGLtyB2mqqUAzvNQp3wE4DuVYolTNsJ57V9uTEzuX1AAeKc8GPubxXxYWB7V/tGurPX6Z11/vnPXloTYMnYS5HUiqcIf3ytk1PQFMY3jt7H/2S/Ro3cC/a1qiPLABxGZJiL5IpJfUZEZQwVNZtpeVkdlfQsPvbsNCLzerm+uvEM0sWvd3sP8/PWNMdcxUonOu/uOYPruzI87X//lw6LEfpHJKrEG/TInZYPz3NHMKgFO9tluJLDfKR8ZoDwgVX1SVcep6rghQ4bEWEWTi/yD9f7q7hObPbtqd499LvrN4p7HCXL8nT6jU2qcyeC+O/PjmEb1RNsqDjZZ3N+X7wiaUgrl/gUFrNlZGfTzWM7pcEN8E+SZ9Is16M8HpjqvpwJv+ZRPEZH+InIa3g7bPCcFVCsi451RO7f47GNMzJpaY5uRMljA+7enujpT/2uWd0x6rAPNXos2DRTkax5cuK1zkZeIDuNznDW7Yu9jMLkpkiGbrwCrgTNFpEREfgg8BEwQke3ABOc9qroFmANsBRYBt6lqR4/brcBTeDt3dwALE3wuxgWueXRlQo4TainGDlv2e29IizXrctfc8CmhQOkoCD0RXWVdM6XVsU3drEFeG/cIe3OWqn4/yEdXBtn+AeCBAOX5wLlR1c6YANbuib/16hvw2oNMBV3X3MZZdy901vBNfIjM332If/nbah7//kV864KTul1cZszbxPcvOSXgfmPv/wAg6M1OIa9nPh9600/dL2lr9xzmyH69OXv4wG7ldltN7rBpGEzaxDoS5Lsze46FDxaTioOMx/cNflv2B79Rq6nVQ0uSFjXZWur9JfGTV9YBPc+hqbU9/qGpBL+oBfLdmR8n7NeUyUwW9E3a/OrNnnfcJlqwAJbu1IbHLxA3trT36J846+5FfMlp1QdTeKCWL/ziXUoOB5/+YsLDywOWp/vfIJftrKhj1PQFYacKTwcL+iZtXlqTuPntVXvOtb/884qAU0BA96CbjuUJX87rfu7/78XAU0qHm/7ilby9tHuU97Z09QH4/2LYeTD6+XL8rdwe/eghVe1xcYvG5n3VfOEX78bcf5Fo2w7URPy30jEp34KNpcmsUkws6JucUFhWy5m/WtStbE5+cZCtu7dyH/lge5JqFZz/r5zlEQzJvMlnZFEwCzeVdruY+me9unXk+rz59dtbeS3Ev9fts9f77Kc9UnPby2o55Hc/xG/eKeD0X7zb+X5DcVXY+oP314vHo7z4yR7aPcqyDJjWe09lPZMeWcmDzv0g4WRyF4gFfZNTvnrGCZ2vQ/1/55vnPliXnjnio40LHxUd7PZ+dt5eVu+o7DxWU2s7t770WcQ5/KLyOn7xxiaKDzXwzKpd/G8Eo43Au47vlX/qnjKa8OcVTHpkRbeyZ1bt6vZ+8l9XsWjzAZpa23l21a6A9dxQXMXER1bw9xWJXw8hHh1/I+uKD0e1Xyam0GxqZZNTfFtYHQExkOYg8+2nUrwLyE/3W2py4eboUgnfeMzb31FUFt1c93XNbQHTZuW13YefivQcSbS7sp4nPqzmsaVFHN2/DzeMO7nb5yWHvamcjSVVHDugLxB8NJKqcvvs9fzbpadw6enHR3UObmYtfZNTfANEQ4jpif/wXmHn63T9Er/v7a0JPd789eFnJd0bYDRTPJPEhRJw7hWFKmda7frmNjaVVDPh4eXUNbfR2NLOkoKyzu2CpUjW7KykrKaJplYP8zfsZ+qz4W9c+2RnJaOmL2C307+xeV81M5ftSOiiP3srG3rMDJuJcxVZS9/kFN8USKR51V450PT5oKCMjwP8smn1+zVRlcJpFCRAU99/orvfLdrG9vI6Lv714m5DYxdtOcCXnda7/z7/+uQnHHdkX1bP8N4qFMkPpnnOhHlrdlUy6oSj+ObjHwHem+Pu/uaY6E4siGsfW0ltcxvXXXBSRt/XkAN/7sYEFqql7ysDG2NBtQUZPRIo4AM8tiR5ndQlhxsYNX1BjxlOf/baBn78cuC+hd8vKuw2f09HQA90L8RqZ94gVfjT+4Xs8FlysdscQHH894t2ZNC6vVWcPmMB9QHSW4FGWmXin5YFfWOyyLx1+9JdhU7r9npH4/hONd3U2s7ctSW8E2Ko4ttRLo5zsK6Zx5cWcfNTa1i6rWtoasfqa5FMqRGMqnexnmAX0w6+DXePdi0+H8ianZVdqS313kF+oLop6PapZkHfmCzyaYomUJv0yIqIlrT0d9bdi8Jv5Lj37a2sKgre2d6hY6x/dWMr//Fc18Lsl/52CdC9Nb15XzU3P70maK5eFVb5pAAXbj7A/7yyjqc+2hVw+2BCZW/+9clPeMPn4s/3VucAAA4DSURBVPzdmau54k/Lojp+MlnQN67XMbFaNoh65s4YbTtQy11zN0S0bbJvQPqgwJs+qg+SrvPtLP3FG5tYuf1g53j6T3cforS6sXPaaqX7TKodVu+o5KU1eyKu04biKjaWBL/vYE+lt8N4wSbvv01DSzu3vfRZxMdPppwN+r0ytx/FmKzge5dvOjWGmX/It+ugI9f+3Me7Abjhb6uZ8PCKzpb5K3mB7wJf/nkFv3wj8mlBps/bxHV/WRX080AduQs2lcZ1h3Ki5GzQP/PEgeE3MsbErGOiuGTbFcU0EjVNXZ2pHSN26prbKHbmJuqYHiGc6obWqObN8Z/YL1ibM9iqZRtLqig8kNg1koPJ2aBvjHGPWR/v5s4566nwuUHszjld6alI+g6ga2nN/355LT96YS1lNaE7YF930m1f+/2H3cr79QkcWoNddK77yyomOnc0Vze0ct4975GXpP4bG6dvjMl698zfkpDj+C+t2bHO8NYg/T4/fW0D3x07skd5sH6iDwpCp8zumruBr44eQm1zG48t2c6L/3lpJNWOiquD/j3fGpPwuyKNMbnjRy+s5azhxzDvs+QNlfW9i3dOfgn5e7zz+8QzFDWUnE3vFJQGH5Hx9o+/ymd3T+Cqs4elsEbGmGyztbQm4QH/3Hveo6m1nQvue58/vV/I//j1jeys8PZh+M9amig5G/RDOW/ksQw+qh8nDz6SvF8GXPWx09H9o/sx1D9ILu/ebyXmVm9jTHara26j5HAD1Y2tPL40cMcueIfNJkPOBv03/vsrHHNEH2ZPGx9yu6HHHBGwvI8z5nP86YMj/s4BfXvzj1vGAd2n+AX44onHRHwcY0z2mPjnFeE38vN0lDeDJVLOBv2LTjmOTfdOZLwzadOkc04Muu3tV47uUTbzprEhjz9sYH/OHdF9WOg3zx9Ou5OH652gGwWiuegYY1KvsCz6FvkrecEXrEm2lAd9EZkkIoUiUiQi01PxnbsfupaZN10c9PPbrxzNb799HoX3T+rxWbC+lNd+9BXe+cnXGNC3d2eZCFx88nEcO6AvP7nijB77DDmmf7f3NwTo9ffXx5kCsl+fXrwUpif/Z1d/kSU//eewxzTGZIdk3MyV0qAvIr2BvwLXAGOA74tISpLdIsK/XzaKS0b1bDn36iXceOkp9O/Tm4e/dwH3XXdOZyv+e1/qWuRh2MCuoH3SoCOc43Yd5+oxJ3LskX3ZcM/VjPP7nvNHDuLqMV0dxycc3Z8/3HBB0Pr+9tvnMejIvvyf0xfw1m2XcdkZJ/DsD77E188cEnCf6y4YweAj+wX87I6rev6aidT3Lzml8zxPHNg9HXbxKYP42ugTAuwVmce+f1Hn63+79JSYj2NMLgo0+2i8Uj1k8xKgSFV3AojIbGAykJJxk/d865yw23zn4q7W9+6HrgXg+gtP4gtDjuYnV46mtLqRXQfr6dPbe73s5UTDjfdezcAj+gY85s7ffoNevYT7rjuHaf90Ov/8h2V8+6KTADj+qH5U1rfwxWFHs6OinjOHHcPxR/fjxktP4UYnCHbUA+DrZw3l62cN5Y/veaeaXbj5AN+64CQe9wmei+74GpMeWdmtDjePP7VzLdiO4+2rauSyh5YC0K93r25/YE9PHccPZ3knt3rwO+fx4HfOA7zznHxQUO5dJxWY6Jc2u+rh5RSV1zH82CNoafPw4f9ezvn3vg94A3xVQws1ja388f3PAThl8JGd+15z7nAWby3rsQKTMW7Vr3fi2+WpDvojAN9kVgnQI2chItOAaQCnnJL+1t8jU7oC6vBjBzD82AGd778w9Gg2FFfRO8BcG9OvOYvHl2ynl5Pf79O7F6cefxSb75vIkU5aaOlPL6e+pY2TBg3osX8oP5t4Js1t7Qyav4U7J5zZ7bOzThzI1l9PxKOwZV81uw7Wc/zR/TmqX28uPvW4zu1GDBrAsp9dzpBj+nNU/z5U1jUz9v4PALjy7GE8OuXCHvUSESaMCT7U9c3bLqO6sZURAc7nugtO6nx96enH09zq6VwSD+D0IUfxyJQLufEf4RcAv/iUQXy2t4qpXz6VWy8/g7KaJr4w9GjOvec9wPvrJNg8K1edPZSaxjaGDuzfYwrgm8afwouf9Nxv/OmD+e/Lz+CJZUV8srPrTskRgwZw7fnDeTKJa7r+7Oovdl4kA7n2vOGdE3uZ3NIrCZOISSqX8xKRG4CJqvqfzvubgUtU9SfB9hk3bpzm5+cH+zjtDte3sKGkisvPHJruqmSsmqZW2tqVwUcFTj2t2VnJ2ScN7Pyl9PGOg9z4jzVM/fKprCw6yBu3Xsa3n1jFfZPP4ean8zj+qH6svXtCwGMVlNYw6Mi+nDjwCGqb27jz1Q38+2WjOFjXzIbiap5ZtavzV1l9cxu/fnsr//m10+jdS9iyv8bbGe9R5q3bR2VdC1efM4zCA7VcdfYw+vXpRbtHafcoa/ccRlX5yhknMDtvL9PnbeKvN17MkoIyfnDZKE49/ijW7KykrrmNg3XN7K9q4s6rv8ihuhYu/+MywDu813et3pvGn8LQY47gnJMG8kre3s7ZJTfeezX7DjdyzaMrO39BVTpjuF/+r0sZf9rx/O69bVx/4QiuebTrF97H06/g87JafvDspwB8adRx3Pb1Mxg4oC8jBw2gvLaZT3ZWcv+CAm69/Avc8uVTeS2/hOOP7sf404/nmkdW0tLu4eTBAyg+1MhFpwxi3d4qjujbi6bW4GmHY/r3oba5jV9dezZXnj2M6x7/KOACI/5++Y2zueDkQSwrLGfbgVqWbisPuN2TN49l2gtrgx7njqtGd/6q9XfLl09l76EGlhVWhK1Puk0YM6xzNGAsRGStqvY8gKqm7AF8GXjP5/0MYEaofcaOHavGdGhpa9fWtvZ0V6Mbj8ej+bsrY9p3aUGZvrtxf9T7NbW2aVu7p0c9/rJ0u24vq9ENxYcjPlZ7u0c9Hk+Pco/HoxW1Tdru9z01jS36/Ord6vF4dFVRhW7ZV60Ha5s6j9He7tGt+6t7HK+splEbW9q0+FC9VtW3aE1ji5bVNGp9c6vWN7cGrFv+7kPa3NqujS1t2tjS1q0uHo9H9x1u0EN1zaqqur+qQQsP1HTWY0d5rT7xYZGqqhaUVmtjS1vnvjsr6nTptjJ9Z8N+bWpt09a2dt1VUac1jS364bYyffSDz7W6sUU/3VWpza3tuu9wgz6y+PNu/wZlNY2673CDvvrpXm1ubdf2do++sHq3Prl8h1bVt6iqakOzt97X//UjrWpo0W2lNfrMRzu1saVNlxeWa2Vds+4+WKdLC8q0rd2jn+05pB6PJyF/40C+BoipqW7p9wE+B64E9gGfAjeqatCJMzK9pW+MMZkoWEs/pTl9VW0TkR8D7wG9gWdCBXxjjDGJlfIJ11T1XeDdVH+vMcaYHL4j1xhjTE8W9I0xxkUs6BtjjItY0DfGGBexoG+MMS5iQd8YY1wkpTdnxUJEKoA9Me5+AnAwgdXJBnbO7uC2c3bb+UL853yqqvaYkjfjg348RCQ/0B1puczO2R3cds5uO19I3jlbescYY1zEgr4xxrhIrgf9J9NdgTSwc3YHt52z284XknTOOZ3TN8YY012ut/SNMcb4sKBvjDEukpNBX0QmiUihiBSJyPR01yceIvKMiJSLyGafssEislhEtjvPx/l8NsM570IRmehTPlZENjmfPSYSYFHfDCEiJ4vIhyJSICJbROR2pzxnz1tEjhCRPBHZ4JzzfU55zp4zgIj0FpF1IvKO8z7Xz3e3U9f1IpLvlKX2nAMtp5XND7yLs+wATgf6ARuAMemuVxzn80/AxcBmn7LfA9Od19OB3zmvxzjn2x84zfl36O18lod3uUoBFgLXpPvcQpzzcOBi5/UxeFdbG5PL5+3U72jndV9gDTA+l8/ZqeudwMvAOy75294NnOBXltJzzsWW/iVAkaruVNUWYDYwOc11ipmqrgAO+RVPBmY5r2cB1/uUz1bVZlXdBRQBl4jIcGCgqq5W71/M8z77ZBxVLVXVz5zXtUABMIIcPm/1qnPe9nUeSg6fs4iMBK4FnvIpztnzDSGl55yLQX8EUOzzvsQpyyXDVLUUvAESGOqUBzv3Ec5r//KMJyKjgIvwtnxz+rydVMd6oBxYrKq5fs6PAHcBHp+yXD5f8F7I3xeRtSIyzSlL6TmnfLnEFAiU23LLuNRg556V/yYicjTwOnCHqtaESFvmxHmrajtwoYgMAt4QkXNDbJ7V5ywi3wTKVXWtiFweyS4ByrLmfH1cpqr7RWQosFhEtoXYNinnnIst/RLgZJ/3I4H9aapLspQ5P/Fwnsud8mDnXuK89i/PWCLSF2/Af0lV5znFOX/eAKpaBSwDJpG753wZcJ2I7Mabgr1CRF4kd88XAFXd7zyXA2/gTUen9JxzMeh/CowWkdNEpB8wBZif5jol2nxgqvN6KvCWT/kUEekvIqcBo4E85ydjrYiMd3r5b/HZJ+M4dXwaKFDVh30+ytnzFpEhTgsfERkAXAVsI0fPWVVnqOpIVR2F9//Rpap6Ezl6vgAicpSIHNPxGrga2EyqzzndvdnJeADfwDviYwfwy3TXJ85zeQUoBVrxXuF/CBwPLAG2O8+Dfbb/pXPehfj06APjnD+wHcBfcO7GzsQH8FW8P1c3Auudxzdy+byB84F1zjlvBv7PKc/Zc/ap7+V0jd7J2fPFO6Jwg/PY0hGbUn3ONg2DMca4SC6md4wxxgRhQd8YY1zEgr4xxriIBX1jjHERC/rGGOMiFvSNMcZFLOgbY4yL/H/pHIX/2nymkAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(losses1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 367,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(0.0947)\n",
      "tensor(0.0035)\n"
     ]
    }
   ],
   "source": [
    "print(tumor_model.sigma)\n",
    "print(1./metacluster_model.gam)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results visualization\n",
    "We compute the conditional distribution of each response variable at a number of new input points to model what our data should have looked like at intermediate points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 378,
   "metadata": {},
   "outputs": [],
   "source": [
    "# new input time points\n",
    "new_y = torch.FloatTensor(np.linspace(-600, 1200, 100))\n",
    "\n",
    "# compute conditional distribution of response variable\n",
    "new_tumor, cov_tumor = tumor_model.conditional_distribution(new_y)\n",
    "new_metacluster, cov_metacluster = metacluster_model.conditional_distribution(new_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 379,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Metacluster proportion')"
      ]
     },
     "execution_count": 379,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAADQCAYAAAD4dzNkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2deZgU1fX3v2eYGWaGfZWZYRVF3BISZsG4oSZxV0AD4pK4JOgvaszrkiC4EdyTmESzuKC4RjGRUYPGPcZopmcYDAguKCAww6AsgoLs0+f941TR1T3d01t1d93q83meerq7qvrWqVv33HPuqbsQM0NRFEVRFLMoyLUAiqIoiqIkjxpwRVEURTEQNeCKoiiKYiBqwBVFURTFQNSAK4qiKIqBqAFXFEVRFANRA64kDREVEhET0dAMXuNmIno4U+krihchopuI6PFcy+EHiOheIro+13JkEjXgHoGItjq2IBFtd/w+J9fyKYofIaKVRLSLiPpG7F+YqJNKRGOJqCVTMiYLEQ21ZC/MtSzZgojOJ6K3nfuY+RJmnpkrmbKBGnCPwMxd7Q3AagCnOvY9kSu5Ml0JEFEBEWk5VHLJpwAm2z+I6FAApbkTJ7e4ofPZdB7yyVGJRCtOQyCix4noJsfv7xLRSsfvFiK6moiWWK32+4loHyJ6mYi+IqJXiKin4/xxRPQ+EW0mojeI6ICItK4hosUAtnUg1qlE9CkRbSCi221DHBn+JqL9iIgdv98moplEVA/gawCDiWhfIvoPEW0hopcB9Im4/8OJKGDJu5CIjopIbwYR/df6/0tE1DuZ/FXymscA/NDx+0cAHnWeQESdieg3RLSaiD63wrOlRNQFwD8BVDgiZhVEVENE9VZ5XUtEfySiYkd6BxPRq0T0hZXetEihorXsrYjBd63vNUTUZOn350R0l3XaW9bnZkuew6zzLySiD4lok1UvDHGky0R0KRF9AuCTKLLYrfopRNRq3dNVjuM3EdHfrXrqKwDnW3n2e+v8Vut7Z+e9EdE0q/5YSY5IIxH1IKJHiWg9Ea0iousc9cv5RPQOEf2OiL4AMAfAvQAOs+53s3Xew0R0syPNnxDRMivPnyeiioj7v4SIPrHy509ERJH54DXUgPuLCQCOBTASwBkAXgDwCwD9AXQGcCkAENGBAB4HcDmAfgBeA/APIipypHUWgBMB9OjgeqcD+DaAKgBnIrwSjMd5AC4E0B1AC4CnAAQA9AVwu3UclryDADwP4EYAvQFMBTCXiJxG/mxIxbsPgC4ArkxCFiW/CQDoTkQHElEnAJMg+uHkDgAjAIwCsB+ASgA3MPPXED1pdUTMWgG0Afh/kPJ8GIDjAPwUAIioG0TnXgJQYaX3egpy/wHAH5i5O4DhAJ629tvObU9LnnoiGgdgGqSO6AfgPwCejEhvHIBaAAd1cM1jAOwP4PsAptrOhMXpAP4OoCeAJwBMBzAGkmffBFAD4DrH+QMg+VMJ0d37KdSQuAdS9+wL4GhI3XKB47+1AFZA6rZzAVwCoN66356IgIiOBXAbgIkAygGsgtQ5Tk4BUG3JOhHA8R3kgydQA+4v/sDM65i5BcDbkAK9iJl3AHgWwLes884C8Dwzv8HMuyEGsztEKZxptTDz9g6udzszb2LmlQDuhiMMmQAPMfOH1vUHQ5T8Rmbeycz/AvCi49wfWvK+zMxBZn4JwCIAJzjOeZCZP2HmbQD+ZqWnKIlit8K/B+AjAGvsA1ZL7CcA/h8zf8HMWwDcCtGjqDDzAmYOMPMeSz/ugxgiQAzFZ8z8W2bewcxbmLkhBZl3A9iPiPoy81ZmDnRw7sUAbrN0bo8l/yhnK9w6/kUcnZ/BzF8z82IAsxGu8/XM/Kylo9sBnAPgV1adtB7ADDgcc4vrLZ3/N6TBMdHhRF1r5c1KAL+N+G8rM99j5W9H8tqcA6lz3mXmnQCuhbTYhzrOuZ2ZNzPzagD/ggF1iBpwf/G54/v2KL+7Wt8rIB4oAICZg5BWcKXj/OYEruc8Z5WVbqI4/1sBYKNlfJ3p2QwBMNkKR262QmRjIq73meP7NoTuVVES4TFIFOd8RITPIS3WMgALHOXvJWt/VIhoBBHNI6LPrJDyrZDWJgAMArDcBZkvgkQFPiKi+UR0SgfnDgHwB4f8XwAguKvzkf8Pq2einL/JimBEHu8LoDjKf5OVNaYszLwVwMaINI2rQ9SAm8PXkErEZkAaabVCFBqAdCQDMBCOVgeARJapG+T4PthKF0hMVmf6awH0ISJnx6HBju/NAGYzc0/H1oWZf52AjIoSF2ZeBenMdhKAuRGHN0Ac4IMd5a+H1eEUiK4rf4G05Pe3QtzTIAYTkPI8PAGxwvTIapnudRqsiNNkSBj5DgB/t97JR5OnGcDFETpUysz/dWZDAjLF0vlo/w+rZ6Kc38uSN/L4Bkh0IfK/HdVP8WSPrPO6QPrZrIn5DwNQA24OCwGcTES9iKgcwM/SSOtpAKdZHUmKAFwDYAuAZMN4vyCinkQ02JJnjkPWo4loEEnHuakdJcLMywG8B+AmIiom6aB2suOUxwCMJ6LvEVEnIiohomOcnVAUxQUuAnBsRKvQjlA9AOB3RNQfAIiokojsd6SfQxxQZ3+RbgC+ArCViEYC+D/HsXkABhDRz62OXt2IyPn6yuZjACVEdLKlp9dB+rLAkuFcIupnybfZ2t0GYD2AIOT9sc29AK4looOt//Ygoh8kmjEOrieiMiudCxDS+Wg8CeA6IupHMkzvBrTvWzDD0vkjIa8W/sbMbZA66hYrb4ZA+rR0ND7+cwADydFRMIK/AriAiEZZHeluBdBgheeNRQ24OTwM4ENIGOgltO+AkTDM/D6k08hfIMp+AoDTrPfRyfAPiLH+H4A6S0ZY8tUBWAygEdIBLR5nATgcEtqbDjHatrwrAYwHcL0l72oAV0HLr+IizLycmZtiHP4lgGUAAlZI/DUAB1j/+whirFZYIeoKAFdDQvJbIMZ/r6Gz3qF/D8CpkLDtJ5DOYZHyfAnp+DYL0lL8GvKqy+YEAO8T0VZIh7azrHfq2wDcAuAdS54xzFwHaaU/Zcm/BNL5Lln+beXD6wB+w8yvdHDuzQCaIM75YgDvWvtsPgOwCdI6fgLAJVZeAtLB9mtIR7W3IQb4oQ6u9QaA9wF8RkQbIg8y8+uQ+uMZSMRvODrow2AKxJxI1ERRFEXJV6zOXp8CKLI6waWb3lgAjzPzwHTTyme0BaMoiqIoBqIGXFEURVEMREPoiqIoimIg2gJXFEVRFAPx9CTwffv25aFDh+ZaDEUxggULFmxg5piTi3gF1WtFSY5Yuu1pAz506FA0NcUa1aEoihMiWhX/rNyjeq0oyRFLtzWEriiKoigGogZcURRFUQxEDbiiKIqiGIgacEVRFEUxEDXgiqIoimIgasAVX8LMaGhpQDAYDPvUiYsUxVxi6XW+6rkacMUXRCp2Q0sDxs8Zj9kLZ+/9HPfUODz0v4fyRrkVxU8wMx7630MYN2dcmF5H+7TrAL/ruqfHgStKojSuacT4OeMx85iZuP5f12PuxLmom1SH6opqHNL/EFRXVAMArvvXdTi438EgItRU1oCIciy5oigdwcxoXNMIZsZ1b1yHm4+9GReMumCvXkf7ZGaMnzMedZPqUDsw2lLr/sDTc6FXVVWxTvigdISt3NUV1ZjfOn/vZzTj7KwIJjw9AXMnzvWVISeiBcxclWs54qF6rSSDHU1LRl8j6wXTdTyWbmsIXTEau+U9v3U+agfWoqCgALUDa6MqKxGhdmAtagfWom5SHQBg/JzxaFzTmG2xFUWJg/1arLqiem9LOpZuR2Lr+vzW+b4OqasBV4wkUrlrKmsS/m+kIa+uqPalciuKyUQ656m0oGsqa3ztrKsBV4zEDeWO9NL9ptyKYiLpOOeROJ31uRPngpl95airAVeMw1bCuRPnpqXcNraX7kZaiqKkhxvOeSREBCLChKcn+MpRVwOuGEfjmkZMeHrCXqVMF9tLB6ChdEXJMZlyqP3oqKsBV4wjU4poe/5+8tAVxRTs0DkA11reTvzoqKsBV4wh0wruRw89HYjoISJaR0RLYhwnIrqbiJYR0XtE9O1sy6j4h2w50H5y1NWAK8aQacXzo4eeJg8DOKGD4ycC2N/apgD4SxZkUnxKthxoPznqasAVY8iW4vnJQ08HZn4LwBcdnHI6gEdZCADoSUTl2ZFO8QuZjqxFYjvqJk/sYqMGXDGGbCmenzz0DFMJoNnxu8Xa1w4imkJETUTUtH79+qwIp5hBrhxm23EwOdKmBlzxPNlWND956BkmWgZFfUjMfD8zVzFzVb9+/TIslmISuXKY/RBpUwOueB710D1LC4BBjt8DAbTmSBbFUHLlMPsh0qYGXPE86qF7lucB/NDqjT4GwJfMvDbXQilmkGsH2Q+RNjXgiudRDz03ENGTAOoBHEBELUR0ERFdQkSXWKe8CGAFgGUAHgDw0xyJqhiIVxzkXDsS6aDrgSuexV4SMFdLATqHleUjzDw5znEGcGmWxFF8hlccZNuRMHHtcG2BK57FKx66oiju45UQtlcciVRwxYDrjE1KJvCKYpkcYlMUL+IlnfKKI5EKbrXAH4bO2KS4jFcUSyMBiuIuqlPu4IoB1xmbcsPHHwO33AL89rfAxo25lsY9vOSdA96JBCiKX/CiTnmt3kmEbL0DT3jGJiUx3nsP+Pa3gRtvBKZPBw49FPiiIxfKILzmnXslEpBPmFiZKonjRZ3yWr2TCNky4AnP2KRTLkZn+3agrS30+447gG3bZN/OncDmzcBTT+VOPjfxoneuZBcTK1PFbEysd7JlwBOesUmnXAxn3TpgzBigWzegpASYOhVgBnbtkk+bYBDYvTt3crqJF71zbRFmFxMrUyU+XtYjL9Y78ciWAdcZm1LkjDOApiZpae/ZA9xzD/Dgg8AVVwBlZXJOQQFQWgqceWZuZfUz2iLMLiZWpkp8VI/cxZWJXKwZm8YC6EtELQBuBFAEAMx8L2TGppMgMzZtA3CBG9f1O21tQH19eOh82zagrg544QXgtdeA++4DunQBrroKqNReBRlDW4SKkj4m6FGuJ5BKBlcMuM7YlBkKCoDiYnn/bUME9O4t3w87TDa/4GXFyfdZ2RTFDUzQI5NmZtOZ2DwMEfCrX4VC5UTS2p4+PbdyZQoNrymRePmdqeJPTIgS2KgB9zhXXw088ggwaRLwk58ACxYAI0fmWqrMYJLiKNlBnTol25jU/0IXMzGAM8/Mjw5qJoTXvBzm9yPq1PkD1ZvMoC1wRUkCbRFmF5NaQ0psTNQbE17fGG3ATchgJT4mPUdtESpK8pioNyY4HUYbcBMyOF127pRJWmy2bgWWLZP9fsGk56gtQkVJHhP1xgSnw2gDbkIGx8JplKOxYwdw/PHSA71bN+CZZ4Dnnwf22QcYNQoYOBBYvDg7smYak5+jkh1MitIo/sAEp8NoA25CBkfyxBNA9+5AURFQWwt89ln08265BXjrLTH027YB550HnHWWfP/6a2DDBmDixOzKnilMfI5KdjEpSqMo2cJoA24ajY0yFGzLFjHM774LnHJK9HOXLJFWuE2nTuFznwNAczMUJS/QKI2ZaOQks6gBzyL//nf4giN79gD/+58sTBLJqaeGT+DSuTPQs6d8B2SGtu98J/MyK9HRiim7aJTGTPwQOfGyrqsBzyK9e4vhdVJcLOF0G2ZGoDmAA78XwE03MQ48EDj6aOCdd8QBGDVKDPlxx5m/fKiXFSMefqiYlOzw+efAtGnAj38M/O1vuZYmu/ghcuJlXScvV55VVVXc1NQU9zxTJgnYvh0YPRpYuVJ6kZeUAHfdBVx8cegemBmnPClx9XmT54GIPH9fqdLQ0mDMnMOReLHMEdECZq7KtRzxSFSv/cC6dcChhwJffCERt7Iy4JprgJtuyrVkSqJkSteTSTeWbvuiBe5lD8lJaaksDXrXXcCMGcA//wlMmSKtUNuYAWK4502eBwAYP2f83uNedrZSwWTvXEO62cfEiM1f/gJs2iTGG5BOqLfc4q9hoH4nU7ruht3yxVSqJhmCsjLgkktCvxta5CHOnTh37z3YBYWZUTepDsxsbEu1I0yYOlXxDiatEmWzfn14vxdAVhn8+mvp16LkL27YLV+0wE1sDdmtieqK6r0VUuQ92PdVO7AWdZPqUF1RbVwLRFHcwiRH3ebEE0OdUQEx3kOHAr165UykrGBitKQjMnE/btgtXxhwk7ALgh0yn986P+5DtB/0/Nb5RrwqUJRMYKKjfvLJ8rqsqAgoLARGjABefjk0msSvmPJaM1G8ej9qwLOMXRAAJN2aMLEFEgu/eOh+uQ8lc1x9tXRg3bQJ+PBDaYH7HT/VVYB370cNeJZIJGRus2OHKPr27eH7ne+MTTcaXvVok8Uv96HExtbdYDAY9pmM/nXqBHTtmkEhPYaJ0ZKOcPN+3HT6fWXAvdgaSiZk/tprwJFHylSrBx0kCl9VBbzxRvh5fjAaXvVok8Uv9xELIjqBiJYS0TIimhrl+Fgi+pKIFlrbDZmWKVt6Hqm7sxfODvv06+gQJbO4Wn8zs2e30aNHczIEmgNc/ptyDjQHkvpfJrFlql9dz4HmAAeDwXbnfPUV86RJzADzsGHMv/gF8yOPMN9wA/Pw4bL/wguZd++W84PBYMy0lPwFQBO7qH8AOgFYDmBfAMUAFgE4KOKcsQDmJZNusnodSTb0PBgM8qwFs3jAbwbs1d22trawz/rV9Z6rb5TM4kbdm0oasXQ750a6oy1ZRfeSYbNlsZU9lkzr1jFXVTF36sQ8Ywbz9u3hx7dtY546VZ7Uuecyt7W1v4YX7lfJPRkw4IcBeNnx+1oA10ack3UDno1yH2gO8IBfD+BZC2bFvE6iOp5P+L1OylUjMS8MuJdI5EFv3Mh84IHMJSXMzz/fcXo33yxP64orkruG1/C7gueSDBjwMwHMcvw+D8AfI84ZC2Cj1Tr/J4CD46Vrgl4nU05N1MNM4fe8yFX9FUu3ffUO3CvYmTt34tyY70Z37QLOOANYvlxmZDv11I7TnDYNuPxy4A9/kGEogJnvX/3w/j4azN7rf+EC0XrsRN7guwCGMPM3AdwD4NmoCRFNIaImImpav369y2K6h/0cASTcacnWQ52nwcw6KRnS6cyWiTpCDXgGaFzTiAlPTwARxXzQP/0p8OabwOzZwNix8dMkAu68Uzq3XXQRsHmzmT09/argPnVMWgAMcvweCKDVeQIzf8XMW63vLwIoIqK+kQkx8/3MXMXMVf369cukzI5rAvfcA1RXAxMnAp99Fv8/qTxHt+ZpWLQIePJJWWbYVEysk7JFRuqIaM1yr2wmhNqcJPpO7NFHJRx+3XXJX6OxUd6X/+Qn7a+rYenc4YVnAPdD6IUAVgAYhlAntoMjzhmA0KJINQBW279jbW7pdbw8f/BB5rIy0bXCQuYRI5jjPZ50nmM678Tvuktk7dZNPm++OenLK1kklXKSTtmKpds5N9IdbaYZ8ETe/yxdytylC/ORR4Z6lTsJBoNcv7qe61fXx6wIfv5z5oIC5g8+SPy6iv9x24BLkjgJwMeQ3ujTrX2XALjE+n4ZgPct4x4A8J14abql1/HK/WmnSQ1nb8XF0u8kGm46YIl0gHOycSNz587hspaUMK9Zk7YoSobIdp2bUQMO4AQASwEsAzA1yvGxAL4EsNDabkgk3VQVPVetoXjX3b2buaaGuVcv5tWro6cRaA5wnzv6cJ87+vCsBbOiDkFbv1489fHjE7uuYi7JPNtMGPBMbNlqgV93nRhC2yj26MG8Z0/0tNyskJ1D0BJJ76OPmLt2DTfg3bszv/tu2qJkjXyrg5K5XzfyJmMGHBkaK8ppKHq2vaNEH9Add0iOP/lk7DTa2tratcCjjTf91a8krfp6t+8mM+SbgrtFMmU53wx4PLZvZz71VAmfDxjA/M47sc91u3wmE07fsYO5b99wA96rF/PWra6IkhXyNQqYSLlxI28yacAzMlaU01D0bBuLRB7Qhx9KmGz8+Ojv4TpKI1plsGULc79+zCed5OadZI58VfB00RZ4ZnGzrti2jTkQYF64MDRfQ6LlftEi5qFDmYmYBw1ibmpKW5yskq8OeiLP1+stcFfHigKYAqAJQNPgwYNTvuFsEu8BtbUxH3GEeNVr10b/byKe+gsLA1x2QzmXVwV4//2ZjztOnuCSJYnJkUu8LJub5PI+89WAp5PnbjmWK1cyV1ZK6LtLF+axY5l37kxOv1euZH71VeYVK9ISRcki8cqeW/VBJg34D6IY8HsizukOoKv1/SQAnySStpc89Wgk+nBmzZKcfvDB9scSrUCam5n79A1ypyEBBoIMSIueSGZoSyYtJXPk8hnkqwFPJ8/dqmCPOkpGh9gh8NJS5ltvTUzGTZuYjz9e3tf36CGfxxzDvGFDWiIpWSSyHNm/3ZpuN6ch9Cj/WQmgb7y0vW7AE6k4Pv9cWt5HHhk9dJ5oBXLeec4KIsioFEPeqZPsb2nJn1aul9EWePYNeC7zfMsW5v/8RzqVOt9hA8zjxrWXMVpL/KijpIe8879FRTLFsmIGti2wOxzbhrujNTCSIZMGPCNjRdkAA55IxfGjH4kyvv9+8v91ElZBVAYYV5VbRlz2TZuWxo1kEHUqUifZvMtXA54K6ZTLTz+V1vXo0RIBizTc9lZQwHz44cy33SYRNOb2Ff3ixUEuLY3+/7Iy5vnz3b1vt1H9FiJb3G4ZbptMDyNzfawou6DouS5c//mP5PDUqe2PJRv2Cx9mEmqB20NO+vcPf+fmFYXSsH7qJJt3+W7AMz1/eUODdEItKBA9rK1lvvFG5nnzmF95RTqVdu0q4fORI5mvuUaMvG3MTz6Z+a23wiv6mbMDUVvvtgGfPTv1/MgGqt/hZKr+zagBz9Tm1WUHE3lIu3czf+Mb0qM02nCQZB/0WWeFKg6nIe+8b4AvuEAM+VNPeU+hvOZQmIS2wJMjmR7BycyW9sknzBMmiM717s187bXMq1a1P2/rVua33pJWs3O8+fLlMibdHio2dixzQ4PI8dZbbVy6X8gZj9z69w91UvUiqt/ZIS8NeKYKVyIVxR/+ILn7zDPuXHPFiijv2SoDXPCLcn71gwDvu6+8S1OF8ga5eA75bsDdHpO7fbu0sDt3lpb1jBnMX32Vunxffy31Qr9+or+TJzNPuXEB09XljMr6sKiacysvjz8FrOJv8tKAZ4p4FcVnn0lY+/jj2yteOhV7bW2EAacgn/EzSevOOzlsSJmSW3IRCcl3A27TkY4lqn8NDcwHHSQ6dfbZzK2t7sn35ZfMv/iF/f48yOixglFRz7iqnGlQoJ0BLypi3rzZvesr5hFLt3U1shSIt+LOL38JbN8O3H23rCLmJJ0VafbsidjBhN7bRY4LLgCKi4EHHkg62YzA7MvlNRPGr6uumUA0HbPLI9DxMqG7dwM33AAcdhjw1Vey1O8TTwDl5anJYl83GAzu/fzwqwZUVQVRPGiRnPTlMKC1FoVzXsSl46tRul8DnKu2du0KdO+e2vUzSb7ruBdQA54EiRTYd94BHnkEuPpqYMSI9sfTqdgvuggoKwv9LisDzjlH5Fq+owHjxzMeewzYsSP3yuXT5TUTxmvrBucTzvW5A80BBJoDaGhpiFseP/4YOPxwYOZM4LzzgCVLgBNOSE0G+xna1529cHbY5z9aZmPnxBOBigag5woAwJ6WUfh8xyqUnD8Opd95CCWljJ49gRdeaN8Q8AL5ruOeIFqz3CtbthY9SJR4YdF4HdfSJRhkvvtu5gMOkOs891y4XHf/9QO251rPdWc2fRefOqk+O2gIPYxAc2hhoI6G9QSDzA88IL2+e/dm/vvf07uuczET+7p2pzn7891327jz8IC8+76qnDFgAReUbmKA+eAxzdzvum/ycwvq+e1PnZODML/2miw1evfd7ob1U71P1fHsEEu3c67MHW3ZWnYwUeIV2I46rmWysNtp79kT5KFDZYpVVS5zSfXZqQEPJxgMLc0bKy/XrWM+/XTR22OPlQmR0iXQnNhyorNnM3cuCXLn4QEu3b+ee91cztfcsoJLS4PcvecuvvXepXvHFL+wMMAjDwxy164yGqWkRLaZM9OXV/E+eW3As2HMWltjd1xjzl6LeOZMearLl2f0MjFRxyGcbOaHGvDkeOYZGaZVXMz829+GFiBJlVSGqG3bJhPD7NgR+u+cNxdyVZX0Rj9l4uf82gcNXDi1nAsGt++pXlbGXFeXntyK98lrA54uiVTCZ50lw00+/jj1NNyguVk89OnTs3tdm1yH7r1GNvNDDXhitLYyT5wotd+3vsX83nvupJvMs37jDZHhuOMkHL5lSyiNAb8ewPcGHuRp04JcUMBcWRnkosoloXB7xJCzb3/bHfkTRZ309tivNysqZDGbiROZ1693L3014GkQTzFfeUVy8qabsixYBLZinXhikCsrZTKJbBtUVe5wksmPdPNODXjH7NzJ/LvfSaSsc2eJVu3alX66yba8p02TlrOzFT1sGPPGje3fn983dwn362e3uIOMioZ2hry0LLu65mUnva1N+h+tWyeNmRUrpFH14YfMH3wg24cfMi9dKhP0fPqpnLd2reT/V1/J+uzJRmPuuy/8mRYVMY8a5d74/Vi6bc9P7kmqqqq4qakp12KAmdG4phE1lTXtehVv3w584xvSS/S994CSksT/6zZ2j9fLu72JaRePwAsvACeemL3rZ4sdO4ClS4EPPgCWLQM+/RRobQXWr5ehPzt2yHmdO8vwm/JyYMgQ4KCDgFGjgKqq9s/JC9jPr25SHWoH1ib9fyJawMxVGRDNVbKt13v2AE89BcyYIeXl+98H/vhHYP/93Uk/mee2bBlw6KGhMmpTUAAUFgKlpcDUqYxjzm0EM2PC0xNwWe/ncdPFVdi92z6bgcpG+TxrAvq8NhcvzKOs6Xg267TQNYHmZhkpsHw5sHo10NICfPYZsG4dsHEjsHkzsGWLe9fs1EmG5hYVtd8KC+W4vX30EbBzZ/j/CwpkJFJRERAMSjncs0eGKu7eDezaJaMennsuviyxdDtvDHimCt20acBttwGvvQYcd1z74zZMCGgAAB2WSURBVOlWyslg3+OofjUYNIhw5JHAM89k9JLtru12/gaDwPvvA/X1QCAANDWJ4W5rC51TUSFb//5Ajx5iuIlEob78Uoz7p5+KggNivL/zHWDcOGDCBKCy0jVxo5Jo3qSbh2rAw9m4UYZ0/vnPUukfeihwxx0yNMyNImo/r+qKasxvnZ/Qc7vrLuDaa6XyjkVZGfDkk8Cpp0r63+xbjX7fmo+tH9UAcKbPKB7WiB/+iDGvy3jcfOzNuPBbF2bMqGbLcDPL83rnHWD+fGDBAhnSt3Vr6JxOnUTny8tF7/v0AXr1Arp1k3HzpaWi57bBLSiQzU6/rS1kVNva5Hk4DWusT9sAt7WF/tvWBrz8cnsD3qkTMGYM0K+fXLuoSPbZMhUXAyNHApdfHj9PYup2tGa5VzY3Q22phn06CmsuWsRcWMh8/vmp/T+TXHmlyPbZZ9m5nns9/SXM9fvfM592mizFaoel+vZlPuEECUE+9RTz4sXSCSjRdFtaZOjdlVeGZtkqKGA+6STZn24npljEyxu3ygh8FEJ/+20JbzrnFI9HMMj80UcSzjzppNASnYcdJh293H6+qZT53/1OwvexVjCzt1/+Mvw6vW4p587DA2HrIRAxjxjBvGVL+2FrmahvMhk6X7+e+bHHZNnkiorQPXbtKlNEX3458733St+BVauSKxfZ4MorZVSA8xn26SNDi90glm7nXJk72tw04KlWkrEK7e7dsl5v//7MGza4JqZrfPCBPN3bb5ffmXYk0kl/927mV19lvvRS5iFDQgqw777MF17I/PDDzMuWuT8f9AcfSGc/u8I44AAZ2uOW0tnEe0fqVsXoFwPe1haqDIuLmffbT4Z4TZ4sFfnUqcw33CCO3M9+JvvHjGHu2TNUdoYOlUp10aK0srRDUinzK1Ywd+rUsfEuK2N+6KH211m8uI1P/EmAhwwN8pAhsjZCp04yLevWrZldzjLV++2INWtk6O2RR4YWaurXTzoE33uvOOiZcqrd5uuvxWksKRGnY8AA5sZG99LPewOeKrEK7S23SO49/XSOBIuBU96jjmIePlyUIJcrs8Xif/9jvuyy0OIOpaUyJve++6T1lS1275bJb771rZAhf+YZ9x2GyGeQyrCjjvCLAd+zh/nNN2VylWuukR69hx0mZblHD+kgBEiEqXt32T92LPP//Z/856OPMrv4R7rOauSqgkRyH0VF4rCcdVZ0w2WXH9s4b9gQ5AsvlDSGDJEoUibWpXbTcG/ZIg75sceG1lI/5BBxyBobzTHYsWhpkfUo3L4PNeAWiRbGjs577z1RtokT07+OTVsb85/+xPyDHzBfdVXqQxCcRuLxx+UJv/Zabldmc7Jnj8x0NWaMyFZSIvc8d654sbkkGBQ57PD62LHSCnAv/XCDbVeybjlVfjHgiZDLQQ7pOMO7donj4TTggwbJq6CNG2WLRTAY5OcXBPj6++q51y3l/PZKuf6bb4bK7AkniBNz+BFBPvj4AN/6aHtD/vHH4jgPGcK8zz7SW3rWrNj6l67zHwwy//e/Ek3r0kXkHD5cVnr78MOUksw71IBbJFoYY523YwfzN78pofOOjGyyhf6880LDEIqKmCsrmTdtSuivYTgN9fbtMjXkD34Q/bgbJO4QSav2gAPkHvfbT95zf/GFK2K4yu7dzH/+s7zD6tRJWoFuTo0b2ZJy61nkkwHPBW5FTC67TAxZUZHo/LPPJva/u+8Wh7dL1yCX7R/gykFt/PwCkWPXLua77pLwrdM5KC0L8q+fFGdxwG/K+exf1HPxsAAXFoUvW9q1q7yCCDiqq3Tvd9Mm5nvuYT70ULlGly5ixN9+W5dHTRY14BbptsCvvFJy7R//cOc6zDKxRGTHltJSaZEHgzI2MdVOG3ZnNnuKSDdD6Yne45IlzEccIfd10EHy2sFrnVCisWED80UXidzDhkkkww0yFQ1RAx4bN/LczY6adXXMv/514u9JV61q30mqYFCAO08LdwTPPTf8HCA0tfIPrw1w5+GRk8G0hU0K07Ur85Il4WH4ZO43GGR+5x3p2FtaKtcfPZr5/vvTW0s931EDHkEshe5I0V96SXLs0kvdlWXp0vaec2Eh8/jxEuIqLJTjv/xl8h2sli+Xd02ZmJktXoW2e7dMblNUJK3Z++93v4NYNnjzTYkYAMwXX+zdikgNeHv27Anyi4sC/Oay8FByMq1Kt/sqpMKsWeEThdgTuxQNCze0V1zR/h17QYF09JNe+UEx2PZkMKNmMa4aIJ9oYwwMcO2Z9bzPneV8x+P1fOcTAV63Lv79Njcz33EH88iRvNcRmDKFecGCLGROHqAGPILIMGa895KrVokROuSQ+EOXkjWSu3ZJyNzu1GGH0SNb5WVlHQ9ZiyXDaafJEKzt21OXMd41Ilm1ivnww0Xus8+WmZFMZts26ZtAxDx4sPSa9xpqwMP55BPmAaMDjKvLuWBIPV9xZ0i/Zy2YFfc1RmSHsFzOPDZ3rvQ6j2xd9+0bvmjLihVB7tFDnP6CAmkFT5zY/r37XkOOtnAjflU506B6LhoW4G7dg9y9u7T8f/xjmcnOSXOzhMjHjg3VXUccIc6GPTWs4g5qwCOIVM6OFHrHDubqaukpunRp/LRTCbV9/LG8HyYS7zWyRW5vnTvLtH/JyPD66/Jf59CUdMKB8Yz/v/8tFUu3bsyPP5508p7mnXdk7C3AfMklzF9+mWuJQqgBD2fEAUHGwHppbSLIZWXMb70VvSNhtGU/szG2OlF27JBonLN1XVbGfNttctzZIHm2KcAzZgT5xhtD9dXo0dHrk73GfEBTu3C6cyspYT75ZOa//lWG7tmd5gDmAw9knjEj9joQSvqoAY9BvPBYMCgdzADxgpNJMxWF37VLeqLak1BEbj16JPYu1ilDMCiRg0MOCQ1vSEfGjoz/Qw+Jt3/AAYk5OyZit8YLCqQH8Qsv5FoiQQ14iF27mDEwYL3vDbA9pvyuu8LPi+XIz1owK6ElQbPJ8uUy/KqwUCY4uvnm9voca/iYPUTSja20lPn732e+804ZsqdkHl8a8H/8Q95JRoZ23OTGGyWXsrnu7q5d0d53hZQnFaWxh5RFOiGpLLYRzdkJBqVCAZi/973UetCbRiAQaomceaY7a0lHI9ExpZkw4ABOALAUwDIAU6McJwB3W8ffA/DteGlmqwXetVswrEVZVMQ8Z070cyPLdi7fd6dDLEM+8Sx5vx2tdZ3oVlgoIzJM7MdiOr404Pba14cfnpkpQ//0J0n//PMTG/bgZgexn/0s1IvT+V58zJjU0tu9m3n//duvkBNoDkRtaWzaJMO8Roxg7t4jyL0PDfCJU+q5723tW97BIPPPfy4ynntuZh0qr7FzpzguJSXidM2Y4d6Qs127pAyOHJnYcDu3DTiATgCWA9gXQDGARQAOijjnJAD/tAz5GAAN8dLNtAG39fBnV4Qbq9JS741++P3vZTnQv/3N3XQjDfnUObOYro7W+zx22DxyKyqSnvNKbLI9wiTnRrqjLRFFf/JJUcyKCub//CfF3InCn/8suXPaaYkbJDeHaO3cKZ1PSkokbF5WJu/hP/88sf8Hg9Jycxaohx+We3ruOed57edRfvPNIHfrJmNInT1WaVA9d943wFMuDlppMzc0MB9/vKR7xRXmz6SUKitWSCsckDkCfv3r1N+P79ghHYHsnu9HHSVTycYjAwb8MAAvO35fC+DaiHPuAzDZ8XspgPKO0s20Abf18JHXA3uHXhUVSU9sL+EcPlpcnJmx0bb+79nTxvuNDXDBYGfv8yi90Dsw5l26MNfXuy+jH8jEDHhOfGvAmWVKzuHDZdKNW29Nz8sOBkNh81NPTa41mQnva/Vq5hdflLHUibBypUzFWFwsHeL2O1oWQgg0B3j3bsmngw4KXwfZWfj63lbOxcPrw4ea7PXag3s7z5x5pvTGtnu3FhbKBDd+fe+dKO+8I68Q7KE0F1/M/K9/xQ87trUxv/su89VXS2clQN5bvvBC4hV7Bgz4mQBmOX6fB+CPEefMA3CE4/frAKqipDUFQBOApsGDByd2Qyni1MPnnmP+7ndlDnU3J+Nxgy1bQs75sGGZv15rK/OgwTL0LKwFbhnx4tpZXHCN9NiPNOSFhe6ub+03Ikc1uW3IM2rAM/GejJP01DdvlhYrIGHmRA0es3QOufBCmYf76KN5b9g8UeOd6YVCEmXVqvDVu+wepsXDAvzAAyLbs8/K/rvukhW9Zs4MdcIKBoNcNS624Y63EcnMb4n0kvc7jY1ShuwWYM+e0vHnggukZd25swwdnDyZedw4abXbFeXpp8swtWSLUwYM+A+iGPB7Is55IYoBH91RuplqgXtFD5NhxQqZGzxbwyy/+ELeY3fvLluPHswlpUHe7+gAz5nTxv9cHOB9j6oPC7d36RrkkSNVrzsisuy5PQQxYwY8U+/JOAVFDwZlmEPv3tIav/RSWfGmI5Yta7860DnnJFd5uhk6T4fzzgsfS+7cunaV3tPBoMyXbE/jWFAgobHp052zPQWTMtzOrXNnWS1KEbZskbnff/zj0JSSkduwYVLmHnoovYo830PoXtFDE9i5k7mpSV47Ll8efqytLcj31AV40lX1XHZDOf/+b4G8fTUWj3hOo1tOZSYNeEaUnNNQ9PXrxXh36iQtmsmTZR7uyJ7Ry5dLB5LICrVv3+Su5xXPP1bPdTvs/eKLct6rr7Y/Xlgo/Qm6d0/OYEfbMhwhNZb772//jIhkxio3yIABLwSwAsAwh3N+cMQ5J0c4543x0tUWuBl4YQY6r5OI0+hGuYyl2wVIn0oAzY7fLda+ZM8BABDRFCJqIqKm9evXpyRQ377AH/8ILF0KXHYZ8OKLwBlnAL16Ab17A0OHAl27AsOHA+++2/7/O3cmdh1mRkNLAwCgdmAtiCgled1i167Yx4LMWLihAcyM/v2B4uLw40TA9u1iVtIl0fzLN778EtizJ3wfM7BxY27kiQcz7wFwGYCXAXwI4Glmfp+ILiGiS6zTXoQY+WUAHgDw05wIC4CIPKGHfsHOz/mt8zF+zng0rmnMtUieo6ayBnWT6lBTWRPznMY1jRnLPzcMeDRtiTQDiZwjO5nvZ+YqZq7q169fWoINHw787nfA+vXAW28Bt90GnH02cNRRwJQpwD33AH/+c/v/TZ6cWPqZfDCpUF0d+1jbgEb8fq3IOnIkUFkpRhsACguB2lrg0EOBYDA9GQoKgCOOSC8Nv3LaaUCnTuH7unQBzjknN/IkAjO/yMwjmHk4M99i7buXme+1vjMzX2odP5SZm3IrseI2tpGqrqhGQ0uDHXnJa5JpvCVi5NMSJJ0NHgyhJ8t998m0nyUl0skoXo9hr4aW3nij/dhxe4jKueeFh3HWrmU+5ZRQJ6s775Se0Pvum174vLRU3q0p0XngAckje47pK690r2cv8nQmNg2dZwftYxAi23kRS7fdMOAZeU/GWTTgyeLlgjxnjvR4LiuTrXPn8B71kZXd119Lr2eAuaaG+ac/bb+ISrT35QMGiAEqKpJ9JSVimB59NIc3bwhbtkgv9Y7Wk0+FfDXgXtZHP6GOUohU8iKd/MuYAZe0cRKAjyG90adb+y4BcIn1nQD8yTq+GFHGiUbbvGLATZtmcdcuGXs8b177iV+iVXZ79jDfe29ire+SEhkGtX699OCfPp35hz9kvv32xCeZUTJDvhpwNSzZRfM7NdJxNDNqwDO15cKARzPW9kxl9kIHJnv6HSnfnj3SMnz+eeYbbmCuqpLWePfuMgytb1/mX/0q+nrYqtS5J18NuJJd8jnikU49l4kWOMkxb1JVVcVNTdnpE8PMaFzTCGbGhKcnYOYxM3H9v67HzGNm4ro3rsPNx96MC0ZdgPmt81FTWZM3PV3XrpWttBQYMSK8ExZzqCMHAEx4egLmTpwLIsqrPPIKRLSAmatyLUc8sqnXivvYdWU+6nhDSwPGzxmPukl1qB1Ym7XrxtJtNeAW9oOxDVB1RTXmt87f++m3wpqOEjqdnVOePAUAMG/yPBDRXgdIDXn2UQOuZJN8NOS5uudYuu3GMDKjsVuR1RXVe72q2oG1KCgoCPv0WwFNZwic/V9ADPe8yfP25lvtwFrUTaoDAE8NsVP8h627Xm6E+BmvDaPNBm7MNeBmuc17A24Xwvmt831pqGORytjOaM7OmEFjMGbQmL35Zhdw25Dr2NH22PkYDAY1b9IgHw2Il8jo+GYf42a5zXsDnq+FMJlZlmyDY79mSMTZ0Vmc2hOZj7MXzta8SYN81V2vkE8z37nZanaz3OatAXd2wMqXQhiNmsoazJ04NzQswUGkwQGQdMHTSjaE89VD3aQ6XDDqAo1SpEE+GRCvki+vMdxsNbtZbvPWgGv4TSAiEBEmPD1hr7G2Q7uRhtsOjSdT8OzCCiAvFD0a8fpZaJRCMZV8qUe92hDJWwPu1QeSC+y8ABAW2gVSN9yR5IuiRyNePwudazo58qXVZwL5Uo9mItrjSjmONjjcK1smJnzQCUdik8k53vM53xO993QnyECeTOSSzxOJKP4hmXIcS7fzrgWezy3BeNheZiaGzuVrKJ2TGDeaL62ZdNF88h6sUZGkcaMc550BV+XPLfnmQCVzv9opKzE0n7yHn/U6U86JG+U47wy4Kn9uyTcHKpX71daMYhp+1msvOyd5Y8C1UvQG+eZApXK/Xq4wFCUaftZrLzsneWPAtVL0Fn53qNK5Py9XGLnE72VG8SZedk7yxoBrpegt/O5QpXN/Xq4wconfy4wfUCcru+SNAddK0Vv43aFy4/60MgzH72XGD/jJyTJB//LCgJvwIPINvztUbtyfnypDN/B7mfEDfnKyTNC/vDDgJjwIxR+46Sz6qTJU8gM/OVkm6F9eGHATHkQ+4sfIiFcXPVAUJTlM0L+8MOAmPIh8xI+REbedRT86Oamg+WAW+ryyQ14YcMWb+DEy4raz6EcnJxU0H8xCn1d28LUBVy/Q2/gpMpKpsuZHJycVNB/MwuTnZZLd8LUBVy/QDExSmFhkqqz5yclJB80HszD5eZlkN3xtwE32AvMJkxQmFlrWFMUfmKTLvjbgJnuB+YRJChMLv5U1IupNRK8S0SfWZ68Y560kosVEtJCImrItp+JtTIyumaTLaRlwVXLFDUxSmDxiKoDXmXl/AK9bv2NxDDOPYuYqt4Uw0QAoIfwQXfMy6bbAPaHkkajSm4mJz81EmRPkdACPWN8fATAuF0KoATAbk6JrJupyugbcE0oeiSq9mZj43EyUOUH2Yea1AGB99o9xHgN4hYgWENGUWIkR0RQiaiKipvXr1ycshEkGQGmPSdE1E3WZ0vE2iGgzM/d0/N7EzO3C6ET0KYBNEGW/j5nv7yDNKQCmAMDgwYNHr1q1Kmm5mBmNaxpRU1ljRMFRBBOfm5dkJqIFyUS4iOg1AAOiHJoO4JEEdbuCmVuJqD+AVwFczsxvdXTdqqoqbmrSN2n5hJf0JBZeljGWbhcm8MeOlDxRDncqORF9FEvJLeN+PyCKnsQ19mJ7fYpZmPjcTJTZhpm/G+sYEX1OROXMvJaIygGsi5FGq/W5jojqANQA6NCAK/mH3bqtm1TnWX0xUZfjhtCZ+bvMfEiU7TkAn1vKjUSVHICt5IpiLCa+L0uS5wH8yPr+IwDPRZ5ARF2IqJv9HcD3ASxx4+J5kL95hZdfhZhc1tJ9B55TJY/E5AehhDDhOZr4vixJbgfwPSL6BMD3rN8gogoietE6Zx8AbxPRIgCNAF5g5pfcuHge5G9e4eV34SaXtXQNeE6VPBKTH4QSwoTn6OUWhRsw80ZmPo6Z97c+v7D2tzLzSdb3Fcz8TWs7mJlvcev6fs/ffMWLzrnJZS2tTmyZJtnOLl7uhKAkjj7H1Ei2E1uu0E5s+UtDS4Pn34V7kVi67auZ2LwcplESx9mZxGveuhdbEIpiCl5q7fpBl31hwP3wIJT2eDGU7kWZFMUUvNTI8oMu+8KA++FBKO3xkrdu40WZFMU0vNDo8oMu+8KA++FBKO3xYijdSy0IRTEVLzS6/KDLvjDgfngQSmy8oOxeaDEoil/IZaPLT7pstAH304NQYuOFCIsXnAhF8QtEhJrKGjSuacx6/e0nXTbagPvpQSix8UIo3QtOhKL4iVzV337SZaMNuJ8ehBKfXCm8jktXFPfJdv1tR2wB+OaVq9EGXN995xe2wldXVGe1Ja6RHkVxn2xH1vyox0YbcCW/sBV+fuv8rCii7bFXV1RrpEdRMkQ2DCszg5kxd+JcX+mxGnDFOGoqazB34ty9Spkp7Iplfut8jfQoSobIRmStcU0jJjw9AUTkKz1WA64Yh62EE56ekBGvXVveipI9IiNrDS0Nrhlyv+uyGnDFSDLhtdvKbi+4oC1vRcketk4DSDukni+6rAZcMZJYXnswGEzaoEcqOwBfeuuK4mVsna4dWBvmnCej0/mmy2rAFaOJ9NpnL5wd1aAHg0EEmgMINAfC9kVTdrsS8Zu3rigmEOmcO3U6lg7nrS7bHYG8uI0ePZoVJRGCwSAHmgPc1tbGgeYA16+u5/LflPOsBbP2fva5ow/3uaNP2L7y35Rz/ep6DjQHOBgM5vo20gJAE3tAb+NtqtdKIkTT6Vg67DddjiSWbhfm2oFQFDdwjimtHVgLZt4bhjuk/yGorqjGwf0OBiCtdnvfIf0P0QlaFMWDRNPpeZPnAWivw/mqy8Qenke8qqqKm5qaci2GohgBES1g5qpcyxEP1WtFSY5Yuq3vwBVFURTFQNSAK4qiKIqBqAFXFEVRFANRA64oiqIoBqIGXFEURVEMxNO90IloPYBVLiXXF8AGl9LyAn66Hz/dC5C7+xnCzP1ycN2kcFmvAS0/XkbvxR2i6ranDbibEFGTCUNsEsVP9+OnewH8dz9ex2/57af70XvJLBpCVxRFURQDUQOuKIqiKAaSTwb8/lwL4DJ+uh8/3Qvgv/vxOn7Lbz/dj95LBsmbd+CKoiiK4ifyqQWuKIqiKL5BDbiiKIqiGIgvDTgR3UREa4hoobWd5Dh2LREtI6KlRHS8Y/9oIlpsHbubPLomHRGdYMm+jIim5lqeRCGilVb+LiSiJmtfbyJ6lYg+sT57Oc6P+pxyARE9RETriGiJY1/SsptSxryKn/UaMFO3TdZrSx6zdTvaIuGmbwBuAnB1lP0HAVgEoDOAYQCWA+hkHWsEcBgAAvBPACfm+j6iyN/JknlfAMXWvRyUa7kSlH0lgL4R++4EMNX6PhXAHfGeU45kPwrAtwEsSUd2E8qYlze/6rUlp5G6bbJeWzIZrdu+bIF3wOkAnmLmncz8KYBlAGqIqBxAd2auZ3kajwIYl0tBY1ADYBkzr2DmXQCegtyTqZwO4BHr+yMI5XnU55QD+QAAzPwWgC8idiclu0FlzET8kOd+0m0j9BowX7f9bMAvI6L3rBCJHQKpBNDsOKfF2ldpfY/c7zViyW8CDOAVIlpARFOsffsw81oAsD77W/tNuM9kZTeljHkdP+o1YEaZj4bf9BowSLcLs3GRTEBErwEYEOXQdAB/ATATUrhmAvgtgAsh4Y1IuIP9XsMUOaNxODO3ElF/AK8S0UcdnGvyfZpexnJKnuo1YJasTvJFrwEPljNjDTgzfzeR84joAQDzrJ8tAAY5Dg8E0GrtHxhlv9eIJb/nYeZW63MdEdVBQmefE1E5M6+1wlDrrNNNuM9kZTeljOWUPNVrwIwy3w4f6jVgkG77MoRuZbrNeAB2D8PnAZxFRJ2JaBiA/QE0WmGSLUQ0xuo9+EMAz2VV6MSYD2B/IhpGRMUAzoLck6choi5E1M3+DuD7kGfyPIAfWaf9CKE8j/qcsit1XJKS3aAy5ll8rNeAgbrtU70GTNLtXPYAzNQG4DEAiwG8Z2V6uePYdEjvwaVw9BQEUAUpfMsB/BHWLHVe2wCcBOBjS87puZYnQZn3hfTeXATgfVtuAH0AvA7gE+uzd7znlCP5nwSwFsBuiLd9USqym1LGvLr5Wa8tWY3SbdP12pLHaN3WqVQVRVEUxUB8GUJXFEVRFL+jBlxRFEVRDEQNuKIoiqIYiBpwRVEURTEQNeCKoiiKYiBqwBVFURTFQNSAK4qiKIqB/H9cqDnFXYvT7AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make a plot\n",
    "sizes = cluster_sizes.numpy()/70\n",
    "\n",
    "# time points as numpy array\n",
    "new_times = new_y.numpy()\n",
    "\n",
    "# means\n",
    "tumor_mean = new_tumor.numpy()\n",
    "metacluster_mean = new_metacluster.numpy()\n",
    "\n",
    "# plus or minus standard deviation\n",
    "tumor_std = torch.diag(cov_tumor).numpy()#**0.5\n",
    "metacluster_std = torch.diag(cov_metacluster).numpy()#**0.5\n",
    "\n",
    "# make plots\n",
    "fig, ax = plt.subplots(1,2, figsize=(7,3), tight_layout=True)\n",
    "\n",
    "# plot tumor\n",
    "ax[0].plot(new_times, tumor_mean, \"b-\")\n",
    "ax[0].scatter(y.numpy(), tumorburden.numpy(), s=sizes, c=\"b\")\n",
    "ax[0].scatter(new_times, tumor_mean + tumor_std, c=\"g\", s=0.3)\n",
    "ax[0].scatter(new_times, tumor_mean - tumor_std, c=\"g\", s=0.3)\n",
    "ax[0].set_title(\"Tumor burden\")\n",
    "\n",
    "# plot tumor\n",
    "ax[1].plot(new_times, metacluster_mean, \"b-\")\n",
    "ax[1].scatter(y.numpy(), metacluster.numpy(), s=sizes, c=\"b\")\n",
    "ax[1].scatter(new_times, metacluster_mean + metacluster_std, c=\"g\", s=0.3)\n",
    "ax[1].scatter(new_times, metacluster_mean - metacluster_std, c=\"g\", s=0.3)\n",
    "ax[1].set_title(\"Metacluster proportion\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 318,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(285.)"
      ]
     },
     "execution_count": 318,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tumor_model.lamb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}