{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IoYdLutn7dUa"
   },
   "source": [
    "# Latent Dirichlet Allocation\n",
    "\n",
    "Recall the following generatize model for LDA. Suppose we have $K$ topics and $N$ documents.\n",
    "\n",
    "For each topic $k =1,\\ldots,K$, draw a topic \n",
    "\n",
    "$$\\theta_k \\sim \\text{Dir}(\\phi)$$\n",
    "\n",
    "Then, for each document $n = 1,\\ldots, N$, draw topic proportions \n",
    "\n",
    "$$\\pi_n \\sim \\text{Dir}(\\alpha)$$\n",
    "\n",
    "Finally, for each word $l$ in document $n$, first draw a topic assignment \n",
    "\n",
    "$$\n",
    "z_{n,d} \\mid \\pi_n \\sim \\text{Cat}(\\pi_n)\n",
    "$$\n",
    "\n",
    "and draw a word\n",
    "\n",
    "$$\n",
    "x_{n,d} \\mid z_{n,d} \\sim \\text{Cat}(\\theta_{z_{n,d}})\n",
    "$$\n",
    "\n",
    "As mentioned in the lecture notes, while this formulation is easier to present, it's more efficient to represent the documents as sparse vectors of _word counts_, $\\mathbf{y}_n \\in \\mathbb{N}^V$ where $y_{n,v} = \\sum_{d=1}^D \\mathbb{I}[x_{n,d} = v]$. \n",
    "\n",
    "This notebook studies Federalist papers in their entirety. We've provided a $N \\times V$ dataframe of the essays represented as word counts. The rows of the data frame correspond to the 85 individual essays and the columns correspond to the 5320 words in the vocabulary. We have already preprocessed the raw essays to remove very common and very infrequent words.\n",
    "\n",
    "Using this data, we will fit a topic model and do some analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "GFN1-wvIue3K"
   },
   "outputs": [],
   "source": [
    "import torch\n",
    "from torch.distributions import Dirichlet, Multinomial, Categorical\n",
    "import pandas as pd \n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm.auto import trange"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DP4qzh_Kue3L"
   },
   "source": [
    "## Load the data\n",
    "\n",
    "We've already tokenized the text and created a bag-of-words representation of the corpus. We removed words from the vocabulary that occur in more than 95% of the essays or only appear in 1 essay."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "2kDPPoqpvQ-x",
    "outputId": "4a14bcf1-05a1-441b-838b-68c69cb8a241",
    "tags": [
     "hide-output"
    ]
   },
   "outputs": [],
   "source": [
    "# Download the data\n",
    "!wget -nc https://www.dropbox.com/s/p6jb2cw5w5626pl/tokenized_fed.csv\n",
    "!wget -nc https://www.dropbox.com/s/ftedra0jyk1j3hx/authorship.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 487
    },
    "id": "Ig2BB1-Eue3M",
    "outputId": "8d8cf501-7c6e-42ce-875d-0e04e337c164"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "  <div id=\"df-e6fdb306-54fd-4327-bd14-55de9f371f9b\">\n",
       "    <div class=\"colab-df-container\">\n",
       "      <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>unequivocal</th>\n",
       "      <th>experience</th>\n",
       "      <th>inefficacy</th>\n",
       "      <th>subsisting</th>\n",
       "      <th>federal</th>\n",
       "      <th>called</th>\n",
       "      <th>deliberate</th>\n",
       "      <th>new</th>\n",
       "      <th>constitution</th>\n",
       "      <th>united</th>\n",
       "      <th>...</th>\n",
       "      <th>chancery</th>\n",
       "      <th>jurisprudence</th>\n",
       "      <th>reexamination</th>\n",
       "      <th>writ</th>\n",
       "      <th>commonlaw</th>\n",
       "      <th>intent</th>\n",
       "      <th>refutation</th>\n",
       "      <th>habeas</th>\n",
       "      <th>corpus</th>\n",
       "      <th>clerks</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>80</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>81</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>82</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>13.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>...</td>\n",
       "      <td>7.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>83</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>26.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>84</th>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>13.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>85 rows × 5320 columns</p>\n",
       "</div>\n",
       "      <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-e6fdb306-54fd-4327-bd14-55de9f371f9b')\"\n",
       "              title=\"Convert this dataframe to an interactive table.\"\n",
       "              style=\"display:none;\">\n",
       "        \n",
       "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
       "       width=\"24px\">\n",
       "    <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n",
       "    <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n",
       "  </svg>\n",
       "      </button>\n",
       "      \n",
       "  <style>\n",
       "    .colab-df-container {\n",
       "      display:flex;\n",
       "      flex-wrap:wrap;\n",
       "      gap: 12px;\n",
       "    }\n",
       "\n",
       "    .colab-df-convert {\n",
       "      background-color: #E8F0FE;\n",
       "      border: none;\n",
       "      border-radius: 50%;\n",
       "      cursor: pointer;\n",
       "      display: none;\n",
       "      fill: #1967D2;\n",
       "      height: 32px;\n",
       "      padding: 0 0 0 0;\n",
       "      width: 32px;\n",
       "    }\n",
       "\n",
       "    .colab-df-convert:hover {\n",
       "      background-color: #E2EBFA;\n",
       "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
       "      fill: #174EA6;\n",
       "    }\n",
       "\n",
       "    [theme=dark] .colab-df-convert {\n",
       "      background-color: #3B4455;\n",
       "      fill: #D2E3FC;\n",
       "    }\n",
       "\n",
       "    [theme=dark] .colab-df-convert:hover {\n",
       "      background-color: #434B5C;\n",
       "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
       "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
       "      fill: #FFFFFF;\n",
       "    }\n",
       "  </style>\n",
       "\n",
       "      <script>\n",
       "        const buttonEl =\n",
       "          document.querySelector('#df-e6fdb306-54fd-4327-bd14-55de9f371f9b button.colab-df-convert');\n",
       "        buttonEl.style.display =\n",
       "          google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
       "\n",
       "        async function convertToInteractive(key) {\n",
       "          const element = document.querySelector('#df-e6fdb306-54fd-4327-bd14-55de9f371f9b');\n",
       "          const dataTable =\n",
       "            await google.colab.kernel.invokeFunction('convertToInteractive',\n",
       "                                                     [key], {});\n",
       "          if (!dataTable) return;\n",
       "\n",
       "          const docLinkHtml = 'Like what you see? Visit the ' +\n",
       "            '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
       "            + ' to learn more about interactive tables.';\n",
       "          element.innerHTML = '';\n",
       "          dataTable['output_type'] = 'display_data';\n",
       "          await google.colab.output.renderOutput(dataTable, element);\n",
       "          const docLink = document.createElement('div');\n",
       "          docLink.innerHTML = docLinkHtml;\n",
       "          element.appendChild(docLink);\n",
       "        }\n",
       "      </script>\n",
       "    </div>\n",
       "  </div>\n",
       "  "
      ],
      "text/plain": [
       "    unequivocal  experience  inefficacy  subsisting  federal  called  \\\n",
       "0           1.0         1.0         1.0         1.0      1.0     1.0   \n",
       "1           0.0         2.0         0.0         0.0      2.0     1.0   \n",
       "2           0.0         1.0         0.0         0.0      2.0     0.0   \n",
       "3           0.0         2.0         0.0         0.0      0.0     1.0   \n",
       "4           0.0         1.0         0.0         0.0      0.0     0.0   \n",
       "..          ...         ...         ...         ...      ...     ...   \n",
       "80          0.0         0.0         0.0         0.0      6.0     1.0   \n",
       "81          0.0         0.0         0.0         0.0     12.0     0.0   \n",
       "82          0.0         1.0         0.0         0.0      7.0     2.0   \n",
       "83          0.0         0.0         0.0         0.0      2.0     1.0   \n",
       "84          0.0         2.0         0.0         0.0      0.0     1.0   \n",
       "\n",
       "    deliberate  new  constitution  united  ...  chancery  jurisprudence  \\\n",
       "0          1.0  5.0           7.0     1.0  ...       0.0            0.0   \n",
       "1          0.0  2.0           0.0     3.0  ...       0.0            0.0   \n",
       "2          0.0  1.0           0.0     4.0  ...       0.0            0.0   \n",
       "3          0.0  0.0           0.0     2.0  ...       0.0            0.0   \n",
       "4          0.0  0.0           0.0     0.0  ...       0.0            0.0   \n",
       "..         ...  ...           ...     ...  ...       ...            ...   \n",
       "80         1.0  8.0          12.0     8.0  ...       1.0            1.0   \n",
       "81         0.0  2.0           4.0     5.0  ...       0.0            0.0   \n",
       "82         2.0  9.0          13.0     6.0  ...       7.0            1.0   \n",
       "83         0.0  9.0          26.0    12.0  ...       0.0            0.0   \n",
       "84         0.0  4.0          13.0     2.0  ...       0.0            0.0   \n",
       "\n",
       "    reexamination  writ  commonlaw  intent  refutation  habeas  corpus  clerks  \n",
       "0             0.0   0.0        0.0     0.0         0.0     0.0     0.0     0.0  \n",
       "1             0.0   0.0        0.0     0.0         0.0     0.0     0.0     0.0  \n",
       "2             0.0   0.0        0.0     0.0         0.0     0.0     0.0     0.0  \n",
       "3             0.0   0.0        0.0     0.0         0.0     0.0     0.0     0.0  \n",
       "4             0.0   0.0        0.0     0.0         0.0     0.0     0.0     0.0  \n",
       "..            ...   ...        ...     ...         ...     ...     ...     ...  \n",
       "80            5.0   1.0        1.0     0.0         0.0     0.0     0.0     0.0  \n",
       "81            0.0   0.0        0.0     1.0         0.0     0.0     0.0     0.0  \n",
       "82            1.0   0.0        5.0     2.0         2.0     1.0     1.0     1.0  \n",
       "83            0.0   2.0        0.0     0.0         0.0     3.0     3.0     1.0  \n",
       "84            0.0   0.0        0.0     0.0         1.0     0.0     0.0     0.0  \n",
       "\n",
       "[85 rows x 5320 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load the data\n",
    "df = pd.read_csv('tokenized_fed.csv', index_col = 0)\n",
    "docs = torch.tensor(df.to_numpy()).type(torch.int)\n",
    "vocab = df.columns.to_list()\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7Nyg-s9Yx1JN"
   },
   "source": [
    "## Write some helper fucntions for Dirichlet distributions\n",
    "\n",
    "Specifically, we need the expected log of a Dirichlet random vector and the KL divergence between two Dirichlet random variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "id": "UoqrOOfvue3N"
   },
   "outputs": [],
   "source": [
    "def dirichlet_expected_log(dirichlet):\n",
    "    \"\"\"Helper function to compute expected log under Dirichlet distribution.\n",
    "\n",
    "    Args:\n",
    "        dirichlet: A torch.distributions.Dirichlet object with a batch shape of\n",
    "            (...,) and a event shape of (K,).\n",
    "\n",
    "    Returns:\n",
    "        (...,K) tensor of expected logs, E[\\log \\pi], under the Dirichlet.\n",
    "    \"\"\"\n",
    "    concentration = dirichlet.concentration\n",
    "    return torch.special.digamma(concentration) - \\\n",
    "           torch.special.digamma(concentration.sum(dim=-1, keepdims=True))\n",
    "\n",
    "\n",
    "def dirichlet_log_normalizer(concentration):\n",
    "    \"\"\"Compute the log normalizing constant of a Dirichlet distribution with\n",
    "    the specificed concentration.\n",
    "\n",
    "    Args:\n",
    "        concentration: (...,K) tensor of concentration parameters\n",
    "\n",
    "    Returns:\n",
    "        (...,) batch of log normalizers\n",
    "    \"\"\"\n",
    "    return torch.special.gammaln(concentration).sum(dim=-1) - \\\n",
    "        torch.special.gammaln(concentration.sum(dim=-1))\n",
    "\n",
    "def dirichlet_kl(q, p):\n",
    "    \"\"\"Compute the KL divergence between two Dirichlet disdtributions\n",
    "\n",
    "    Args:\n",
    "        q: A torch.distributions.Dirichlet object\n",
    "        p: A torch.distributions.Dirichlet object over the same domain\n",
    "\n",
    "    Returns:\n",
    "        A (batch of) KL divergence(s) between q and p.\n",
    "    \"\"\"\n",
    "    kl = -dirichlet_log_normalizer(q.concentration)\n",
    "    kl += dirichlet_log_normalizer(p.concentration)\n",
    "    kl += torch.sum((q.concentration - p.concentration) * \\\n",
    "                    dirichlet_expected_log(q), dim=-1)\n",
    "    return kl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-bCID4jque3N"
   },
   "source": [
    "## Implement Coordinate Ascent Variational Inference (CAVI)\n",
    "\n",
    "_Note: The `torch.distributions.Multinomial` object doesn't work well when you have a batch with different numbers of counts. We hijack this object by not giving it a count so that it defaults to 1, which is equivalent to a categorical distribution. Then we multiply by the total counts to get the necessary expectations under the multinomial posterior._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "id": "x-TIiIlxue3O"
   },
   "outputs": [],
   "source": [
    "def cavi(docs, \n",
    "         num_topics=10, \n",
    "         num_iters=200, \n",
    "         tol=1e-5, \n",
    "         alpha=20.0, \n",
    "         phi=1.0,\n",
    "         seed=305 + ord('c'),\n",
    "        ):\n",
    "    \"\"\"Run coordinate ascent VI for LDA.\n",
    "    \n",
    "    \"\"\"\n",
    "    docs = docs.type(torch.float)\n",
    "    N, V = docs.shape\n",
    "    K = num_topics      # short hand\n",
    "    \n",
    "    def cavi_step(q_c, q_pi, q_theta):\n",
    "        \"\"\"One step of CAVI.\n",
    "        \"\"\"\n",
    "        # Update the topic assignment counts\n",
    "        E_logpi = dirichlet_expected_log(q_pi)\n",
    "        E_logtheta = dirichlet_expected_log(q_theta)\n",
    "        q_c = Multinomial(logits=E_logpi[:, None, :] + E_logtheta.T)\n",
    "        \n",
    "        # Compute the mean of q(c) since we'll use it twice below\n",
    "        E_c = docs.unsqueeze(2) * q_c.mean\n",
    "        \n",
    "        # Update the topic proportions\n",
    "        q_pi = Dirichlet(alpha + E_c.sum(axis=1))\n",
    "        \n",
    "        # Update the topic word probabilities\n",
    "        q_theta = Dirichlet(phi + E_c.sum(axis=0).T)\n",
    "        \n",
    "        return q_c, q_pi, q_theta\n",
    "    \n",
    "    def elbo(q_c, q_pi, q_theta):\n",
    "        \"\"\"Compute the evidence lower bound.\n",
    "        \"\"\"\n",
    "        elbo = 0\n",
    "        \n",
    "        # KL to prior\n",
    "        elbo -= dirichlet_kl(q_pi, Dirichlet(alpha * torch.ones(K))).sum()\n",
    "        elbo -= dirichlet_kl(q_theta, Dirichlet(phi * torch.ones(V))).sum()\n",
    "        \n",
    "        # Entropy of q(z) [a little different from multinomial entropy]\n",
    "        E_c = docs.unsqueeze(2) * q_c.mean\n",
    "        elbo -= torch.sum(E_c * torch.log(q_c.probs))\n",
    "        \n",
    "        # Expected log p(z | \\theta)\n",
    "        E_logpi = dirichlet_expected_log(q_pi)\n",
    "        E_logtheta = dirichlet_expected_log(q_theta)\n",
    "        elbo += torch.sum(E_c * E_logpi[:, None, :])\n",
    "        elbo += torch.sum(E_c * E_logtheta.T)\n",
    "        \n",
    "        return elbo / torch.sum(docs)\n",
    "    \n",
    "    # Initialize the topics by randomly clustering the documents\n",
    "    # and using their word counts\n",
    "    torch.manual_seed(seed)\n",
    "    clusters = Categorical(logits=torch.zeros(K)).sample((N,))\n",
    "    q_pi = Dirichlet(alpha * torch.ones((N, K)))\n",
    "    q_theta = Dirichlet(phi + torch.row_stack([docs[clusters == k].sum(axis=0) \n",
    "                                               for k in range(K)]))\n",
    "    q_c = Multinomial(logits=torch.zeros((N, V, K)))\n",
    "    \n",
    "    # Run CAVI\n",
    "    elbos = [elbo(q_c, q_pi, q_theta)]\n",
    "    for itr in trange(num_iters):\n",
    "        q_c, q_pi, q_theta = cavi_step(q_c, q_pi, q_theta)\n",
    "        elbos.append(elbo(q_c, q_pi, q_theta))\n",
    "        \n",
    "        if elbos[-1] - elbos[-2] < -1e-4:\n",
    "            raise Exception(\"ELBO is going down!\")\n",
    "        elif elbos[-1] - elbos[-2] < tol:\n",
    "            print(\"Converged!\")\n",
    "            break\n",
    "        \n",
    "    return torch.tensor(elbos), (q_c, q_pi, q_theta)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 347,
     "referenced_widgets": [
      "678444446c6b403fa165ecf9abbb7624",
      "eb674660b30041baa6de6eddc6926cf6",
      "10bfdab627d24288bffa6b00d1214955",
      "02f48be17ec246409841be7c4429dff9",
      "b14a145cf2154b1a9366a34350b0d2c6",
      "eb963e8602844150b388f15247e9c90d",
      "50f0438726d74985bb222d18e2844996",
      "1b30e510e3de423d9427e14306e4a54f",
      "d54823f5694b40f69c97eb47f948d916",
      "4e68636d615c443f9be9a18071ed02d8",
      "6189223a819546c9bcf080bf463f9fe8"
     ]
    },
    "id": "7iZzXB-9ue3O",
    "outputId": "34da3019-f44d-45c6-fe68-1554460b3d46"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "678444446c6b403fa165ecf9abbb7624",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/200 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converged!\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'ELBO per word')"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXzddZ3v8dc7e5ou6crShYJlUVEKDUhRtsp11OEK13FgZlDAASqiI3DHBeUxXJlxHGVQh9G5KoIsCsoAMnB1EEVFQdYApYBby562ULqkW5rkLJ/7x++X5iRN0tMkJydJ38/H4zzOb/n+zvmc/uD7yff7/f1+X0UEZmZmu6ui3AGYmdnY5ARiZmaD4gRiZmaD4gRiZmaD4gRiZmaDUlXuAEbSjBkzYv78+eUOw8xsTHn88cfXRcTM3tv3qAQyf/58mpubyx2GmdmYIumlvra7C8vMzAbFCcTMzAbFCcTMzAbFCcTMzAbFCcTMzAbFCcTMzAbFCcTMzAZlj7oPxMysP/l80JnL05HN05HN0ZnNJ69cvsdyJhdkc3kyuTzZfJBLX9l8EBHk8pCPIICIICJ9h2Q53U66DhBEwXJXuZ5lKDiOHp+182d0fU7hB5x1zHymT6wdxn8xJxAzK7OIpOJuz+Rpz+TY3pmjPZvbsd6R7X7v6HrPdlfqOyr7goq+o2A5k+uVCHolhK7lbH78zo0kwfsWznYCMbPy6czm2daRZWv62taRZVtnLnnvyNLWmWNbZ5a2juR9e2eObZ05tndm2daRoy2To70zR1smy/bONGFkcuSGUHlXVYiaqorkVdnrvWB5Yl0VNZUVVFdVUNvH/sL12urKPstU73gXVRXJe0WFqJSorOherhBIQgKRLkO6nmws3Afs2J8sa8cyhccVlOk+pmtPWq6wQIk5gZjtAfL5YFtnls3tWTZvz7B5e4Yt7Vk2t2fY2pFlS3vXK1nf2p5lS/relSy2tmfpzOWL+j4JGmqqmFBTmb6S5cl1Vew9uZYJNVXUVVdQX51sr6uuoL6mivrqZLmu672qktrqCmqrkvXaqkpqqrq311RWUFExchWm9eQEYjaGdGRzbGrLsLEtQ2tbJxvbMmza3klrW4ZN2zO0bk/eN7VlaN3eyebtWTZtz7ClPcOu/sivqapgUm0VE+uqmFRXxcTaKvZtrGNSXTUNtZU01FYxqbaKhvQ1sWu5pnLHen1NJRNrq6itqhjRv4StPJxAzMogItjSkaV1W4aNbZ1saOtMEkK6vrErObQl663pe1tnrt/PrKoQU+qrmTKhmsb6amZOrGXBzIlMrq9mcl11sq++mkl1VUxO3yfVdb1XUVtVOYL/AjYeOIGYDYNsLs/Gtgzrt3Wwfmsn67d1smFrBxu2Jcsb2zrZsC1JEOu3Jcmiv0FbCRrrq5k6oYbGCdXsNbmOg/eexNQJNUydUE1jur2xPn2fkJSdUFPpv/ptRDmBmPWjq5WwdnM7r23uYO2W5P31LR2s29r9vm5rkiCij3zQlQymNdQwraGG+TMmcMR+jTROqGFamgimNdTQmCaHaQ01TK6rdr++jQllSSCSFgLfAuqALHBBRDzaR7mfAkcDD0TEyQXbrweOBzalm86OiGWljtvGj4hgc3uW1a3bk9emdta0bufVTe2s2dTOa5vbeXVze59dRvXVlcycVMuMiTXsP6OBpvnTmNFQw4xJtUxvqGVaQw3TJ9YwPU0MlU4GNk6VqwVyBXB5RNwt6b3p+gl9lPtXYALwkT72fSoibitdiDaWZXJ5Xt3UniaH7axubWdVV7JoTda3dmR7HFNZIfaaVMs+jfW8cd/JnHDwLPaeUstek+uYNamOWZNrmTWplom1Ve4qMqN8CSSAyenyFGB1n4UifiHphJEKysaGiGBjW4bVrdtZ1bqdNWkLojA5vLalfacupWkNNcxurGf/GQ0c84YZzG6sZ9/GevZprGPfKfXMnFTr1oLZbihXArkIuEfSlSTP4zpmEJ/xz5IuA34BXBIRHX0VkrQUWAowb968QYZrIymfD17f2sHLG9po2dhGy4YkUXS9Vrdupz3T836EmqoK9plSx+zGet5x4Az2nVLH7KlJgti3sZ59p9RTX+OrjMyGU8kSiKR7gb372HUp8E7g4oi4XdJpwLXASbvx8Z8FXgVqgKuBzwD/2FfBiLg6LUNTU9P4fVbBGLSlPcOKtVtZ+dpWnnt9K8+v28aL67bx0oY2OrM9E8SMibXMnlrPIXtPYsnBs7oTQ2Md+zbWM72hxt1KZiOsZAkkIvpNCJJuBC5MV28FrtnNz16TLnZIug745KCCtBGzcVsny1paWf7KJn63ZhPPrt5My8btO/bXVFaw3/QJzJ/RwAkHz2TetAnMnTaBOVMnMGdqPXXVbj2YjTbl6sJaTXIV1X3AEmDF7hwsaZ+IWKPkT85TgWeGPUIbtIhg5dqtPPbiRppf2sATL23kxfVtQHJZ6/7TGzhsbiN/fdQ8DtprEgfOmsjcaRM8/mA2xpQrgZwHXCWpCmgnHaOQ1AScHxHnpuv3A4cAEyW1AOdExD3ATZJmkjxLbBlwfhl+g6UighVrt/LQc+t5+Pn1PPrCBtZv6wRgekMNi/abyl8dNY/D5jTyljlTmFjr24/MxoOy/J8cEQ8Ai/rY3gycW7B+bD/HLylddLYrEcEL67bx4HPreej59Tzy/HrWbU0SxuzGeo4/eCZH7z+dpvlT2X9Gg8cmzMYp/yloRWnZ2JYkjOfW8+Bz63htc3LR296T6zj2wJksfsN0Fh8wnbnTJpQ5UjMbKU4g1qdN2zM89Nw6frNiHb9duY6X0jGM6Q01SbJ4w3SOecMM5k+f4BaG2R7KCcSApFvq92u28Ms/vMYv/7CWZa+0kg+YWFvF0QdM46zF8zlmwXQO3muSE4aZAU4ge7R8Pnji5Y3c/cyr/PSZV1nVmlxW+5bZU/j4iQs49qCZLJzbSHVlRZkjNbPRyAlkD/Tc61v50RMt3PHEKlZvaqemsoJjD5zBJ965gBMPnsWsyXXlDtHMxgAnkD1EJpfnZ8++xg0PvcijL2ygQnDcQTP5zHsOYckhs5hUV13uEM1sjHECGee2dWS5+ZGXufaBF3h1cztzptZzyXsO4f2Hz3ZLw8yGxAlknNrWkeX6B1/kmvufZ2NbhsUHTOcLpx7KiYfM8h3fZjYsnEDGmWwuzy3Nr/C1n69g3dYOlhwyi4+duIBF+00td2hmNs44gYwjDz23nn+48xlWrt1K035TufrMRRwxz4nDzErDCWQcWLe1gy/+5Pf86MlVzJ1Wz7c/tIh3vWkv369hZiXlBDLG/ffTa7j0jqfZ2pHl4ycu4GMnLvDESWY2IpxAxqhN2zN8/q5nuePJVRw2ZwpfOe0wFsyaVO6wzGwP4gQyBi1vaeWCm55gzaZ2LjrpQD524gLfLW5mI84JZAyJCG5+9GUuv+t3zJhYw63nL/YguZmVjRPIGNGZzXPZnc/ww8de4biDZvJvpy9kWkNNucMysz2YE8gYsKktw0dvepwHn1vPx09cwMX/4yDfDGhmZVeWjnNJCyU9LGmZpGZJR/VT5iFJz0paLun0gn37S3pE0kpJt0gat3+Kv7Khjfd/87c89uIGvvKXh/HJPzvYycPMRoVyjbxeAVweEQuBy9L13tqAMyPizcC7gX+T1Jju+zLwtYhYAGwEzhmBmEfcite28IFvPci6rZ1875y38ReL5pQ7JDOzHcqVQAKYnC5PAVbvVCDiTxGxIl1eDawFZiq5O24JcFta9Abg1JJHPMKWt7Ry2rcfIgL+8yOLOfqA6eUOycysh3KNgVwE3CPpSpIkdsxAhdMurhrgOWA60BoR2XR3CzC7hLGOuOYXN3D2dY8xtaGa75/zNvab3lDukMzMdlKyBCLpXmDvPnZdCrwTuDgibpd0GnAtcFI/n7MP8D3grIjI7+7jOSQtBZYCzJs3b7eOLYcnXt7I2dc9xqxJtdx83tHsPcWPXDez0UkRMfJfKm0CGiMi0i6pTRExuY9yk4H7gC9GxG3pNgGvA3tHRFbSYuDzEfFnu/repqamaG5uHs6fMqyWt7RyxnceYdrEGm5ZutjJw8xGBUmPR0RT7+3lGgNZDRyfLi8BVvQukF5ZdQdwY1fyAIgk4/0K+EC66SzgzpJGOwL++OoWPnTtozQ2VPMDtzzMbAwoVwI5D/iKpKeAL5J2MUlqknRNWuY04Djg7PRy32WSFqb7PgP8b0krScZErh3Z8IdXy8Y2zvzuI9RVV3DzuUezb2N9uUMyM9ulsnRhlcto7MLasK2TD3zzQdZt7eDW84/h4L39QEQzG13668LynehltL0zx4evf4xVrdv5/rlvc/IwszHFCaRMcvngolueZHlLK9/+4CKOnD+t3CGZme0WPwO8TL509++559nXuOzkN/GuN/d1tbOZ2ejmBFIG33v4Jb5z/wuctXg/Pvz2/csdjpnZoDiBjLDf/Ol1Pn/Xs5x48Ez+4eQ3lTscM7NBcwIZQSvXbuFjNz3BgbMm8vW/OYIqzyJoZmOYa7ARsmFbJ397fTO11ZVcc1YTE2t9/YKZjW2uxUZAeybHR77XzKub2/nh0qOZM3VCuUMyMxsyt0BKLJ8PPn3bch57cSNfPe0wz2FuZuOGE0iJffXnf+Kup1bz6XcfzMlv3bfc4ZiZDRsnkBL64aMv841freSvj5rLR49/Q7nDMTMbVk4gJfLTZ9bwuTue5viDZvKPpxzK7s5jYmY22jmBlMCDK9fxiR8sY+HcRr75wSOo9uW6ZjYOuWYbZstbWjnvxmb2n9HAd88+kgk1vtDNzMYnJ5BhtHLt1nQu8xpuPOcoGifUlDskM7OScQIZJqtat3PmtY9QIfj+OW9jr8meUdDMxjcnkGGwpT3Dmdc+wpb2LDf87VHMn9FQ7pDMzErOHfTD4Ms//QPPr9vGzecezZv3nVLucMzMRkRZWiCSFkp6OJ3nvFnSUf2UeUjSs5KWSzq9YN/1kl7oY670Edf84ga+//DLnH3MfBa/YXq5wjAzG3HlaoFcAVweEXdLem+6fkKvMm3AmRGxQtK+wOOS7omI1nT/pyLitpELeWcd2Ryf/dHTzG6s55PvOricoZiZjbhyJZAAJqfLU4DVOxWI+FPB8mpJa4GZQGvvsuXyrfueZ8XarVx39pE0+Om6ZraHKdcg+kXAv0p6BbgS+OxAhdMurhrguYLN/5x2bX1NUu0Axy5Nu8maX3/99eGIHYC2ziz/cd9K/vwt+3DiIbOG7XPNzMaKkiUQSfdKeqaP1ynAR4GLI2IucDFw7QCfsw/wPeDDEZFPN38WOAQ4EpgGfKa/4yPi6ohoioimmTNnDtOvg2dXb6Yzm+d/HT572D7TzGws6bffRdLXSbqa+hQRnxjogyPipAE++0bgwnT1VuCafspNBn4CXBoRDxd89pp0sUPSdcAnB4qlFJ56JelJe+tcX3VlZnumgVogzcDjQB1wBLAifS0k6U4aitXA8enykvRze5BUA9wB3Nh7sDxtlaDkCYWnAs8MMZ7d9lTLJvaZUsesSb5h0Mz2TP22QCLiBgBJHwXeERHZdP1bwP1D/N7zgKskVQHtwNL0s5uA8yPiXOA04DhguqSz0+POjohlwE2SZgIClgHnDzGe3ba8pZW3znHrw8z2XMVcOjSV5IqpDen6xHTboEXEA8CiPrY3A+emy98Hvt/P8UuG8v1D1drWyUvr2zj9yLnlDMPMrKyKSSBfAp6U9CuSv/iPAz5fyqBGu+UtmwA4bE5jmSMxMyufAROIpArgj8Db0hfAZyLi1VIHNpotb0kG0A+d7S4sM9tzDZhAIiIv6T8i4nDgzhGKadR7qmUTB8xoYEp9dblDMTMrm2LuA/mFpL+Q52Td4alXPIBuZlZMAvkIyb0anZK2pK/NJY5r1Hp1Uztrt3TwVo9/mNkebpeD6BExaSQCGSueSsc/DvMNhGa2hyvqCYCS3kdy9RXAfRHx49KFNLotb2mlskKe98PM9ni77MKS9CWSx478Ln1dKOlfSh3YaLW8ZRMH7zWJuurKcodiZlZWxbRA3gss7HqQoaQbgCfZxRN0x6sX1m1j0X5Duo/SzGxcKPZpvIUjxnt0300ml6ferQ8zs6JaIP/CzneiX1LSqEaxTC6orizXNCpmZqNHMVdh/UDSfSRzb8Aefid6Jpt3AjEzo4gEIun7wK+B+yPiD6UPaXTrzOWprvI9lWZmxfwpfS2wD/B1Sc9Lul3Shbs6aLzK5PLUuAViZlZUF9avJP2GpAvrRJK5N94MXFXi2EadXD7IB+7CMjOjuC6sXwANwEMkE0kdGRFrSx3YaJTJJVOyO4GYmRXXhbUc6AQOBd4KHCqpvqRRjVKdOxKIx0DMzHaZQCLi4og4Dng/sB64Dmgd6hdLWijpYUnLJDVLOqqPMvtJeiIt86yk8wv2LZL0tKSVkv59JJ4WnMkmCaSmyi0QM7NiHmXycUm3kNx9fgrwXeA9w/DdVwCXR8RC4LJ0vbc1wOK0zNuASyTtm+77Jsnc6gemr3cPQ0wDyuQCcBeWmRkUdyNhHfBV4PGIyA7jdwfJXOuQ3N2+eqcCEZ0Fq7WkCU/SPsDkiHg4Xb8ROBW4exjj24nHQMzMuhVzFdaVJfrui4B7JF1JkhiO6auQpLnAT4AFwKciYrWkJqCloFgLMLuf45cCSwHmzZs3pIA9BmJm1q2kf0pLulfSM328TgE+ClwcEXOBi0nuN9lJRLwSEW8lSSBnSdprd2KIiKsjoikimmbOnDmk39PVAvF9IGZmRc4HMlgRcVJ/+9Jup64bEm8FrtnFZ62W9AxwLPBbYE7B7jnAqqFFu2udWXdhmZl1GbAmlFSZPkSxFFYDx6fLS4AVfXz/nK5LhiVNBd4B/DEi1gCbJR2dXn11JnBnieLcYccYiK/CMjMbuAUSETlJeUlTImLTMH/3ecBVkqqAdtJxinR84/yIOBd4I/AVSUHyJOArI+Lp9PgLgOuBepLB85IOoAN0ZruuwvIYiJlZMV1YW4GnJf0c2Na1MSI+MZQvjogHgEV9bG8Gzk2Xf05y82JfxzeT3Nw4YjwGYmbWrZgE8qP0tcfzZbxmZt2KuYz3hnQcYl5E/HEEYhq1nEDMzLoVcyf6/wSWAT9N1xdKuqvUgY1Gnemd6DWeD8TMrKj7QD4PHEX6/KuIWAYcUMKYRq2ML+M1M9uhmJow08cVWPlSBDPauQvLzKxbMYPoz0r6G6BS0oHAJ4AHSxvW6OQEYmbWrZia8O9IZiDsAH4AbCZ5jtUeZ8cYiBOImVlRV2G1AZdK+nKyGltKH9bo1H0nugfRzcyKuQrrSElPk8xM+LSkpyTtdAPgnmDHhFJugZiZFTUGci1wQUTcDyDpHSSzEvZ5h/h4lsnlkaCywi0QM7Ni/pTOdSUP2PEIkuGcWGrM6MwF1ZUVjMDsuWZmo14xLZBfS/o2yQB6AKcD90k6AiAinihhfKNKJpd395WZWaqYBHJY+v5/em0/nCShLBnWiEaxTC7vJ/GamaWKuQrrxJEIZCxIEohbIGZmUOIpbcebzmw4gZiZpVwb7oZMLk+NZyM0MwOcQHaLx0DMzLoNOAYiaRbwMZJHmQA8C/zfiHhtKF8qaSHwLaCO5JLgCyLi0V5l9gPuIEly1cDXI+Jb6b77gH2A7Wnxd0XE2qHEVAyPgZiZdeu3NpT0duCxdPXG9AXwSLpvKK4ALo+IhcBl6Xpva4DFaZm3AZdI2rdg/xkRsTB9lTx5QPd9IGZmNnAL5CvAqRHxZMG2uyTdAXybpFIfrAAmp8tTgNU7FYjoLFitZRR0t2Wyvg/EzKzLQAlkcq/kASQTSkmaNMTvvQi4R9KVJInhmL4KSZoL/ARYAHwqIgoTzXWScsDtwBciIvr5jKXAUoB58+YNKehMLk9ttROImRkM/Fe9JE3tY+O0XRzXVe5eSc/08ToF+ChwcUTMBS4med7WTiLilYh4K0kCOUvSXumuMyLiLcCx6etD/cUREVdHRFNENM2cOXNXYQ/IYyBmZt0Gqg2/BvxM0vGSJqWvE4C7030DioiTIuLQPl53AmcBP0qL3koyZe5An7UaeIYkWRARq9L3LcDNuzp+uHgMxMysW79dWBFxtaTVwD/R8yqsL0TE/xvi964GjgfuI3kUyoreBSTNAdZHxPa0JfQO4GuSqoDGiFgnqRo4Gbh3iPEUxc/CMjPrNuBlvBHxY+DHJfje84Cr0mTQTjpGIakJOD8izgXeCHxFUgACroyIpyU1kIyfVAOVJMnjOyWIcSe+D8TMrFu/CUTSDJJ7QDaQzP9xBXAc8Bzw9xGxcrBfmj4SfqdJqSKiGTg3Xf45fcw5EhHb+jp2JGSyHgMxM+syUG14M8nlswcBjwIvAh8gaZFcU/LIRqHOXFDtR5mYmQEDd2HtFRGfUzJ70ksR0XWz3x8kfWwEYht1PAZiZtZtoNowB5DeX7Gu1758ySIaxTwGYmbWbaAWyAGS7iIZwO5aJl3fv+SRjUK+D8TMrNtACeSUguUre+3rvT7u5fNBxveBmJntMNB9IL/ub5+kW4B+949HmXzSa+f5QMzMEoOtDRcPaxRjQCaXPGrLYyBmZgn/OV2kTDZtgbgLy8wMGPhGwiP620UywdMeJZNLEojvAzEzS+xqPpD+/GG4AxntOrsSiFsgZmbAwIPoJ45kIKNd1xiIu7DMzBIDTWn76YLlv+y174ulDGo0yrgFYmbWw0C14V8VLH+21753lyCWUa0z25VAfBWWmRnsYkbCfpb7Wh/3PIhuZtbTQLVh9LPc1/q45zEQM7OeBroK6zBJm0laG/XpMul6XckjG2U8BmJm1tNAV2FVjmQgo133Zbx7XO+dmVmfyvbntKSFkh6WtExSs6SjBig7WVKLpG8UbFsk6WlJKyX9ezpvSclksm6BmJkVKmdteAVweUQsBC5L1/vzT8Bvem37Jsnc6gemr5JeGbZjDMSD6GZmQHkTSACT0+UpwOq+CklaBOwF/Kxg2z7A5Ih4OJ3w6kbg1FIG6zEQM7OeBhpEL7WLgHskXUmSyI7pXUBSBckjVT4InFSwazbQUrDekm7biaSlwFKAefPmDTpYj4GYmfVU0gQi6V5g7z52XQq8E7g4Im6XdBpwLT2TBMAFwH9HRMtghzgi4mrgaoCmpqZBX37c1QLxZbxmZomSJpCI6J0QdpB0I3BhunorcE0fxRYDx0q6AJgI1EjaClwFzCkoNwdYNSxB98OD6GZmPZWzC2s1cDxwH7AEWNG7QESc0bUs6WygKSIuSdc3SzoaeAQ4E/h6KYPdMaGUB9HNzIDyJpDzgKskVQHtpOMUkpqA8yPi3F0cfwFwPVAP3J2+SsZjIGZmPZUtgUTEA8CiPrY3Azslj4i4niRhFJY7tHQR9rTjKqwKt0DMzMBT2hYtk8tTVSEqKtwCMTMDJ5CiZXLhAXQzswKuEYvUmc17/MPMrIATSJEyubwfY2JmVsA1YpEyuby7sMzMCrhGLFImF26BmJkVcI1YpE63QMzMenCNWKRM1gnEzKyQa8QiZXJ5anwVlpnZDk4gRfJ9IGZmPblGLJLHQMzMenKNWKRMLu8n8ZqZFXCNWCSPgZiZ9eQEUqROX4VlZtaDa8QieRDdzKwn14hFcgvEzKwn14hFSh6m6DEQM7MuZUkgkhZKeljSMknNko4aoOxkSS2SvlGw7T5Jf0yPXyZpVqlj9sMUzcx6KteUtlcAl0fE3ZLem66f0E/ZfwJ+08f2M9JpbUeEx0DMzHoqV40YwOR0eQqwuq9CkhYBewE/G6G4+uUbCc3MeipXC+Qi4B5JV5IksWN6F5BUAXwF+CBwUh+fcZ2kHHA78IWIiFIFGxG+D8TMrJeSJRBJ9wJ797HrUuCdwMURcbuk04Br2TlJXAD8d0S0SDtV3GdExCpJk0gSyIeAG/uJYymwFGDevHmD+i25fBCBWyBmZgVKlkAioq9WAwCSbgQuTFdvBa7po9hi4FhJFwATgRpJWyPikohYlX7HFkk3A0fRTwKJiKuBqwGampoG1UrJ5JLD/CgTM7Nu5erCWg0cD9wHLAFW9C4QEWd0LUs6G2iKiEskVQGNEbFOUjVwMnBvKYPtzOUBt0DMzAqVK4GcB1yVJoN20i4mSU3A+RFx7gDH1pKMn1QDlSTJ4zulDDaTJhCPgZiZdStLAomIB4BFfWxvBnZKHhFxPXB9urytr2NLKeMWiJnZTlwjFiGTTcdAnEDMzHZwjViEHWMgHkQ3M9vBNWIRPAZiZrYzJ5AieAzEzGxnrhGLsKMF4i4sM7MdXCMWodOD6GZmO3GNWAR3YZmZ7cw1YhG6B9H9z2Vm1sU1YhF2tEA8I6GZ2Q5OIEXozHkMxMysN9eIRchk3YVlZtaba8QieBDdzGxnrhGL0J1APAZiZtbFCaQInZ5QysxsJ64Ri+DLeM3MduYasQhdg+geAzEz6+YasQiZXJ4KQWWFx0DMzLo4gRShMxdufZiZ9VKWWlHSQkkPS1omqVnSUf2Uy6Vllkm6q2D7/pIekbRS0i2SakoZbyaX9/iHmVkv5aoVrwAuj4iFwGXpel+2R8TC9PW+gu1fBr4WEQuAjcA5pQw2k8v7Ciwzs17KVSsGMDldngKsLvZASQKWALelm24ATh3W6HrJ5PK+B8TMrJeqMn3vRcA9kq4kSWLH9FOuTlIzkAW+FBH/BUwHWiMim5ZpAWb390WSlgJLAebNmzeoYDuyeY+BmJn1UrIEIuleYO8+dl0KvBO4OCJul3QacC1wUh9l94uIVZIOAH4p6Wlg0+7EERFXA1cDNDU1xe4c2yWTC4+BmJn1UrIEEhF9JQQAJN0IXJiu3gpc089nrErfn5d0H3A4cDvQKKkqbYXMAVYNY+g7ybgFYma2k3LViquB49PlJcCK3gUkTZVUmy7PAN4O/C4iAvgV8IG06FnAnaUMNhlE9xiImVmhco2BnAdcJakKaCcdo5DUBJwfEecCbwS+LSlPkui+FBG/S4//DPBDSV8AniTpAiuZI/abytaO7K4LmpntQZT8Qb9naGpqiubm5nKHYWY2pkh6PCKaem93x76ZmQ2KE4iZmQ2KE4iZmQ2KE4iZmQ2KE4iZmQ2KE+mJeLcAAAXtSURBVIiZmQ2KE4iZmQ2KE4iZmQ3KHnUjoaTXgZcGefgMYN0whjMa+TeOD/6NY99o+337RcTM3hv3qAQyFJKa+7oTczzxbxwf/BvHvrHy+9yFZWZmg+IEYmZmg+IEUryryx3ACPBvHB/8G8e+MfH7PAZiZmaD4haImZkNihOImZkNihNIESS9W9IfJa2UdEm54xkqSXMl/UrS7yQ9K+nCdPs0ST+XtCJ9n1ruWIdKUqWkJyX9OF3fX9Ij6bm8RVJNuWMcCkmNkm6T9AdJv5e0eLydR0kXp/+dPiPpB5Lqxvp5lPRdSWslPVOwrc/zpsS/p791uaQjyhd5T04guyCpEvgP4D3Am4C/lvSm8kY1ZFng7yPiTcDRwMfS33QJ8IuIOBD4Rbo+1l0I/L5g/cvA1yJiAbAROKcsUQ2fq4CfRsQhwGEkv3XcnEdJs4FPAE0RcShQCfwVY/88Xg+8u9e2/s7be4AD09dS4JsjFOMuOYHs2lHAyoh4PiI6gR8Cp5Q5piGJiDUR8US6vIWk0plN8rtuSIvdAJxangiHh6Q5wJ8D16TrApYAt6VFxvRvlDQFOA64FiAiOiOilXF2HoEqoF5SFTABWMMYP48R8RtgQ6/N/Z23U4AbI/Ew0Chpn5GJdGBOILs2G3ilYL0l3TYuSJoPHA48AuwVEWvSXa8Ce5UprOHyb8CngXy6Ph1ojYhsuj7Wz+X+wOvAdWk33TWSGhhH5zEiVgFXAi+TJI5NwOOMr/PYpb/zNmrrICeQPZikicDtwEURsblwXyTXd4/Za7wlnQysjYjHyx1LCVUBRwDfjIjDgW306q4aB+dxKslf4PsD+wIN7Nz1M+6MlfPmBLJrq4C5Betz0m1jmqRqkuRxU0T8KN38WlfTOH1fW674hsHbgfdJepGk23EJyXhBY9oVAmP/XLYALRHxSLp+G0lCGU/n8STghYh4PSIywI9Izu14Oo9d+jtvo7YOcgLZtceAA9OrPmpIBvDuKnNMQ5KOBVwL/D4ivlqw6y7grHT5LODOkY5tuETEZyNiTkTMJzlnv4yIM4BfAR9Ii4313/gq8Iqkg9NN7wR+xzg6jyRdV0dLmpD+d9v1G8fNeSzQ33m7CzgzvRrraGBTQVdXWflO9CJIei9Jf3ol8N2I+OcyhzQkkt4B3A88Tff4wOdIxkH+E5hH8tj70yKi90DfmCPpBOCTEXGypANIWiTTgCeBD0ZERznjGwpJC0kuEqgBngc+TPKH4bg5j5IuB04nuXrwSeBckjGAMXseJf0AOIHkse2vAf8H+C/6OG9p4vwGSdddG/DhiGguR9y9OYGYmdmguAvLzMwGxQnEzMwGxQnEzMwGxQnEzMwGxQnEzMwGxQnEbBAkbU3f50v6m2H+7M/1Wn9wOD/fbLg4gZgNzXxgtxJIwR3U/emRQCLimN2MyWxEOIGYDc2XgGMlLUvnraiU9K+SHkvnbvgIJDczSrpf0l0kd1Ij6b8kPZ7OdbE03fYlkifPLpN0U7qtq7Wj9LOfkfS0pNMLPvu+gnlBbkpvPjMrqV39JWRmA7uE9C53gDQRbIqIIyXVAr+V9LO07BHAoRHxQrr+t+mdxvXAY5Juj4hLJH08Ihb28V3vBxaSzPsxIz3mN+m+w4E3A6uB35I8L+qB4f+5Zt3cAjEbXu8ieW7RMpJHw0wnmQgI4NGC5AHwCUlPAQ+TPCzvQAb2DuAHEZGLiNeAXwNHFnx2S0TkgWUkXWtmJeUWiNnwEvB3EXFPj43J87i29Vo/CVgcEW2S7gPqhvC9hc+ByuH/t20EuAViNjRbgEkF6/cAH00fl4+kg9JJnnqbAmxMk8chJFMLd8l0Hd/L/cDp6TjLTJLZCB8dll9hNgj+K8VsaJYDubQr6nqSOUfmA0+kA9mv0/d0qz8Fzpf0e+CPJN1YXa4Glkt6In0EfZc7gMXAUySTDX06Il5NE5DZiPPTeM3MbFDchWVmZoPiBGJmZoPiBGJmZoPiBGJmZoPiBGJmZoPiBGJmZoPiBGJmZoPy/wHAwWvLunKiEgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "elbos, (q_c, q_pi, q_theta) = cavi(docs)\n",
    "\n",
    "plt.plot(elbos)\n",
    "plt.xlabel(\"Iteration\")\n",
    "plt.ylabel(\"ELBO per word\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 301
    },
    "id": "VKgIDiPBue3O",
    "outputId": "f38e196a-fb59-4b0f-f3e3-e188f37601dc"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fd4c70b60d0>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEKCAYAAAAyx7/DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9f7Aty1Ue9q3u2eec+5709AOBLCQRESzsEJIKlJDiiEphMLZsU1YluFSIRIGE8kuqguME20G4UvyKUyVXOThKRQU8g0CKHQthh0jYAkwREwrjYAnbASNCEALBkwVCv5Ce3j1n7+le+WN196zp6ZnpmT37nnPvPatq19579kxPT+/pb1av9a21iJlxK7dyK7dyK/ePmOvuwK3cyq3cyq0sk1vgvpVbuZVbuc/kFrhv5VZu5VbuM7kF7lu5lVu5lftMboH7Vm7lVm7lPpNb4L6VW7mVW7nP5Ba4b+VWbuVWTiRE9GYi+jAR/cuR34mI/mcieh8R/SIRfXFNu9cC3ET0KiL61dDZ119HH27lVm7lVu6B/CCAV038/icBvDS8Hgfw3TWN3nPgJiIL4E2QDn8BgNcS0Rfc637cyq3cyq2cWpj5ZwB8bGKXVwN4K4v83wCeTUQvmGu32aqDC+TlAN7HzO8HACJ6G6Tz7x07oLnzKO+e9dz1Z2SAwnv6zv3t5Fl+Z4A4fuZuOxhQ21NDmwpt1M69jIYt9Hn0Mja6vq2GqSSDoeO6326EVAzMmv9m8/E+5R8IfHL/ux9h5s9ce/yf+KOP8kc/5qr2/YVfvPplAJdq0xPM/MSC070QwG+r70+GbR+aOug6gLvU0VdMHbB71nPxea/7pnVnY4B8eDmAHMM4gFrAtIA9MMyBYffybvYeZu9g9g508DD7Fji0oNYBhxbwHty2gHMdgHsGTOXNSCP70YaLH/YTv1WAjQ/7jF2Tvgbd77A/9X5Xn03hGsfGY0qmjqn9H7TE69Vjk/5bnz6z/r+BbpxLYzo3hlOyZkyA8XtI9aHqv8nPv7Y/U9e+ts0K+fHf+K4PHHP8Rz/m8E9/4nOq9rUv+LVLZn7ZMedbI9cB3FVCRI9DbD5oHnsO/G5lQxq4DWCI4IlhggLtPUCe4BvRutkAsAS2BuQ8mKh/s/sMFH02mbXkN+4Y4BX2pYU3dj/njO33RwM50RBoSn0vbR9cj+ltJ6LuGiMQ6G15G1PXOPWAKfU37u9Gjpt7YOn/Ve07CtZ5m7VjWJJ8XEt9rbkfSg9sMr0+MJSSwdzdZ85Nn6f00J0T3Z0TAvXWwgA8JpSfbeWDAF6svr8obJuU6wDuqo6G5cYTAHDx2S/mVcAdTR9eTCNkAFAAmDA5yBO8Z5AD2BDYErgluc+Iyl4ADdb5ZNGAPKaJlwAv/TajoZbE+66NqB3qQ334UprYObCMTf64rwljlwA6A21jOrCO167BO7zzjCZKvQdPphHHa/M8/G0EgHMZTa42NR5zQD210tGSP7THwF2P0ZpkcETlPrmuD5ybe8b+l1xhmTzt/QPSJWEwDlxnKtlA3gngG4PJ+BUAfp+ZJ80kwPUA97sBvJSIPhcC2F8D4GunDiAW08ZSIQ3cyVQCmFbMJaZVphMvYEEuvDyL9uYRbNv5hB7T7Pw8eMd9wm+sNR8Ninqy5CBemkj5kj6eP54z73cNaJeEKC5fun5b2wds6l5sMuCumNhsMQTl6HOI20h9jtuzB9gY4OQPunTe3gOPu/+SfX+1YqhOmy5JbwU08XDW7a8xueh7aex8Y2aUbuPy88716T6QrTRuIvo7AL4MwPOI6EkA3wZgBwDM/D0A3gXgTwF4H4CnAfynNe3ec+Bm5paIvhHATwCwAN7MzL88dQy1wMVHl//hEbg1eBvHffA+MKhlmIOHPfhg13agg7zQOtChBVoHtC3Ye5l4JW07XWRuTomTRS9TnUwKl34pT2JTAPFe2zOa3xZAHR80CZADaBvTAbY1gDEC0sZ034kAGwA8rmAUeHOOK8nvq0E5OIwRHkp6W/quPqfVVMF23buuoVZPuRYfbNwcO6r/d4Phg3Gp1ByXmTuqZMKsNTBp6YeqXiXp9xqZUm7G/ouB2e6emShGhcFwGz1gmPm1M78zgP9yabvXYuNm5ndBnjRVQh7YPb1Wuwkae9S6fQfc5FjeA3BT62FaeafWC2i3Tm6m6KTSExuos23HyZlPwAFOjyxrF1/zjA27xmQSpaTpxckfJ30E6fAOa8DGiK8gvQsoaO07B+0oPfBmAWj2HB7EHXBroCbXB214Bps6sOiBNakHASAPJmZQBHADeRAn8C5o3ktAtkaTXmuCKZyraNKK/VCrJECZs6rs6901Dx6YzPLwZgYwsVJhXmdLP4H4G8cY6suNdU5qIc84+9S6JzF5pMkPbQ7xgHEe1MokNQcHOBZN27k+aLcOcA7sMm17qVMvXyIP9lt1ifPnz889td+YhIlNTSNg1lgB610DNBbcGPDOgm14bwi+MWBL8I0ANjfBj0DBJKEAXIN138QlD9r4v5lWzCXUMsiFh6xnoJWHK7UecN2DlpzvA8jUGPSAvzMtsQZy9r2HQ29s15hOah78tVJSENI95QEHcNK8lV8kAnn8HN5TL9b2J5fStY75d65RGIC7Be7jxRwYF793tb6BBNwBFLwPIO4FrH032SlO+vTuZOJG+p8C7bnqQWliFCfntNNo/bUuMKnEQyaugzyUOcQCTQDqXQOc7eDPGvDOwJ838GcG7tzAnxHcjuB38R3wDQTMLQADeafwSh1BB9gO4ntwLL6ISN1s0VE3Dwx75ZR5i0TzBrr/MTx0EwDXjM2UP0CPb3bckmpSkw68o23bqo/pnoptdiuFQW+XAGeNrXzMHKP3X2OSuQdyq3FvJNSut33lS2yKtssE3soU0rObZprUlF17K8mdm0uOG5OVoJ1Lb0IqmzXbvmbtw2d5KdCOn23UuPvmkp6mbQCYeE4O5mRZbpOVfTgwgUjZZllriuq/7GnNxYHIxm8OtI+UnkP61FK6p0qrwiV23eijSe117XP836Jz1PsOvNM27jOUSo7UaxIGcLjhTtT7A7hZLaVXHa/sb9p5NWXrLImhnn2TiAbAt8lkXAveYzLBgMj7q69HA2Lqj7aDGsiEDLbrZA6xBA5aNVsEUO+DNgeNm6PWzeoF2eaBwLcncODYy3ki8LMcbwJgrx2yKdDuDZYZ7qt/jqafwvGL7outTBO1A7LF+TIzTeKLxwdUzhM35sbyvBl8ayq5dsnYCUl89ruWBbbKVUC91cRces6K65q8njyQRjsZw4sp2rA7gE1AazrQTt9zU0n8X1zAybFXOBcMdVp2epgsYULMaNFTYzYypkc9vJf0fYuH+zH3YklLLtjZ86AfYKg0ycYbAt48Hsd1U+Q+Ae4ZTXj28BEHlH4fYR5wD/j7TslaM8NgIlcFXByxLK8A6Oq+ewhpc/BDNJfIZ7bUgbXVr6iFBw08fWYB+qjIx7+AAeOo07gRh54kstVH8Cew5aTp6wcJEEwofgIIRqmcleaUI+zbwApw33IFpmUuOnZOxgB4lBIbDosa7bFBRicQrT/cVLlPgHulJO+/vJH63n0WxEhmEm3fzsE6d0pWck4ZqKM5qcf8cSHvFVLLl7VWrjsiZgRKBO062LfZiDMy2bZ3gD+Td3cWzCQ7wO+4M5cYHjgnyYtZhBxALYEP4bQGiDYV8gRiBjmCaZUJhSDccvJDp9ecf2KKQjkC2IMxrxnTcB8MTFJjMpMaYbHMPdTXcP0HD5XSOXy57zcSIQnuxImwjpUHG7ij08Ogf4MYDPjRbDpGwkkkOmhOIItAe2GAAzOD9E1cpHQpk0fmcOybUcJLbe85J9Fvh4iTVj9oM51L2eG30t5OBdpxv6n7YC5h2dzvW4sG5dqIzzHZkvp4QhHn5M3rl5b7A7iJ4M9K6/UZCfcJJc0ZnUZtRKsmIsAT4IQA1Pu7nJMby6FzTC6NiMkm6UC7GotwWyC9FgZL1GxJXxM6P3e+sOqgQJ+k1sAcfHJK2sDRlgCYaNsMirCPJg6EYBz0nJNdegKSDI4OMHtJTyAUwI4WaA7C6zYHnzjdkcopKQsUnVPzr0syp2luAdr5/gXtGwjmqSTqfov3ypJbcCoaN+1TAKlT0VKBG2MSGRMGbjXuTYQgUXkrhUNQB0eTiA80FUcyS1ynUzIC0IegBGKWoIWYdMaQUNFqbj4FkpQ79oD+5JjLmjem+fQvFMkgnQIw+hrjgPlgTB3wBEpdYuWEXC7k4wsp74tPJgwCE8MQgRiS4ZGRAD53TlIwLvZ43C1gDiFFwSEAeBu+txyiXLvcMuSCmUvTO716n73OmQffsaA9d6wx4/dWxYpwqBhMaMmU3YdzSc+mhPNxmvLrVLJ4rlH8rca9gTAkoOKY42MATrRphwmeQqV1AI6ycydtDRjauOdELYt7vN3EYS1EugHlCMqabGWnWuKnB13H3SVpTGzNrQcdLMzewl9Z2DMWO/dZtHlz4HZD2bfRmTqCxAIWlIG3pCfowNrsVW6ZvYc5qNwyMXd6rm3fL3KkSW3AD9emFZ3cDOhHVQIpWKfnOBy51wdBNjpnDQm3nsPnSBkd5K4J78nJrXwnR8lHjjv8VuPeSMgz7NOH9Q1o7rbOaeG5D9YxGIe5HDGpgXGJTTN2A+gmZc8RmU2ODTJKjkcJLl/WszQIcgTeH2SCpWyABsYaGKLB5E0TN/G+9UTVKxD0lLDkOM415vBKjmTv+7+lhFDxmCwC0qgH5wrZVNuekg3AOwrFYJck3cpx9D6rMJMIk8eCrJFo2l0j0bS7Bv6sgb9owGcG7tzCnRm4c4I7k5ffITithXGUHui9h/oRWvjPrT8UABgEd8PrqN8XwA3guOWU5nCXwCBvP4+YBLbT2gqT8tiAjaoVwLEg41UQMBuxI0ezCzA0/5QSFClAoBqn1AjAyoNkBZ0vSn5ufZ48yEZxtQcBV7VmplqZAetNIy199jDTUnWvB7YR1MPZmpCzxnSgfW7Q3jFozwnuHALg54A/QwBwls8W4vuwEP/TNePmralkA/FnBndf9Ix1B0enFyOZSGKCqWgTJdc5t2JWQDjXLbljOtegfZMyGYyB5uwkGwOurZ2TXYdGj68B/mqn6pitfo2za2SMpC8z7S3MaCdfw/fSAzuu2PLw94nI0353MuaLHi9ta9Y5zvPkTyU/yZwMnK48/Jw/fMbGRUkvy2Bs49DK6rj1MFctbEg85gOY+52BtxTSH1Dg85PStFVcwBG4+YH1hwIQjXvPK8gQ91DuC+B2O8JTn72iqwq0O7sp95IXmYOAt92HlK4qHzdSPm4rIB4BPNIMPYPmtJMarz5Q7wSaFVPWOseWnuz7dL8ameMWHwvYqZmRh0CpPuIUqI2NbW+Vxclun37zIfqIuQ/qUdNkFl+wAvPRscyBWuc01+lxoykptw/nJiZTe10F01O6jv6qc5CLPLz3zjKWAItZ8tYfpOJJdVGGmxItqYQB+OtW+WfkvgBuboC7z1v+B5MGbhdfFJxcJCyFA8MeAL8nmL2BvfJga2AasdvG5bE4zqywTIB+YM6YaAeklmIWwIml6xIZzXU88oBZau8tXZPPgVyfd6SW4ZQdPy+rlTnJEhCNgbkvnG9KSmaysG2Q0lXvX5ExsOt/NL+ZMD4e8AYM37GY4rs13TbyXQ5zUtp4aWUzFkRUAuPSb4PrKlxL9lB4UOXWObmFPOrQvOLjiw8T6wjBOYO2NXCthd9bYG9AVwb2imAv5dXcNbCXjOZpYPe0R3OX0TztYC9b2KctcHcPwr5LlsMefGiX2d5vinaxBf1q6loWONbW2G3n6lUu1vIrcpVzCdz7O8yfZ2Buig+obBx6oDxCGR1rt/K/Hey10mE7foJa5/3G591AmAmOH1KNm4jeDOCrAHyYmb8wbHsugB8C8BIAvwngNcw8i8h3dgf8m5/5O4v74EFovcGl2+HKNXj6sMPTV2e43O+wv2zQ3m3g7lrYM8A3hKYhAAySMicgZ2CcOFyoscDBKJueYjI8rDIG0FNjUuGYnRV3D8Z9rl9bnN91y46TwdfcQ3TWF1MJYJWxAPeL+IdY4/5BAP8LgLeqba8H8FPM/AYien34/s1zDX36qQu8+x//4XW98NpUEswjDnh0LxF5ds8we6C58rBXjObSw156mCvRtumyBV3txX53tQcfDmLrXnsTXndppqWAYzpbK1mTiilQKKgAa8GNDXQwqYTjwzs31DmkohPK9ivgDPNxh2AeRnAcR79EKDHnxA9BwR8xqFh0aIXK2TrAuxAo5Kb/r9KYzGmx2il3E+UYwF5zj6455oaOnTgnb7Yx4mS9Y+afIaKXZJtfDal4DABvAfDTqADu8495fN4PP7WmE713zeEuBt6MVb5xHqyDb2JU3pTU3Mg1ZoLJvBaVN35OYattJ4A2WQFm2u2A8zPw+Rn8xU5edxq0dxqhfN0htBfychcQ2tc54M450L4YfMZC92oYZMLn1A8CewJakmjWA8EcCGYvL3sJ2D3k/ZLRXDJ2dxn20qO562AuHezdA+hKXtgfwoOWBNzHgqfmxmTMhqyPPQUIneIhvzVgD5pfoKlae2Nzldw6J/vyfGb+UPj8OwCeP7YjET0O4HEAuNg9Brq7PACnl/N3LJgjBdx0gTdg7upLxsg7Ddpalt7sYzf2llpO3sda7XH0PCqIpkdVo14NyVLebR/zklgIAyOCtvUgy/2uGQYchXo3CEAe2+IUKh/PkdLCEoRtEd4HUX29oaij7xUBuQbAb6IcAdabcMdrwPlUaWtXirvlcZeFmZkGIYO9358A8AQAPPbYi/jwvEdWnScWCZbwbAQKXxby3vquhJlzwunWRYKNBBuQKhZcBPExKUyO4oTYrPLJsJ0iQC1IQxoakcCb1gH7VqIlgZCrxMLuxeHrLgnNWQi02JGkdd3FKDmbAi1S5sDY7R4DiFKtSWH/aLMWi+YdzFv2SrRts29h9q38bzFfSRwPMijV+aSpquPAMCnXkrE7lRzxkBgF4rX33hTg1jhTl7Z5D+Q2cnIov0tEL2DmDxHRCwB8uOYgNkD7yEpCPCvbabSbek6JiST4hsGtBzkHag042kuBQMFy4BaiLXoBAPJhUusJPDGhZgNYgPEbdmwCjBZkGCZKKmqacyYCY7rjuMsImDLwHRyICAZyIzkGyBOMMwFsCX7HsDuk+pOp0nsEbgApOyAUcPsA2iGhFLViJjEtw+45gLiHvRLefcxX0kswlctY6bFinuhgXhvT0NeasE6omVdrx1sFf02BfZGzXRnTcAPEP6yskhF5J4CvA/CG8P6O2gNnKWAFES5smMMcNKsQmZWquoTtZAB4FegQubSAgL4NrGLyIc1rCLZYqsEsvXmnJqOlMkDlGQL1My/WyxzLsW2zB6QJ2qruR9C6GQgVasQERS3D74xw45t+YQW2nEXJAZPFgn0E8BAo5WIK1+CkPIR0rq1KMuW4S+kaxo6MCf9bYdynwub1MKjxobHjpsawJKc2H9RquxuBaXXJu955biY4MvDwatxE9HcgjsjnEdGTAL4NAthvJ6JvgESmvqamLaZgI10oDEnBGrMDcrTRegZa1XjT/5MIEPsqROtDg6R5gw2YlP27pl+TN3XFDTKhEU6ylgwKIJOZPoDpa9Dh2EBHYdsfxOF3aAFrYI0BGiurEEvqPUQBxnBmS8keDWCYHRBI/1cycYVc6uR8MHf50eyO6ZpUJCLpay3JGEdbbZ8Mhy+1MyVHaKOLbM4jvPHi7zUrwiXnOna/axQG4fCwhrwz82tHfvqKxY2RaG/rOhI1OQqmEaTKVuQAQwZoswnHDOJM847hx7G4gheud3lJvpE2NTV59G9TQQxTlcnzftYCT0yU5aL2jmT3J2OEEWKkpBlRAPEYoh0pgKr/TNR3JEczRTBxIYB4ygXuVF5wnYY3gfzIfzJ2ffo3bf5S21OCKV0gOB/bKefl2D1RUZpsFKyrTSNHAPYxQFs1D24ekMst9pBq3FuKOwM++a8dYePWS++QqyQl5w/5Suw+LL/3HvZKlt4xZwnTQbTtFh1lkH0CruJJxyTlRp6I+Y4389L0rrNBI5XA3G9U3oyAJytHHxFJMYpIGwxMkzzHRkpIBHS/Rxl7eOQ5NeI2Dc45UyjuozTvWb59bYqAonN3yXiWmCjUvw+I+v95HvqvjxvrU0nCg6hKW18LpCWzx8Asc6RWf8+EHuoAnM3EN8Ddz1wDOsGk7QB4wLQEaqkrh3Wg4OSiwFRg2EvJXmavBIQMgpklOuaiRJogMJzAUzdkYWIOpaK9JaBxTPhxXi6LwkrDqIRHI1nt+t/V5yglkMgpnLpfOTiHz8MEUIW21oZWL0nnO5e3JpdBzhru7zeX52WhDXzwAOj9uJAlMzi3yvMddxn7rxPVEvPgfg3CuNW4NxGCRDsuluTsoo6lEMphdUmnwraolftAFYxlsCIFLhZSULxuAOXJOpZBDZhfGtfeuHZkv1KioamyZ3OANpW3uSRR2+59V6A9B9zabBE/Ryes7rriWcdWZK5ZdU3qAD1eNZqqNoek4+Zzmsy2B2TgHM0yKpXCdfCeaxyqMzLQ6EtgXfptats1yUPrnNxUPGCfXuHpjsAdzSWxhmGsWRi5weHd7j1MfLWB7hYL0IYXcx+8Z/Mv64kYNNOwQ7cNSNungkdGo/SmnGjUaaYJxAsrhTmTQipgazxS+atot047dWaTBNhZ4A7nppJ4XOn80b6cHh6ZSQQQX0PcL2ricc4prWmRyYT9kLGTj112zJIUCJTjf8weWAveS0C7xmY+58SslfvUpp0Lg24LKWwhxCs1bmBg4465uJPGHbVsHxxggckwyFccQWHlZC1ONgXaRcAe01LGKrUzD51oGmy0Y03J4rwrGmSWHDaZt7ugEQdAZyOO5Z7zL4K9ocQ+6Wn7MW962H/KvpvGC+iP2ZRDMz9+gQzqQqYfKsa1dtyXAnaNRlwjx3LJb4AwgMPDmqtkU2GxSS+VnBcsdu5YeFY5J+PnXtVyBdhRTIi+8wDIYCLwM+w/c3OqHM0c7ZuKxaCzx4WNo0DCuXao28+PicAUtfxa82apOj0zoHJLw/vOGamrt0ABxxiAL3UkTu0/pzWqVQvpcbWxPJnteP89DVyZYozixC+xoc+lbc1XY/kxKIDwnCzlUuuH55K2B9eQrbaybcVjr13oNh/3FkIesFcrD+agUXPfpt2BNicwp5ZTZGWyqUYQ73WIAkf6CFtkT6vsALzrdsEmWuNsq+UW55r4GvHqYcMujQena3bjzqpu4/LzHpuqFOhsugUOd99mDvW/xPHKtwOlcPrx/lWAdY2TL7VxQnvsFIVy0I++aYxzU1leyUdRRAH0q/tcozBuIyc3ESZhliyVqHHDSLKiEEMn048FcMhDovkcugRJidoWjjUE8p3WyGTkgCmONLDcnDC2VK7V5tYyA5YcN6YFjk3k9POGzqi5KugabMbOsdQ8VJIpTveYTIF2afdJzvU1gkvuQCbq/B5ByHEfwLkDcAK6SGWEOQfI8WvNohvKrca9hRigXZFjKi9d5h0FMwmDW8C3km3OWIIlTqqWhYC7iYEfkVESgzNMuAkNMCjbNXkdJbsm9yfkVtXkY9ulPgxWEEdcw8Sy/lo1xShTAD1RyWbSbq1NKzl410jpgTc1fvcKsNcmzjKdgzgReVSfKfsux2Smn/Th+gGTmTbVuInoVQDeCIGW72PmN2S/fw4kzfWzwz6vZ+Z3TbV5XwB3Sry/+EB5MfqrL2YCB+2ArXxOVaYtgVuS+ydq3VEMoRdZt0SmzBE5eG8hk/bfhTbZMZnRskf32xp4tsjUd4wGXjueYyuUKSd1qTBy/nlOxq4t37723r6hBRHWijgntwl5JyIL4E0AvhLAkwDeTUTvZOb3qt3+OwBvZ+bvJqIvAPAuSJWwUbkvgNteAc96/7qbo08JZJW4KDgpD5LvQt4DBTDSAA9tqO4uWee4bVNhhdlCwVoMDZd/Az731EWsWO4PHHlHcJC1xOuY4qeH31LrRTPJPV7m16YIqAngWcqF14m6rC1XD2oMeGeH1YMaUSbkHZKoywBQyky0DUsagO6eNy675w/iwxG/jhefzsF1933ru0pCMa0xc1dU5AbWhzyNbFpz8uUA3sfM7wcAInobpKCMBm4G8Fj4/CwA/2qu0fsCuJunWzznF5cXCx5E4eVFFHTdyPg5K6SQbthSHu5ax41DT3MSrSqLNCtVfq8NbZ4C6TH64gQXfFao8CDScsrk/FulJC3JBCAPTCe140cEgpUFsLVSRehsB744gz+XCkLuooG7MFJFKFYQOldVhHaAP2f4BvA7BjfyCmG9ypsKgAnk5P8xBwK1IcXunkIOc8BcAc0lS6m+K5ZqQnddv1zf/iDNHkKZvqiwrJUtTYAnFnFOVt+jzyOi96jvT4RaAlFeCOC31fcnAbwia+PbAfxDIvrzAB4F8MfmTnpfALesXdrZ3caP5/Q+muciADInEPfq9wWgrX87BUCdQhbzuGeub8Z5OBu0NCVjwSkrueWr5JiHnpbAqkhmuvAuJjx0KXAbhrdIgM07lmAk23cGhppbIY2tgI8BJ9OgdwSyAMW2LYED8LPV7A/V5hYmvPsItKMsiJz8CDO/7MjTvRbADzLz/0hEfwTA/0pEX8g8PnD3B3ATJEXnWlG5LpgybZuipi1gQyz2bonIK6iVtXbAkRu+B04lXnSNjTiKdpABZYZDjf21MtBktj8ZWN+Iaiu1UkoNoMY1BTYBw3ugdvx0zhWOVZgkfsA7hnEE0zJ8SzAHTmAu482iBXoCu07jZtPXuMmLtk0u5OWJFYT2BBsqCcXCFKblfpqHpMCo64JiUlmsfEjd7BSpuWwcOflBAC9W318Utmn5BgCvAgBm/idEdAHgeZgoNHN/ADdwpPYa3rX50qALPNGReHoSRu2jZBZYAnYl7bMENEvzluSBKzP0tB74bCk1gL0EWI9NeSudqDvXgoffJHhPnoI71kRiX8h7skunwC9JgsbBzEyBdUqOQCbW6Aw+BJMNVTKVyP4mRge31I8UDu8xOC3l9NG3Rh5AE/00WzhF7wPZsFjwuwG8lIg+F0Ea8VIAACAASURBVALYXwPga7N9fguS7voHiejfAHAB4PemGr0vgJsNwZ/vVh+vK9kkTYeD6SNWeffBDOK87O+cVByPodbOg6JD0tpuAs8UDx5EC+YBFrWBFxMyHjUZaFoWCihiUiMaHlu6nig1dTPHIgLz3+Zk7LqPyV09I4McL3r8gIRqafzCfdA7FijeD70+pdJvDrRvEzyQs8FJznAHA3clVe2lXidCzc5Q/s2KQzIVTU4XgQTGxOglVZOKQQhad1b+LeXoCY75Qvk3IgJbO796A+6tme8EwgwcltB8J9viloi+EcBPQGbim5n5l4noOwG8h5nfCeAvAvibRPTfQP7Fr+cZDeu+AG4Ygnv0GOBGAG30CgULR9tLscSsWDCsURXgXb/qjeceiOd97U5cAOkx6lcWHl4VChzBGOjZ8eW3buL1GAG6co/nYSmuKeCrYJIM2pgzA2UyG6hTQ48bdWCOOHvTOFSMoU4jkD8M837ovsREUmyEuRHOi9bBHhxMY2ADo8Q3BrwzwiqxgaLaoPscGSWEfsHlqDWH984MEswirmOV9JhUrQ+MEqm7ijbMgSjWBvOin2Y/pWse2WnLIKwTiphKtvOXBE72u7Jt36o+vxfAK5e0ecrSZS8G8FYAz4c8RZ5g5jcS0XMB/BCEp/ibAF7DzJOUESaCPz+mWHAH2uwBsgz4AMyeQORDTiFK56PoC3W+Wyrmdm+9vB4DtSnQ1lzdAX+3AOC56HtLK3rM3VKewxLfQJkBgplnLBHV2DXl16b3130d+X119ORcAqQa0Aa2WbqX/Adz6QPUeKS87k7+f2plchDboCRY2GDzpoOHaUyiA5rouCShBubl33ql31LgWGfDJg8Y51PKYimUHcq/tT5o2twfp5jEC6jj30/9n0uUgmuWhzlysgXwF5n5nxHRMwH8AhH9JICvB/BTzPwGIno9gNcD+OaphvwZ4anPPlvVic6G2N3MMUOgOGkUr/XgRQs5OKn6fhDeNkVua9sGbrILID4yWWtBOwdsnfJ0LGpOSy+0O0446iafssMSoJyuvgPvMZkD7ZFrk00zD6Wp9qOUluQ5vVN/HqkVOSXF1egYX7uW0108kQ8RtsFERyQc6YPcCxTuBYpVg4yBJRUqrqoJJRNXGt/C+fQCoDRm2iyUb4srDEOg2PiYzjQ3zqXfx465QRxxxiI64LXIKWtOfgjAh8LnTxHRr0A4ja+GFBEGJMzzpzED3GyBq2evG0hi0RjERILktNFly0zLMHsDc/CwVx7UGJh9CxiADqFoMHMwLjJSuPvSnAo1oK22jaZBTTc/9RNgJU0afS1Jla7qgXfoUzHXdHEwR0A7B+zStaX9R8wbUZLzDwWgQQagndmibKPGPLhomQOPY2ht8djsnmFgaLuPMvFwK/5SMlGV/pu478x9x9YM2x2TkYdquj+1wz//X2+cbGsqOYXcExs3Eb0EwBcB+HkAzw+gDgC/AzGllI55HMDjALB77Dlwd1aePGnblAopsOuoVj0bIYxkIox5SWJVYemQajPYN0cnulNeeBO0Ig84CK3KUFeAGJBzsZo8YTIl1sscHXAAWtxpoSUbbU/j8v1jpyRmwIuOG3VdICMsB0CCNTRIuHit3AfxGi78VN/iQ4lomNFPp2WdkzgGsUpOwQYufS483LbQFEuZBWOxisG+FSswIvW/yoCk/yY+7HT5sPi/hP+E0r0o86QWYHtKxNTqKMqSNL33WB76mpNE9AwAfw/Af83Mn9R2TmZmGklqHaKPngCAixe9mK+eve5PTelcPaWSZaalUHMSsHvA7wn2imFN/xyGWeyBRqnXEbCdA8d82RqE0wV0HETqaZ0C3KJp9bfJppJZYeTpXzAPFOsv5s5Jn+9bAI7SNaUI0FgoGEjh3LoOZdL2/PB6xpyLczI2qXMWT32L6rptefwAhLSS/W3q8+gYahAdE/1gg/rvZWMZpGvYDr1dOgBn9eAk3T9j0Kv/qBQIWpOpb87MVfreO/56A3aEVXKzuecnBW4i2kFA+28z8/8eNv8uEb2AmT9ERC/ABMk8CkNrUgvOD31/MGBFi000quiZN0he+vhKJq658+Zab7zp8+jBlLvao6exUn9bvM7elB27kdXNP2kqqAXtqcmVXVeviouOWmSG2KZEYxxcjz7fKbMDzj0Uog8g9ik6+BDGpseJL2udi8cw71em/Q4q45QcxDURovq4+ADm7j+Boe5c2j8wqKxUaDvnu9dICaTzFWzvHNdrpnioS5eR3IHfD+BXmPm71E/vBPB1AN4Q3t8x15ZxwPknVg5kMpVEbRuBChU07oNwZ+2BYa4Y9kr4rInTGmtO9m5wQXaC4nP3r111fqhVT9ICgb4dsuoauVsGB9sws5pgqe4kYcBHjs6u2uXw2LUN+jQECiB7IMUVyxaUsHyy1y6755yaAxaJWk1ER7At3wdzsiiydG3q3RF64uw9N+dIHqV/Do/vVb5Zco5rlIfZVPJKAK8D8EtE9C/Ctr8CAey3E9E3APgAgNfMNUQO2H1yeQd0trSudBmkRJlyUNq9ZEqz+8As2TvQwXfZ0lynqaZc3GRk8uZayRiPe8wskjuM9PvY5Bi1wXZaFkVAogDi0d5L0cYdbMAWAAcQ0m2VpJbHnYsauyLAZdsWBdHolcAxMge8yfEbx7IDcADHlS8bPecKzXNJEq450K5Nv5A7Nkm91Pde1Zu4f9yWt73Fw3ylMB5uVsnPYtzk+BVL2iIPnD210nHBAuAdHTBo3C6muQTMwQ/ogAigrSu8J7MDc8/s0JPeMlVrnSiaRQgYOvOiJKdeuJDi9WWOnywJVjH4JvYzbtN9n5J86T5XmTxf1k9xw3uX1PVjFsRLtvilUjI5VR23YOxyyfcvjcsYi2UK0HW7Pcfi0FmJaC6J9x+AVBSBKNne++ObtZ+DdmQUWQM2Rmp4xvfGBKCmLve9oc40aSiQBa4fNG9ZJRuI3TOe+Vvrik4mJ0wMSGDuJdWJppAUTZk07ADcznUZA73rnH6eRyd5XwvPGCYAEJ2gTCpntUHJITUHXMOQ9YKNW2/Xv2XgMQdaRFQAhqjZd2aRIrtBA7w+fkZmTVGyU77DbLulY6sKLufb/cjxS8V1x86vOCo9hrnWGJ3IQBp77mm5Gcd1io5IJBGVREBjhTporXzeWfiz+G7gziS/uDszITWtBBT5Bl0IvyV4C2DA9FohP3bEsQCYCe0tcG8gzqN5akWZ9yg98I6as/oc85W4/veUQD7k6K4BbTldsOemZTz3Qa5Hp4uieL460q5UNHhMuysFjIwFi6wAnEnnWWbTTn0oOTALfejJDKAP+jHcoQ68x677XoP2oFvHtdW776Lo+w9Q7BT9f7n+uEUML2j4bCDh8U3Tf/ilEH0J2XdnJoA3wZ0R/I56+Vd8gxDOL22miNCohV+TPLSmkk2F5EZYLCmXRFzOMRgmaOFePkdGQXxpbXGOawxMsiR6ADMGcsC4lgoMTQ25TFV1n9GyRwFiInFW8aGUriF7MAHlh9OcTABG6sfITx3TpR78Nte0FyTqWiKrC1Cs9U3k+w47VD5P0Q6ODoypvz1hJOFGgPZDbePeWgZRhDViCLHmJEUtG0Ehg+mAHLKRDQnIG+otXRGruueST9ACiPeALge5kn04pxXWRuvNAXa2vRps8t9K16XbHtXsgJ52VyM8wTop8ZnHnJ81UhqzpaA95yCdySRZks2rBU2Bdn6uWufomFMx2KtLD1kmKgK2/AagHN5xz+QWuLcQWgnc6P5/ZoSbyoOYwPASXGBY7hSrnDKK7UHMfXMFKpeyGY+7qH3n9mFgCOBAGbjmlvnANqA9aH94XdLFAoDrPpRs3GMycHLOmD4W+gVySbzt2L9IASwxSbas5lLBac//q0XXlv8XwLjCICfrj/PYyk/vr98Hn+U7+Zk+5y4Kvm6N+yHmcW8pTAR3sS6SqUtzKao3OROyppmQGY0EvFvhUzO5LjNgvJEjBRDcTYaaupMZILKenE7baVX9yTWRaloKwFgF1LWaauG6gDD3suvrDiFUX5geAyqYXDIOfGo/p7bFffO+pw6r/ulAmhIPPkZPRvA2BUZRTUm7MRkZ05LMtl56ALjcJ6Cc4INapyOrIr3iJAK86ezcJDlz2HlxUDpOcy7Nv8joEitlKhrBjC4/mkfSumfB/sTyMPO4txPCOhs3gqYb5h6xtAWPAA5GfgtaeMwXkkwmNVS0ud96zry+hjWgvdVyqGe01nLWuyPAeuqYsesDBmaVORmyVjJtfewYDdo1XGTZIfSZO8065J2G91kEZWBjRPDO5ZQc8qWmkpn/YNK0Vcv8SShrs20IK1SkFMrwHB6McQUWTsXdK4J3LARxvUYS6U+7USGFU8l9AdxsCfvH1ubj7j/1pc5eyBTYegnGaYUWKFxuF9K6OimmYI2kdDUGTKKpUMhVIs2P32aTFWKAISDVmoOmhsJzX1dICZTGTRVL7cIni/jL29BBS+q3qgIUU9GDvfOFd413oc3IOkoZFeElEVUegaraXzKWVWaPtRGFtTnT1+TIifvF1agxXSZBH+ZQoNqag4G/NGhUXnFuRDnyDYWCyJTSTmjnJS99aG0ot6aSDcQ3wOVzVg4kk1qeQeXiBowzEj15iC8POtgucvLgpLp8E3Jy25CTO1TIiZzwoixx/Kh9j3FGjXK6ARSpgvGcS4NH1ji1KgFocP0lICmCTGHbkrG0Qev2GGjgwjxiAKYLaDLAoLwZFoxldVTiylB3FMB4Lv97+Mx6+1i/xiQUZ8A+4W9fz9iQNnkqubVxbyir+PDKNELhSZ4sIOFz9G4zSYJ7AoIWEdrI7c5baElrQplH7LP9XSpZFTnNUSdUmpI1YdojxxzFloi+B+1My4N/lrYX29Dfg9mk220D0FmiQR9R1X642quz/a8lAQxkbKzuA+AG0Pk4bqjcH8DNgFkZf5M7R1KekvgeCykcGCaUcYo1+EhHUoZ6lOy9aFrOjzMMekCvAh2iZOljtQu990stwJf6UbRz5+ddaJsd5Vcrh1YpiKMUwDFmySwxSqJk9ttq8K9k5RT53AAG9MCx32uk5KMt3SNyknGAnjol5RG54cRTCabCZ+ptnxjfTFPv5SixJmxTIe8GwaRCvbD3WIoN1L0DR7JKfu2IY4PcOic3EGIxcayS6ATx3GUIVEmmElD7YOsONfnkFSIp9SvmKokyp6mWkuSXPPlpf61Zr7jomqjKuX2nJOdpAxgk/deXlgcU6WNKMtcn1dbmutvah92acdRSukeiFJWAOdGOR9VAQEPOzXTe9x+CU9q6FgXUBBPSJQfQbowAdWPgd1aAuiF4a8J7AO1g607plak7/3UpvTLNb4H7aCEPNHeXT47I4ZYcJVHrFg0bQfs2B58A2xy8FE09OKE5uVhr0nXh7zG5VJbbelT7K3FpAVTlVd5KFlZtmUzyVEqOlPN9e/twWXObeygtMaVsNo4FP8CcGWnkt0WJsmrZREs0+7HYgPggHRuzGuek1raTdk1drhJL4J3kL+FYsb4xyTEpTknJTSKJpiBArsEbwPUpvQR3yyo5XqhlXHz0sO7YHLxjcqlY5dpzqHodcm8Hkwg532UHdE4ckhNJpmbZJSWK1b0A74kMdjU229l9BtVbCoC8RnuZspX3tMNtnLxynbYbI0I3dol5wugF4VTmKtkq+RSw5LoUVzvlxsn8AnmbORd+UgsP2nZklDQWbC3QGEkutRPAdjvJVeJ3lJJLuTxHSUMpT8mwnOD1yK2NeyNZFQEb52AMvgnmDgrcUrF/cwqHl2M4ZRFMzUyxNZaKGQGd/DegXuusdTTWOiHvF1kK2hUO3jVjtGWCqWPOdZTDdyNh9RCItmz5HskA8lnnLWGlaevfr0sYt6aSTYQbwuVzl3c1gb0yk0j9SQ6UwM7GbaJDsvXAoTOPUOsQcxZzC4hGaQH2oSbBjEZUAuqcunak1jhwqiVfV0ZX85wKJwAoJv8fTVVb6svUtclBWSNHcLoXyCi4TQHsHHWyBOYB5LM6qlV93Cxh1ORJCvdV4by9MmY6ilObSnJNPeaKD4oPc1zBevDBwB4cTLBvp2yBTcjJbZHMJdFMEksHssG18rcBADx9q9wEOWXpsgsAPwPgPJzn7zLztxHR5wJ4G4DPAPALAF7HzJOcEbbA/pkrTAqcs0pYeNwxECcUUzCtAR9CIYVWHDVkCXSgxCgEkKLCABeWn245oBUCIFYX1I3RaoCabLK9t/SPy/7IVdYmAdNf9k8CysLrGrRXG9wRpaLWpnwdmWVLNOepLItTkmno9wSQgfqHYGUgUg+84741SbGiqYxlkrH3UlzbGdDBgPZtclJK8QSTGCW9QgqWboRjMsrDzCq5AvDlzPxUKBr8s0T0YwC+CcDfYOa3EdH3APgGAN891RBDgnDWSA+4nTzZjRM4jjEUngF4gmED9gxqgpZKIaES+c6DrvOW5MvqscjINaA9FYiSh9Mzd8eFNLVxb9Y5NqYSJ02F1NeASu74CrzqXqj1BLVvst21UtPvZNPW2qXvj21sK99X16LM2zu2X7kcW8ZsYnU3GzVZY3rSwiy+IsgsY+8F2I3kB0q0wDaAtg74uQEaNz/MzkmWGftU+LoLLwbw5QC+Nmx/C4BvxxxwG6B9ZOWfmbTtzlRiWpLIyRYwDcEcgnlB8auNXETKWSH2cStaNxnZ15s+lWvMZq0nzZKwbd1GTzKg0CCj8m2kPW0A0HQz5k63vN8Fut+YlFLSRi63ZqLNtTMnSwNOasWWzBwjTsq0YsFwP7Vpkv1xL5hEtSHuY8rCGI97aowjD9wHTQmQmIeoRKj3XttTDtBrlIfWVAIARGQh5pA/COBNAH4dwCeYOebfexLAC0eOfRzA4wCwe+ZzjnNWEBAJ/gROThBNO5IIymBWKN1A+U1NQeueCu8eA+wI1j3wppFziOQRbSm02kIBpjLOedN9jzxd4iGAa+07Nb4QXPL0oLEvvQ7TNKCNycqSZ73zTklG2+MxDTslYsrAu+TMvFc0z5JMjU3JfKVBewpUS99L0gv/B1Bg23MejXwD5aFmlTCzA/DvENGzAfwIgD+84NgnADwBAM94zov4Wb/Rzhwx1lC4gRhdEI6PdMAu+CZxuCMtsA0c7rz2ZF7GTEsEaGNBqgYfrAU38tnvbOC2WvjGTDptUvKd8ODpXRPn16QCioLjNRZApvjuXCiErPjpbTt/XTUyW7VmhfkFQLGYQuFcrINOescvUJ1qcrscG316aunVNqW+CcubZL7qsgT6DrwLVMEkCx6ASQaFRpRPZuyYGyCi7zzEwB2FmT9BRP8IwB8B8GwiaoLW/SIAH5w73l46PONXP77u5L3JyOmVKH95RCRzCrbhCGI5oKnAmyhSPBWyJLRGavHtGvDZDnyxA581cOcW/tzCXRi4c4P2PNThO0N69ztkRVRDdro80JIRGDIGpgUomn4OgNkD9oph91Jo2e4Z9tLDXnnYy1aSaF22oP0BtDegQyuw17bJ6L+IXdITlVu8YrdJmYsunQL1wfEjMpcuYCq0/ZSZFo+R5HCOAB5/yMBbOlDTyenf9YoP6Od8Sdt8+syD/QopBa5ZHlo6IBF9JoBDAO07AL4SwF8D8I8A/FkIs+TrALyjssF1HekFT2jTQ1zKpS29c4khZOGNFOh3aVKkB0EwV3jTURIdwzgCOwY7Cto/idUkdIc9wD5q3+ohwQCYAnAj2O0DaLeieZsQ2m8i7dHHhxUk4EhPqiWXOBUlqmUsYvQYOUVWwty+nzsmUxCLMottqGUP8mSfSojSinBgrtOmuvDigbmv0D99/8Tr8Cz+nx5oR60/zAOE657izF9zvMENeoYU5ZQa9wsAvCXYuQ2AtzPz3yei9wJ4GxH9VQD/HMD3zzXExsA/er6+J0ojiJGTGlDhww0Vq7sHrziRA9hI1CR5sVl4RqyAUszGF7XzVkw7cUISR/qhh9l72CupgM0NSSXsRswkMUexhAOj47cSdWtMDu35DrhjDpYE3DFV7d7L970L6WqDuWR/SOYSblsxlcxMllXgsgWA1zrb9O8l6lrOFOn1M2p9JuySaYXkFRDFcw8jKGM/Ts7nnm3YKDCWvPIgEvNd/ByjHk008Rm5z2wA76hFlMB7oGUjFE/w6bOeWwg877iNwspW/C4RwNUYb/nAXygMgn+IWSW/COCLCtvfD+DlixozgD9bX0ihF0FpFGiTBxyBDAOtgCO1SJQ6GBNsgOPOlGIuD/JgpnTzUut6fPDEWPEM3xiQC7bthuXdhDwOwc4twN0Hgpg4C4zETe9s25IwK+UYbz3MQYF2lvFwTjYBl5IDc/KkE8A7x87p7TvhaAY6oIhsHCBphQD6dMrcmauplCs0xHsS6agBXIE2N1bGKwA3R7AOGf1A1GXyA7qbN3/4cQTpMM88CWnJ+7AqZNnmvMwBojDvfChz5kQDN0AqVqHbvya54Qr3fRI5aQj7x3arj08glzsnfYj2aj2otcERaQeOSTGbIIBc95cOoiY9RAtzAIKNnJwDrAUdDLA/AI0FGQPT2C7tZdNNHrZGnJFx0uTadjo5VB3Nzm4frw0hLS21vrsO130uOltHqvpMJkwa466XfpMGRv6lgsykBBg+NAsBI3MAoI+pscHm2yrzlZTkpPbw6IyEA5jkXmQGHCk+tTKXBECnzGySgHtKAZ1waPLOQJdSECvfxHVsoeh+5MjjGZs6J4noVQDeCBmI72PmNxT2eQ2EGs0A/h9m/tp8Hy33BXD7hnD1nLUadwRuqNSuMdwdopHGLIEuli4zwMEI2EVTByDOSsOj2jdHTR4xTWZw+iV7oVreJ16r6Yo3oLDk16KDbPon7vcB6ANMBKQR0Fmda0NfV20E5VZaZslElUthNVEVZTkV8r4BYK+V6oyDOio2sUw82NPoqiW1l3G6J8vUaZt4NL8UFBIfypal6MkmRk6ieycV9n7sLfKeI48HNlO5g6n4TRAf35MA3k1E72Tm96p9XgrgWwC8kpk/TkSfNdduFXAT0TkzX81tO5WwAQ6Prvw3A2iLSYGSUzCVL2sIpmWwIZiD3MDxoR+dKBQ50SbYuSloMkAfRNVSvTfJ0i6ub3YhtTQMkygdtRTgcvAYA6IawBm5psUyB9pzba8swrsYSJeC9ty5pvp9zHjOnRcZ0EbzVC+aNf2YQJ1V6t1U8b7E9dbi0QG3EbMgWwO0AFnufAWNOMy9VFoGWZLPhmUuBZJArDmZVpjXLBtq3C8H8L5gIgYRvQ3AqwG8V+3z5wC8iZk/LufmD881Wqtx/xMAX1yx7STiLoBPfP66Y6O2DR+ZF0aocwd5JdrcldDmmkuGvdrB3nUwVy3MVQvsW9DVIXnjmQ9y07kcLCMI92/0HnOgF0iiATWfHJVOmpJtdYa2VgXYelv2QCpqYSUbtgp77yUxGml7tA/dyQubKjToORnLUzLz8FsM2jW/57IQ6IvmOwBlpQFArLITrUH6/9HVc4p96z+Ye1p7Znaxmq2iH+Rxn/g9nfs6nZOA99Xnfx4RaR3/iRCDEuWFAH5bfX8SwCuyNj4fAIjoH0PMKd/OzD8+ddJJ4CaiPxBOfIeIvgidpfUxAI9MHbupNAz/WSuUeyaZg47ECdka0EGSR9k9YK4Idk+wl4C9BJpLSrxpJhmcmPpViqA6aWvupioBEgqgpxkXpWAevc/stZZAtxLoasBmI01x8bmj1ID21FjVUPhq+NtT51i5QpiUtW2ajB0DDeJAL193ryKTvidnHNcR7Nf18OaKGOJr9/4IM7/syDM2AF4K4MsgsS0/Q0T/FjN/YuqAKfkTAL4+NPZdavunAPyVY3q6SIiF+bFYxKPNTOJwMLJUJIKqd4eeba3LB9z9vurZX9C+q1K+ln6bPRcPj9eJo7YAG3U9aQURD+0FfHA36anT8Ablso6RpQyOrXjXY2amUwD2sTL3sB279yqzCS6Wm06MzmTD7n4QwIvV91LQ4ZMAfp6ZDwB+g4j+PwiQv3us0UngZua3QLjYX83Mf29VtzcSY5ePJAMBdT0YRjQEG5whbQDyyNwwFL7H5eKU9haX0n74D9fc7FPpUUttjEUA6uRHul9AZ3rJ6Go97vmEk7MoRL1jOAJDqoITjncYgvOa4Jg5uRfAMueMXDqGp5bcFAUMATx3lve2FX7L90ntz1xrr6Rd4R7N5ZqDbnqyXVfeDeClIZ31BwF8Dboke1H+DwCvBfADRPQ8iOnk/VON1tq4/z4RfS2Al+hjmPk7K48/WlaNo17ucPZeOMHiKjvGVPGgR2WKqzy1j/4tnxxjE6PENY489Zrji3btCft3bvNOSZrU9SyZqCWQz23mJZkav9qxi+e/ScBSkikm0pRMgXZF0qqe5MnKBibAwjjfuHGlzZyTzNwS0TcC+AmI/frNzPzLRPSdAN7DzO8Mv/3xEJzoAPxlZv7oVLu1wP0OAL8PyfR3T5gkPfEEd3cFczHaqpw4X6glmD1JWtdg445OSglcQYo8TLUoGUC+Eo7RaMRSZ29C+pnYKpang2vQ4DITcq0nhA4SicdEZ5ImThxhu54tItHbeeF5ajngC0xMw6hG9d957mxisegE0D38PaeUrb3KQdbeMzrgEined6X88DW0zal7JF9xKO1nUMezxg9zU5J2bfiXMvO7ALwr2/at6jND6hR8U22btWj4ImZ+VW2jmwsD2K8DGAqRXOQoJGIi0AEwB5J8HiFBUypppoJYkmMyNUbBY97RmajmPqsxG0SZ0iKPuam19qNs35tE720F1pNtzT/0itdS2NbbL0VOUh/QE6DEMYtfTdplqgTctcuYOa7As19VLLi0PTPBpeegjj7VxTy03BTABsQfVs8quRapBe6fC17OXzppb8bEEZpPrgjAiUqUR8jlEUDbQTTug2TPM3uguYpZ9BjmSvKJxBBxcioNqnNyk4XcHqM1J03QZlS+P/7MBgAAIABJREFUiEGCH6N/pz5dqqj5FJyNOpkPq6yHMU8Eq+/MknclBuUQD2y4q6L5xkArL8hQI74AylNBPWNgo/+HyfOpsQvvsj2Ml9a42aOXtCxr45Sa9+IHrDZVjfG4AcArHnc0dUXz35IgMKAHxoM6qLEfY5JXILp2eTCA+0sBfD0R/QbEVEIQDf/fPlnPlNgD8OgHVwxktFszVA7uzhxiWsAesoRMIRmT2Tupl3doQQd557aV5FHOCQBq52RchpuY0McA1oCappiXW/JvW6Ax8DtJ9iMh710tvmKCqd61sTyU4vXFLIBtiAZtfZdjPM/D7X2XXMp5AJLKtuhwLUm83snxX6ZxCziptLBeR2SuePBpjrDq84D2px98zncPO/3dh8hDDeAOfdBOiaoqAbzKkT2k9dXKgP4H9P0+OgAMEM79mj7OyZzvYCAro6S3lBu2gMqlFrj/5El7MScMmP26kdQBOOSDth3D3UMiJknI5GGcD+AeQNn5LhlTynQWXxnA5SaOqHEn80oIC7YmgTbvwuedgbeSKdAn8EYG3v1rEuAmGZs2hPIbkqIQJDt5EnxhxOtnsO1SuhIVJmqtlCZjPskXRhAWU5x67pUXS8eq4I5eyLUGbKWR59F4XYRgXJEEALamA2xAcnjE+olExVS/q0C7tO+YAzi/9urmR8YTGMYORI08lxPR92+8PAjAzcwfIKIvBfBSZv6BkGv7GaftWif+DHjqc5YflxcKpmTTpp6ZJBYcaK48zBWjubQwVw5mZ2EuJaMa7UOObo6mBlMGbx+yoHkGjAN7CxwOkmjqysh7ysimgN2afma2wCPvmU1ySWYQoEujCZDvHjb6IYSQVEpyqgRNWxdOqNW2x2TJsROgzsYEbTEEiTABoQiA5FEPbBbnBVgRHkoaZDwQK96DqAPqUp+1SSmam/LfMLL831Jyf0pJajnjYwA/5ahM517BuZ9yMpYiUfN9bpKPgIHjE6acVmpzlXwbgJcB+EMAfgBS+PdvAXjl6brWCe8Y7Wfvlx/HgKSZpBA5KVGTZh+A+4pgrwBzBTSXQHtpsbvL4IZgdwY2aG0y/8OEjpo4jZRZCKDIgEyGQ9std5HZwIG+kygeA6QEVHUXOj5JeMyGXciKd09lks/rE9uFPOQJbIAI3iDuqIjRbuu9VBTXQ+bCKmhsGDVAa9BOD8S+GaWq71tJDad+SkoBOGNJwdbkNtfnAbqHJOIKZIQBpe7HG3MvFuSGdGNUak0l/wEkt/Y/AwBm/ldE9MyT9SqXlmA+uiKta3hykgsmkzawSWKlGF3i64rFQXnJsJdOynxdKTt3qM8YbZysnX5A3zlmDFKeBmtlm5V0riknsk5iX8qJrLTt0sM/sa6SLTb0x/Gkxk0FjXtgr98qEvCYMPleLcQA2AG8OQT4JC2aKTEVqJcBj/qFaQsMCCqA96xTN9q59UpFtXm0VlwjU8iSuPQyholfr5NORf+BaDdh/6BgxFs5tldyVmbnH6xG8nJkettN0q7H5AFhleyZmSmExhHRoyfs00B2nwb+wM+ttXEzOuckutSurYdxwaHX+s6Rd3DijHRZQd22TYUHOOWwzjmsBLIW2DWg3Q7Y7aTu5MUZ+KxBe2cHPjNoLyz8mUF7ITUn3TnBnQF+R1JzcodQOBjB1s2IleoBdOYfHx5KbbDdR076XtWbvAqv8DCKibPo8gDaH2RFcLXvrnHL8O21WQajvyCYnQKfbMDY4fgZUGCdmQAwzcgYpMEFOoCO26e0xB5gLRy7U4TK5znInaQeSFsCKhNx90AbiaQczLgJ4E5yTFbK/BzXKIuD8e6x1AL324noeyGFfv8cgP8MwN+sOTDko30PgA8y81eF0M+3AfgMSEDP65h50g5i9h6PfOiysqtKwuDHCh1dwqiglbKA9kA7zQsNlDTTQQVrL9q0rjbSWPD5Dv58B3/RwN1p4C4M2kekUHB7J4D2BeDO5eXPGH7H8DsAjTgTYTksGdR1BRMQtZ0JyO4hZqA9obkEzBWhuSsriWZH8DuPxspKwEDGQkqqWTB7kOsma3Xu51OJdvaGsdYh9im8fowloYC7ag6WAKfWvFQA4BI43ZNxLD0sNfMpH0to9kmWeArIzHWVaFai803FDswEsd1zYTwwzsm/TkRfCeCTEDv3tzLzT1ae4y8A+BVIRkFACgb/DWZ+GxF9D4BvAPDdk+c3BHexruZDXE6LfTq+h8o3noVrHCrFsPGplBKIQG3b/X8ORbt2z3MfltSEVra3QiW0l3uYXQN7tgPvLPjcwu9E6/Y7A78L1d4bClp2KF1mFaskmz+9dLW9mpOK3njwsHsPOvhpiuP+0F9JjFxjGtM5AFprAlgAbL3lf5IC8CyRmnwkE4A9R9m753Uoczt3/jBUVMNeqgKgzwEvdrIfo0A2mAGjoz2WRsupr40cx6lEHykTIbBJPu4fPe7wAY3rBko1GjLzTxLRz8djiOi5zPyxqWOI6EUA/jSA/wHAN5HcHV+OLsnKWyDleiaBGwZoH1n+VO5soAHkgi3YOCmOQC2DmlDeywRTSVgiEgCwHVa+UflJ9ESMdmPyHnxANymUzZsMgQK/29ZU3JZB7L93Jxy+l2yzcbt3Pdusz1kk3lcDy2Sl96kivcfIBH2wJ3mO9CUyE1RS2ucUQTdr2yznSc98MAXtW2d8TO1M2aE108RaiVXYNeBdA5zt4M8a+IumZxZ054T2PCgowSzozgBuEAplB7Og5Q431942RwM3HgyNm4j+cwDfAeASshYlyKX96zOH/k8A/lsA0ZH5GQA+wcxt+P4kJN936ZyPA3gcAM4vni3OtxUSzSPE6OpMhjwk8Mqp532oyef7zALuwG72XNoLHwFbBePAWtlHaSY9WmAs5Jo5J4GhFkIasLVzMlxDNP1Q4p9H+7wHvAtjoIJujJGHTgVolHOUVAL2qfJ6byF50q0K6WVbPJEcpYHn492LNVDmkvwcUzRArW1HUQoDefEZeUOwe5+06Zjf3nlKwXC+IZgGYYUJVRx7/SVvIidwP2wptRr3XwLwhcxcXYaTiL4KwIeZ+ReI6MuWdixUkXgCAB57xgv57BPL6YAJ3LR9O4Jcznl2rs/AKNm4M9u2nrQUgTpGSp7t5PPZDny2A5834LMG/ryBO7dw5wb+nODODNwZgqmkc1Am5+SEqSRVeA+5VswhRoIGx2RwTpqDOCfNXqr60JUDXe1Be3FQ8uEgzslYLCID8CrtejDxJ0Ah7XPNs7NUaGKqIg9z/5pjTo7s+qaAfDUIr3ngFUB6io46OGaKjqqpg+GeodYBl/ugrEjVm16sgpH3FKOgP2cry2OA+1fXHyrCR3bgHkgtcP86gKcXtv1KAH+GiP4UgAuIjfuNEAdnE7TuUlLxgZBn2KeXA3eSBOBQgSpKQ4hhzaUglRmNe6hlE9A0HWhfnMFfnMHfUc7JO8IoaS8I7oLQXkTHJODPo3OSwU1wTBqIqSZ3TjqSwKJQ1UcAm2CuAkf9EmguObwM7F2L5q6F3bUwwa5IRGn5RP6QNPtJgFmjxa1JRHXqnBWF/NKDDII6N0e8Hv1bQUM/2kZ9zKpk5AE6CdgjjJwl18Hed/NKNoT3TnlKSbmU05eAog/humHzQWGVfAsk0dTPQ6V1Zeb/auwAZv6WcByCxv2XmPk/IqIfBvBnIcySr4OkjJ0Ud27w6Zeso43HQsEpr4fnlLaV2uCkbH33rnN6OPkcc3qQpgOOLaujPdkZoDXAoYWJkXsMmNaAWgt7ZWCvCO4uoTkjAe4dpRfHpWOgAxadk44CtbFAB7wS7XuSDnho5friKmMLyZxewIzNdE0iqlPJWMGJyWOWm1cGx0/JmodA1uYoaEeJTkudjArZKq9G5pJJxeyAYZzTQ3JJ3dF7JQ8IcH8vgP8TwC/heOvPNwN4GxH9VQD/HMD3zx3gzgkf//x1rJKO84yUnySaFUwbTAoHDtxnyQporhzMPpgUgjkBrQsmhaChjoA3M0uQi3wRU8yhBV1a0F1x4NhdyFOys8IqaUgSTTXCLGEL8bYrM4m2cWumTHdtwekaE2gdxM5IBw9z6DIdpmCiQxu4265LnBXtlFNgVAIa7fAChmMyBU4zvotpzX8G1GoiT3OQGQNvPSYlzRuoA5yxsVjoD1ik1U+VyIuiwTvKoID1hJTAuvTgS1q9MjlG/VrlOL+VaalFwx0zVyf5zoWZfxrAT4fP74eUrL82icug0fwVJxCpcRltexBKVPScR4dMokOF7Rq0tR+IKNg2OO1HJNuJONXO5GAKkfY608isTAHQVC3DHMDH2js2UVJql8vgvSQPeKk4xRrN+1jRVL0oE2M9OS655Clege66e79l/99aU9WAiXPD1deCPCimkh8LLI8fRd9UMkkH3ErsFeM5v9bO71gSVuaSpKGGjIDOd+aSgxNTSes7rTQzlcSQ92LkZAjAiayRGDnJ5zvg/Az+QhyT/tyivbBwFyZETFIIvhGnpDsLjsmdUKOSuQTo2biFxx0jJ1UY/wG9UP7miiV8/9KH9xbmsgVdHeS1P8hDZQ+5Pj10Y7S/Oc1yDpinyr2NaJ5FoFpTUWhMtHaZFVpO4K217soAnMH4LQ2HnymNN4ov+XkjlTX9zNlv2f5rimCkTo2YEEu73kPlqVqiD+kGSy1wvza8f4vaVkMH3ETM3uORJ5f6RoPkjskYORmBt5TTIwBzpM71iiYED7o0XWBeGAKM5N5G5LVeNPAXO7SPNnDnBodHgnPyDsQ5eQfw54A7Z4mcPGPwjoHGgxoGjFS51yt09gR2JHb0mDzrykio+yUl56S7C9hLi90umGMMwRqSjMdeTDnkHNiaYuTk3MTaHNj1+eaCRdJ+wzD3xbbh2CVdqaVG1LWWxmo1MCnA3szZqXKO8OC3/EGj8qLXyoRmPToON8GeXZIb+DzRUhs5+bmn7siU+HOLT33uuvQoXe7qqGkjFVWQYgpB6z64TvOODsrgvEtUubaVmz4AO5XCi4NzktuOsWFCaD0dHOx5A3vXwJ8btBdGNO4QlOCiY1IHJJhgQsmxTtvtY76SFqEghMpVsu80brN3sJctaN/PVcKtUB5ZLdeX2VBP62Cb1bKXUtlyieW0AAxKlJnIhAi7Rq1bs4kAefifKNBosSwIjlrE356TqXqfYzwRmzsmbwZiPhCmEiL6T0rbmfmt23anLIdHgd99+bobmnxkYASQO4TakweE3B4T2QGTSaEVkNsfwG34nIeIaxrZXhgbfHWVsgOStWhS1ZaR4JvIa81zcqv2w0nkXa8iEmumT22cojqmjHd5lrslDrQoNaDQ278SIMZMIDmFrdBmLZD2ypMRkIoFMyPVn9SFgnPKYMxT0zU4cT1HgOrgmCNAtubBdizPfq7/x5hjTi0PAnAD+BL1+QLAV0BSvN4T4IZluMfW2rgJscJ7zxa8F5uyDSwOHWJLbIKGbmG8VMVh50AhdJ1jcv9YTEGLTuLjVFJ+04Jz7myhYCvlmtzcza/BQ32Poe6s+bSarjWVQKlWc6wFm6UAUJrQC80gq7Tf6JTMi9oqZ2X5OLNu/OKxmWwyhr0GFwDkqcB67j+9afIgADcz/3n9nYieDeFh3xsxDPuMdcDNHmBnwI7ABwNuCX5PsDG8FghaFinN3ATqoE+1ICmvbzgnizLaoZ+iFFhuXug5y9TnQg7k0QRKS6WgoRft/rV5L1IjhT45qMmv2jP6WtKJl8+7ufFT2zcbv/z4ifqSdNRplINzDixLvtBFwJ8fO8JmmZFryUgZz80PiKmkIJ8GcM/s3sYwLu6si5z0nuCcgXOSmYwPoURYyOdMAUTJiR/StwTfcuBSm1RmbBKwx2hw0oFut/hBT45T5+5YAywlB2HNObJrqQqbX2TTdEPgmSJdrOFxR8n6NQDTuXFdOoZjbU6A+SpxEwnCxg8Kfam5FzLT1eTcGNI5rxOwe/IgsEqI6EfR4Y4B8AUA3n6qTuXi9xZ3P7AicjI+OQNtbhdt3Cn4BhIeHiINmytGc9en6jfm7rDoQHTkjRYMnpIaBsaxQL4EUGqkZv+RwJteM5XXlSoH6cRceeWgpkvSxU3cJowZqSoUuOuWehz4QRAToxdR20v36ySQihyLozrm4tD52kuso/yeqB3DfCyjnCAqc5TmucZ0MfVwLHHEOaumkz14b4qi+6Bo3H9dfW4BfICZnzxBf4py/jGPP/hDK+mAAAZFFFKOkmlKoE4yxSOFFMY0oUW5PrqDtqFHraBerdHoFmlHM3zktKrZNZLrZbcDLs7BZzv4O2dwd3b9QhQXUohC8r3EQhQs+V7OGHzG4MYDDYMsg6zvZiOLQ5IdpdQEtCcpa7cXKmWsQ2rvMpqY8+Wulwf7063w4e8egKs96HIPPgROfIxCnbvepeO3RAr0v8Fvvf3nHb2DdqbuF71aycehphLOTaAIPiDA/VsAPsTMlwBARHeI6CXM/Jsn65kSutrDvm82F1VZpqK4MvuvH3H0zXF1y6edMBOMRcSVoueWSql/IxPh2OX35nmjIzODRMPmxoLPm8CDt2gfsWgfMdg/SmgfIbSPAu0doL3DcI948LkH3XHYnbc4Pz/gYtfiomlx3rTYGQdD8X8meCZcuQZ7Z3F3v8Plfof9vsHhssHhysA+bdDeJTRPE9xdwH1aTGdsVQ1QxzBBA+8V3agEntUpdOeklIN7lM2SgXatnyX/be5alharvk7wfoBs3D8M4N9T313Y9iXl3TeWpgE+4znLjyuZMTir5B2pcpGXzdOFYSkH9TmZ0HYms+ctpUrlttoKZyTJj8vOUyMredtdrvLod/BS3OLKwRhCEwoEkzewB4LdE+xdQnNBcE+FVLlnDdqzMxzOGJ+yDOw8yLIEManZKAFMBnDUBTDtCbs9pQyLzSUCPRTY3Q1mtEuP5mknhaQD3z8l6ArsILZ2u3QKxzzI57jwJQ17LRc+0igBoVICHZ0SSEgo95zv+NuniKi8u/7QrgMbtHFCqQXuRteFZOY9EZ2dqE8DYWvQfsYR9YmTLTOYTHq5uIdV0fV3zXcmxXsGpm+u6qAR2bn/W/57jZRWEioZ1eC32uKtlVKlGdYkhYr5VKI5KmRXlI4CpjUwBwt7adDcDTU7dzqACSmIiQ26lAGEvhrFlAKYdBCT3YcgppTTHLB7LwB+5SQB2ZVKQNaGKNuUcIpAsEeP6TbBPBPKQA1QL7Z52/59GA+PvPiu4bCdE67nct2uweNYPKeXWuD+PSL6M8z8TgAgolcDqC6qcKywJRweW/GciMDFCOHdCBpbdD4FZ5QL9u7Wp/wkFOpPCpA4cPQyO0g03QQbYtGkuxde9BIPOUukBCwH8E3AuiApTzPL/8MhLN8cCB4QIPcsfHtvYFrAtYDZQZXBUqlxDbpEW+kkmM4aGep22n14v3KSOfIQ0wAHR2X0gxwpm9q71wK2Om5tf4qpYPM8MN1JTrPiewikFrj/CwB/m4jeFL7/NoDXnaZLQ3HnhE983m75gVHTZvTDw8NElRzWXWFdCRcPRXUPHkYV1qUY+h55gySsg1kZA67ocdd0qa1SatZmZ8sA/SRaXq2Q0sJCWlxyTlglIfiJGwujC9A2Enma3kPhWSmynBWfLXWLw2okvHf1SLNc7W0w2cSH+hi7ZOuxHJOlUac5WMfSY1D9JCoHfdVeR1SS1OcuDa5eoSqtXAc55XLdoe83/HlSG4Dz6wD+XSJ6Rvj+1El7lZ/fAvvH5vcrC2WaFVReD8AehCJo91IfTyY7wRgn2h0yDZBjNCJXRNQtmLw53/VepdScu4Ya2TqqL4I3ma5CirUCmNYAB6ECUkwdQKFS+CBlAAI1sNy/nvks+DtS7c6CCU0zkOB9qPrC6BVi3lK2GNc50M4BO0XuLgBvZkDn1E4Rp10eF7DOspjdczl43wDQfiCck0T0LADfBuDfD9//LwDfycy/f8K+ded3wNknVx7M0TzSt2fGpbFNGnfQuvceJiSaSsvh6MRMHTJim3MYtdElGSxbK5axmYxpb0XTRn7T58a6vJZezMOxtazWvCmBSwKWkO+ll98l8riNEb52LLBsjGQjsEawJJqBiho3AYaFdBFAmymANjwYplctrtcEhbRJxCHlru9OsvWsP2YsgQGI97TsfN8oOQd7LvgsHpOfPwJ4rM8ZDzNAT/setYdfkzwIwA3gzQD+JYDXhO+vA/ADAP7DqYOI6DcBfAoCcS0zv4yIngvghwC8BMBvAngNM398qh3jgDsfXRMBGPqhlsQmVJeWiu8AtV62HXzIGBjLl6mlcUzUJBcVNAYzDXgzjh4amzRTgSqkJhIKq3/v+yHHzIgbeGyCGWyvKR6TPEgDduR2FwBbvlP/PQK2UcUjzIiphAMYhwAceAabYCYhliExAFz4r3zQ3H3oX8wQGLMCskn5YVZ5trZKuFQA+dl7rcYcUkn3K5x8qH3rPg0echtHiq6VBwS4P4+Zv1p9/w4i+heVx/7RrDr86wH8FDO/gYheH75/81QDTMIYWCpdiS8kFgHHCHYjWjdgxPTRGBB7AQEv+8JTf/m4REqVRoA0sXrRa1qjKaWKTW2O3E055zXbbxS0T1WMd67dGZBKpqn+RpAPGrEPa1lPvXeCgISUTiQw5DcygtKpEHKywYbuJHoo0BWRVqCObnvaN/QpvifQXm3iyo5bC+TZfSbdm7jXjMn8LPpeWXjuuYjROe42bgBgIyzUHhBWyV0i+lJm/lkAIKJXYj1b8tUAvix8fgukpNkkcLtz4JOftwY8qcAeQC93daR/dbmrFe1r70D7FrgKy+I0qbkX6jyQyormPGdGmTh2EiDm7NyV+TkmZcnyfQBCE32PDzYisScbArkuXwzFNLjKvh23sTFyKoo2btG4Z5f5EaSzyNqY1oB0xKwutBFs8emeAFCiW05d56wcuXqJknozYJOMJ6Ba42CtMt+lnW8oOm5s4yaiVwF4I2SN/n3M/IaR/b4awN8F8CXM/J6pNpewSt4abN0A8HFIhfY5YQD/kCTy4XuZ+QkAz2fmD4XffwfA82cbaYCrz1xBuWJ5kVMlvlqE0ObglLwCzB5oLgm+YSkXRoAl6kwr1oOdSgXPit9divDKykR1kiXbj5dUBfRTQF0PxJuEGYeUtXMyqKYy2y6pLIAMGAaYwJ4B48UZGTVG8iAfANsr/nc0ZVGstxmdbdm5EuFBac8qGCtSEXvbdGWkkQAtabtiLCtKqMv4HUE1LBn283tK33sePWAf/MNLlIfefhP3ynU7Isdko24RkQXwJgBfCeBJAO8moncy83uz/Z4J4C8A+PmadieBm4h0geC3AohRMJ8G8McA/OJM+1/KzB8kos8C8JNE9P/qH5mZicrPtlDj8nEAsM99NtisGEkGiMOSOS6jmeA9y1LaILwIbFjxfdVLaXAprSsAkBFbaEXypKoqI8fkKy5RBzmflDzoSw/EN8pSuJoGN3atzAA84A0YXoBWMyFCNGXSshP4k3IgTvRJO54nomqjxs3MfdDOrj3VpjxSNqETloopix1pfJ+c3ZQfWyu1Jo8tWE2nkO269HIA7wtF0kFEb4NYHd6b7fffA/hrAP5yTaNzGndMyfeHIOHt74DA2X8M4J/ONc7MHwzvHyaiHwkX8btE9AJm/hARvQDAh0eOfQLAEwBw/pIXMey6kUxOJpnBCMqbRNZ5CdYghxCsEQBccYEjxazvMDMY04Crsq5NMUlKxxcpc+r8tjDxouc0TgqbHePLGeJW1ZickmMobT2uuwccVDEKtXpJWrXpA/VCDvKAf6wjIr0C7PgdGIDZ0qrwJ09jOgbewND/EkXHF9wLuYEFFRaYSp5HRNqs8UTArigvhMS9RHkSwCt65yL6YgAvZuZ/QETHAzczf0do+GcAfDEzfyp8/3YA/2DqWCJ6FIBh5k+Fz38cwHcCeCfEzPKG8P6Omo5ijcYdxQsgy/JZbmQ23LENKJjDIxMBmdaxkcN/VmpBW28fOLUqJ90KTeeegnYUDd4xSIlMYIGEfsV9Q6SrnNssA59SBaF4fmAatI/QGEerwW8pJfAGhtr33P6nkvtb4/4IM79s7WmIyAD4LgBfv+S4Whv38wHoSgZ7zNumnw/gR8IN2QD435j5x+n/b+/cY635yvr+edbsfc77vvzwh+VnKQWqWI3RmnipxbtS0QatkcZaQaO1isFetNbGtGgbbU2b0EvamrSpJaBFqyCiRqJ4qfe2SREEL4DSKIL8KBehXJTfe87ZM+vpH+sya2avue09+5y933e+yc4+e5+ZWTOzZ77zrOfyfUReBbxURJ4NvIU6xbAbRiluT/f1uVmtyxLR0rlGbKiooxW4LAVTqi+ZdhKfGGeFD063EzW2rLXVFq83Cdl2dnYJWQAjj3soADnUGKCxbGvdEWL+ecW/mQigXWUaZjv+3GnuvLXTzjowmHEzgqzzAbkRboWOzkFzw6UvZs5Rdp+64i0D7qYhjCHnYwlW6qxZJW8DnpR8fqL/LuDRwMcDv+yv1T8DvNxLjHQGKMcS9/cDv+bdHQB/DfivfSt4n84nZL5/D65n5WiIwGo9vXWZqnj3pMEa43zY3ny2FhebKcVlRbW6qgffd9yBABsCU74renqDdt2sLULfvig6OruMzSjIXfBjiHpsQLLnuNrb35awnZG8YXsar1Xz9wnPudgSboAw+uQCcu6QOduXHVq6tPXAbbSSi4ZBYnVPyXLqWm5sVskQSd+0FT7f8K8CPlpEnowj7GcBXxmHcUWMD4XPIvLLwLfOklWiqv9SRH4a+Gz/1deq6msn7f6emDqTdDNrF4x06+dv4TTwngvCb5VLTw089Um5wvhAZR9ybpOMDskWebeb3B4KU8m7IUTUOheZPOVsQG0X622o7LpT8+WazuMUDF138R89gcqp2Ie0b5qoW5grHVBVSxH5RuBncZGm71XV14vIdwGvDsJ9UzG656SqvgbX2f01Sba/AAAgAElEQVQkkLp6E6+1+0V6ngLO7600KkAyGSV9AcpRnbuH9Ld3kXVNt5HJKsmqAKb7OpZ8eh5enQ0kdrkxjXQTcC6TZt9YRJ+7aTAve0bynrkP6ShjIYe+DJPccinG1A4cQot7Lsy4C6r6CuAVre++o2PZp47Z5q7Ngk8GPh078w/YLs+buu1MN/Fdbrq9dZA7ttlBelMzHxqYcHydvQ3HImddx41nAmt9y48dax/sSt4HbBi9M2mnyM5o9jxfx07aR7AbfbjnibsBUfZm6wS9fuOOjucNP2OAad3sU/zb0GOVbu/f3gU4Hcd3EAy5kaYIeHUh27XFo9GxvBmUzD785vR1z3R+e687mHjt7chmueuzS1DtxlsoeJtuIe6bw9jMuFFINSpGLd9P4AFbEf+I7lLkxj4NYNCC2YdsRpDLzlb3Pt1bhoSWkn1rKIVZ7yILRGPS703tfuoj77kws+Weve6g+9o7VDrgXGJaB8ZC3HNAmv0Cx0LTaGPfdTjmGj3CIoFTwEHzk1PBpNR10iG01Im+NMBjxo5Wee9vcp8TdsSR//ynQdx7QCLpByuNJlGH+z7j897KKJlcgNJ9sQ4S2pgbKFPSnhunl7R28cvuUua/L9Jqv8axpqltrZz3lvZGfrsDgbWOAOW1+WJndEkdJKNpLhxLDnfAQtw3hzQd0H2uf4+ssP7k7WeCkwH73HC7WD09lZAHndJ3jLcX+o6/Xa4N+eNOmwAMYUzK2nWmqx0gfnBwvZyAqefp2AgbuGc64Nw4VLB2OhmoCmoF69+1ErACvloyVkiFH0qJPSqBWjmuDzvcZIPEtks+7T7EMjNRdKYEjkVXdkhf9s0ugcm4HZoEkvq+c1ovY1IrjwSTtXO6lhmLWLnacT2OIepjyOk+gl3ow2kQN2Cr6TeFqrOItWqRdiXulZC32G0STzY028U0msjGpraN3K/rtroDdhkzW+GXIg1ODnUoH3O+VaHdeSW21Ep83iHLJMnEqP99HHf6JP/12DTUfcvd4SSKblIceyOF4zMR5kajNLL9v8x37f/fJKb6uQdwcCW6Q6Kr+W17samkPWW5OPaIAqsbwMH2oYucj+RhdQiIjnvdFE7G4t4Z7cBk43+Z79r/v27s6t8OOGIrpg+TZiIdaoGYelYh7VzQrL60tj62rO19haauGZNTL3Pa3NBR8DQlFTazbK4gbMitclNYCnCOB1uBybHG2Iw35LWkYR2rMD07WIS7+PlzPT3dF72rbpH2iaKXvMc0VgjLwfVVoB4jgR/RruRwEsStlbD54HqHFT1bW+/TLgWzcYqA5tK1LysucX0nr1z/yWIDZqOYsu78TuWLb/YVyN9Hn6S9vb59ac/hEktHyBP7Xr7oXTCGFPqKblquk1xThe3tdWeaSEiJ9OcuFudYJWrSxIKclq+7K1g5AQdxc1xXcLexTf+eDUIOaKQXHAV5C0tWyTyohNX/24G4wQUjQwDS95yUEsymRdiXrlnw6kIpLq1rGFxaZFO5mzo0i40VkULoHJ9F46YYWfXXFWAbwHbpdWtKn94wVutiQd0Osk3GVAKeus0MUTdIOhD0FnGP8HXnytuDqyX81uICl40u7oUMBit3xlwzr65K0qEHX7oMTM+UiTUBSUVqhqyPwbXUBzmCB0gfToO49/A5SbK+JO/hiZpmkAw+ZQ/xY85Q8HBTGSODJLPLsfUR7txVfbNqIsyAY6jO3Ye0wzrtgq7Meb6xa3YMFh/3TCiU8sEdu12rS/+jEu/+8J3er5zbxK6EYoXrihOMEStgDVIqujJQGtdhvDDeyqqygaluGdNE+6El9F9Pw6WhEtu4bsZYjClmltTs9ZkGzEE6W/5XTdwbNlq4GqolVZsxiDksxZlal03uFtSxnZ0x1GUptAL0WVdx5IalvquOTVcWyojtHQmZ39euEhF5DPACXGseBb4OeCPww8BHAG8GvlxV39u7oUIxj9pMHl/V5W6rFSiFqhJkI+iVb1+WuiYsVCpIBVIZ10B4ZdGNIIXkydNq4wbtvVlzXW/aH+eyuOYQn+pZdnSbsq6ehm30VkqG6ZA0b/yUxMM2Wl1pxribBt1MfR1xprSD22G5XdAtHtUSjopk3nr4prIBYyQDUszVwuwYcOS7eWiL+7uBn1HVLxORM+AO8O3AL6jq80TkucBzgX/cuxUr2Ls77qqVaHXLxvWSlBKKKxecNBv/Kv2rwlnm1gUlnSslE5gc8nHn0Bt8nDGlvt3RPUWibDcLOos2Bo5n7EOqL2jWFxMYeXzttEGNmy961AK3fdtpgPJocN3Byfbh567BInOOjsTKTnHfWtwi8iDwOfjuxap6BVyJyDOAp/rFXgT8MkPEXQmr9+24qyEwmQlOmk0ISkIRg5QWc+VfpYXSIqVtkrcnbcQMTyfHlBPPcdNs3SSZCH4amOxcbyL2tag7t9sf3O3MItmlACecJ1XHPVsZJokVmwYnY9YJ/Q/L60LXb9FF0GOC4rs+jJJq1OFFj5Alj3CXUhzS4n4y8EfA94nIJwC/Dnwz8DhVfbtf5h10dIsXkecAzwFYPeZDd4vWh2BkIO7Kv6y3rCsa5e5i1f1fib7TnfO45yTsMXrc7bZl6fbTB07XeruQzqFIe2Bbg6Td8HUPjG91+3cIATZvjTcCaaGQJFjqad58T9ehyZhrBrYrae87exhR8JQd9xigM2QHHRiHJO4V8MnAN6nqK0Xku3FukQhVVekQ2lbV5wPPBzj/8Cdp+cCugRIXbJTKW9qls7ztyudsF4IarSVcNSzfCk5CfTGGpgpqe3x23lfY9itCJlAUxk78i5JZrw8z9/0bhwHB/S0J1hFoq/plqiPzY40g7VxBTlfBUpId0bAIJwYpd8augcE2coFJqIOTaUC8i0DHPkT2eWgdkdUt3MeuEuBh4GFVfaX//DIccb9TRB6vqm8XkccD7xrcklF4oJy8A6rUwUlfgGM3zretBa5EWgARRDUGJ+3aFewYI6gRd922L2qrW8HJLKqc77MrUNTY+8nHW686jqD3naJ2d++B3ka/cQMZQmhX7aWZJYSvRviTBzRNGrORQN596YHBTZL7bFuZKHui97zuhEw3pR2TtGbHsQYrj+hBksPBiFtV3yEibxWRj1HVNwJPA97gX18DPM+//8TgxkQx692e5mpBK4MaQb3PzSqICtYqUglagC0EKRQtcMsaabroorVt4w2raUFOG0k62qh2ZXNZWDAp1c8tv1sjhd7Uwb4bMpJyz7g2mcq7J3D9vdckacyCMgTf+334X9dna/M3b1spMkfaY89nR8ri2AfAXpK5Q9jVxXWsRDwR97PFDfBNwA/6jJI3AV+Lo8OXisizgbcAXz60kaKwfMijH5k8uKpQqVBVhrIsqMqC6sqga4OuDKYwqHg3ibofq/IBTLMRbGGQldm+wVShqrpvbnD/79ov2CaTPguoKyd5CuHOaUH0HdsYMulaPTnO6EELgUExdbqaDwyr8YkMVUWjihJqF0A7TbCN6O5I3v3vqilJBxdJsLTT2VYrB3w0es7jGB/zzr/omBz36siZ65DwsbFjxkGJW1V/A/iUzL+eNmU7AqyK3SxSY+u+0arOstbCZ1cYP002Lg3MWdpai1AN3Tt9FvfWjmTIfwipld+3D1Oxa4PgHPY5rjaS/VJj6plJGlTNdm5pkXa77L1X3CuMney7928LPj1wjBW5i8Ud9yFDpPucxyH07d8RNoO4CdzPwcnZUF6seM8bHzt5vVjibsWnA8J5yOPe4PK4vU7J6lLde9AquVthLkvMZYlcbJCrDWxKtEyqJtuE0bszHTfZPn0pPfob4eZyaYt5LPCJxzTeL+2rVMUghXH7awysCjR8LtysyRYFFIKuTO3iKiT+rUIdy0j3RXG51/EaUaRS/+7FxUqLVJX/u4KqQsoq/o217r1y79F1Nif5zZZvP1JD/FBl9yfmQlmIewas7sJjf2vXfFKf6qc+1c8qpvLqf6UrtjEb6wWlLGZjkavS5W5vSmRTuhu1DKRdgdokPWzEfrVumsH+fxNTweLauT6JRZ24vV0luMfVOeWY+o6nna6WkvWqgKJA1yv/KtC1wZ659+rMYNcGu3ZZQnYl/oWPVfh3acUroCbrJFXUhOshvF/5a2OjTnTsqnKiY1dl49rQzQYqg5QliKB9LpDrRMP11JFtMyRGtQtyRkHOaXzTee9dUO7f4OScEIXVxY4n0vuu65vUEbZUjsClVEyVkPamcv69MvFhh/fe9L8Mxlg5e5J2//jNIOiWsM8uHd7DeglGNaLtwkjtZ9FY3I4oyd9K26clmtx34d/pYu2fMOfTTD5PzuXf9bzOgRxZpzOZdtu3NDe+T22xD410yW2/f3eDCsl/fwS434OTs0Aq5ewDe4hM+Seo+BJ2R94d0+LSQumnxNZGAs9lhujIKXGniyDbaWSkvsdYZMg7DqW697S+U7ckjA3jbsiQEmmUNhWL31ejitoCtery7Dcu0Gw26mIX3uKuM4NC7EK3gqahKXRafGUq/0Av/QN94yxuKS1m463tcH1UNrpI6o16ATK4MV+xiLjzHkhZBIrCi6Q5FxPGOJdTEE4zBrsybrZQiP+/uHPm02XbriZIyE29RIT17/5ew9r4TmX9Z1t/9kF+DcF+bWVb3SSRL8Q9AywUF7v/iPEmVXV51VajTzNeXKV1utuppR0urnYKWNhujzTlaH9uDruQd99F3lFkskXiEzDK/54tbBk4bmtALKrifo/Kk0/IxU/GNbiYolGQQrC+YYb6AHQk7yTQrFITjiOZhLjLxJ22se7hXno3WmndTCxos6cZRS2VRLHz5XRPwRZpm8Ltz2oVSVq9+4lV4eICq5q4tXDxAVu4dB0tkjhBOO1h9gJ5o6hqxQr8eQtGkfqHnRh/rwFSVWhROPJOg8FzGjATICwW9yxwJ3KHMxl+f9WYNRBu1lj9aHHkHQg6Wg+1myTehK0bte/mHJRHHZJE3dfa6HHpzEEqY7bRa43n0C6/DgSzXrn3M+fnjj7utcGeGeyZNPzbdiXYIvFzm2A5hp0nPshTKYQgMJbzcZuryr+cxK+aMh6fqPqi10yl5TVC/fFgbJ1yqeJmKIUBWziLt6hg44g8yBWLMY40vbWNd53Uwd38A1e8URPF2Cp1D9zE6m5Y28EgsmlAV13sCLYLnW4CgQOOGCdB3Kibuu66LqRZBMkFZamzBvwFFqdytibyeGElF1Vv8U0bHQUrvcUqU6P7hy66mQpjplvxFpemKW7aHkn7bI09W2Fvr7BnBdWtgupcqM4N5blQnUF1Ltg12DNccHLl3wtFjbqUz3Qsb2ljwZQS5RDMRlqdkYzvjGQoLgrMZUVxWSGXjtgE6uvECGxaPt5rhhoTc7Aj+RhBNzTdJ9SfRZKZjDF5/3buemy4NdrGTT0jCeSswJZkwJ4SuQfDcfP2qRA3jlD3WT+1DIJ7JLHEG0QdCXvbtz0Ju6Zg7dQ5Jl/m3ukOCft2CILZ1b9rpD4XSZAsZoWIJL7rxIdd4Ii5aL8UXTn3iGuUkRy/ddakCFgUo+I4pQBbgPh3k+b4e1eCWy8N8rUI7qaKVzoLtfxDJTaj8J9xlrmmglOpjnlOfCtuM3PdpEQdxoUkGNnxfQY3LTy1uEpmgFSW1fvu7reRRhVc4hZp+7ET/+XWNC7dJz+Vz3a/CUij+Y3P299nb5axBDimewuJuyn9vsgLQY1qnjAngm82+mllmzhUG/GKmJ9fBd+qs5yj8YhgVaHw1nVDZIqoFBnlfi1Rkz24TtLgpbsWWvuM266EB4wYV/WZXBspro2Qcil/uetvSHExbqO13yHNNJNiGjXOixAD8stG33UrqwSOiymDK/WIcRLEjVXkkYv9ttGobNPm1K6dutQucw7rBLQ657S/c//I3CxD6Vdp+lZ6owxpbbT2MfrrwT+ETPO4JByX15POaEmPaoI8N9okkoOGgBjeF0kkcFM5i9lULpZJ5chbVRut6QBI1osEXuFjHglha80pMRCXegiM1JopPismEFJvEcchzyOQzc/OGAzu6/Tay1Sgdn1WbT0MdWsmlzaqiE0qgg6NUF93hWw/BG4SR7QrOZwGcRcG+yF39tuGpRGATF0jqdXdsMDT4Ekgu9gN3F9wY/SeU6Ekn/YW+iZGKzi9ccINP9Y6azyUMtb3LlPWQxNLG4FMGrriGq1sDTn4lc/0KHzKWgVaOlIQUYyIOzzvPtWqdqukwUnxP0cg6LTBRrFxAcqYGljV2SZim9dQbxDrus9hijQzqXENAqlgF55QYyyodT1OHrfHV92X6ndMpM28EwAReTquG1gBvEBVn9f6/z8Evh4ocT0Mvk5V39K3zZMgbrs2XDz+gcnrxZPfUdpMzOPOpwRKyNcNEfDKNgm8L/qd6Y9YP8YtoTnw1vWx71Q6dwPkyOXQkft9U7laMyKpKigFMYIpxVm2uNlEIF9bCrZ0Co91Pnf9cvtFJO4QnHTE7cg5tLFz5O1kEFxWieuKJGX9aqQFHhnxAPVvvKWLHog8ZHL4r0X6Dc32b7rLNXSM5ymDubJKRKQA/hPwBTip61eJyMtV9Q3JYq8FPkVVHxGRvwP8a+CZfds9CeIGZzVNXse/S5ABrZR0vuwyAlzkSq31vVPTqPtMFlP7Au9tMKC7W2pDF9t1plnNXUjUBdn+W9vfSWZZmj91dntTcdNpbF1ICTwg5/prk+qWa2TC8Z0IQWcRXHHz4CnA76nqmwBE5CXAM3Dy1m441V9Klv/fwFcNbfQkiFtKy613TJd1jeun1lulTZdJyDENgcmqtqQ0pAWmqYA5Oc8c+srC0wyQLq2IesX+g9PWDdgm77b4f1xtO6g0GqOFsTIC/p0rtnyrxlf2rVfo2Qp77l7VuYnpgOUtcWmBaTpgfPmMEoNTgxTqQpJgcVdu1mMqlwYowdq+AnMprg/pJV6ATLBrYXXhgqcuScVfP+E4g3stXEP7Yobqy8Hrbqi12dY++WWGjIToluo4D0dM7ELCGcN4SERenXx+vu/eFfAE4K3J54eBT+3Z3rOBnx4a9ESIu6J49/t3Wzn9AfpyTFP/7xjN5SHSDqXQ/mbQ8D2tgGaXZboVGOqYpkrLWkrJu6tDS5d+9BiySY8ts6/pcTYWGdq08alqYVUjaGHQoohFN9XtwuVu3zaUt4TyNo68b0F1DvZMXR73mWLPLKwUVooUFjHasKatFbQSV5lZGmQjmCv3Ki4FcwnFpbC6q9iVsE6OaWWd/1s2FVqYuoLWaqwM3Jm403O3q1hVx+/RGUhP1xnTQBi2urV31wsMXLdu5fy6N4nxP9+7VTUnXT0ZIvJVOBnszx1a9iSIG8Slre2CXDZJkDVtaC5nfNa7li+Hm9aYenyR+L1u3Uj+j1x/wDiN7biJuyzsnHWdkkkuoDkG6bG1txULipLtjemWE/Y3ZrcU7QXroUzzFbXUC3zRjaIrhcKT9spiCkVEXZk1gAo2yr66jBt3OfjWdStFNmm+uCYVmBIfLiqyl2cli8y5m4zc7+ELomqN82aAsuH3Fmn0R9UpfoOchd1J6kdI2B47NwnfxtuAJyWfn+i/a44n8vnAPwE+V1UvhzZ6GsQtgp6t99tGyFBoF9r4akoncuPT5lLBG2MdeaeKb6ngfx8h5QgxsVgba86d39t5s3ST9NADqkm0HcfWRmvm0Rivta5UVZ390L6pE41t64tj7Cpxj5ypc4+cK3pmkbMKs7asVhVFYSkK26x4j52RDNWqwF4VqECFEzcJKYLW53Tbjdc/CQJWRdKL1OfhN8J/U85lF/ZVbkz3wXcJit9UmeWn7NueuPHKyD7M6+N+FfDRIvJkHGE/C/jKdAER+STgvwBPV9XhHrwckLhF5GOAH06++kjgO4Dv999/BPBm4MtV9b39G8NpM++KhAjUeAIPaX/i/ZRe5EbAV5JZdv31eotyugjvEBfyyBt/7E002KQ3Z43XK29/134QFEV0T3WO0go0pt2KVEBD4rWAEfXucre9wluZquImN/7/4pyayUvqd9JtJ4eTKxA6BvQ1ckh/gx7jIi51HQqHNyV/24v5tEpUtRSRbwR+FjeV/F5Vfb2IfBfwalV9OfBvgAeAH/H31x+q6pf0bfeQzYLfCHwixJSYtwE/juv0/guq+jwRea7//I97t1UYygdv77wvWSGcmGbWDEpK5VMCQ1cTnwoYrfKW33vSrTtDt5sxiPtXdDzsWjfL3iMP3OCjji1UTnqfMeClU50gkrnrL1YtkMpQbAzFFawu8LolQnUmvrGComtnlVcr5bJQXzLvx/IWVQxOboS1D1CaK98d6dJ3RrpIOiTdtawuKsxFhbkonXZ7CGD7Y5DCOFnamyT1EYS7l3rlnNjVBXpozGhIqeorgFe0vvuO5O/Pn7rN63KVPA34fVV9i4g8A3iq//5FwC8zSNzC5sGz3UYOlW420eK2Wsu7lja+O/3tnjxuU4AGlTPbWS4O9F/4uWDkjDdKp7UQq9Qy44+1MKZkh0xdP1Qfhmo7n1MvIhj/O0qpFFeu+83qEacMWK2dy6RaSyIuFRQCk+KbdGil0f1GSlxnJJ9ZYjbqcrk3TiWwuLSJQmDSBaesmj5lMUjBYazxOcl0bFei68SxzGB0RDD9hnFdxP0s4MX+78ep6tv93+8AHpdbQUSeAzwH4PzWY3afuqTFN8HHbWsCj1a4KlFwqu0Lz0m7zo1Q3TbHdnLoy8GdK0g0lELWu640l01dEf7v4KqIbookxU+Tzyotsk5cKXE4mm4WSdwhmll+a1pyEyQz9DuNvX6uk7SPhYyn4ph98FwDcYvIGfAlwLe1/6eqKpIvLvW5kM8H+JAHnqCrP9nstgPBGEpLlEMxTnCbdHXoSPO5q6r2lY9pYdZbCt1BokkkfycMFUjMQdBDOuIptvRbhq31qMNd+C4tieC/Swt0GtzBNeLcIy6X23rLO8i6xmyTIpBxve9iJVZNRkGp0km6uupLcRMTEW9+GcS67BP1Of+SSqB2naPrxNjy+/Q6aV9ridtnZ/T9zmOu7ZuUCgg4bt6+Fov7C4HXqOo7/ed3isjjVfXtIvJ4YDiKqk7IfmdEnzRR5L3x95giHKhJO252+Ncd1UygfaHOVYE3gkCmRPe3jqVvlrADaWfFtYK1bZJXktmhRdoo2JG2XdeEbVfOv51a1w2Bqkrq1OVwSOoySGwhSKFIIRijSTeY7dlAG4fImtjLbz7XjG4IXft4YPfg3JCjDJrWuA7i/gpqNwnAy4GvAZ7n339i1Fb20Dhuiksln0NmSZu0g381tbA1Ie00QDmAUfKoc1g57W2N3J8pyB5LY7yqPob0Oeu0BPo37lX1FJwl29b/MOIIuKjJ26UE1mQdmij0Vk6GXbJ4dUCnURKaDcRlVDy5i1MbXAmyEqQyaOmlSw/RuagHu/x2W7nyU2dJzY1NHj8id14adQtHBGVKAc6N4KDELSKPwomrfEPy9fOAl4rIs4G3AF8+uB3F5fjuglColfqoU9JOUwOTBqaor4LT5Lu0AnHqE3moCw4MW+Kdx9h9Qx+yA042bSzzgO3whjVRUWfB+JJ3aT0g1VvedZuy0PVGsGc4d8mZI227Bl2pI+0Clwaank5P3FROi9sWYFbApSTLOWKvvGVuvIVf53GbWtYV4oNmp+tjTiS/x1Zq6thZ0iF807lrsauw7AYh6JwFOAfBQYlbVT8IPLb13XtwWSbTsMd9sEXaCdENynKmfmfj3CZTOrzHcbpuhK6baIolktHTjpvJFrLpQXJ0R6f97QJfKBWa+tqV80vbUpACjFcC9OLYgKKVyzLRAv+bhetAvDEAeItaNl6vxEu7iu87mW2oECyy9s0d0hh3uD7mxGCXpS59nDFaJXvtGAedkcyK+5m4Z4Pu8QRM18tU5KkvawdcmlxC8mJtvJijrnXoJtPedhu5m2eEn28fP2bWus6UxGdHmHJ+dw0+DZB2FPQvisTlYpFN5dIBrSJXluLSpQJWZwZd+eDkWuo0QE/Wdcl6koESxvLk67KLQoBSm82CQ8PgMsi6Vq77+1Xp0kdDuihAUdTFW2qnXa9zW7cD11njGmuIULXWO9WMkDmwEPdM2OdExuDkQEpfuFBDfnbSvcORuALTmuBuEXH784w3S/Zh0p6ya3O8fXzenQ+ZvmPos0TT85+KZ5VVPDZTqROFWhnMpbhsk8K5TeqAZSBsasLO7VMIUKon8JKo1x4bKJTWf3Y5/lJVkbRTsbIQpJTCoLrHbzgHWWYt7p7rrEHeE4PKu+DYHwj3u497LmhhqB4432nd6F7NBSd9WqAmmSXR5x1K4aH2ebes7rjdrUHrCzP+t20FRXlMO0zuY9Daj23Z1m3Le2gbW2jtV2PpsTncHb7frZ6H6ezHhIbBBlYGu3apgq5KUqLVbUOWScPaTl0ojsTrvH6J7g/xAlPB6rbB6t4IprRY8dpSQSahrW4XqmzTYPaEcxk31b/Wbi4M3wS4HrrHl517sB470R4AS1bJDNCVcPnQrR1XJhbZSJVUT5Z+6p12v/HWlJaVq4gDZ1nhibCqfADTDmaWTPEz6lxBoT6y7tDldqvtma0QsIe/PirvifFuhlV9LkScyJSXd7Vnzk1i18brcNeukpAK2O6AE1MCSYg75nF70i7V5XJv0spJobg0mML5UwJ518JkCWmX5bYswr7ndBcM/A6NPZtiYd83HXD06PfzJIgbWuQ2AaKuaCLetZW3lgq89eRMMKcrlVjKqe8baDQYHtFMoR2c2kqdS4twUs3ifWJCI5soQM4in3ATmry7KFsmHDTHhxC1uF2xSy7wFysbQz53Utoeim7sKi17Jxbf5FqXRfKuXJk7IhjUnaoQwLQ4VcBgwRfSndLpX5OzSnqkb3dBd6Pn1j9ynZi6rj/p+B3H7vOYoqRjCVx6Q++YcRLELRbWHyynrxjcj0k5e8gKcGJStUaJs7prjZgwYhsAABkTSURBVJKY0x2sqVAt2SDD8RdaNi0rIL2hdq2ezF30mZslH8Cc54bJqiKO1XJR68i7SKzsWHxTV07q2lvb50J57iso10kHnLOkctIX4ESL2w/l/OU1aUslLpOkguLK/W2uXPWke1gYViQzt3Xhr5miVpa01j+kZHrNQfv875mRkv7Go/O4d/Vtd8QOBnEdVb774EieIV04CeLGKsXdHYjbrwv1TZcqA8bqybHl7pPGtdkbMHtT5S7SrrL47Fj9+9Zrye1C2n3SofQQR2MbmX1qixnGLJNEi3vtM0lWQVjKk/a5z+NOW5cZfC538HOH4HQ7HdCnF/rmCYXB+dTRuGxlTez0bivFlAVSWNcBp9zW494L+zxIW7/LoBRvG1vxlpEPkUYZvfTfL23SvmmSzuC+zuM+Foyp/7hOTM7pHsJWALB5wLG1VnZdsx9RXBPib5i5yV2wWZKgo3spECsgFRoq33E2Vq+TG69+17jeUd7Uh8obH9P0+VhcHHPiGH/jBPcFccdMAnA+PKeiTzuG3/Brjyom6SC9jptoi7Czvf/2vAHVNrcbZhxppkv7ogz7O5bAR5LE3sG2rptH65fo9mdV3GwqVNpbvEZ2sm8WF9PwAUpihkn3ttP9CZlJcwnuNzAjCU/+DULc55Dwue4RYxsQXxdCltkR494n7kzncxVpWk3BNxkuWrODn3IOTJVFzV3ouZtizA0xte/kCHKZPE1P12t+EV/iLV5JiDYEEQPpBh0SrLqHdoxIpscgcd2oEtjYlsbU0fhdsLa3gsAzXysjz+9kRMlc0/wMSQZPR1VlL7yPq3FeWimzbXJuX6fHhsXingEC9nyP1mUhMKnE0mkqby2XjqRFBEpvmUUCd6lprhrORw3bWQO5sueOYFMnkXXleG8tlxSl5JYbcVM0rO0cWY/NPx7bFm0CATX0aEJBS/hsoQgpnVXhGxwYqivBXnhZ1yAwtaqrKNMKysZYGqzuIOlKXTUZ0gGvfErgpWukUFzWTRTkqkQuN7ApoSzRskwyS0bkcWdPQDpT2s/PHUWwUtXFoiDK5hr/vfHyuSFPPnwnPnG9nV+fIpWSaNVHBN0fSUXbEi2gRk1ESLFNjIwb70l50+MP4CSIW41Q7UrcmUYKWiliXNEFIohY58oDKHFKdbA9bWyTYe7HTb9LSW4MgY3N5+66iVKrJkPenaQ95SLdJUg7dtOhYTCglLFaUvx2RNexitFeGcxZQXFlfKCyFp4KkqxRj1vYTnML10Qg7koxIZe7guLKImXogFNhSotcVsimQq42LvMoJe2QgbRLIDvu0wxk0dYbKYwjbGNgtXI656vCkXfQOveFTe5BZ+oK1CAV4DN92gWhoztLlb4dYKiPCJlbVYWGFNyKeL3ePGlzPG6bDpwEcbv0vV1vBvcm0dJOiDw0U4D+9J/gSikKoKq7vsMoH3e0snMFOCZj0QxpmmRdJNREHVd3zXeDCJWQuE2KYr9c7rhf06f0vcHZaCGaOqskIGYDiat0LC1iJB5ucGdY7/ZQP1FKu+PEfYjErS3y9qXvnsRNSBv1DTcagmVxY4YgSysU0/O454L/LaK13efqiJaxIqKotbgzWefRS+EqVsWqe4fa3ZQGbZN7KdtZKoyHjyOpxMrYGGkSA7YC8Q/qG3Wj6HG7cTgR4hZViosd5R/b2QDhomqkAyZdcFSbFxc4a8UYlxoY+k6m287udE9a1S6CP32IJFE0j5m2lZ383Rab2tXKmZRqNm7WEafzcbofrDIvDyDutysqL1lQCEUh2MJ4K9tbiyYh7bZUR/CHB+Lxfm2pEo2SYDlW1r2HvqSpq8qLS6Gm7kl6kxZj6qsO5zK4QUKcoKzAeLIWgY1Bci4S8VID0F2Y4y89aZN06jpJJJTjPhoDa6nXHdIRmoL37rm+sgQnZ4FSq7DtgkDWUFsZqcWRXlRQZwqkFp8x3uJQoojQ0HSqT/mvi6R3iegX24Tt9i/RQVGNxS2OzBPXk9Vm4G4qRgaxRgUqw02dFOAATVKoFMGipeMTVUF9PjaFoGVLETB0rmkP1TXV99ebVFoXZ4XZmW1eK3F/ra0Fpm56mu1/jwZpB1h/HajGWYl7GLn9FiMu26bt2+575qbkDS0Cb50L8eJfXus8LmeYh7TnwjHtSwYnQtzeV7bjusB28CR8DsJAoQCnHUhJ25cF8aB4wXXsU0ew0U0JZXsZX2zTK/4z5hjjx3DMTQt7a9kxwlNd2KGbe++tkHMVJRZjmPWk38WgmkitKBisxKHAGnRbiF3BtUDaOdGxMdfFWOyaEhpmJ9Yg3reN92frqoD1ylefFq6QaW2wXrQrxAmizouhjhN4ka4tVxNhphJmK02Xk9loHSeoLLKxThZ3U0HpYgKp37vRyOSmcT8Tt4h8C/D1uHv2t4GvBR4PvATXYOHXga9W1ave7Sjux94VuRsz3IQpYed6Tdoq3pwarP6xWQOJte4+tnzdLT+35sgrYKDIZouY2yTSjtbvGpzMIUeMA66esWmCcc+2SrK7z09jyaFxtmYpzYedJt/NEtwdgx3zqCXomK9WID74eLZGz8/QWyuqWyuqWwXVrYLytlDeMlS3XOVpFSpPz6E6c23fXPs330UoVJ4Gp7SKnwmL6yK0AVMK5grMlaG4AnMJqwuluIDVpbK6sBQXzu1ZXJSYCx+Arhxp62azt1DXPND7l7hF5AnA3wc+TlXvishLgWcBXwT8e1V9iYh8D/Bs4D8PbnDH4GS2+01C4FttyzpIu24YnFhgfQg+xbAOxIBmLQLkezTG3ntaE5BWzRu477mVs6KnEvauwbT0OFMMbG9yh5iQsll/keyD9J+fLgt2QN9lNFHPFYgM52QX0vDVsZI+tqKrBOc6KiRK4boXVGculbI6dy97ro64zxxps7LISsH4TCzxu6egVtDKQOnUGfXKZ6QYJTw+KxvavwmmFGQFJoiESUawKwh13aTVrdxMcHkCDu0qWQG3RWQD3AHeDnwe8JX+/y8C/hlDxG0t5pGL/fakj8C73CEx3zSxZq0dZwkkF168OEOq21DzhF2LL7aaJvQTzCwWTccNNmhRD92YmXOkmf9F1cUhizy7k8kyMSNnwEedeQDOZhnuQVbiGVUBKRKLUQQtXNpftfbiXLeEzW2hui2Ud6C8BdUdpTpX7G0L5xXFecV6XbFel6yLinVhKRL9nMoaKmu4KgvKyrC5WlFeFZSXBXJpsGcGexHcLYKKK5pyVrrBVAY2Bt149xfUM9vkvrwx3K8Wt6q+TUT+LfCHwF3g53CukfepalCMehh4Qm59EXkO8ByAW8WjXc7sfjtUvydEnPVTpnrbLQs76z8eGrrNITlS7/g/kCfysQUwXRfggQNoSoZMp6CqmucmnINQ4NQ+Jyl5jy1oSs9BXyVfu9iqi7RvMCiphrrVXkN5sra4dWWwoVuQL1aq1t7K9qQttypWtzacn5fcWpfcXm84X5WsTcXKWIx/fJZq2FQFl9WKi3LF3cJyWazYGMXKyk2Q1Qlz2RJs6fXSoxSvwRSSf74Gd+VNukru16wSEflQ4BnAk4H3AT8CPH3s+qr6fOD5AA+eP27/XzBVLOvSGAlBLWvA2Dpf21oX6NG6GIQiL9zUnaOcuUJHWtpd2xzTwd2lqmUs72LiNkfszywYOk9taxvy+fFx+Z59bVvXgbz992H7Wzn76T6k11Eja2c69j6vxtR+7lj85eM7PrXRVBIbRpgNFBvQS/GJSQYLbBRsVVCeFVxsVpytSgqjFMZifFQytbg3ZcHmakW1Mai3uIsLQ3GJ83n7StTQx1NCkU6lNOonYh5/neN9I1DQ+ziP+/OBP1DVPwIQkR8DPhN4jIisvNX9ROBto7a2i/WWLVRJCNwXFWgga7dSk7ypbwJN1oORGXRTiGZiJ5y4xECArf1dX9n7wW+WuQp2+si6/f/e/UkILt1W4haom0lnipYy+3ZjhBNK1lv7JNblpJvKNT8uCid16zJFvPvCCrIBuymozg3VmaVarWBlMSuLFI5IxRO39SmYWglaOreHbITiUlzXoAuhuMQFJu+q+/sKiit1mSVBBz+9F8W5VcRnCN1oBeVNp3QO4JDE/YfAp4nIHZyr5GnAq4FfAr4Ml1nyNcBPjNraPicyvQByrpL4v+5skc6ONy0yH4U+Qft0u/taYPtYDe1c5bkwQU2x3hX1u5HsR3oOG1rQpv7/LujRip5UaZq7jg45U+lDI1uJmrCF5svgAotG/WcfkDQ4S1u0TpSy3tZRv4643yP29wzpg61319RiIE2TGeMGu+Kmxx/AIX3crxSRlwGvwSmAvBbn+vgp4CUi8i/8dy8c3JhV9Ko3Y3Bw/XrHetLm0gyMDr/mPtPgrTUP0Zh1ZIbHjd8YbAdth6DQT/JzEOMJnb82pCjq4CQQ8tt1ZbDrgvJWweaOca9H+cDkHReYLO9Y9JbF3C45v1Vy62zDnfMrbq1Kbq82nJnS+bi9xV1aw5VdcVmuuFuuubtZc/dqzcXFms3FivJugbkrVI8Iq3PBPoJv4mx856HCuUw2FbJxv2no66pe9+XGgpM3HRgdgYNmlajqdwLf2fr6TcBTJm4IrjYz7dRAmtxQIHJHKzZ7m++rvT1q4OO9ADupr++8tING+wRA+9CbWXKc51RJsmwgpgNqKLI5E5f+dw7VbfcqH6VUj7LwQMnZ7Q0P3LngMbcvePDsLh96dpdHry94oLjkjrlibUoK/6td2DWXuuKR6owPlLd4/+Y277u8zfsubvMnF+fcXZ85V4uEql5xLeIqqDaClAZz6drSxV8wqCsGsa6bdFcc4YM5xUlUTuqtM8qP/XM7ruzeRLUWlYraJIlGSVXVrcuCclllXS53/Ns2JSkhP00PSEuPoS7l9n8PFeSMtiB7im4aD6au2cPUCspDPXAyJJzNLHH/aHzXXG5kauAYt0h25tVqlrxPBWoXdj3H6rKitPQqhjhxAxeYXFNcrljfNWw+6NICyzsSi3KqW2e8//wO7z1XdK1wZhHv4zbGYkx6fgRbGax1Pm69cvnc5sJgroSzC3EFOHe9j/tCWT9iWd1VVo+4ApzikQ1yuUGufOFNct4klsbfBG44j3wEToK4rz7E8PBfftTk9RqtrCxec9m9FxunAuc0l52MZ0N3+TJoLpfuBrjaoOECi9VdFk3qgLcIOpBziPQb9y7G1CJKhRew8iI/DWGltHw7h1w1qCdoSUu20wdOkB2tKqf61kp53NIxGXOe53BR5HLb22p3qbZ0TkQp/p08+MZWTrZy+8Wfk2y6aFW5oHZnuuge2vF7QHyqpF5tnDFyeemusw8YTFFQGMNZTtq1cLnesbfnygcHo1CXATGNyt5mKzdv+CguG6eytdZLIukaq5RjibsziNQbR9HQ8PfGTWaV3M/BydmgBWwe3OFEpsRduZLcQN7VxpXn2kuwV2ALQxFKez0MuHQla31k3qJqEQ1KgaaZDgbbpB3kNdukHQjbv0d/pDFsqbP1EHezClRd9KeyriotuBU82UlVuQlIVTUlaqHe5319e/t0bmlZ19IgYdM8j12ELfVLh2YuIegZCDsIHalCJTWBi/o+Gj7jCK8CmKaLxqrYI3GjVJW7dqoqnkvtePiZ9rmE7QdgH7ZmfFo/zDSNG6lLs0uXz2Q9AftdR3PgSN1hASdB3AB2NZ24Y8/AykXQLYpRcWJ41j0QtHBBE5N0StGVT3PaCFL4Czq8isLJixaFm071TN/rnVegAhU0CFqVnkwT0nabaBKPtLcdkJJO+NxRFdqwHHNl/IF8wuygy783IUd9MIWvD7kUyUA0KcLxRcnSpoUtQ+mVXYQD2+etJ+OoURfQd/66MHe2SVcKaof7JbaKS/PUw/kcs3/trK0UPv22dn0UiQ/euZ3S9W7Myk6g4O7TmSAiTwe+GzcVe4GqPq/1/3Pg+4G/CLwHeKaqvrlvm6dB3OIIdjICcYu3oAr/g6g4nirw00TXMcUUoIV6vYXQDcQprWlhkGC5BJKlO7sg3jBWXe5UcJn5m2fLGhzo9dcmwjBuY/Sc373H/93AlBZjW4VDKUH2HMcUghrIzY5t4IL7p2v7fcfVoT+iOYswFw8w0p/jPQP2ckPlKknbed4+rbLR5zM9L11aNH1ozzradQ+x+CvJJjkmqDKXxS0iBfCfgC/AVYq/SkRerqpvSBZ7NvBeVf0oEXkW8K+AZ/Zt94bnI8cHndv6GRwwQ7aZp7166y+8GgjR+K5tH7m/rgHvT50NXa6LY3Fp9GAW0h76bkEW6jNbhl4j8BTg91T1TV4F9SW4ivIUz8DpNgG8DHiaDPz4J2FxX7314Xe/5e9/61smrvYQ8O5D7M8y1smOs4x1OuPsM9aH7zPoH/Pen/15fdlDIxe/JSKvTj4/38t1BDwBeGvy+WHgU1vbiMuoaiki78fJXnce+0kQt6p+2NR1ROTVqvoph9ifZazTHGcZ63TGue6xUqjqaE2lm8LiKlmwYMGCw+FtwJOSzzl9priMiKyAB3FByk4sxL1gwYIFh8OrgI8WkSeLyBmumczLW8u8HKfbBE7H6Rd1IGJ7Eq6SHfH84UWWsY5krHvxmO7Vse7FYzoYvM/6G4GfxaUDfq+qvl5Evgt4taq+HKfX9AMi8nvA/8ORey/k6FJxFixYsGBBLxZXyYIFCxacGBbiXrBgwYITwz1J3CLydBF5o4j8nog894DjfK+IvEtEXneoMfw4TxKRXxKRN4jI60Xkmw841i0R+TUR+U0/1j8/1FjJmIWIvFZEfvLA47xZRH5bRH6jlXs79ziPEZGXicjvisjviMinH2icj/HHEl4fEJF/cIix/Hjf4q+J14nIi0Xk1gHH+mY/zusPeUwni3ZF3qm/cAGA3wc+EjgDfhP4uAON9TnAJwOvO/AxPR74ZP/3o4H/c8BjEuAB//caeCXwaQc+vn8I/BDwkwce583AQ4ccw4/zIuDr/d9nwGOuYcwCeAfw4Qfa/hOAPwBu+88vBf7Wgcb6eOB1wB1cAsXPAx916HN4Sq970eIeU2I6C1T1V3FR4INCVd+uqq/xf/8x8Du4G+kQY6mq/on/uPavg0WwReSJwF8FXnCoMa4TIvIg7oH+QgBVvVLV913D0E8Dfl9Vp1YYT8EKuO1zje8A//dA43ws8EpVfURdb9pfAb70QGOdJO5F4s6VmB6E5G4CIvIRwCfhLOFDjVGIyG8A7wL+u6oebCzgPwD/iGa/70NBgZ8TkV8XkeccaIwnA38EfJ93/7xARKaLyU/Hs4AXH2rjqvo24N/iesm+HXi/qv7cgYZ7HfDZIvJY37P2i2gWsdz3uBeJ+56FiDwA/CjwD1T1A4caR1UrVf1EXJXXU0Tk4w8xjoh8MfAuVf31Q2w/g89S1U8GvhD4eyLyOQcYY4Vzn/1nVf0k4IPAweIsAL6w40uAHzngGB+Km7k+GfizwKNE5KsOMZaq/g5OIe/ngJ8BfoNaX3MB9yZxjykxPTmIyBpH2j+oqj92HWP6Kf4vAYfSbvhM4EtE5M04l9bnich/O9BYwWpEVd8F/DhTe5+Ow8PAw8ks5WU4Ij8kvhB4jaq+84BjfD7wB6r6R6q6AX4M+IxDDaaqL1TVv6iqnwO8FxfXWeBxLxL3mBLTk4KXeHwh8Duq+u8OPNaHichj/N+3cTrCv3uIsVT121T1iar6Ebjf6RdV9SBWnIg8SkQeHf4G/gpuSj4rVPUdwFtF5GP8V08D3tCzyhz4Cg7oJvH4Q+DTROSOvx6fhou1HAQi8qf9+5/D+bd/6FBjnSLuuZJ37SgxPcRYIvJi4KnAQyLyMPCdqvrCAwz1mcBXA7/tfc8A366qrzjAWI8HXuQF4A3wUlU9aJreNeFxwI97meMV8EOq+jMHGuubgB/0hsObgK890DjhIfQFwDccagwAVX2liLwMeA1QAq/lsCXpPyoijwU2wN+7pgDvyWApeV+wYMGCE8O96CpZsGDBgnsaC3EvWLBgwYlhIe4FCxYsODEsxL1gwYIFJ4aFuBcsWLDgxLAQ94KjgVfV+7t7rP+KkIO+YMG9jCUdcMHRwOuw/KSqHqTEfsGCewWLxb3gmPA84M97bel/41+v8xrazwQQkaeKyK+KyE95zfXvERHj//dmEXnI//03ReS3vK74D9zgMS1YMDvuucrJBSeN5wIfr6qfKCJ/HfjbwCcADwGvEpFf9cs9Bfg44C04EaIvxWmCACAifwH4p8BnqOq7ReRPXeMxLFhwcCwW94JjxWcBL/ZKhe/EaTL/Jf+/X/N66xVOo+OzWut+HvAjqvpuAFU9uGb6ggXXiYW4F5wi2oGZJVCz4L7CQtwLjgl/jGvNBvA/gGf6pg4fhusq82v+f0/x6o8GeCbwP1vb+UXgb3iRIhZXyYJ7DQtxLzgaqOp7gP/lmy9/OvBbuJ6hvwj8Iy+ZCk669z/iZEX/AKetnW7n9cC/BH5FRH4TOKgU7oIF140lHXDBSUFEngp8q6p+8U3vy4IFN4XF4l6wYMGCE8NicS9YsGDBiWGxuBcsWLDgxLAQ94IFCxacGBbiXrBgwYITw0LcCxYsWHBiWIh7wYIFC04M/x/p8g/0YartzAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "topic_usage = q_pi.mean\n",
    "plt.imshow(topic_usage, aspect=\"auto\", vmin=0, vmax=1)\n",
    "plt.xticks(torch.arange(topic_usage.shape[1]))\n",
    "plt.xlabel(\"topic\")\n",
    "plt.ylabel(\"document\")\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "GICzRqWWue3O",
    "outputId": "589371f5-6eaa-439f-88e9-9a7f23b5f149"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "topic  0 usage :  tensor([0.3847, 0.0888, 0.0834, 0.0931, 0.0613, 0.0777, 0.0585, 0.0508, 0.0508,\n",
      "        0.0508])\n",
      "\t people : tensor(0.0149)\n",
      "\t federal : tensor(0.0077)\n",
      "\t national : tensor(0.0070)\n",
      "\t union : tensor(0.0062)\n",
      "\t representatives : tensor(0.0058)\n",
      "\n",
      "topic  1 usage :  tensor([0.4278, 0.0656, 0.0646, 0.0685, 0.0777, 0.1051, 0.0534, 0.0457, 0.0457,\n",
      "        0.0457])\n",
      "\t executive : tensor(0.0141)\n",
      "\t legislative : tensor(0.0097)\n",
      "\t senate : tensor(0.0074)\n",
      "\t body : tensor(0.0074)\n",
      "\t president : tensor(0.0068)\n",
      "\n",
      "topic  2 usage :  tensor([0.2959, 0.1060, 0.1093, 0.0754, 0.0586, 0.1825, 0.0463, 0.0420, 0.0420,\n",
      "        0.0420])\n",
      "\t cases : tensor(0.0093)\n",
      "\t courts : tensor(0.0088)\n",
      "\t jurisdiction : tensor(0.0065)\n",
      "\t trial : tensor(0.0060)\n",
      "\t court : tensor(0.0059)\n",
      "\n",
      "topic  3 usage :  tensor([0.1641, 0.0610, 0.0538, 0.0989, 0.0484, 0.4167, 0.0448, 0.0374, 0.0374,\n",
      "        0.0374])\n",
      "\t peace : tensor(0.0058)\n",
      "\t military : tensor(0.0055)\n",
      "\t time : tensor(0.0049)\n",
      "\t militia : tensor(0.0047)\n",
      "\t war : tensor(0.0045)\n",
      "\n",
      "topic  4 usage :  tensor([0.1576, 0.0593, 0.0542, 0.0715, 0.0478, 0.4437, 0.0471, 0.0396, 0.0396,\n",
      "        0.0396])\n",
      "\t powers : tensor(0.0105)\n",
      "\t constitution : tensor(0.0075)\n",
      "\t congress : tensor(0.0068)\n",
      "\t authority : tensor(0.0065)\n",
      "\t confederation : tensor(0.0062)\n",
      "\n",
      "topic  5 usage :  tensor([0.1520, 0.0574, 0.0452, 0.0517, 0.0385, 0.5027, 0.0535, 0.0330, 0.0330,\n",
      "        0.0330])\n",
      "\t nations : tensor(0.0073)\n",
      "\t us : tensor(0.0062)\n",
      "\t commerce : tensor(0.0044)\n",
      "\t war : tensor(0.0041)\n",
      "\t foreign : tensor(0.0034)\n",
      "\n",
      "topic  6 usage :  tensor([0.1873, 0.0498, 0.0726, 0.0520, 0.0590, 0.4437, 0.0382, 0.0325, 0.0325,\n",
      "        0.0325])\n",
      "\t confederacy : tensor(0.0045)\n",
      "\t members : tensor(0.0044)\n",
      "\t empire : tensor(0.0034)\n",
      "\t cities : tensor(0.0027)\n",
      "\t among : tensor(0.0021)\n",
      "\n",
      "topic  7 usage :  tensor([0.1334, 0.0510, 0.0439, 0.5379, 0.0393, 0.0635, 0.0378, 0.0310, 0.0310,\n",
      "        0.0310])\n",
      "\t kind : tensor(0.0003)\n",
      "\t reasons : tensor(0.0003)\n",
      "\t equally : tensor(0.0003)\n",
      "\t union : tensor(0.0003)\n",
      "\t connected : tensor(0.0003)\n",
      "\n",
      "topic  8 usage :  tensor([0.3348, 0.0920, 0.0871, 0.0862, 0.0592, 0.0634, 0.1420, 0.0451, 0.0451,\n",
      "        0.0451])\n",
      "\t kind : tensor(0.0003)\n",
      "\t reasons : tensor(0.0003)\n",
      "\t equally : tensor(0.0003)\n",
      "\t union : tensor(0.0003)\n",
      "\t connected : tensor(0.0003)\n",
      "\n",
      "topic  9 usage :  tensor([0.7449, 0.0341, 0.0359, 0.0272, 0.0263, 0.0312, 0.0260, 0.0248, 0.0248,\n",
      "        0.0248])\n",
      "\t kind : tensor(0.0003)\n",
      "\t reasons : tensor(0.0003)\n",
      "\t equally : tensor(0.0003)\n",
      "\t union : tensor(0.0003)\n",
      "\t connected : tensor(0.0003)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Analyze the topics\n",
    "usage = q_pi.mean\n",
    "topics = q_theta.mean\n",
    "\n",
    "# Sort the topics by usage\n",
    "topic_perm = torch.argsort(usage.sum(0), descending=True)\n",
    "usage = usage[:, topic_perm]\n",
    "topics = topics[topic_perm]\n",
    "\n",
    "keywords = []\n",
    "for k, topic in enumerate(topics):\n",
    "    # if torch.allclose(usage[k], usage.min()):\n",
    "    #     continue\n",
    "    print(\"topic \", k, \"usage : \", usage[k])\n",
    "    \n",
    "    inds = torch.argsort(topic, descending=True)\n",
    "    keywords.append(vocab[inds[0]])\n",
    "    for i, ind in enumerate(inds[:5]):\n",
    "        print(\"\\t\", vocab[ind], \":\", topic[ind])\n",
    "    print(\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "a4BTtABuue3P"
   },
   "source": [
    "## Exploring topics by author usage\n",
    "\n",
    "Using the model, plot the total topic usage for each author"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "id": "MC1TVd-hue3P"
   },
   "outputs": [],
   "source": [
    "# load authorship and convert \n",
    "authors_df = pd.read_csv('authorship.csv', index_col = 0)\n",
    "to_tensor = lambda df: torch.tensor(df.to_numpy())\n",
    "author_names = [\"HAMILTON\", \"JAY\", \"MADISON\", \"DISPUTED\"]\n",
    "authors = torch.column_stack([to_tensor(authors_df == name) \n",
    "                              for name in author_names])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "id": "z2mzQVVD43c5"
   },
   "outputs": [],
   "source": [
    "author_usage = torch.row_stack([usage[authors[:, i]].mean(dim=0) for i in range(4)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 347
    },
    "id": "wG4j2LI05JNW",
    "outputId": "270bbce3-2432-4561-94c7-f4a01131d83a"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAFKCAYAAAAXGTnxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5gsZZn38e/tAUQFA3JEBfSAYjgioBIVEyZYBXQVAQOwoqwrpjXia2Z1FeMaMK0BBV0FET0uGFFEUMIhCiqCwIKYiIqKku73j6rh9PSZma6u6emumvp+rquvmap+nq5n5vymzt3VT1VFZiJJkiR1ze0mPQBJkiRpEiyEJUmS1EkWwpIkSeokC2FJkiR1koWwJEmSOslCWJIkSZ1kISxJkqROshCuICIujYgn9q3bLyJO6lt3QkRcGxG371t/WERkROzet/6D5fr9ZnrN/u1GxF96HrdGxA09y88t2yyPiBUR8aeIuD4ifhgRj+x5jWXlNo/rG8sREfG2ur8jdUuZ9Rf2LO9X5mrPcvkOEXFhROzb1+8tEXFyRLjvUSXlfvDGiFi/b/1ZZeaW9ax7W7luu762+0XELT37y0si4nMR8YCeNlP7xjXK5Y0i4uiIuKrcn543ta8un799RLwrIi4r98UXRsRrIyJ62pwQEX+PiI171j0xIi4d3W9ITVNm9oby/+DrIuInEfHiqf1eWRO8o6f9/hHxy7L9HyLiuIhYt6ftjWVur4mI70XEg2Z6nXLdbTkeVDOUfy839bW7rue1MiL+Wq6/OiKOn9rHLyb+ZzQi5c740UACu83Q5FfAPj3t1wCeDfy66jYyc52pB3AZsGvPui9GxP2Ak4GfAZsA9waOAb4bETv0vdx2vQWyNE/7AtdQZjwzbwD2B94XERsARMSDgVcD+2fmrZMaqFrpEmDvqYWIeChwx94GZQG6Dz057PPTct95F+CJwA3AGRGx+SzbPBy4HLgvcHfg+cAfep4/CngC8E/AuuXzBwAf6nudvwJvHvgTarHZNTPXpcjPu4HXA5/pbxQRjwX+E9i7bP9g4Ct9zd5TZncj4I/AYVUGMKhmKJt9pbddZt6172W2LPs/sNzuRyPirVW23xYWwqOzD3AKRVD2neH5bwI7RsTdyuWdgXOB349wDG+j2Nm/MTOvyczrM/PDFDv0Q/ravgd45wi3rY6KiPsCj6UoAp4SEfcEyMwTgSMpdpwBfBp4V2b+cmKDVVsdzvTidl/gC31tHg3cC3g5sFdErDXTC2XmLZn568x8CfAjiv3mTLYBDsvMv2bmzZl5VmZ+CyAingA8GXhmZp5XPn8K8DzgwIi4f8/rfBjYuzxQoY7JzD9l5gpgT2DfGd54bUPx//ZZZftrMvPzmXn9DK/1N+BLwGxv3hZMZl6VmYcD/wa8ISLuPu4xLBQL4dHZB/hi+XjK1FGwHn8HvgHs1dO+f0c+X0+iOErR70jgURFxh551HwMeEH1TPqQa9gFWZubRwC+A5/Y893qKHf3RwO2B945/eFoETgHuHBEPjoglFPvRI/ra7EtxwOHIcnnXCq/7NYoCerZtHhoRe0XEffqeexJwamZe3rsyM08FfkNxpHjKFcB/A2+vMB4tUpl5GkU2+vN2KkXN8PaIeFT0Ta3sFRHrUOxfz1q4kQ70DWANYNsJjmGkLISr+3o51+e6cg7Nx6aeiIgdKT7+ODIzz6CY7vCcGV7jC8A+EXFXiiNoXx/xGNcHfjfD+t9R/Fuv17PuBoojwu+Yob00jH0ojlJQfr3tyF1m/gU4EHgGxZSIW8Y/PC0SU0eFn0TxhuuKqSci4o7AHsCXMvMm4KvMPD2i32+Zvl/stQfwY4ppDZdExNkRsU353Gz7Wsr16/etexewa0Q8pMKYtHitlrfM/DHwz8DDgWOBqyPiA+UbvimvKeuOi4B1gP1GOKZn99Y2EfHDuRqXf19X9f8cbWYhXN3TM/OuUw/gJT3P7Qt8NzOvKpe/xAzTIzLzJGAp8Ebgf8t5lKN0FcVHg/3uBdwKXNu3/tPABhFR5ciJtJqIeBTFfPQvl6u+BDw0IrbqaXZ+31epjsMpDjDsx+qfpj0DuBmYOgn4i8AuEbF0wGtuSDGneDWZeW1mHpSZDwE2AM6mOCASzL6vpVx/Ve+KzLwS+Chw8IDxaHGbMW+Z+a3M3JWiuNydIuMv7GnyvrL2uGdm7paZU+cW3Qys2fdya1L8f1/1PIwje2ubzHz8XI0jYk2KOmbGv5s2shCep3K6wbOBx0bE7yPi98C/A1tGxJYzdDmC4oShUU+LAPg+xVGMfs+mmIP0t96VmXkjxcd1/wHEDP2kQfalyM7ZZfZP7VkvjUxm/h/FSXP/RDGlode+FEfKLitzeBRFQTDTJ3O9nkFx1HfQtq8C3kdxAvJ6FPva7XqvBgEQxdUqNgZ+MMPLvBd4PPCIQdvT4lN+mrAhcNJsbTLz1sw8niI/VeYBXwYs61u3CXD5Ap6QvDtFAX7aAr3+2FkIz9/TgVuA5cBW5ePBFDvXmT6a+zDFR3snVnz9NSNi7Z7HGnO0fTvwyIh4Z0SsFxHrRsTLynG8fpY+hwNrU5y8Jw1jbYo3WQewKvtbAS8DnjMgq1Id+wM7ZeZfe9ZtSDEn92msyuCWFCcIr7YPjoglEbFJRHwEeByzzN2NiEMiYvPyMlTrUpwkdFFmXp2Z3weOB46OiIeUr7k9xYGOj2fmhf2vl5nXAe8HXlf3h1f7RMSdI+JpFJ+aHZGZP+t7fvdyHvrdorAtxdTJUyq8/NHAUyPiyWUG7w28iVWf0I3y51gvisu0HgockplXj3obk2IhPH/7Ap/LzMsy8/dTD4qPwZ7bXwyUZ4Qen5lZ8fWPo5jPO/V422wNy53vjhT/CVxKMVftmcBTMvPkWfrcAryFRTTfR2NzF4pMfqEv+5+lOJnCN1caqfJqDyv7Vj8aODszv9uXww8DW/Scpb9DRPwF+DNwAnBnYJv+wqTHHSkuP3kdcDHFeSC9l8Z8JvBD4NvAXyiK4M9QvBGczYcoDpxo8ftmRFxPcQm+NwIfAP5lhnbXAi8CLqTI5hHAe3subzarzDyf4rKC76KYqvBTik/lhjkxc8++6wj/JSLu0fP8OeXfzUUU0zX+PTPfMsTrN15Ur8ckqRARZwIHZ+aoT/iUJGlsPCIsaSjlme8PZrKX8JEkad4shCVVFhGHAN8FXl+evCRJUms5NUKSJEmd5BFhSZIkdZKFsCRJkjppYtf5XH/99XPZsmWT2rwWiTPOOOOqzBx096iRMruaL3OrtjK7aqvZsjuxQnjZsmWsXNl/OUhpOBEx9hO2zK7my9yqrcyu2mq27Do1QpIkSZ1kISxJkqROshCWJElSJ1kIS5IkqZMshCVJktRJFsKSJEnqJAthSZIkdZKFsCRJkjppYjfUkCRJmo9lBx1bue2l737qAo5EbeURYUmSJHWShbAkSZI6yUJYkiRJnWQhLEmSpE6yEJYkSVInWQhLkiSpkyyEJUmS1EkWwpIkSeokC2FJkiR1koWwJEmSOslCWJIkSZ1kISxJkqROshCWJElSJ1kIS5IkqZMshCVJktRJFsKSJEnqJAthSZIkdZKFsCRJkjrJQliSJEmdZCEsSZKkTrIQliRJUidZCEuSJKmTLIQlSZLUSRbCkiRJ6qRKhXBE7BwRF0TERRFx0BztnhkRGRFbj26IkiRJ0ugNLIQjYglwKLALsBzYOyKWz9BuXeAVwKmjHqQkSZI0alWOCG8LXJSZF2fmjcCXgd1naPcfwCHA30c4PkmSJGlBVCmENwQu71n+TbnuNhHxcGDjzDx2rheKiAMiYmVErLzyyiuHHqw0KWZXbWRu1VZmV+My75PlIuJ2wAeAVw9qm5mfysytM3PrpUuXznfT0tiYXbWRuVVbmV2NS5VC+Apg457ljcp1U9YFNgdOiIhLge2BFZ4wJ0mSpCarUgifDmwWEZtExFrAXsCKqScz80+ZuX5mLsvMZcApwG6ZuXJBRixJkiSNwMBCODNvBl4KfAf4BXBkZp4fEQdHxG4LPUBJkiRpIaxRpVFmHgcc17fuLbO0fdz8hyVJkiQtLO8sJ0mSpE6yEJYkSVInWQhLkiSpkyyEJUmS1EkWwpIkSeokC2FJkiR1koWwJEmSOqnSdYQlSYvXsoOOrdz20nc/dQFHIknj5RFhSZIkdZKFsCRJkjrJQliSJEmdZCEsSZKkTrIQliRJUidZCEuSJKmTLIQlSZLUSRbCkiRJ6iQLYUmSJHWShbAkSZI6yUJYkiRJnWQhLEmSpE6yEJYkSVInWQhLkiSpkyyEJUmS1EkWwpIkSeokC2FJkiR1koWwJEmSOslCWJIkSZ1kISxJkqROshCWJElSJ1kIS5IkqZMshCVJktRJFsKSJEnqJAthSZIkdZKFsCRJkjrJQliSJEmdVKkQjoidI+KCiLgoIg6a4fkXR8TPIuLsiDgpIpaPfqiSJEnS6AwshCNiCXAosAuwHNh7hkL3S5n50MzcCngP8IGRj1SSJEkaoSpHhLcFLsrMizPzRuDLwO69DTLzzz2LdwJydEOUJEmSRm+NCm02BC7vWf4NsF1/o4g4EHgVsBaw00hGJ0mSJC2QkZ0sl5mHZub9gNcDb5qpTUQcEBErI2LllVdeOapNSwvO7KqNzK3ayuxqXKoUwlcAG/csb1Sum82XgafP9ERmfiozt87MrZcuXVp9lNKEmV21kblVW5ldjUuVQvh0YLOI2CQi1gL2Alb0NoiIzXoWnwpcOLohSpIkSaM3cI5wZt4cES8FvgMsAT6bmedHxMHAysxcAbw0Ip4I3ARcC+y7kIOWJEmS5qvKyXJk5nHAcX3r3tLz/StGPC5JkiRpQXlnOUmSJHVSpSPCkjSsZQcdW7ntpe9+6gKORJKkmXlEWJIkSZ1kISxJkqROshCWJElSJ1kIS5IkqZMshCVJktRJFsKSJEnqJAthSZIkdZKFsCRJkjrJQliSJEmd1Mg7y1W9I5V3o5IkSVJdjSyEpTbylsKSJLWLUyMkSZLUSRbCkiRJ6iQLYUmSJHWShbAkSZI6yUJYkiRJnWQhLEmSpE6yEJYkSVInWQhLkiSpkyyEJUmS1EkWwpIkSeokC2FJkiR1koWwJEmSOslCWJIkSZ1kISxJkqROshCWJElSJ1kIS5IkqZMshCVJktRJFsKSJEnqJAthSZIkdZKFsCRJkjrJQliSJEmdZCEsSZKkTrIQliRJUidVKoQjYueIuCAiLoqIg2Z4/lUR8fOIODcijo+I+45+qJIkSdLoDCyEI2IJcCiwC7Ac2Dsilvc1OwvYOjO3AL4KvGfUA5UkSZJGqcoR4W2BizLz4sy8EfgysHtvg8z8YWb+rVw8BdhotMOUJEmSRqtKIbwhcHnP8m/KdbPZH/jWfAYlSZIkLbSRniwXEc8DtgbeO8vzB0TEyohYeeWVV45y09KCMrtqI3OrtjK7Gpc1KrS5Ati4Z3mjct00EfFE4I3AYzPzHzO9UGZ+CvgUwNZbb51Dj1aaELOrNup6bpcddGzltpe++6kLOBINq+vZ1fhUOSJ8OrBZRGwSEWsBewErehtExMOATwK7ZeYfRz9MSZIkabQGFsKZeTPwUuA7wC+AIzPz/Ig4OCJ2K5u9F1gHOCoizo6IFbO8nCRJktQIVaZGkJnHAcf1rXtLz/dPHPG4JEmSpAXlneUkSZLUSRbCkiRJ6iQLYUmSJHWShbAkSZI6yUJYkiRJnVTpqhFt4IXTJUmSNAyPCEuSJKmTLIQlSZLUSRbCkiRJ6iQLYUmSJHWShbAkSZI6yUJYkiRJnWQhLEmSpE6yEJYkSVInWQhLkiSpkyyEJUmS1EkWwpIkSeokC2FJkiR1koWwJEmSOslCWJIkSZ1kISxJkqROshCWJElSJ60x6QFIktpp2UHHVmp36bufusAjkaR6PCIsSZKkTrIQliRJUidZCEuSJKmTLIQlSZLUSZ0/Wa7qyR7gCR+SJEmLiUeEJUmS1EkWwpIkSeokC2FJkiR1koWwJEmSOqnzJ8tJkiR1VdcvGuARYUmSJHWSR4Rr6Pq7J0mSpMXAI8KSJEnqpEpHhCNiZ+BDwBLg05n57r7nHwP8F7AFsFdmfnXUA5UkSZovP9VVr4FHhCNiCXAosAuwHNg7Ipb3NbsM2A/40qgHKEmSJC2EKkeEtwUuysyLASLiy8DuwM+nGmTmpeVzty7AGCVJkqSRqzJHeEPg8p7l35TrhhYRB0TEyohYeeWVV9Z5CWkizK7ayNyqrcyuxmWsJ8tl5qcyc+vM3Hrp0qXj3LQ0L2ZXbWRu1VZmV+NSpRC+Ati4Z3mjcp0kSZLUWlUK4dOBzSJik4hYC9gLWLGww5IkSZIW1sCT5TLz5oh4KfAdisunfTYzz4+Ig4GVmbkiIrYBjgHuBuwaEW/PzIcs6Mg7wsu8SJIkLYxK1xHOzOOA4/rWvaXn+9MppkxIkiRJreCd5SRJktRJFsKSJEnqJAthSZIkdZKFsCRJkjqp0slykhYHr0IiSdIqFsKSpMar+ibON3CShmEhPEYejdNMzIUkSZPhHGFJkiR1kkeEpRbyKLIkSfPnEWFJkiR1koWwJEmSOslCWJIkSZ1kISxJkqRO8mQ5SXPyxDxJcl+4WFkIS5LGxmJCUpNYCEuSJDWIbxjHxznCkiRJ6iSPCEuSJC0Aj+w2n4WwpMbwPw1J0jhZCC9SFhSSJHWH/+/X4xxhSZIkdZJHhCW1Xp0jIXWPnnjUZfGr+m/sv6/Ufh4RliRJUid5RFi3GecRMo+qSVpM3KepSxZTveARYUmSJHWShbAkSZI6yakRkiRNSJ0T8xb6Y2mnbqhLPCIsSZKkTrIQliRJUidZCEuSJKmTLIQlSZLUSRbCkiRJ6iQLYUmSJHWShbAkSZI6yUJYkiRJnWQhLEmSpE6qVAhHxM4RcUFEXBQRB83w/O0j4ivl86dGxLJRD1SSJEkapYGFcEQsAQ4FdgGWA3tHxPK+ZvsD12bm/YEPAoeMeqCSJEnSKFU5IrwtcFFmXpyZNwJfBnbva7M78Pny+68CT4iIGN0wJUmSpNGKzJy7QcSzgJ0z84Xl8vOB7TLzpT1tzivb/KZc/nXZ5qq+1zoAOKBcfCBwwRBjXR+4amCr0fRzW+3Z1n0zc2mN7QxlAtldjP9WbmuVxZrbuv3cVnu2ZXbdVlu3NXN2M3POB/As4NM9y88HPtrX5jxgo57lXwPrD3rtYR7AynH1c1vt2lbTH03//bmtyW2ryY82/P7c1mS21fTHYv29u63R9Ot/VJkacQWwcc/yRuW6GdtExBrAXYCrK7y2JEmSNBFVCuHTgc0iYpOIWAvYC1jR12YFsG/5/bOAH2RZrkuSJElNtMagBpl5c0S8FPgOsAT4bGaeHxEHUxyWXgF8Bjg8Ii4CrqEolkftU2Ps57bata2ma/rvz21NbltN1obfn9uazLaabrH+3t3WaPpNM/BkOUmSJGkx8s5ykiRJ6iQLYUmSJHWShbAkSZI6qRWFcETccdJjaKuI2DUiWvHvvBiZ3XrM7WSZ2/rM7mSZ3fq6mt1GnywXEY8EPg2sk5n3iYgtgX/NzJfM0WcD4D+Be2fmLhGxHNghMz8zYFt7AN/OzOsj4k3Aw4F3ZOaZc/QJ4LnAppl5cETcB7hnZp42YFt3Am7IzFsj4gHAg4BvZeZNo+xT9jsC2AE4muKKH7+co+2r5nqtzPzAqPrV3VZbjCu7dXJb9mt0dofJbdne7I5A0/e5Zb+hszuP/WcjszvOvLeF2Z1/v65mt+mV/weBp1DenCMzzwEeM6DPYRSXert3ufwr4JUVtvXmMtQ7Ak+kuCTcxwf0+RhFaPYul68HDq2wrROBtSNiQ+C7FHfrO2wB+pCZzwMeRnG3v8Mi4qcRcUBErDtD83XLx9bAvwEblo8XU/yhz6ZOv7rbaotxZbdObqHh2R0yt2B2R6Xp+1yol91a+886/caU3XHmvS3M7jz7dTa7o7g93UI9gFPLr2f1rDtnQJ/TZ+hzdoVtnVV+fRfwnP7XmKXPmcOOr6/fy4DXVRljnT59/e9O8Qd+KfAt4ELgZbO0PRFYt2d5XeDECtsYul/dbTX9Ma7s1sltX54and1hcjufPJnd8ea2t/04sls3g03P7jjz3vSH2TW7dR9NPyJ8eflxR0bEmhHxGuAXA/r8NSLuDiRARGwP/KnCtq6IiE8CewLHRcTtGXzE/KaIWNKzraXArRW2FRGxA8VHJMeW65YsQB8iYreIOAY4AVgT2DYzdwG2BF49S7cNgBt7lm8s1w1Sp1/dbTXduLJbJ7fQ8OzWzC2Y3flq+j4X6mW31v6zTr8xZ3eceW86szvPfl3N7sA7y03Yi4EPURz+voLi8P6BA/q8muKWz/eLiJOBpRS3fR7k2cDOwPsy87qIuBfw2gF9PgwcA9wjIt5ZbudNFbb1CuANwDFZ3KVvU+CHC9AH4JnABzPzxN6Vmfm3iNh/lj5fAE4r/yAAng58vsK26vSru62mG1d26+QWmp/dOrkFsztfTd/nQr3s1t1/Nj2748x705nd+ffrZHYbfbJcXRGxBvBAIIALcsCk8p5+OwKbZebnyndq62TmJQP6PAh4Qrmt4zNzzneg5bvBQzLzNVXGNF8RsQnwu8z8e7l8B2CDzLx0QL9HADuWiydm5lkVtzd0v7rbWozqZLdObst+jc1u3dyWbc3umI1zn1v2q5zdtuxzy7Zjy6DZLZjdadvrZHYbWQhHxEcoPzqYSWa+fI6+5wJfBr6Smb8eYptvpZiE/cDMfEBE3Bs4KjMfNUef7YHzM/P6cvnOwIMz89QB2zolM7evOrayzwOA1wDL6DmSn5k7Dei3EnhkZt5YLq8FnJyZ2wzot4Ti44bebV1WYZxD96u7rSYad3br5Lbs1+js1s1t2dbsDqkt+9yy39DZrZPbsl/jszvOvDeR2Z21n9mtqKlTI1bOo++uFPN2joyIW4GvAEdW+CU9g+JsyTMBMvO3MfuZklM+zvSzFf8yw7qZnBURK4CjgL9OrczMr83R5yjgExSXh7llwOv3WmMq1OU2bizDPauIeBnwVuAP5baCYkezxaj71d1Wg407u3VyC83P7tC5BbM7D23Z50K97NbJLTQ8u+PMe4OZ3ZmZ3YoaWQhn5rT5HuW7ppx6FzWg7/8B7wHeExGbAW8GDmHw5PIbMzMjYmoS+50qDDWy55B6Ftfrq/I7XZviEi+978wSmCvYN2dmlcuz9LsyInbLzBUAEbE7cNWAPq+geKd79ZDbqtOv7rYaaQLZrZNbaH526+QWzG4tLdrnQr3s1sktND+748x7I5ndWZndihpZCE+JiK2Bz1FcHiMi4jrgBZl5xoB+96V4l7cnxTuG11XY3JFRnAV614h4EfAC4L8H9Lk4Il7OqusHvgS4eNCGMvNfKoyn3zcj4iUUE+3/0fNa1wzo92LgixHxUYp3TpcD+wzocznVzpwdRb+622q0MWa3Tm6h+dmtk1swu/PSgn0u1MhuzdxC87M7zrw3mtldjdmtqJFzhKeU83cOzMwfl8s7Ah/LzLkOt59KcdmPoyjm/Qz8z72n75OAJ1ME4DuZ+b0B7e9BcRboThTv0I4HXpmZfxzQ7wEUfwgbZObmEbEFsFtmvmOOPjNNws/M3HSubfX0X6fs8JcKbT9DcfLAsUz/A5rzzi11+tXdVtONM7vD5rbs04rsDpPbsr3ZnYem73PLPkNnt05uy36Nzu448950Zne1fma3okYfEQZumQo1QGaeFBE3D+izT2ZeMOyGyo82fpCZ34uIBwIPjIg1c44zSMvw7jXstijeOb4W+GT5OudGxJeAWYOdmZvU2A4AEfFU4CEUd5mZer2D5+hyWflYq3xUVadf3W013ViyWye35Xgan90auQWzO1+N3ueWY6qT3aFzW7ZrenbHmfemM7vTt2V2K2r6EeH/Au4A/A/Fu6c9gb8DRwBkz329I+J5mXlEzHI/6grvMs4AHg3cDTiJYgL+jZn53Dn6rA3sTxmanm29YMC2Ts/MbSLirMx8WLnu7Mzcao4+a1LcVnDqlpEnAJ8c9IcXEZ8A7gg8nmLS/LOA0zJzrmsCap7Gld06uS37NTq75nYymr7PLfsNnd06uS3bmN2WMLur9TO7FTX9iPCW5de39q1/GEXQeyePT01Wn+nMzSrVfuSqi0Z/PDPfExFnD+hzOPBLivubH0xxB5dBd7IBuCoi7jc1roh4FvC7AX0+TvERzsfK5eeX6144oN8jM3OLiDg3M98eEe+nuGXiaiLivzLzlRHxTWb4nWXmbqPqV3dbLTKu7NbJLTQ/u5VzW47D7I5G0/e5UC+7dXILDc3uOPPeImZ3OrNbUaML4cx8/BBtP1l++/3MPLn3uYiY89p+q5rddjvCqXc/g84cvX9m7hERu2fm58uPK348oA8Ud7v5FPCgiLgCuKTc7ly2ycwte5Z/EBHnVNjW38uvf4viWodXA/eape3h5dcfAaf3PTfXpWHq9Ku7rVYYY3br5Baan91hcgtmdyRasM+Fetmtk1tobnbHmfdWMLurMbtVZWZjH8BdgA9QfOywEng/cJcBfc6ssm6GNo+huNXi68vlTYEPD+hzWvn1RGBzYH3g4iF+vjsB61ZseyZwv57lTSv+XG8G7kpx68TfU7yTPLjCtjbvWd4bOLXiGIfqV3dbTX+MK7t1clu2a3R26+R2Pnkyu+PNbdlm7NkdJrc9/8aNze448970h9k1u7WzM86gDj04OBp4e/kPuCnFRx5fm6XtDhT3Db8ceFXP423AOQs0vhdSzBF6LMUlUP4IvLhCv18DX6S4VMlDKm7rCRQTxE+geEd0KfD4AX1uR/FRx9Ty7QftGMp2m5ZhexDwIop3rAvSr+62mv4wu/WzWze3Znfx57Zuduvktg3ZHWfem/4wu2a37qPpJ8utNiF8tkniEfFY4HEUYflEz1PXA9/MzAsHbGspxfUD+yexz3kL4zoi4vbAdhST7R9FcSmQczPzGRX6PbBcvCAz/zFX+7LPbRPshxzjA4CvU/whPSMzb1iofnW31WTjyu44c1tubyzZrZvbsq/Zrcl97qx9G5vdcea9yczurH3N7gCNnk+G8x0AAA1dSURBVCMM3BARO2bmSXDb3J0Zf+DM/BHwo4g4LIu7xQzrixS3V3waxR/HvsCVc3WIiA2A/wTunZm7RMRyYIfM/MyAbd0C3FR+vZXiXeGM1xGMiJ0y8wcR8c99T90/IsjBt1k8PiKeSfHOeM53PRHxM6ZPQl+PYt7TqeW2ZrweY51+dbfVIuPK7tC5LcfT9OxWzm25LbM7Go3e55ZjqpPdyrktt9Ho7I4z7y1idjG7dTT9iPBWwOcp5v4AXAvsl5mzTviOiB8y8xmFc75Ti4gzMvMRUZwtuUW57vTM3GaOPt+iuJPNGzNzyyhulXhWZj50wLb+BvyMYj7T93OOWwVGxNsz860R8bkZns4cfLmr6ynmFt1CsVOIst+dZ2h737lea7YdRp1+dbfVFuPKbp3clm0and1hclu2N7sj0PR9btlm6OwOk9uyfaOzO868t4XZva292R1SowvhKVHcO5zM/HOFto/oWVybYtL3zZk5520TI+KUzNw+Ir5DceeX3wJfzcz7zdGn7vX9dgd2BLYFbgR+ApyYmcfP0WeTzLxk0Do1y0Jnt05uy35mV7Nq6j637FfnWtZD57bsZ3Zbxuze1s/sVpUTmthe5QFsAHwG+Fa5vBzYv8brnFahzdMo3kluDvwQOAPYdUCfE4C7U56JCWwP/GiIcT0I+Hfg/4AbBrSd6ezWMypsI4DnAW8ulzcGtp30v+1if4wru3VyW/ZrdHbN7WQeTd/nlv1qZ3eY3JbtzW5LHmZ3tfZmt+Kj6XOED6P8GKFc/hXFvJxZ59NExHo9i7cDHsGqj0rmsgdwUmaeBzy+fJ33Ad+co8+rKC6hsmlEnAwspbgTy5wi4miKi3//muKMx32AU2dp+yCKCfl36Zvzc2d6JunP4WMU84p2Av4D+AtwKDDnRziat8MYT3br5Baan11zOxmH0ex9LtTI7jC5Ldub3fY5DLNrdmtoeiG8fmYeGRFvAMjMmyPilgF9zqCY8xPAzRQXn65ye8AtMvO6qYXMvCYiBp09+XPgGOBvFGebfp3ij2+Qd1HMCxr0s0BxxufTKK7tt2vP+uspLh0yyHaZ+fCIOAsgM6+NiMVyb/kmG1d26+QWmp9dczsZTd/nQr3sDpNbMLttZHYLZndITS+E/xoRd2fVrQW3B/40V4fM3KTmtm4XEXfLzGvLba3H4N/PF4A/U5wFCvAciruf7DGg3znAgRExdQ/wHwGfyBnuAZ6Z3wC+ERE7ZOZPq/0o09wUEUtY9TtcSvGOTwtrXNmtk1tofnbN7WQ0fZ8L9bJbObdgdlvK7GJ262h6IVznY4QDgS9OvVuLiLsBe2fmx+bqR3EXmp9GxFHl8h7AOwf02Twzl/cs/zAifj6gD9S7B/hZ5c/Wf93COa8aQTGR/xjgHhHxTorf35sqjFHzM67s1sktND+75nYymr7PhXrZrZNbMLttYnanM7sVNfqqERGxNvBS4CkUh/V/CnwkM/8+R5+ZLqpd6SLRUVzTb+qyKT/IzDkDGhFHAB/NzFPK5e2AAzNznwH9zsnp9wCfcV3f80cBv6R4B3kwxb3Gf5GZrxjwY03NGXoCxcc/x2fmLwb10fyMM7vD5rbs0/jsmtvxa/o+t+wzdHbr5LZsY3Zbwuyu1s/sVtT0QvhIio8Rvliueg5w18yc9WOEKC66vEWWP1h5mP/czHzIAozvFxTzcS4rV90HuIBirlHm7BeGPhPYIzN/XS5vSnHplYfPsa2zMvNhUV63MCLWBH6cmdvP0n69mdZPycxrBvx4mgezO61P5eya28lqem7L1x86u3VyW7Yzuy1hdlfrZ3YravrUiDofI3wb+EpEfLJc/tdy3ULYuWa/11L8LBeXy8uAfxnQZ2o+0HURsTnwe+Aec7TvPQngPhQXFw+KCfSXAXXnRqkas7vKMNk1t5PV9NxCvezWyS2Y3TYxu9OZ3YqaXgifGRHb932MsHJAn9dThPnfyuXvAZ9eiMFl/TuZnAx8kuLjh+uA71B8jDOXT5Xzl95EMQ9qHeDNc4xtE4CI+G/gmMw8rlzeBXh6zXGrOrO7SuXsmtuJa3RuoXZ26+QWzG6bmN3pzG5FTZ8aUffj2zsA98nMC8Yy0CHV/Ajn9hR3vVlGMXEeit/BwQO29bPsu33jTOs0WmZ3Wp+hs2tuJ8PcrtbP7LaE2V2tn9mtqOlHhIf+GCEidgPeC6wFbBLF/ccPzszdRj24eajzEc43KC4FcwbwjyG29duIeBNwRLn8XIrbQWphmd1V6mTX3E6GuZ3O7LaH2Z3O7FbU6EK45scIb6W4J/cJ5WucHRFNm99S5yOcjTKzzvyivSl+J8dQzAE6sVynBWR2p6mTXXM7AeZ2NWa3JczuasxuRY0uhGu6KTP/FBG965o2/+MRwE8iYtpHOOUZrLN9hPOTiHhoZv5smA2VZ3u+IiLulJl/nd+wtcDMbsnctspizS2Y3cXO7PboanYXYyF8fkQ8B1gSEZsBLwd+MuEx9atzZHdHYL+IuITiY45g7j8CACLikRST/9cB7hMRWwL/mpkvqTEGLSyzWzK3rbJYcwtmd7Ezuz26mt1GnyxXR0TcEXgj8ORy1Xcp5vwMM6+2cSLivjOtH/RxUEScSnF3mBVZXiQ8Is7LzM1HP0rNh9md1sfctsRizS2Y3cXO7K7Wp5PZXYxHhJdl5hspwg1ARDyOcg5QW83jcldk5uV9H/3cMv8RaQGY3en9zG07LMrcgtntALO7er/OZfd2kx7AAjgyIl4XhTtExEeAd016UBN0eflxR0bEmhHxGmDR3zKxpczuKua2PcztdGa3PczudJ3M7mIshLejmEz+E+B0ikt/PGqiI5qsFwMHAhsCVwBblctqHrO7irltD3M7ndltD7M7XSezuxinRtwE3ADcAVgbuCQzb53skMYvIg7JzNcDj8/M5056PKqk89k1t63U+dyC2W0ps4vZXYxHhE+nCPbWwKOBvSPiqMkOaSL+KYqJPm+Y9EBUmdk1t21kbgtmt33MbqHT2V2MhfCLgAuB/5eZvwNeBpwz2SFNxLeBa4EtIuLPEXF979dJD04zMrvmto3MbcHsto/ZLXQ6u4vx8mkfB24FdsrMB0fE3YDvZuY2Ex7aRETENzJz90mPQ4OZ3VXMbXuY2+nMbnuY3em6mt3FOEd4u8x8eEScBZCZ10bEmpMe1KR0MdQtZnZL5rZVzG0Ps9sqZrdHV7O7GKdG3BQRSyhvkxgRS2neLRPHJiL+OSIujIg/demjjpYyuyVz2yrmtofZbRWz26Or2V2MUyOeC+wJPBz4PMVdUt6UmV2cAE9EXATsmpmL/lqAbWd2VzG37WFupzO77WF2p+tqdhddIQwQEQ8CnkBxb+3ju/aP2isiTs7MLl8XsVXMbsHctou5XcXstovZXaWr2V2UhbBWiYgPAfcEvg7cdv/0zPzaxAYlDWBu1VZmV23V1ewuxpPlNN2dgb8BT+5Zl8CiDrZaz9yqrcyu2qqT2fWIsCRJkjppMV41Qj0iYqOIOCYi/lg+jo6IjSY9Lmku5lZtZXbVVl3NroXw4vc5YAVw7/LxzXKd1GTmVm1ldtVWncyuUyMWuYg4OzO3GrROahJzq7Yyu2qrrmbXI8KL39UR8byIWFI+ngdcPelBSQOYW7WV2VVbdTK7HhFe5CLivsBHgB0ozv78CfCyzLx8ogOT5mBu1VZmV23V1exaCC9yEfF54JWZeW25vB7wvsx8wWRHJs3O3KqtzK7aqqvZdWrE4rfFVKgBMvMa4GETHI9UhblVW5ldtVUns2shvPjdLiLuNrVQvsPzRipqOnOrtjK7aqtOZnfR/4Di/cBPI+KocnkP4J0THI9UhblVW5ldtVUns+sc4Q6IiOXATuXiDzLz55Mcj1SFuVVbmV21VRezayEsSZKkTnKOsCRJkjrJQliSJEmdZCE8RhFx14h4yTz6HxcRdx3lmKQqzK7ayNyqrczu+DhHeIwiYhnwv5m5+YSHIg3F7KqNzK3ayuyOj0eEx+vdwP0i4uyIeG/5OC8ifhYRewJExOMi4sSIODYiLoiIT0TE7crnLo2I9cvv94mIcyPinIg4fII/k7rB7KqNzK3ayuyOS2b6GNMDWAacV37/TOB7wBJgA+Ay4F7A44C/A5uWz30PeFbZ51JgfeAhwK+A9cv16036Z/OxuB9m10cbH+bWR1sfZnd8D48IT86OwP9k5i2Z+QfgR8A25XOnZebFmXkL8D9l2147AUdl5lVw220QpXExu2ojc6u2MrsLyEK4mfonbjuRW21hdtVG5lZtZXbnyUJ4vK4H1i2//zGwZ0QsiYilwGOA08rnto2ITcq5PnsCJ/W9zg+APSLi7nDb/cClhWR21UbmVm1ldsfEQniMMvNq4OSIOA/YATgXOIciqK/LzN+XTU8HPgr8ArgEOKbvdc6nuP/3jyLiHOAD4/kJ1FVmV21kbtVWZnd8vHxaw0TE44DXZObTJj0WaRhmV21kbtVWZnc0PCIsSZKkTvKIsCRJkjrJI8KSJEnqJAthSZIkdZKFsCRJkjrJQliSJEmdZCEsSZKkTrIQliRJUif9fw6ouXp/VUXZAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 864x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 4, figsize=(12, 4), sharey=True)\n",
    "for ax, name, usage in zip(axs, author_names, author_usage):\n",
    "    ax.bar(torch.arange(10), usage)\n",
    "    ax.set_title(name)\n",
    "    ax.set_xlabel(\"topic\")\n",
    "    ax.set_xticks(torch.arange(10))\n",
    "    ax.set_xticklabels(keywords, rotation=90)\n",
    "    # ax.set_ylabel(\"usage\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9cDAAFhl7bto"
   },
   "source": [
    "## Conclusion\n",
    "\n",
    "This notebook demonstrates a very naive implementation of CAVI for LDA and applies it to the Federalist papers -- a collection of documents authored by Hamilton, Jay, and Madison that were influential in forming the US Constitution. \n",
    "\n",
    "There are many ways this implementation could be improved. For example,\n",
    "- While it does work with word counts, it does not take advantage of the sparsity of the data matrix. It explicitly instantiates parameters for the posterior over $c_{n,v}$ even when $y_{n,v}=0$. We could improve performance by leveraging this sparsity.\n",
    "\n",
    "- It operates in \"batch mode,\" which is fine for small datasets like this one, but can become intractable for massive corpora, like all the pages of Wikipedia. For those regimes, it is better to work with _stochastic variational inference_ (Hoffman et al, 2011)."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Lecture 9: Latent Dirichlet Allocation and CAVI",
   "provenance": []
  },
  "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.11"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}