diff --git a/examples/ode/SEIR model with ODESolver.ipynb b/examples/ode/SEIR model with ODESolver.ipynb new file mode 100644 index 00000000..2eee65c9 --- /dev/null +++ b/examples/ode/SEIR model with ODESolver.ipynb @@ -0,0 +1,489 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## SEIR Epidemiologic Model with TorchTS ODE Solver\n", + "\n", + "This example solves a compartmental model used in epidemiology, known as \n", + "[SEIR](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SEIR_model) \n", + "model, using the TorchTS ODE Solver.\n", + "\n", + "The ODE (ordinary differential equation) system, in generic form, is\n", + "\n", + "\\begin{equation*}\n", + " \\frac{d \\mathbf{A}}{d t} = \\mathbf{F}(a_n)\n", + "\\end{equation*}\n", + "\n", + "or in uncollapsed form\n", + "\n", + "\\begin{align*}\n", + "\\frac{d a_1}{d t} =& f_1(a_1, a_2, \\dots a_n) \\\\\n", + "\\frac{d a_2}{d t} =& f_2(a_1, a_2, \\dots a_n) \\\\\n", + "\\vdots \\\\\n", + "\\frac{d a_n}{d t} =& f_n(a_1, a_2, \\dots a_n) \\\\\n", + "\\end{align*}\n", + "\n", + "In the case of the SEIR model, the ODEs are:\n", + "
\n", + "\\begin{align}\n", + "\\frac{d S_t}{dt} &= - \\frac{\\beta_t I_t S_t}{N}, \\\\\n", + "\\frac{d E_t}{dt} &= \\frac{\\beta_t I_t S_t}{N} - \\sigma_t E_t \\\\\n", + "\\frac{d I_t}{dt} &= \\sigma_t E_t - \\gamma_t I_t \\\\\n", + "\\frac{d R_t}{dt} &= \\gamma_t I_t\n", + "\\end{align}\n", + "\n", + "Here, the the compartment $S$ (susceptible population) represents the first variable $a_1$, and $f_1$ is denoted by the right-hand-side of the top equation.\n", + "The coefficients $\\beta$, $\\sigma$ and $\\gamma$ (either constant or time-dependent, still to be implemented) are optimized using PyTorch." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In order to solve various different systems of ODEs, the following quantities must be somehow parameterized and passed to the solver:\n", + "- The equations $f_n$. \n", + "- The coefficients (and a flag denoting whether they are time-dependent)\n", + "- The data used to train the model. \n", + "- An optional output modifier which takes the numerical solution and brings it into a shape consistent with the training data such that the loss can be calculated, as explained below.\n", + "- Other user-controllable parameters including the temporal discretization (time step),\n", + "optimizer, scheduler learning rate and loss function." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Specification of the variables, parameters and functions\n", + "\n", + "Variables, ODEs and initial values for coefficients are passed to the function during initialization:\n", + "\n", + " ODESolver(inivals, cfuncs, inicoeffs, dt, time-dependent=False, \n", + " solver='Euler', outvar=None)\n", + "\n", + "##### Variables\n", + "For working with the solver it's easier and more intuitive to assign actual variable names to each quantity.\n", + "They are provided as the keys in the dictionary passed to \n", + "the positional argument `inivals`. The values of the dictionary provide the initial values assigned to each quantity.\n", + "For the SEIR model, one could use, for example:\n", + "\n", + " inivals = {\"S\": 0.95, \"E\": 0.,\"I\" 0.05 : ,\"R\": 0}\n", + " \n", + "(Here, the population in each compartments is normalized by the total population size).\n", + "\n", + "##### Functions\n", + "\n", + "A function specifying the right-hand term in each of the system of ODE's is passed to the solver as a dictionary in the\n", + "positional argument `cfuncs`. The equation pertaining to each variable is stored under the key representing the respective variable.\n", + "Each function receives two positional arguments, `cvars` and `coeffs`. These will be dictionaries containing the \n", + "system's current variables and coefficients. As an example, the function describing the ODE for quantity $S$ would be defined as:\n", + "\n", + " def fS(cvar, coeffs):\n", + " return (-coeffs[\"beta\"] * cvar[\"I\"] * cvar[\"S\"])\n", + "\n", + "##### Inicoeffs\n", + "\n", + "Initial values for the coefficients are provided in the dictionary inicoeffs. Each coefficient must be present,\n", + "and the keys in the dictionary passed to the solver must represent the names of the coefficients that will be optimized through data.\n", + "\n", + "In the SEIR example, one could use\n", + " \n", + " inicoeffs={\"beta\": 0.50, \"gamma\": 0.20, \"sigma\": 0.20}\n", + "\n", + "##### Output quantities (and time skip, still ToDo):\n", + "\n", + "By default, the network returns a time-dependent value for every variable and every discrete time step resolved \n", + "during numerical integration. Depending on the model and data, a training value may not be available for quantity.\n", + "For example, only data on the currently infected and susceptible population was typically be available during the Covid-19 pandemic, but not on the exposed population.\n", + "(Alternatively, one might only have data on cumulative reported cases (`cumsum(I)`), not currently infectious cases.\n", + "Handling such cases will require functionaly that is not yet implemented.)\n", + "\n", + "The keyword variable `outvar` designates the names of the output quantities that are present in the data and used \n", + "for computation of the loss. In addition, it indicates the order in which they are present in the training dataset\n", + "(format described below).\n", + "By default, `outvar` is the same as `variables`. In the case of the compartmental model, one would use\n", + "`outvar = [\"S\",\"I\",\"R\"]`, as no data on the exposed population $E$ is available." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Training the network\n", + "\n", + "The solver is trained using \n", + "\n", + " ODESolver.fit(train_data, num_epochs=100, lr=0.001, optimizer= None, \n", + " scheduler = None, loss_fun=torch.nn.MSELoss()):\n", + "\n", + "The PyTorch tensor `train_data` is assumed to be of the shape `(nt,nvar)`, where `nt` is the number of time steps used for training and `nvar` is the number of output variables (consistent with `len(outvar)`). The sampling interval of the data is expected to be the same as the timestep `dt` passed to the solver during initialization.\n", + "\n", + "By default, the value `None` is passed for the optimizer and scheduler, \n", + "and the network uses \n", + "\n", + " optimizer = torch.optim.Adam(self.coeffs.values(), 0.001) \n", + " scheduler=torch.optim.lr_scheduler.StepLR(optimizer, step_size= 1, gamma=0.95)\n", + " \n", + "To learning rate can be changed using the keyword argument `lr`. If a custom optimizer is provided,\n", + "the optimizer's coded learning rate is used. A warning is issued if the user tries to set both.\n", + "\n", + "### Predicting\n", + "\n", + "Predictions are made using \n", + "\n", + " ODESolver.predict(nt)\n", + " \n", + "Where `nt` represents the total number of time steps in the prediction (starting from the same origin time as used in the training data)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "cpu\n" + ] + } + ], + "source": [ + "import torch\n", + "import numpy as np\n", + "import torch.nn as nn\n", + "from scipy.integrate import odeint\n", + "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n", + "import os\n", + "os.environ[\"CUDA_DEVICE_ORDER\"]=\"PCI_BUS_ID\"\n", + "os.environ[\"CUDA_VISIBLE_DEVICES\"]=\"4\"\n", + "#device = torch.device(\"cpu\")\n", + "print(device)\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from tqdm import tqdm" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If torchts is properly installed, the module should be in the path already. This code adds the\n", + "module `ode` to the current path assuming the example is called directly from the repository's tree:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " from torchts.nn.models import ode\n", + "except ModuleNotFoundError:\n", + " import sys\n", + " sys.path.append(\"../..\")\n", + " from torchts.nn.models import ode" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Defining all the functions. Here, the coefficient $\\beta$ is assumed to be normalized by the total population $N$ already. Population sizes are assumed normalized by the total population as well." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def fS(cvar, coeffs):\n", + " return -coeffs[\"beta\"] * cvar[\"I\"] * cvar[\"S\"]\n", + "\n", + "def fE(cvar, coeffs):\n", + " return coeffs[\"beta\"] * cvar[\"I\"] * cvar[\"S\"] - coeffs[\"sigma\"] * cvar[\"E\"]\n", + "\n", + "def fI(cvar, coeffs):\n", + " return coeffs[\"sigma\"] * cvar[\"E\"] - coeffs[\"gamma\"] * cvar[\"I\"]\n", + "\n", + "def fR(cvar, coeffs):\n", + " return coeffs[\"gamma\"] * cvar[\"I\"]\n", + "\n", + "cfuncs={\"S\": fS, \"E\": fE, \"I\": fI, \"R\": fR}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The training and validation data is loaded from the provided `PyTorch` file and normalized:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total population: 3338330\n" + ] + } + ], + "source": [ + "SIR=torch.load(\"SIR_data_SD_county.pt\")\n", + "\n", + "npop=SIR[0,0].numpy().copy()\n", + "print(\"Total population: %d\" % npop)\n", + "SIR[:,:] = SIR[:,:] / torch.tensor(npop)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Preparing the training data. Here, just a short excerpt of the full time range is used.\n", + "Time-dependent coefficients are needed to fit longer time windows. This is not yet implemented in the ODE library, but in the specific solver for the SEIR model:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "nt=37\n" + ] + } + ], + "source": [ + "training_data=SIR.float()[350:380,:]\n", + "nt_train=training_data.shape[0]\n", + "test_data=SIR.float()[350:410,:]\n", + "nt=test_data.shape[0]\n", + "print(\"nt=%d\" % nt)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'S': array(0.94708073, dtype=float32), 'I': array(0.01257096, dtype=float32), 'R': array(0.04034832, dtype=float32), 'E': 0.025141911581158638}\n" + ] + } + ], + "source": [ + "#The values at the beginning of the observation are taken as initial values\n", + "inivals={}\n", + "for n,var in enumerate([\"S\",\"I\",\"R\"]):\n", + " inivals[var] = training_data[0,n].numpy()\n", + "\n", + "#The fraction of the initial exposed population is assumed twice the infected fraction\n", + "inivals[\"E\"] = inivals[\"I\"] * 2.\n", + "print(inivals)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "inicoeffs={\"beta\": 0.50, \"gamma\": 0.20, \"sigma\": 0.20}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The function is initialized using the initial values, initial coefficients and functions (right-hand-sides of ODEs) provided above. Also specified are the output variables given in the training data, in the order in which they are present:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "ename": "NameError", + "evalue": "name 'seir' is not defined", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmyopt\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptim\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSGD\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseir\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcoeffs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalues\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.005\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mseir\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mODESolver\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minivals\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcfuncs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minicoeffs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdt\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moutvar\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"S\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\"I\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\"R\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moptimizer\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmyopt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mNameError\u001b[0m: name 'seir' is not defined" + ] + } + ], + "source": [ + "myopt=torch.optim.SGD(seir.coeffs.values(), 0.005)\n", + "seir=ode.ODESolver(inivals, cfuncs, inicoeffs, dt=1, outvar=[\"S\",\"I\",\"R\"], optimizer=myopt)" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 1000/1000 [00:13<00:00, 72.02it/s, loss=4.66e-5]\n" + ] + } + ], + "source": [ + "seir.fit(training_data, num_epochs=1000)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The values of the optimized coefficients can be retrieved like this:" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'beta': 0.3824509084224701,\n", + " 'gamma': 0.3044351637363434,\n", + " 'sigma': 0.0070380899123847485}" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "seir.get_coeffs()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Training using custom optimizer:" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [], + "source": [ + "#myopt=torch.optim.SGD(seir.coeffs.values(), 0.005)\n", + "#seir.fit(training_data, optimizer=myopt)" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "y_predict=seir.predict(nt)" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9rklEQVR4nO3deZgU1dX48e+ZfQYYwAFEQQEjbiAMsigGAUURgSCIiuKrgIoxRiKJBs0bf6IS12BelyBGjUYMKEnEBAnGBUHEqGyyyaKILIPKKsswa/ec3x+3e6anp2fohu7paeZ8nqee7qpbVX1qljp97626JaqKMcaY+isp3gEYY4yJL0sExhhTz1kiMMaYes4SgTHG1HOWCIwxpp6zRGCMMfVczBKBiLwkIjtFZE015SIiT4vIRhFZJSLnxCoWY4wx1YtljeAvwIAayi8D2vumW4CpMYzFGGNMNWKWCFR1IbC3hlUuB6ap8ynQREROiFU8xhhjQkuJ42e3ArYFzOf5ln0XvKKI3IKrNdCgQYOuZ5xxRq0EaIwxx4ply5btVtXmocrimQgkxLKQ412o6vPA8wDdunXTpUuXxjIuY47K7NmzARgyZEicIzGmgohsqa4snokgDzgpYL418G2cYjEmaq644goAPB5PnCMxJjzxTASzgdtF5HXgXGC/qlZpFjIm0YwZMybeIRgTkZglAhF5DegLNBORPGAikAqgqs8Bc4GBwEagALD/HnNMeOGFF+IdgjERiVkiUNVrD1OuwM9j9fnGxMuMGTMAGDlyZJwjMSY8kmjPI7DOYlPXpaS471fWR2DqEhFZpqrdQpXFs4/AmGPSz372s3iHYExELBEYE2XPPPNMvEMwJiI26JwxUfbCCy9Yh7FJKNZHYEyUWR+BqYusjwCYOxeuuw48HsjOhpwcaNkSmjeHxo0rpuxs99qoETRo4KaGDSveN2gAmZmQZHUpU41f/epX8Q7BmIjUm0TQrBmceips2ADffuum1ashJcWd+Pfvd0kiXJmZLilkZbn3/lf/+8DkEZxIsrIqb5uZCRkZbvK/z8yE9HSQUANxmDrt8ccfj3cIxkSk3iSCHj1gyRL3/tAhWLUKPv8c9u2D//1fUHXrBLc6deoE994L+fkwY4bbNiUFkpPdSTo93dUeCgth7163vx073HxBgdvu0CEoKzuyuP1JwZ9g/JM/YaSnQ1paxZSeXpFUAhNLqCktDVJTK179k//4UlLclJ7uPjM19Wh+A/XHk08+CcD48ePjGocx4bI+ggAlJbBlC2zaBNu3w3ffuSakW2915X37wpo18MMPFSf2oUPhzTfd+xYtYNeuiv1lZMANN8Bzz0FREVxwgdsu8GTbrZvbb0EBvPhixbaqbt3jj3fNV/n58MUXrtZSWupi9XrdfkTc/vPzK8r85dGUnFy59hP8GphcghNMSkrNr+G+D9534Ks/KaanV/xc4sH6CExdZH0EYUpLg/bt3RTKggXutawMDhxwNYDAvoJnnnFJ4sCBihpBbq47IWVmwkknuWVFRW46eNA1WQ0Z4ta//XYoLnYn8+Jit8///V946CHYudMlhWCPPAL33OOS149+VLV88mQYNcrVfvr3r1p+550uEa1aBb/9rTuewGnUKDjtNJeE/vEPt01RkYvvhx+gY0f3c8vLg61bXQILnBo0cK/+bfzLa0NgohBxkz8B+2s8WVnuVdWV+5f7X9PT3fukJPcaqgYVXJs6//z3SU4u45FHKmpV1U3+WIKn4N+DfwqMJfA4Al+D91PdfFKSNT0axxLBEUhKgiZN3BRoxIiat/PXHELJzITduyvmVd03ev9JMyfH1VJKSyu+9Xs8FcnhhBNg/nxX5vW6V4/HNW01awZdu8JLL1XUKDwet96ll8JZZ8HZZ7sTtdfryvzlY8a4faxa5WIJLPN4YOJEV/7hhzBpklvuL/N63Wd26OCSyL33Vi4rLYV//xvatIHnn4ff/a7y9gBvveXif+EFt69gzz3nTsCvvQbvv1+5rLQUfv5zl7jnzXPJLFBysqvReTyuSXD79srlqanQvbuLZ8MG1+wHFSfPlBT3eyktdX1MFb+vvkDFF4e6zJ8g/VNw0klJqfiyE1zmn/zb+RNpuEmopoQXGEtwjIGxBr9Wt89Q+wtW03bBP5/gn9XhYgwVc6jPCN5f8Pt27dwXs6j/HVjTkKmrvN6Kf4TiYldr8ic6f7I46SS3zo4drlkuOFGdf77bfu1a1+znL/OftK+80n3WggXw1VeVE2FGRkWz4Ouvu30EJrKmTV0tCuCJJ1y5xwOrVn2B15vMxRefwaRJbtmvf+2SiT8Ber2u5ulPjrff7mpVgfF16eJqe2Vlrma2Z09FWVkZ9OwJ48e7+dGjXa0r0PnnwzXXuH2GupDp3HNdLbGoCH7/+6rlPXq46eBBeOWVquXnnOO+ROzbB3PmuGWBJ70OHaBVK1e+bFnFydf/2qaNu0ovP9/9boI1a+ZqZIWFrvYZLDvbnfiLi10/HFSucaamVv7yElzuj9PfDOsvr8vuvhseffTItq2pacgSgTFRVlt9BP7mLHAnwsCaltfrEtlxx7n1/EkoMNG0aOFOxh6PS4SB23q9cPrp7mReUAAzZ1ZOkh6P6/Pq2tUlqGefrZyEvV644gqXrLZsgYcfrpqkb78devVytc177qlc5vHAY4+58gULYNy4qrXVmTPhvPNcbfPmmyt/ttcLixe7+P70p4qEHmjDBvftevJkl6iDbd7sLjGfNMk1zwbbssU1fd57r6uZBtu61b3efbersfqbJpOSXBL76COXgO6/3x1jYE2pWTMXtyr84Q+ub/KUU+Dpp93v7EhYIjCmFj322GMA3H333XGOxIBLEEVFVRNJ8+au1rB3r6tRBiZKj8fVyNLSXP/b119XTWRDh7oT++LF7lL04ETrv2hszhxXIwrcPjW1Irm88AJ89lnl+Bo3rkgu993nyk8/3SWCI2WJwBhj6rmaEoHdH2tMlE2cOJGJEyfGOwxjwmY1AmOizO4jMHWR3UdgTC2aPHlyvEMwJiKWCIyJMhtawiQa6yMwJsomTJjAhAkT4h2GMWGzPgJjosz6CExdZH0ExtSiqVOnxjsEYyJiicCYKBs7dmy8QzAmItZHYEyUjRs3jnHjxsU7DGPCZn0ExkSZ9RGYusj6CIypRdOmTYt3CMZExBKBMVE2cuTIeIdgTESsj8CYKBs7dqx1GJuEYn0ExkSZ9RGYusj6CIypRbNmzYp3CMZExBKBMVE2ZMiQeIdgTESsj8CYKBs1ahSjRo2KdxjGhM36CIyJMusjMHWR9REYU4vefffdeIdgTEQsERgTZRdddFG8QzAmIjHtIxCRASKyQUQ2isg9Icobi8hbIrJSRL4QkTGxjMeY2jBixAhGjBgR7zCMCVvMagQikgxMAS4B8oAlIjJbVdcGrPZzYK2q/kREmgMbRGS6qpbEKi5jYu2NN96IdwjGRCSWTUM9gI2quglARF4HLgcCE4ECjUREgIbAXsB62ExCW7RoUbxDMCYisUwErYBtAfN5wLlB6/wRmA18CzQCRqhqWfCOROQW4BaAk08+OSbBGhMt5513XrxDMCYisewjkBDLgq9VvRRYAZwI5AJ/FJHsKhupPq+q3VS1W/PmzaMdpzFRNXToUIYOHRrvMIwJWyxrBHnASQHzrXHf/AONAR5VdzPDRhH5BjgDWBzDuIyJqTlz5sQ7BGMiEstEsARoLyLtgO3ANUDw+LxbgX7ARyJyPHA6sCmGMRkTcytWrIh3CMZEJGaJQFU9InI78A6QDLykql+IyK2+8ueAScBfRGQ1rinpblXdHauYjKkNHTt2jHcIxkQkpjeUqepcYG7QsucC3n8L9I9lDMbUtssuuwyAt99+O86RGBMeu7PYmCh777334h2CMRGxRGBMlG3aZN1cJrFYIjAmyuxeF5No7HkExkRZv3796NevX7zDMCZsViMwJso+/PDDeIdgTEQsERgTZXv27Il3CMZExBKBMVHWuHHjeIdgTESsj8CYKOvduze9e/eOdxjGhM1qBMZE2X//+994h2BMRCwRGBNl9tB6k2isacgYY+o5SwTGRFnPnj3p2bNnvMMwJmzWNGRMlC1dujTeIRgTkWMiEZSWlpKXl0dRUVG8Q6lzMjIyaN26NampqfEOpd4oLS2NdwjGROSwiUBE0oHhQNvA9VX1wdiFFZm8vDwaNWpE27ZtEQn1hMz6SVXZs2cPeXl5tGvXLt7hGGPqqHD6CP4FXA54gEMBU51RVFRETk6OJYEgIkJOTo7VlGpZ165d6dq1a7zDMCZs4TQNtVbVATGP5ChZEgjNfi61b/Xq1fEOwZiIhFMj+K+InB3zSIw5RpSUlFBSUhLvMIwJWziJoBewTEQ2iMgqEVktIqtiHViieeihh+jQoQOdOnUiNzeXzz77LOR648ePZ+HChTXuq23btuzeHd6jm1evXs3o0aMjDdcYY8qF0zR0WcyjSHCffPIJc+bMYfny5aSnp7N79+6Q3wj37t3Lp59+ypNPPhm1zz777LPJy8tj69at9kCUOqJTp04ArFpl35dMYjhsIlDVLSLSGbjAt+gjVV0Z27CO3PjxsGJFdPeZmws1nbu/++47mjVrRnp6OgDNmjULud4//vEPBgyo6G6ZN28ed911Fx6Ph+7duzN16tTyffz+979n/vz5AMyYMYNTTz2Vv//97zzwwAMkJyfTuHHj8prFT37yE15//XUmTJhw9AdrjtqGDRviHYIxETls05CI3AFMB1r4pr+KyLhYB5ZI+vfvz7Zt2zjttNO47bbbqn0wyccff1x+NUlRURGjR49m5syZrF69Go/Hw9SpU8vXzc7OZvHixdx+++2MHz8egAcffJB33nmHlStXMnv27PJ1u3XrxkcffRS7AzQRKS4upri4ON5hGBM+Va1xAlYBDQLmGwCrDrddrKauXbtqsLVr11ZZVts8Ho/Onz9f77vvPj3++OP15ZdfrrLOJZdcop988omqqq5YsUIvuOCC8rL3339fhw0bpqqqbdq00a+//lpVVUtKSvS4445TVdWf/vSnevHFF+vzzz+vu3fvLt/2yy+/1HPOOafa2OrCz8cYE1/AUq3mvBpOZ7EA3oB5r2+ZCZCcnEzfvn154IEH+OMf/8gbb7xRZZ3MzMzya/rd76V6gZd9+t8/99xz/O53v2Pbtm3k5uaWPwmrqKiIzMzMaB2KOUpnnnkmZ555ZrzDMCZs4SSCl4HPROR+Ebkf+BT4c0yjSjAbNmzgq6++Kp9fsWIFbdq0qbLemWeeycaNGwE444wz2Lx5c/n8q6++Sp8+fcrXnTlzZvmrfwCzr7/+mnPPPZcHH3yQZs2asW3bNgC+/PJLOnbsGJuDMxHbvHkzmzdvjncYxoQtnM7iP4jIAtxlpAKMUdXPYx1YIsnPz2fcuHHs27ePlJQUTj31VJ5//vkq6w0aNIg//elP3HzzzWRkZPDyyy9z1VVXlXcW33rrreXrFhcXc+6551JWVsZrr70GwK9//Wu++uorVJV+/frRuXNnAObPn8+gQYNq52DNYRUWFsY7BGMiItU1UYhItqoeEJHjQpWr6t6YRlaNbt26afDojuvWrUuYqnivXr2YM2cOTZo0icr+iouL6dOnD4sWLSIlJXReT6SfjzEmNkRkmap2C1VWU41gBjAYWAYEZgvxzZ8StQjrkSeeeIKtW7dGLRFs3bqVRx99tNokYGrfaaedBrgmO2MSQbVnD1Ud7Hu1YSuj6Nxzz43q/tq3b0/79u2juk9zdPLy8uIdgjERCec+gnnhLDPGOAUFBRQUFMQ7DGPCVm2NQEQygCygmYg0peKS0WzgxFqIzRhjTC2oqWH5p8B43El/GRWJ4AAwJbZhGZO4TjnFdZ9t2rQpzpEYE56a+gieAp4SkXGq+kwtxmRMQtu5c2e8QzAmIuHcUFYmIk38MyLSVERui11IiSlaw1A/99xzTJs2LezPLSkpoXfv3ng8nohjNrGRn59Pfn5+vMMwJmzhJIKxqrrPP6OqPwBjw9m5iAzwPcdgo4jcU806fUVkhYh8ISKhR2ur4wKHoV61ahXvv/8+J510UpX1/MNQ9+7du9p93Xrrrdxwww1hf3ZaWhr9+vUrvxPZGGMiFc7F50kiIr5BixCRZCDtcBv51psCXALkAUtEZLaqrg1YpwnwLDBAVbeKSIsjOIYq+vatuuzqq+G226CgAAYOrFo+erSbdu+GK6+sXLZgQc2fd6TDUN9zzz3Mnj2blJQU+vfvz+TJk7n//vtp2LAhd911F3379qVLly4sW7aMXbt2MW3aNB555BFWr17NiBEj+N3vfgfA0KFD+c1vfsN1111Xc6CmVvifC7F169Y4R2JMeMJJBO8AfxOR53A3kt0K/CeM7XoAG1V1E4CIvA5cDqwNWGckMEtVtwKoakI2rvbv358HH3yQ0047jYsvvpgRI0ZUGjfI7+OPP+ZKX5bZu3cvb775JuvXr0dE2LdvX8h9p6WlsXDhQp566ikuv/xyli1bxnHHHcePfvQjfvnLX5KTk0PHjh1ZsmRJLA/RRKC636UxdVU4ieBu3BVEP8NdOfQu8GIY27UCtgXM5wHBd1OdBqT6xjJqBDylqlUayEXkFuAWIKyncNX0DT4rq+byZs0OXwMI1rBhQ5YtW8ZHH33E/PnzGTFiBI8++miVR0h+9913NG/eHHDPG8jIyODmm29m0KBBDB48OOS+hwwZArgnkXXo0IETTjgBcFembNu2jZycHJKTk0lLS+PgwYM0atQosuBN1B04cCDeIRgTkcP2EahqmapOVdUrVXW4qv5JVb2H247QQ1UHD2yUAnQFBgGXAv9PRE4LEcPzqtpNVbv5T6R1TaTDUKekpLB48WKGDx/OP//5z0pNRoH8zU1JSUnl7/3zgR3ExcXFZGRkRPOQjDH1xGFrBCLSHngEOAsoP9Oo6uHGGsoDAntMWwPfhlhnt6oeAg6JyEKgM5BQg7Rs2LCBpKSk8qEeDjcMdd++fcnPz6egoICBAwdy3nnnceqppx7x5+/Zs4fmzZuTmpp6xPsw0dOqVSsAtm/fHudIjAlPuM8jmAp4gAuBacCrYWy3BGgvIu1EJA24BpgdtM6/gAtEJEVEsnBNR+vCDb6uyM/PZ9SoUZx11ll06tSJtWvXcv/991dZb9CgQSzwtTsdPHiQwYMH06lTJ/r06cP//d//HfHnz58/n4GhesBNXBw6dIhDhw7FOwxjwlfdo8v8E7DM97o6YNlHh9vOt95A3Lf7r4Hf+pbdCtwasM6vcR3Ia4Dxh9tnXX1UZbh+/OMf6w8//BDVfQ4bNkzXr19fbXki/XyMMbFBDY+qDKezuEhEkoCvROR2YDvuIfbhJJm5wNygZc8Fzf8e+H04+zsWRHsY6pKSEoYOHcrpp58elf0ZY+qfcBLBeNzgc78AJuGah0bFMKZjWrSHoU5LS4voBjQTey1btgTg+++/j3MkxoSnxkTguynsalX9NZAPjKmVqIxJYKWlpfEOwZiI1JgIVNUrIl0D7yw2xtRsz5498Q7BmIiE0zT0OfAvEfk7UH4phKrOillUxhhjak04ieA4YA9wUcAyBSwRGBOC/6bHXbt2xTkSY8JT7X0EIvIwgKqOAWao6piA6cZaizBBNGzYMOTywsJC+vTpg9db/c3YCxYsqHaIiVDuuusuPvjgg4hjNMaYUGq6oSxwzIPHYh3Iseqll17iiiuuIDk5OWr7HDduHI8++mjU9meia9euXVYbMAklnKahxDJ+PKxYEd195ubCk08e0abTp09nxowZgLt5b8KECbz99tuICPfeey8jRowA3EBlw4YNY8OGDfTu3Ztnn30WVeWmm25i6dKliAg33ngjv/zlL2nTpg179uzh+++/L79U0RhjjlRNiaCFiPwKN3ic/305Vf1DTCM7BpSUlLBp0ybatm0LwKxZs1ixYgUrV65k9+7ddO/evfwhNYsXL2bt2rW0adOGAQMGMGvWLNq1a8f27dtZs2YNUHl443POOYePP/6Y4cOH1/ZhmcPIyckB7OohkzhqSgQv4IaGDn5ftx3hN/dY2L17d6U7iBctWsS1115LcnIyxx9/PH369GHJkiVkZ2fTo0eP8oeeX3vttSxatIh+/fqxadMmxo0bx6BBg+jfv3/5vlq0aMG33waP4WfqAhv8zySamh5e/0BtBnIsChx2GqCmWzFEpMp806ZNWblyJe+88w5Tpkzhb3/7Gy+99BIARUVFZGZmxiZwc1TsjmKTaMIZfdQcoaZNm+L1esuTQe/evZk5cyZer5ddu3axcOFCevToAbimoW+++YaysjJmzpxJr1692L17N2VlZQwfPpxJkyaxfPny8n1/+eWXdOzYMS7HZYw5thx7ncV1TP/+/Vm0aBEXX3wxw4YN45NPPqFz586ICI8//jgtW7Zk/fr19OzZk3vuuYfVq1fTu3dvhg0bxurVqxkzZgxlZWUAPPLII4AbwmDjxo1069YtnodmquFvDrRHVppEIYk2ckS3bt106dKllZatW7eOM888M04R1ezzzz/nD3/4A6++Gs4jHMLz5ptvsnz5ciZNmhTW+nX553MssgfTmLpIRJapashvj+E8oSwdGA60DVxfVR+MVoDHsi5dunDhhRfi9Xqjdi+Bx+PhzjvvjMq+TPRZAjCJJpymoX8B+4FlQHFswzk23XhjdG/Evuqqq6K6P2NM/RZOImitqqGfrG6MqSI7OxtwNwkakwjCuWrovyJydswjMeYY0aRJk6g9gc6Y2hBOjaAXMFpEvsE1DQmgqtopppEZk6C2bt0a7xCMiUg4ieCymEdhjDEmbg7bNKSqW4AmwE98UxPfMhPgaIahHjhwYETXnM+ZM4eJEydGGqKpJQ0bNqz278GYuuiwiUBE7gCmAy18019FZFysAztWhDMM9dy5cyNqUx40aBCzZ8+moKAgChGaaGvRogUtWrSIdxjGhC2cpqGbgHNV9RCAiDwGfAI8E8vAjkrfvlWXXX013HYbFBTAwIFVy0ePdtPu3XDllZXLFiw44lACh6H+7rvvGDFiBAcOHMDj8TB16lQuuOAC2rZty9KlS8nPz2fAgAH06tWLTz/9lM6dOzNmzBgmTpzIzp07mT59Oj169EBE6Nu3L3PmzOHqq68+4thMbGzatCneIRgTkXCuGhIgsF3D61tmDiN4GOoZM2Zw6aWXlg9FnZubW2WbjRs3cscdd7Bq1SrWr1/PjBkzWLRoEZMnT+bhhx8uX69bt2589NFHtXQkxphjWTg1gpeBz0TkTd/8UODPMYsoGmr6Bp+VVXN5s2ZHVQMIFDwMdffu3bnxxhspLS1l6NChIRNBu3btOPtsd7Vuhw4d6NevHyLC2WefzebNm8vXs2Go666srCwAa7ozCSOczuI/AGOAvcAPwBhVfTLGcR0Tgoeh7t27NwsXLqRVq1Zcf/31TJs2rco26enp5e+TkpLK55OSkvB4POVlNgx13dW6dWtat24d7zCMCVu1NQIRyVbVAyJyHLDZN/nLjlPVvbEPL7EFDkOdkZHBli1baNWqFWPHjuXQoUMsX76cG2644Yj2bcNQ111ffvllvEMwJiI11Qhm+F6XAUsDJv+8CYN/GGqABQsWkJubS5cuXXjjjTe44447jni/8+fPZ9CgQdEK0xhTj9kw1DEWi2God+zYwciRI5k3b15Y69fln8+xyN9kV1hYGOdIjKlQ0zDU4dxHUOVsE2qZCS1wGOpo2bp1K0888UTU9meiq23btuVXihmTCGrqI8gAsoBmItKUiktGs4ETayG2Y0a0h6Hu3r17VPdnomvdunXxDsGYiNR0+ehPgfG4k/4yKhLBAWBKbMMyxhhTW6pNBKr6FPCUiIxT1bp7F7ExdYz/kt/iYnuOk0kMh72hTFWfEZGOwFlARsDyqhfBG2M4/fTT4x2CMREJ55nFE4G+uEQwFzcs9SLAEoExIaxatSreIRgTkXDGGroS6Ad8r6pjgM5Aes2bOCIyQEQ2iMhGEbmnhvW6i4hXRK6sbp26Ljk5mdzcXDp27MhPfvKTaoeVDmdY6gULFjB48OCwP/uuu+7igw8+iDRkY4wBwksEhapaBnhEJBvYCZxyuI1EJBnXqXwZrjZxrYicVc16jwHvRBJ4XZOZmcmKFStYs2YNxx13HFOmhO5PD2dY6kiNGzeORx99NGr7M0cnLS2NtLS0eIdhTNjCGXRuqYg0AV7AXT2UDywOY7sewEZV3QQgIq8DlwNrg9YbB7wBROWayPH/Gc+K71dEY1flclvm8uSAJ8Nev2fPntU2DwQOS62qTJgwgbfffhsR4d5772XEiBGAe/D5sGHD2LBhA7179+bZZ59FVbnppptYunQpIsKNN97IL3/5S9q0acOePXv4/vvvadmy5VEfrzk6/kEDjUkU4XQW3+Z7+5yI/AfIVtVwGkFbAdsC5vOAcwNXEJFWwDDgImpIBCJyC3ALwMknnxzGR8eP1+tl3rx53HTTTVXKgoelnjVrVvmQ1Lt376Z79+707t0bgMWLF7N27VratGnDgAEDmDVrFu3atWP79u2sWbMGoFLz0znnnMPHH3/M8OHDY36MpmbLli2LdwjGRKSmG8rOqalMVZcfZt+hnlkQPJ7Fk8DdquoVqf4RB6r6PPA8uCEmavrQSL65R1NhYSG5ubls3ryZrl27cskll1RZJ3hY6kWLFnHttdeSnJzM8ccfT58+fViyZAnZ2dn06NGDU05xLXDXXnstixYtol+/fmzatIlx48YxaNAg+vfvX74vG5baGHOkaqoR1DSGgeK+xdckDzgpYL41EHym6ga87ksCzYCBIuJR1X8eZt91jr+PYP/+/QwePJgpU6bwi1/8oso6gcNS1zTOU3BiFBGaNm3KypUreeedd5gyZQp/+9vfeOmllwAblrouSU1NBaC0tDTOkRgTnmo7i1X1whqmwyUBgCVAexFpJyJpwDXA7KDPaKeqbVW1LfAP4LZETAKBGjduzNNPP83kyZOrnAgCh6UG93yCmTNn4vV62bVrFwsXLqRHjx6Aaxr65ptvKCsrY+bMmfTq1Yvdu3dTVlbG8OHDmTRpEsuXV1TKbFjquqNbt2506xZybC9j6qRw7iMIOWD+4W4oU1WPiNyOuxooGXhJVb8QkVt95c8dQbwJoUuXLnTu3JnXX3+d66+/vlKZf1jqiy++mGHDhvHJJ5/QuXNnRITHH3+cli1bsn79enr27Mk999zD6tWr6d27N8OGDWP16tWMGTOGsrIyAB555BHAffPcuHGjnXzqiE8++STeIRgTkcMOQy0igcNLZODuKViuqnG55j/RhqEOFothqd98802WL1/OpEmTQpYn0s/HGBMbNQ1DHc5VQ+OCdtYYiN5ZrJ4JHJY6WvcSeDwe7rzzzqjsyxy9lBT3bxX4aFFj6rJw7iMIVgC0j3Yg9Um0h6W+6qqroro/c3TOP//8eIdgTETC6SN4i4rLPpOBM4G/xTIoYxLZwoUL4x2CMREJp0YwOeC9B9iiqnkxiseYhLd//37AXUFmTCIIp4/gQxFpiRsyQoGvYx6VMQksJycHsD4CkzjCaRq6GbgP+AB3t/AzIvKgqr4U6+CMSUR9+vSJdwjGRCSc0Ud/DXRR1dGqOgroCtwd27ASTzSHoR44cGC124cyZ84cJk6cGGHEJlbmzZvHvHnz4h2GMWELJxHkAQcD5g9SeTA5Q3SHoZ47d26lMYkOZ9CgQcyePZuCgoJIwzYxsHXrVrZu3RrvMIwJWzidxduBz0TkX7g+gsuBxSLyKwBV/UMM4zsiff/St8qyqztczW3db6OgtICB0wdWKR+dO5rRuaPZXbCbK/9W+V65BaMXRPT54Q5D/d133zFixAgOHDiAx+Nh6tSpXHDBBbRt25alS5eSn5/PgAED6NWrF59++imdO3dmzJgxTJw4kZ07dzJ9+nR69OiBiNC3b1/mzJnD1VdfHVGsJvr8gwVaH4FJFOHUCL4G/knFJaT/Ar4DGvkmE8A/DPWQIUOqlAUPQz1jxgwuvfTS8qGoc3Nzq2yzceNG7rjjDlatWsX69euZMWMGixYtYvLkyTz88MPl63Xr1o2PPvooVodlInDJJZeEHH3WmLoqnKuGHgAQkUZuVvNjHtVRqukbfFZqVo3lzbKaRVwDgCMbhrp79+7ceOONlJaWMnTo0JCJoF27duUPOunQoQP9+vVDRDj77LPZvHlz+Xo2DHXd8fbbb8c7BGMictgagYh0FJHPgTXAFyKyTEQ6xD60xOLvI9iyZQslJSUh+wiCh6Hu3bs3CxcupFWrVlx//fVMm1Z1HL/09IrHQyclJZXPJyUlVWp6sGGo6441a9aUPzzImEQQTtPQ88CvVLWNqrYB7sQ9ttKEEMkw1Fu2bKFFixaMHTuWm266qdKw0pGyYajrjtzc3JC1O2PqqnASQQNVne+fUdUFQIOYRXQMCByGOph/GGqABQsWkJubS5cuXXjjjTe44447jvgz58+fz6BBg454exM9gwcPZvDgwfEOw5iwhTMM9ZvAcipGHP0foJuqDo1taKHZMNRV7dixg5EjR1Z77Xoi/XyMMbFR0zDU4dQIbgSaA7N8UzNgTPTCq18Ch6GOlq1bt/LEEzU9WdTUpk8//ZRPP/003mEYE7aaHl6fAdwKnAqsBu5UVXsIaxREexjq7t27R3V/5uj06tULsPsITOKo6fLRV4BS4CPgMtzw0+NrIaYjoqpVHvhu3M/F1K7hw4fHOwRjIlJTIjhLVc8GEJE/A4trJ6TIZWRksGfPHnJyciwZBFBV9uzZQ0ZGRrxDqVdmzpwZ7xCMiUhNiaC8Gcj3IPpaCOfItG7dmry8PHbt2hXvUOqcjIwMWrduHe8w6pUPPvgAgIsuuijOkRgTnpoSQWcROeB7L0Cmb15wdxhnxzy6MKWmptKuXbt4h2EM4C4RBusjMImj2kSgqtF5srox9cx1110X7xCMiciRPLzeGFODV155Jd4hGBORcO4jMMZEYPbs2cyePTveYRgTNqsRGBNlV1xxBWB9BCZxWCIwJsrGjLEb701isURgTJS98IINzmsSi/URGBNlM2bMKH8cqTGJ4LCjj9Y1oUYfNaYuSUlxFW3rIzB1SU2jj1rTkDFR9rOf/SzeIRgTEUsExkTZM888E+8QjImI9REYE2UvvPCCdRibhGJ9BMZEmfURmLrI+giMqUW/+tWv4h2CMRGJaSIQkQHAU0Ay8KKqPhpUfh1wt282H/iZqq6MZUzGxNrjjz8e7xCMiUjM+ghEJBmYgnu62VnAtSJyVtBq3wB9VLUTMAl4PlbxGFNbnnzySZ588sl4h2FM2GLWRyAiPYH7VfVS3/xvAFT1kWrWbwqsUdVWNe3X+ghMXWd9BKYuilcfQStgW8B8HnBuDevfBLwdqkBEbgFuATj55JOjFZ8xMfHb3/423iEYE5FYJoJQz7YMWf0QkQtxiaBXqHJVfR5fs1G3bt0S6zInU+888MAD8Q7BmIjEMhHkAScFzLcGvg1eSUQ6AS8Cl6nqnhjGY0yteOyxxwC4++67D7OmMXVDLPsIUoAvgX7AdmAJMFJVvwhY52TgA+AGVf1vOPu1PgJT11kfgamL4tJHoKoeEbkdeAd3+ehLqvqFiNzqK38OuA/IAZ4VEQBPdYEakygeeuiheIdgTETszmJjjKkHaqoR2FhDxkTZxIkTmThxYrzDMCZsViMwJsqsj8DURTbWkDG1aPLkyfEOwZiIWCIwJsrGjx8f7xCMiYj1ERgTZRMmTGDChAnxDsOYsFkfgTFRZn0Epi6yPgJjatHUqVPjHYIxEbFEYEyUjR07Nt4hGBMR6yMwJsrGjRvHuHHj4h2GMWGzPgJjosz6CExdZH0ExtSiadOmxTsEYyJiicCYKBs5cmS8QzAmItZHYEyUjR071jqMTUKxPgJjosz6CExdZH0ExtSiWbNmxTsEYyJiicCYKBsyZEi8QzAmItZHYEyUjRo1ilGjRsU7DGPCZn0ExkSZ9RGYusj6CIypRe+++268QzAmIpYIjImyiy66KN4hGBMRSwTGRNmIESMAmDlzZpwjSXxlWoa3zItXvZVey7Qs5BRYHriNp8yDt8z3GrC8uv37ywInVa30OeXvA+KpaTtFq3xedTFXt68hpw9h5NnRv2HREoExUfbGG2/E/DNUldKyUkq9pXjKPOXvA19LvCWUen2vZaUUe4op9haXv5Z4S/CUeaqd/CfO4JNnldegE27gZwbGFBhnTZ8bOCmJ1Yd5OIKQnJRMsiSTJEkkJ/leJbnSe3+Z/71/6n5i99jEZZ3FxlQo0zKKPEUUeYoo8ZaUn0z9J9ZQJ6sSbwkFpQUUlhZS6Clk3VfrKC4rptkJzSgsLaSgtICC0gKKvNXvM/BEGTwf6r1XvbXy8wg8SYXzmpKUQmpyKmnJaaQmpZKanFrp1V+ekpRSPu+fkiW5vCxwSiIJEUEQgPL3mamZpCWnUVpWSkFpAYoiKiDuhNuiQQsapDWgsLSQnYd2uu0FUEDgrGZnkZ2Rzc5DO9m4ZyOIS7AAinJR24toktmEL/d8yfLvllf5Zj86dzTZ6dl8su0TFm1dVL7cP026cBIN0hrw5ro3eXfTu5SVleFRT3nSnHHFDESEqUumMnfj3PLk61UvaclpvH3d2wDcN/8+/v3Vv+nYoiOvDH3liH+X1llsEoqqUuwtrnRyLfIUHfa9f/3AE2+xp7j85FviLaHYW1x+og+cCkvdvkrLSqN3IF9CRkoGWalZpCenk56cXnFCTE4hLTmNBqkNyErNoiSphCQqf0NMT0mneVZzUpJS2FOwhzLKKn07zE7L5pTjTiElKYX1u9fjKXNXKSWJO3Ge2PBEup7YlbTkNOZ+NRdvmRcRceUIHVt05LL2l5GSlMLDHz1c0ZSBoqpc2O5CRnQYQbGnmJvfurlKTeCqs65i5Nkj+aHwB67+x9VVym/rdhvXd76erfu3MnD6wCq1iUkXTuL6ztezesdqLnzlwko1jzIt489D/sx1na5j4ZaF9PlLnyo/3jdHvMnQM4by7y//zeDXBlcpf//69+l3Sj9eX/M6E96v+ujQz27+jB6tevDi8he5e97dVcq/uO0LTss5jX9/+W8mLphYpfynXX9KiwYtWLljJZM/mVwlMT544YOkJKWwbvc65nw5p3y5P+mVaRnJkszewr18e/DbyttLcvnnNEprxAkNT6BFVouI/vwiYTUCA7iTb6hmgKzULDJSMij2FLPj0A68ZV5KvaUUego5WHKQxumNSZIkdhzawfpd6zlUeohDpYcoLC3kUOkhcjJzEBF2HtrJlv1bKPZUnIgLPYVkpmRS7C3mh8If2F+8v/zb8pE2Cfir3ilJKSRJEic2OpGMlAzyi/PZX7zfnQj93zBFGHDqADJTMlm7ay0b926stK+UpBQm9plIWnIaM7+YyWfbP6vU7tswrSHTr5hOWnIa9y+4n/mb55e36wK0a9KOTXdsAuDiaRcz75t5lfbf6fhOrLx1JQDnvngui7cvrlR+/knn8/GNHwNw1pSzWLd7XaXyS390Kf/5n/8AcPL/ncy2A9sqlV951pX8/aq/A9D0sabsK9pXqXx07mhevvxlABo90ghVrXQiu+WcW3io30MUe4o5c8qZVb7939L1Fm7vcTv7ivYxcPrAKuU3dbmJER1HsCN/B7fNva3KifKGTjfQ75R+5B3I49FFj1Ypv7rD1Zxzwjls27+NV1a+Un4C9a8z+LTB/Oi4H7Ft/zbe/frdSifZ5KRk+rTpw/ENj+fbg9+y4vsVVfaf2zKX7PRsdhfsZtv+bVXib9OkDWnJaeSX5JNfkl9l+6zULJIkCVVFRML/I42TmmoElgiiqKC0gPySfA6VHKLQU0hhaSGN0htxWs5pAMzeMJv8kvzyk2GJt4TTm53OgFMHAHDvB/e6b6UB7ap92/bluk7XUeot5Zo3rqnShDCy40jGdh3LvqJ9XPTKReUncv9J/bbutzGy40i++eEbBr82uMq3sjG5Y7igzQWs372ehz56qMoxdT6+M40zGvPtgW/Z+MPGKuXhSklKKf/GGqj7id1pld2K7w9+z6fbP61Sfl/v+2if054PvvmAl1e8XP5t2N92+u7/vMupOafy8ucv89RnT1X5Z14ydglNM5vy1KdP8fKKl6uUzx81n9TkVKYsnsJbX75VqTw9OZ3Xr3wdgBeWvcB/8/5b/m0tOSmZRmmNeOySxwB4deWrrN21luSkZB5+6GEog6cfe5rbe9wOwKx1s9i2f5s7Ufk+o1lWM4adOQyA975+jx+KfqgUW05WDuefdD4Ai7cvpshTVKm8SUYT2ue0B2Dj3o1VTuRZqVkcl3kcAPuL9ldqd/a/JsIJzESHJYIIlHpL2X5wOzvyd7CrYBcpSSnlJ+r75t/H2l1r2V+8n/1F+9lXtI+uJ3blteGvAdD2ybZs2b+l0v6GnjGUN0e8CUDz3zdnd8HuSuXXd7qeacOmUaZlNHm0Cd4yr/sn9f2jXtTuIoacPoS9BXt5eNHD5dV2f3tki6wWNM5ozIHiA2zet7nSif5IpSalkp6STkZKBs2ymtGiQQsykjPYX7zfNWWkZdEgpQGZqZl0aNGBk7JPwlvmZeuBrWSlZpGVmkXDtIY0SG1AbstcTmx0Ynl8wSei1tmtyUzNpKC0gAPFB6q0STdIa1D+rQtIiBPXmjVrAOjYsWOcIzGmgiWCIHsK9rBu9zrW7VrH/uL93HX+XQAMfX0oszfMrtQscVrOacy4YgaFnkJ+/e6v2X5wOxkpGaSnpJOWlEbLhi05/6TzKfGWsOTbJRR7ixGkvGMpOSmZzNRMijxF7C3cW6m9usRbQlGpayIp9BSGFXtqUiqNMxqTnZ5Nw7SG5SfchmkNaZDWoLzNOSs1i8yUzPKmnfSUdNKS0ypN/u0bpjWkUVqj8n0kiY08YsyxxhIB8OHmDxnzrzHkHcir1CHor4IXlBZwqOTQUV2u5r86Ii05jcyUTDJSMshMda/pyemV3mekZJR3JAZP2enZZKdn0zi9ccV738k/PTk9Ib4V12eXXXYZAG+//XacIzGmgl01hLssUETIycqhcXpjmmY0JScrh6YZTWmQ1qDSN+jM1Ipv0pkpmeUn8MCTtf+bt/+bdmpSqp2gDQDvvfdevEMwJiL1JhFc2O5Cvv7F1/EOw9QDmzZtincIxkSk3iQCY2rLySefHO8QjImI9QoaE2X9+vWjX79+8Q7DmLBZjSCaysrA6wWPx70Gv/d6K9YpK3OTiJuSkipeU1MhLa1iSk11y01C+PDDD+MdgjERqT+J4OOPYfx42LsXMjMrTrCqUFpadQo+cXu9bl3/VVb+17Iyd7L3eCqWxYI/QaSkuCk1teqUlubKghOLSNUEpOq2SU9326Wnuykjw/18/FNGhlvu33+oz/FPycmV1ws1paRUJLjAz/YvOwY63Pfs2RPvEIyJSEwTgYgMAJ4CkoEXVfXRoHLxlQ8ECoDRqro8JsEcOgSbN8OePZVP2JmZcPHF7iSVlwfFxdCoEWRlualBA2ja1J3k/N/K/Scr/7f51FRX7j9JJydXzPvfB07+k6Z/P/4Ts/8k7fFASYmbSksrXktLXVmoxBU4+RNW4H4DP9d/HKWl7nhLSiA/H4qKKqbCwoqprCwmv5KQgmtDgT/TcH++wT/r4CkwkQYn1cAkFvwz8//cgn++/r8B374a+/cZmDj9x+JPzMHJ2j8Fzof6mwk+Fv++Aqfgbf3bHQNJ1sRGzBKBiCQDU4BLgDxgiYjMVtW1AatdBrT3TecCU32v0de/P+za5U6k33wD69fDhg2wfz9MmuTWGTIE3nqr8nYdOoDvTlEuvdS9b9SoYuraFR5+2JU/+KD7jIyMin/+00+Hq65y5dOnuxOv/x81JQXatIFzfYf87ruVT9pJSXDiiW4fAP/9b9V/+pYt4eST3bf9VavceoGJ6vjj3Tqlpe54A/nLmzVzcX39ddWTRcuWkJ0NBw/Cxo0VScqftFq0cMny0CHYtq2iduRPWk2auOM8cAB27KicyEpKoGFDdzLdvx9++KFyQispcYkaoKDAJavgmlpamtu+qMitE1jzCSz3J9PAGpx/Pf8x+T/3WBX4d+d/9Tc7lpVVTS7+Gl5SUsXPJfDvw7+9iPv5BSYhkYpaZlKS+/sKVZ6R4fZVWFg1KWZlud+/qvv7CiwDV56e7n5/+fkV8fljbNiworygoGqibdLEfX5pqdt/4LYi7gtgWpqL7eDBymUAOTnuZ1BY6Lb38++neXNXfuhQ6PLjj3c/i4MH3T4Cy0TghBNcnPv3u/K0NLjwQujSJSp/DpWoakwmoCfwTsD8b4DfBK3zJ+DagPkNwAk17bdr164aUwUFql9/rfrZZ6pz56q+805F2WOPqd54o+pVV6ledplqr16qY8dWlJ93nmp2tmpamv87o+rgwRXlJ5xQsdw/XXNNRXnDhlXLA/cfXAaqv/qVKzt4MHT5ffe58m+/DV3++9+78g0bQpc/95wrX7YsdPn06a58wYLQ5f/6lyt/663Q5R984MpnzAhdvmSJK//Tn0KXr1/vyp94InT59u2u/P77Q5fv3+/K77wzdPm+fao//KB6ww1Vyxo0cGV796pecUXV8mbNVJcuVf3kE9Xzz69a3qqV6r//rTpnjurZZ1ctb9tW9a9/VZ02zb0PLj/1VNVnnlF9+mnVFi2qlp9xhupDD6k++GDov6327VV/9jPVn/5UNTU19P5HjnR/o6F+Nu3bqw4frjpkSOjyU05R7ddPtWfP0OWtW6t27ap65pmhy5s3V23XLvT/Dbj/tRYtVJs0CV3eqJFq06bu9xSqPCkp9PK6PN19dyRns0qApaqhz6sxu7NYRK4EBqjqzb7564FzVfX2gHXmAI+q6iLf/DzgblVdGrSvW4BbfLOn+xLGkWgG7D7sWonDjqfuOpaOBY6t4zmWjgXCP542qto8VEEs+whCNUgGZ51w1kFVnweeP+qARJZqNbdYJyI7nrrrWDoWOLaO51g6FojO8cTymsQ84KSA+dbAt0ewjjHGmBiKZSJYArQXkXYikgZcA8wOWmc2cIM45wH7VfW7GMZkjDEmSMyahlTVIyK3A+/gLh99SVW/EJFbfeXPAXNxl45uxF0+OiZW8fgcdfNSHWPHU3cdS8cCx9bxHEvHAtFoNo9VZ7ExxpjEYOMWGGNMPWeJwBhj6rl6kwhEZICIbBCRjSJyT7zjiZSIvCQiO0VkTcCy40TkPRH5yvfaNJ4xhktEThKR+SKyTkS+EJE7fMsT9XgyRGSxiKz0Hc8DvuUJeTzgRgYQkc999/ok+rFsFpHVIrJCRJb6liXk8YhIExH5h4is9/3/9IzGsdSLRBAw3MVlwFnAtSJyVnyjithfgAFBy+4B5qlqe2Cebz4ReIA7VfVM4Dzg577fR6IeTzFwkap2BnKBAb6r4BL1eADuANYFzCfysQBcqKq5AdfbJ+rxPAX8R1XPADrjfkdHfyzV3XJ8LE2EMdxFIkxAW2BNwHz5kBzACcCGeMd4hMf1L9yYVAl/PEAWsBw3ZlZCHg/ufp55wEXAHN+yhDwWX7ybgWZByxLueIBs4Bt8F/lE81jqRY0AaAVsC5jP8y1LdMer774L32uLOMcTMRFpC3QBPiOBj8fXlLIC2Am8p6qJfDxPAhOAwGFnE/VYwI1W8K6ILPMNVwOJeTynALuAl33Ndi+KSAOicCz1JRGENZSFqV0i0hB4AxivqgfiHc/RUFWvqubivk33EJGOcQ7piIjIYGCnqi6LdyxR9GNVPQfXNPxzEekd74COUApwDjBVVbsAh4hSk1Z9SQTH6lAWO0TkBADf6844xxM2EUnFJYHpqjrLtzhhj8dPVfcBC3D9OYl4PD8GhojIZuB14CIR+SuJeSwAqOq3vtedwJtADxLzePKAPF9tE+AfuMRw1MdSXxJBOMNdJKLZwCjf+1G4tvY6z/dAoj8D61T1DwFFiXo8zUWkie99JnAxsJ4EPB5V/Y2qtlbVtrj/kw9U9X9IwGMBEJEGItLI/x7oD6whAY9HVb8HtomI7wEl9APWEo1jiXcHSC12tAwEvgS+Bn4b73iOIP7XgO+AUtw3g5uAHFyn3le+1+PiHWeYx9IL1zS3CljhmwYm8PF0Aj73Hc8a4D7f8oQ8noDj6ktFZ3FCHguuXX2lb/rC/7+fwMeTCyz1/a39E2gajWOxISaMMaaeqy9NQ8YYY6phicAYY+o5SwTGGFPPWSIwxph6zhKBMcbUc5YIzDFJRHJ8o02uEJHvRWS7732+iDwbo88cLyI3hFjeNnDU2Ch8TpqILBSRmD1h0NQv9odkjkmqugd3zTUicj+Qr6qTY/V5vpPyjbg7PWNKVUtEZB4wApge688zxz6rEZh6RUT6Boyxf7+IvCIi7/rGrL9CRB73jV3/H98wGIhIVxH50Ddo2Tv+2/mDXAQsV1VPwDYrReQT4OcBn99WRD4SkeW+6Xzf8ldF5PKA9aaLyBAR6SDuWQcrRGSViLT3rfJP4LpY/IxM/WOJwNR3PwIGAZcDfwXmq+rZQCEwyJcMngGuVNWuwEvAQyH282MgcKC2l4FfqGrPoPV2ApeoGwRtBPC0b/mLwBgAEWkMnA/MBW4FnlI3oF033F3l4O5g7n6Ex2xMJdY0ZOq7t1W1VERWA8nAf3zLV+Oe/3A60BF4zw2RRDJuqI9gJ+B7kIvvRN5EVT/0lb2KG/kSIBX4o4jkAl7gNABV/VBEpohIC+AK4A1V9fhqFL8VkdbALFX9yre+V0RKRKSRqh6M0s/C1FOWCEx9VwygqmUiUqoVY66U4f4/BPgixDf7YIVAhu+9UP0w578EduCeLpUEFAWUvYpr7rkG19+Aqs4Qkc9wtZZ3RORmVf3At3560PbGHBFrGjKmZhuA5iLSE9zw2SLSIcR664BToXwo6v0i0stXFtiW3xj4TlXLgOtxNQy/vwDjffv4wvd5pwCbVPVp3CiTnXzLc4Bdqlp69Ido6jtLBMbUQFVLgCuBx0RkJW6k1PNDrPo2EPjAkzHAFF/TTmHA8meBUSLyKa5Z6FDAZ+3AJZSXA9YfAazxPf3sDGCab/mFuD4EY46ajT5qTJSIyJvABH87/hFsn4XrmzhHVfcfZt1ZuOdubziSzzImkNUIjImee3CdxhETEf/DbJ4JIwmkAf+0JGCixWoExhhTz1mNwBhj6jlLBMYYU89ZIjDGmHrOEoExxtRzlgiMMaae+/+Qz/pHl4AdigAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig,ax=plt.subplots()\n", + "clr=[\"b\",\"r\",\"g\"]\n", + "\n", + "for n,var in enumerate(seir.outvar):\n", + " ax.plot(test_data[:nt,n], label=\"%s (obs)\" % var, color=clr[n])\n", + " ax.plot(y_predict.detach()[:,n], label=\"%s (sim)\" % var, \n", + " linestyle=\"dashed\", color=clr[n])\n", + " \n", + " ax.set_ylim(0, 1.00)\n", + " ax.vlines(nt_train, 0, 1.00, color='k', linestyle='dotted')\n", + " ax.legend(loc='center left')\n", + " ax.set_xlabel(\"Time (days)\")\n", + " ax.set_ylabel(\"Population Fraction\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Figure above shows the observed (solid) against the simulated (dashed) compartment fractions. The vertical, dotted line marks the end of the training window; values beyond 30 days are predicted from the minimum misfit model." + ] + }, + { + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/ode/SIR_data_SD_county.pt b/examples/ode/SIR_data_SD_county.pt new file mode 100644 index 00000000..634eb8c6 Binary files /dev/null and b/examples/ode/SIR_data_SD_county.pt differ diff --git a/torchts/nn/models/ode.py b/torchts/nn/models/ode.py new file mode 100644 index 00000000..76933707 --- /dev/null +++ b/torchts/nn/models/ode.py @@ -0,0 +1,58 @@ +import torch +from torch import nn + +from torchts.nn.model import TimeSeriesModel + + +class ODESolver(TimeSeriesModel): + def __init__( + self, ode, init_vars, init_coeffs, dt, solver="euler", outvar=None, **kwargs + ): + super().__init__(**kwargs) + + if ode.keys() != init_vars.keys(): + raise ValueError("Inconsistent keys in ode and init_vars") + + if solver == "euler": + self.solver = self.euler + else: + raise ValueError(f"Unrecognized solver {solver}") + + for name, value in init_coeffs.items(): + self.register_parameter(name, nn.Parameter(torch.tensor(value))) + + self.ode = ode + self.var_names = ode.keys() + self.init_vars = { + name: torch.tensor(value, device=self.device) + for name, value in init_vars.items() + } + self.coeffs = {name: param for name, param in self.named_parameters()} + self.outvar = self.var_names if outvar is None else outvar + self.dt = dt + + def euler(self, nt): + pred = {name: value.unsqueeze(0) for name, value in self.init_vars.items()} + + for n in range(nt - 1): + # create dictionary containing values from previous time step + prev_val = {var: pred[var][[n]] for var in self.var_names} + + for var in self.var_names: + new_val = prev_val[var] + self.ode[var](prev_val, self.coeffs) * self.dt + pred[var] = torch.cat([pred[var], new_val]) + + # reformat output to contain desired (observed) variables + return torch.stack([pred[var] for var in self.outvar], dim=1) + + def forward(self, nt): + return self.solver(nt) + + def get_coeffs(self): + return {name: param.item() for name, param in self.named_parameters()} + + def _step(self, batch, batch_idx, num_batches): + (x,) = batch + nt = x.shape[0] + pred = self(nt) + return self.criterion(pred, x)