GMM_EM.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Copy of GMM-EM.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9OD12qixL5Gt"
      },
      "source": [
        "# EM-algorithms in Gaussian mixture model"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "E8pSoiW4HixU"
      },
      "source": [
        "import numpy as np\n",
        "import tensorflow as tf\n",
        "import tensorflow_probability as tfp\n",
        "import matplotlib.pyplot as plt\n",
        "from tqdm import tqdm\n",
        "\n",
        "tfd = tfp.distributions"
      ],
      "execution_count": 117,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-U3QrU8oHrZX"
      },
      "source": [
        "# for reproducibility\n",
        "np.random.seed(52)\n",
        "tf.random.set_seed(1234)\n",
        "# set style\n",
        "plt.style.use(\"ggplot\")"
      ],
      "execution_count": 118,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F6-H3n01I4Rc"
      },
      "source": [
        "## [Gaussian mixture model](https://en.wikipedia.org/wiki/Mixture_model#Gaussian_mixture_model)\n",
        "We have a mixture of Gaussian distrubtions, with observed data $\\textbf X =(x_1,...x_n)$ and latend variable $\\textbf Z =(z_1,...,z_m)$  (class of the observed data).<br/>\n",
        "\\begin{align}\n",
        "\\textbf{Z} &\\sim Cat(\\pi) \\\\\n",
        "\\textbf{X} &\\sim \\mathcal{N} (\\mu,\\sigma)\n",
        "\\end{align}\n",
        "<br/>\n",
        "With the exmaple here,  $\\textbf Z$ is two-class mixture of Gaussian distribution. Small vertical lines are observed data $\\textbf X$. Color represents the class, which is only for visualziation but unknow in practice.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FO_N_AhjHsDU"
      },
      "source": [
        "# we take 1D two-class Gaussian mixture model as an example\n",
        "# true values for generating data\n",
        "true_pis = [0.25, 0.75]\n",
        "true_mus = [2.0, 5.0]\n",
        "true_sigmas = [1.0, 1.8]\n",
        "sample_size = 1000"
      ],
      "execution_count": 119,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6nUIj9l9HsF7"
      },
      "source": [
        "# two-class categorical distributions\n",
        "states = tfd.Categorical(probs=true_pis).sample(sample_size).numpy()\n",
        "# Bernoullie also works for 2-class situation\n",
        "# states = tfd.Bernoulli(probs=0.75).sample(sample_size).numpy()\n",
        "Gaussian_models = [tfd.Normal(loc=2.0, scale=1.0), tfd.Normal(loc=5.0, scale=1.8)]\n",
        "# sample with known class for plotting\n",
        "x = np.array([Gaussian_models[g].sample().numpy() for g in states])"
      ],
      "execution_count": 120,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "VvBSlrFfH0Ay"
      },
      "source": [
        "# # a simpler way to generate distirbutin, but not able to show the latent class of sampled data.\n",
        "# univariate_gmm = tfd.MixtureSameFamily(\n",
        "#     mixture_distribution=tfd.Categorical(probs=true_pis),\n",
        "#     components_distribution=tfd.Normal(loc=true_mus, scale=true_sigmas),\n",
        "# )\n",
        "# # x1 = univariate_gmm.sample(sample_size).numpy()\n",
        "# plt.plot(\n",
        "# xx, univariate_gmm.prob(xx), label='GMM',ls=\"dashed\", alpha=0.5\n",
        "# )\n",
        "# plt.legend()\n",
        "# plt.show()"
      ],
      "execution_count": 138,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 391
        },
        "id": "15ajf3tVH1d5",
        "outputId": "1e7c8d45-5b73-43c1-9a43-7bf7089a377b"
      },
      "source": [
        "# plot the sample data\n",
        "xx = np.linspace(-5, 12, 1000)\n",
        "plt.figure(figsize=(16, 6))\n",
        "plt.title(\"Gaussian Mixture Model\")\n",
        "plt.plot(xx, Gaussian_models[0].prob(xx) * true_pis[0], label=\"Gaussian distribution 1\")\n",
        "plt.plot(xx, Gaussian_models[1].prob(xx) * true_pis[1], label=\"Gaussian distribution 2\")\n",
        "plt.plot(\n",
        "    xx,\n",
        "    (\n",
        "        Gaussian_models[0].prob(xx) * true_pis[0]\n",
        "        + Gaussian_models[1].prob(xx) * true_pis[1]\n",
        "    ),\n",
        "    label=\"Gaussian Mixture Model\",\n",
        "    ls=\"dashed\",\n",
        "    color=\"black\",\n",
        "    alpha=0.5,\n",
        ")\n",
        "# plt.plot(\n",
        "#     xx, univariate_gmm.prob(xx), label=\"mixturemodel generation\", ls=\"dashed\", alpha=0.5\n",
        "# )\n",
        "plt.fill_betweenx(Gaussian_models[0].prob(xx) * true_pis[0], xx, alpha=0.2)\n",
        "plt.fill_betweenx(Gaussian_models[1].prob(xx) * true_pis[1], xx, alpha=0.2)\n",
        "plt.vlines(x, 0, 0.005, color=np.array([\"C0\", \"C1\"])[states])\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "execution_count": 122,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA64AAAF2CAYAAABu9MUZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeVxVdf748dflXnZkxw3FBcQNQQ1RlEE2QUtyt8w0JUcrp1HHnHHJGae0Mq2cGls0R00dUykzt9zSHBWX3AVMUDRc2WTfvPee3x/+vN8IEFDgIryfj4ePh2f7fN7nfC5c3ufzOZ+jUhRFQQghhBBCCCGEqKNMjB2AEEIIIYQQQgjxMJK4CiGEEEIIIYSo0yRxFUIIIYQQQghRp0niKoQQQgghhBCiTpPEVQghhBBCCCFEnSaJqxBCCCGEEEKIOk0SVyGEEA3SvHnz8PDwMHYYZarLsdUVrVu3Zv78+VU6RqVSsXbt2hqKSAghRE2SxFUIIUS1ycjIYNasWXTq1AkrKyscHBzo2rUrc+bMITk52djhlfDGG29w9OjRWqmrdevWqFQqPv7441Lbpk2bhkqlIiws7LFi02g0rFq16nFDfSxXr15FpVKhUqk4f/58qe3dunVDpVJVOeEUQgghJHEVQghRLZKTk+nWrRsbN25k1qxZHD16lDNnzrBkyRLS09NZvHixsUMswcbGBmdn51qrz83NjS+//LLEusLCQr766itatWpl1Nh+7969e491vJubG8uXLy+x7vjx4yQkJODk5PRYZQshhGiYJHEVQghRLV577TWKi4s5ffo0Y8aMwdvbm1atWhEUFMTnn3/OkiVLDPvu2bOHoKAgHB0dsbOzo2/fvhw/frxEeWUN6wwLC2PcuHGG5S1bttCtWzesrKywt7fHz8+P06dPA/eTr7/85S+0aNECc3NzmjVrxvPPP2849vfDcZOSkhg6dCjNmzfHysqKLl26sGbNmhL1BwUFMWHCBN5++22aNm2Ko6MjY8eOJTc3t8Lr8/zzz3PlyhWOHTtmWBcdHY2DgwN9+/Ytse9vY1MUhWeeeYYePXoYEkq9Xk9YWBiBgYHodDpat26NTqdj/Pjxhh5PgFWrVqHRaEqUff36dVQqFQcOHADgwIEDqFQqtm/fTkBAABYWFoYE+5NPPqFDhw5YWFjQrl07FixYgFarrfBcX375ZdauXUthYaFh3bJly3juueewsbEpsW9OTg6TJk3CxcUFc3NzfH192b17d4l9zp49S+/evTE3N6ddu3Zs3LixVJ25ublMmTIFV1dXrKys6NatG99++22FsQohhHgySOIqhBDisWVkZLBjxw5ef/11bG1ty9znQTIF95OM1157jZiYGI4cOUK7du3o378/6enpla7z9u3bjBgxglGjRhEbG0tMTAxTp041JGqffPIJGzduZO3atSQkJPD999/Tq1evcsvLzc0lJCSEnTt3cv78eSZOnMj48ePZv39/if2io6PJyMjgwIEDfP3112zbto2FCxdWGG+jRo14/vnnS/RELlu2jAkTJpS4Nr+nUqlYtWoVN27cYNasWQC8++67nD59mnXr1qFWqzlx4gRqtZolS5Zw69Ytbt26VWE8vzd9+nT+9re/ER8fT2RkJPPmzWPx4sW8++67xMfH869//YsvvviCf/7znxWWFRwcjLOzM9HR0cD95PTrr7/mj3/8Y6l9o6Ki2LVrF2vXruXMmTP06dOHgQMHcvHiRQAKCgp4+umnsbe35/jx43z11VcsWrSIlJQUQxmKohAZGcnZs2fZsGEDFy5c4NVXX+X5559n3759Vb4WQggh6iBFCCGEeEzHjh1TAOXbb78tsd7f31+xtrZWrK2tlU6dOpV7vE6nU+zt7ZW1a9ca1gHKmjVrSuwXGhqqvPTSS4qiKMqpU6cUQElKSiqzzD//+c9KcHCwotfry9z+j3/8Q3F3d3/oeT377LPKhAkTDMt9+/ZVvL29S+zzyiuvKL169XpoOa1atVLefvtt5dixY4q1tbWSnZ2txMfHK6ampsrt27eVl156SQkNDX1obD/++KOiVquVefPmKRqNRtm8eXOJ7Wq1Wlm5cmWJdStXrlTUanWJdcnJyQqg7N+/X1EURdm/f78CKF999ZVhn7y8PMXS0lLZuXNniWNXr16t2NnZlXueSUlJCqD873//UxYuXKgEBgYqiqIon332mdKlS5cS10JRFCUhIUEBlO3bt5cop1u3bsr48eMVRVGU5cuXK9bW1kpGRoZh+/nz5xXAUM7+/fsVc3NzJTMzs0Q548ePVwYNGmRYLuszJYQQ4smgKTejFUIIIapIUZQSyxs2bKCoqIhPP/20xLDNpKQk/v73vxMTE0NKSgp6vZ78/HyuXbtW6bq8vb2JiIjAy8uLfv36ERQUxNChQ2nZsiUA48ePp1+/fnh4eNCvXz/69etHZGQkZmZmZZaXn5/PW2+9xdatW7l16xbFxcUUFRURHBxcYj8fH58Sy82bN2fXrl2VitnPzw8PDw/Wr1/PxYsXiYyMpEmTJpU6Njg4mOnTpzNv3jxeeeUVBg8eXKnjKsvPz8/w/9jYWAoKChg2bFiJ3mCdTkdhYSGpqam4uLg8tLxx48Yxd+5cfvnlF5YvX15mb2tcXBwAgYGBJdYHBgYSExNj2Kdjx444ODgYtnt5eWFnZ2dYPnHiBMXFxbi6upYop7i4mHbt2lV06kIIIZ4AkrgKIYR4bB4eHpiYmBAfH19i/YMk0tHRscT6gQMH4uzszNKlS2nZsiVmZmYEBARQXFxs2EelUpVKhH87aZBarWbnzp2cOHGCvXv38s033zBz5kw2bdrEwIED6dq1K0lJSezZs4f9+/czZcoU5s6dy9GjR8sczjxjxgy2bNnChx9+SPv27bG2tmb69OlkZWWV2O/3ia9KpUKv11f6Wk2cOJHPPvuM5ORk1q1bV+njdDodhw8fRq1Wc/nyZRRFeegQYwATk9JPBJU38ZK1tbXh/w/OZ9OmTXh6epba9/ftWZbGjRszaNAgJk+eTHx8PGPGjKnwmEel1+uxs7PjxIkTpbaVd6NCCCHEk0WecRVCCPHYHB0dGTBgAJ988kmpRO/30tPTiYuLY+bMmURERNCpUycsLCxKPLMI9xOfmzdvGpaLiooMPXQPqFQq/Pz8mD17NgcPHqRv376sXLnSsN3GxoYhQ4bw8ccf8/PPPxMfH89PP/1UZlwHDx5k9OjRjBw5Eh8fH9q2bculS5eqeikq9OKLL5KQkECjRo3o169fpY+bN28eiYmJHD58mOPHj/P++++X2G5mZoZOpyuxrnHjxuh0Ou7cuWNYd+rUqQrr6ty5MxYWFly5cgUPD49S/9RqdaVinjRpEvv27WP48OHY29uXWQ/cv/a/dfDgQby8vADo1KkT8fHxZGZmGrbHxsaW+Jz5+vqSmZlJYWFhqVjd3NwqFasQQoi6TXpchRBCVItPP/2UPn360K1bN+bNm0fXrl2xsbHhl19+Ydu2bYZkx8HBARcXF5YvX467uzvp6en89a9/xdLSskR5YWFhfP755wQGBtKoUSMWLFhQokf2yJEj7Nu3j/DwcJo1a0ZCQgLnzp3j5ZdfBmDRokU0b96crl27YmVlxfr161Gr1WX2IAK0b9+eLVu2MGzYMGxsbPjwww+5efNmpYfyVpatrS03btzAxMSkzB7Rsvz000+89957bNu2jZ49e7Js2TJefPFFgoODDUN827Rpw/79+xkwYABmZmY4Ozvj5+dHo0aNmDlzJrNnz+by5cu89dZbFdZnY2PD7NmzmT17tuEds1qtlvPnz3P69OlKTUYFEBoaSmpqaqmZhB9wd3dnxIgRvPbaa3zxxRe0atWKzz77jAsXLvDf//4XgBdeeIG5c+fy4osvsmDBAgoKCpgyZUqJz0tISAhhYWEMHTqU999/H29vb+7evcuRI0ewsLAoc5iyEEKIJ4v0uAohhKgWbm5unD59mhEjRvDuu+/Ss2dPOnfuzPTp0/H39zfM7mpiYsKmTZu4fPky3t7ejBs3jqlTp9KsWbMS5S1evBgvLy8iIiIYMGAAgYGB9OjRw7Ddzs6OmJgYBg0aRLt27YiKimL06NHMnTsXuJ8gfvjhh/j7+9OlSxc2b97MN998Q/v27cuM/6OPPqJVq1YEBwcTGhqKq6srw4cPr5FrZWdnR6NGjSq1b0ZGBi+++CJTpkwhIiICgJEjRzJu3DhGjRpFTk4OAB988AEnT56kdevWhudPHR0dWb9+PUePHsXb25u33367VE9teebOncuHH37I8uXL8fHxISAggI8++ojWrVtX6VydnZ2xsLAod/uXX35JREQEL774Ij4+Phw+fJht27bRoUMHAKysrNixYwfp6en4+fkxevRopk2bRuPGjQ1lqFQqvv/+e4YOHcq0adPo0KEDzzzzDNu3b8fd3b1K8QohhKibVMrvHyASQgghhBBCCCHqEOlxFUIIIYQQQghRp0niKoQQQgghhBCiTpPEVQghhBBCCCFEnVapWYXPnDnDypUr0ev1hIaGlnrpeVxcHKtXr+batWtMnTqVXr16AXDhwgVWr15t2O/mzZtMmTIFPz8/li5dSlxcHFZWVgBMnjy5yhM+CCGEEEIIIYSo/ypMXPV6PStWrODNN9/EycmJWbNm4evrS4sWLQz7ODs789prr7F169YSx3p5ebFo0SIAcnNzef311/Hx8TFsHzNmjCHJFUIIIYQQQgghylJh4pqYmEjTpk0N77Hr3bs3J06cKJG4PpiSXqVSlVvO0aNH6datG+bm5o8bsxBCCCGEEEKIBqTCxDUjIwMnJyfDspOTEwkJCVWu6PDhwwwcOLDEuvXr1xMdHY2XlxejR4/G1NS01HF79+5l7969ALz33nslXj5f12g0GrRarbHDEDVE2rf+kratv6Rt6zdp3/pL2rZ+k/atv6qjbc3MzMou+7FKraS7d+/y66+/lhgm/MILL2Bvb49Wq+WLL75gy5YtZb7oPSwsjLCwMMNyWlpabYT8SJydnet0fOLxSPvWX9K29Ze0bf0m7Vt/SdvWb9K+9Vd1tG3z5s3LXF/hrMKOjo6kp6cbltPT03F0dKxS5TExMfj5+aHR/F+e7ODggEqlwtTUlODgYBITE6tUphBCCCGEEEKIhqHCxNXd3Z1bt26RkpKCVqvlyJEj+Pr6VqmSw4cP06dPnxLr7t69C4CiKJw4cYKWLVtWqUwhhBBCCCGEEA1DhUOF1Wo1UVFRLFiwAL1eT3BwMC1btmTDhg24u7vj6+tLYmIiixcvJi8vj5MnT7Jx40Y+/PBDAFJSUkhLS6NTp04lyv3444/Jzs4GoFWrVkycOLEGTk8IIYQQQghRWYqiUFhYiF6vf+jEq4/jzp07FBUV1UjZwrgq27aKomBiYoKFhUWlP2cqRVGUxw2wNt28edPYIZRLxuvXb9K+9Ze0bf0lbVu/SfvWX9K2xlNQUICpqWmJR/yqm0zOVH9VpW21Wi337t3D0tKyxPpHfsZVCCGEEEII0TDo9foaTVqFeECj0aDX6yu9vySuQgghhBBCCIAaGx4sRFmq8nmTxFUIIYQQQghRZ6SmpjJ58mT8/f3p378/kZGR7Ny5s8brPXv2LHPnzq32co8cOcLYsWMB2L17N//+97/L3ffChQvs27evUjF+8MEHfP7551WKZfny5RQUFBiWx4wZQ1ZWVpXKKMvRo0eJiIigefPmbNu27bHLK4uMAxBCCCGEEELUCYqiEBUVxYgRI1i6dCkA169fZ/fu3TVet4+PDz4+PjVaR3h4OOHh4eVuj42N5dy5c4SGhpbaptVqHzvGL7/8kmHDhhmeK12zZs0jl/Vbrq6ufPTRRyxbtqxayiuL9LgKIYQQQggh6oRDhw5hZmZm6KEEaNGiBVFRUQAkJyczZMgQIiIiiIiI4MSJE0DJXk2AOXPmsGHDBgDeeecdgoKCCAsL46233gJg69athISEEBYWxtChQ0uVcfr0aSIjIwkPD+fZZ58lMTERgA0bNjBhwgRGjx5Nnz59mD9/fpnnsX//fgIDA4mIiCjRW7xhwwbmzJlTZgzFxcUsXryY77//nn79+rFlyxY++OADXn/9dQYNGsSf//znUucZGxtLZGQkffr0Yd26dQ+9FitWrODOnTuMGDGC4cOHA9CzZ08yMjIA+OKLLwgJCSEkJITly5cbrnffvn2ZMWMGwcHBjBo1qkSP7QMtW7akU6dOmJjUXHopPa5CCCGEEEKIUvRfL0dJTqrWMlUt28CLr5a7/dKlS3h5eZW73dnZmfXr12NhYcGVK1eYPHnyQ4cRZ2RksHPnTg4ePIhKpTIMi12yZAnr1q2jWbNmZQ6V9fDwYPPmzWg0Gg4ePMjChQsNyVxsbCy7du3CzMyMwMBAxo8fj6urq+HYwsJCZsyYwcaNG2nTpg2vvPJKmbH9PgYzMzPeeOMNzp07x4IFC4D7w4ETEhLYvHkzlpaWHDlypEQZ8fHxbN26lYKCAsLDw8vsqX3g5ZdfZtmyZWzatAlHR8cS286dO8fGjRvZtm0biqIwcOBA/P39sbOzIykpiaVLl7Jo0SImTZrEjh07GDZsWLn11BTpcRVCCCFEnaMoCnl5eaSlpXHjxg3DHf6CggLS0tIoKCjgCXujnxDiEcyePZuwsDCefvppAO7du8eMGTMIDQ1l0qRJXLp06aHH29raYm5uzvTp09mxY4dhiKyvry/Tpk1j3bp16HS6UsdlZ2czadIkQkJC+Oc//8kvv/xi2BYQEICtrS0WFhZ4enpy48aNEscmJibi5uZG27ZtUalU5SZ5FcXwQHh4eKlXxjwQERGBpaUljo6O9O7dmzNnzjz0epTn+PHj9O/fHysrK6ytrRkwYADHjh0D7vemPriZ4O3tTXJy8iPV8bikx1UIIYQQRqXX69HpdJiamnLjxg327t1LZmZmiXcBDhs2jNatW5OcnMzWrVsBsLOzw9LSEmdnZ/z8/HBwcDDWKQhRL5k8/8dar9PT05MdO3YYlt955x0yMjIYMGAAcH9yIRcXF/bs2YNer6dt27bA/Ver/PZmVlFRkWH99u3bOXToENu3b2flypVs2rSJhQsXcurUKfbt28eAAQNK9douWrSI3r17s2LFCpKTkw1DawHMzMwM/zcxMXnkd9JWFMMDVlZW5Zbx+1l5VSpVudfiUZmbmxv+r1arKSwsfKzyHpX0uAohhBCi1hUWFnLx4kV27NjBZ599Rnx8PHD/DzRbW1u8vb0JDg5m4MCBDBs2jKZNmwL3X0wf1v9pegX0pb2XN2qNGZcvXza8dzIxMZHDhw+Tnp5utHMTQjy6gIAAioqKWL16tWHdb5+pzM7OpnHjxpiYmPDNN98YeipdXV25dOkSRUVFZGVlcejQIQDy8vLIyckhNDSUefPmERcXB8DVq1fp3r07M2bMwMnJiZs3b5aIIycnx/B7Z+PGjVU6Bw8PD5KTk7l69SoA3333XZn7lRWDjY0Nubm5la5r165dFBYWkpGRQUxMDD4+PuVeC6Dc8nv27MmuXbsoKCggPz+fH374gZ49e1bpvGua9LgKIYQQotZotVp27dpFQkICOp0OS0tL2rZti7OzMwAODg4MGTIEgPx7Os7fyScxvZDLSancyC7mboGWIp0KsPn/JXbB2lnFxZ9SaeWQQ1HSL+RevcCRmBjcWraka9eueHh4oFarjXPCQogqUalUrFixgnnz5vHZZ5/h5OSEpaUls2fPBuCll15i4sSJREdHExwcbOiNdHV1JTIykpCQENzc3AxDW3Nzc4mKiqKoqAhFUfjHP/4BwPz580lKSkJRFAICAujcuTMxMTGGOF599VWmTp3Kv/71r4c+N1oWCwsL3n//fcaOHYulpSU9e/YsM1ksKwZXV1eWLl1Kv379+NOf/lRhXR07dmTEiBFkZGQwdepUQ7Jd1rUAGD16NKNHj6ZJkyZER0cb1nfp0oURI0bwzDPPADBq1Ci8vLwqPSz4zJkzvPzyy2RlZbFr1y4++OAD9u/fX6ljK0ulPGEPiPz+bkhd4uzsTFpamrHDEDVE2rf+kratv6Rt64aioiJu3bpF69atAfjmm2+ws7OjY8eONGvWrMQslPn3dBy6lsOha9lcSMlHpwcTFTSxMaWJtSl2FhpszNSoTcDK0pLMnHxyinXcLdByM6eYzEIduqICVGlJ2GUl0cysGL9OHowYMcJIZy8ehfzsGk9+fv5Dh6ZWB41G88jDa0XdVtW2Levz1rx587LLfqzIhBBCCCHKkZeXx88//8y5c+fQ6XRMmjQJS0vLMicquZVTzOa4DA5czaJIq+BireEPbrZ0dLGklb05purSTzfZ29uTmZlZYl1OkY6E9AIupjkTe7sT524nc/2uKZq4dMLaNCLt9g3atGlTY+cshBCiZkjiKoQQQohqVVBQwIkTJzh9+jR6vZ527drx1FNPlTkr5q2cYtafS+Pg1Ww0Jiq6NrOmVwsbWtmbl5p0pDIamavp3tyG7s1tuNfZifN3nDl6PZdVp1P5z/b/4ZJ6nv5+nYjoF4adnV11nK4QQohaIImrEEIIIapVQUEBJ0+epH379vj7+5c522/+PR2bLqSz5WIGapWKvm1s6dvaFlvz6vvTxFRtYkhir2UWsjvBjJ9P6vh0zzliYi8zbnAE3t7ej5QgCyGEqF2SuAohhBDisSiKQnx8PLdu3SI0NBRHR0cmTpyItbV1mfv/fCOXpcduk1GgpYerDU972ldrwlqWVvYW/LFHc0LdHdl0yoPjRw+QsOwbop6+w9Bnwmu0biGEEI9PElchhBBCPLLbt2+zf/9+bt68SdOmTSkuLsbMzKzMpDX/no4vf05h35UsmtmYMqVXM9zszcsotea0dbRgRqgnMe2a8vXeGNbcssU6IZNwDzvpeRVCiDpMElchhBBCVFlBQQH/+9//uHDhAlZWVkRERNC5c+dyk7+rdwt57383uJ1zj9C2doR72KMxMU6iaKJS0aeVHR1fCGPThTSWHrvFlh0/8FKAJz2f6m6UmIQQQjxc6Sn6hBBCCCEqoCgKV65c4amnnmL8+PF4eXmVm7QeSMpixq5r5BXredWvKU97Ohgtaf0tR0sNf/RtQn/3Rpz7NZWZX3zDpu27ecLeFChEvZOamsrkyZPx9/enf//+REZGsnPnzhqv9+zZs8ydO7fayz1y5Ahjx44FYPfu3fz73/8ud98LFy6wb9++SsX4wQcf8Pnnn1cpluXLl1NQUGBYHjNmDFlZWVUqoyxffPEFQUFBBAUFMXLkSK5fv/7YZf6e9LgKIYQQolKysrI4c+YMgYGBWFlZERUVhZmZWbn7K4rCf8+lsfFCOu6OFozxcaGRuboWI66YiUpFuKczbSaO5l9fb+ejbw6QcjebSc8PRqORP5OEqG2KohAVFcWIESNYunQpANevX2f37t01XrePjw8+Pj41Wkd4eDjh4eU/Vx8bG8u5c+cIDQ0ttU2r1T52jF9++SXDhg0zzPK+Zs2aRy7rt7y8vNi5cyeNGjVixYoVzJ8/v8pJdUWkx1UIIYQQD6XX6/n5559ZtWoVZ8+eJS0tDeChSes9nZ6Pjtxk44V0/FrYMMm3SZ1LWn+rnbMV/xw/iOZefqzee5J3v/za2CEJ0SAdOnQIMzMzQw8lQIsWLYiKigIgOTmZIUOGEBERQUREBCdOnABK9moCzJkzhw0bNgDwzjvvEBQURFhYGG+99RYAW7duJSQkhLCwMIYOHVqqjNOnTxMZGUl4eDjPPvssiYmJAGzYsIEJEyYwevRo+vTpw/z588s8j/379xMYGEhERESJ3uINGzYwZ86cMmMoLi5m8eLFfP/99/Tr148tW7bwwQcf8PrrrzNo0CD+/Oc/lzrP2NhYIiMj6dOnD+vWrXvotVixYgV37txhxIgRDB8+HICePXuSkZEB3O81DQkJISQkhOXLlxuud9++fZkxYwbBwcGMGjWqRI/tA3369DEkw0899RS3bt0qu4Efg9xKFEIIIUS5UlJS2L17N3fu3KFt27aEhoZia2v70GPy7+mYf+A6sSkFDGhnT2jbJ2PiIydrM/7xQj/+ZduIQ3mw8UIaI72cjR2WEEbz5c93SLpbWK1ltnGw4JVeruVuv3TpEl5eXuVud3Z2Zv369VhYWHDlyhUmT5780GHEGRkZ7Ny5k4MHD6JSqQzDYpcsWcK6deto1qxZmUNlPTw82Lx5MxqNhoMHD7Jw4UJDMhcbG8uuXbswMzMjMDCQ8ePH4+r6f+dUWFjIjBkz2LhxI23atOGVV14pM7bfx2BmZsYbb7zBuXPnWLBgAXB/OHBCQgKbN2/G0tKSI0eOlCgjPj6erVu3UlBQQHh4eJk9tQ+8/PLLLFu2jE2bNuHo6Fhi27lz59i4cSPbtm1DURQGDhyIv78/dnZ2JCUlsXTpUhYtWsSkSZPYsWMHw4YNK7ee9evXExwcXO72RyWJqxBCCCHKpCgK27dvp7CwkIEDB+Lp6VlhAppXrGPej8kkZhTygrczTzW3qaVoq4elqZoZg3qz4UIa686mcT3pMq+EdsHKysrYoQnRIM2ePZvjx49jZmbGjh07uHfvHnPmzCEuLg4TExOuXLny0ONtbW0xNzdn+vTphIWFERYWBoCvry/Tpk0jMjKSAQMGlDouOzubqVOnkpSUhEql4t69e4ZtAQEBhht4np6e3Lhxo0TimpiYiJubG23btgVg2LBhrF27tlQdFcXwQHh4uKE38/ciIiKwtLTE0tKS3r17c+bMmQpvLpbl+PHj9O/f3/C7bsCAARw7dozw8HBatmxpuJng7e1NcnJyueVER0dz9uxZvvnmmyrHUBFJXIUQQghRwq1bt3B2dsbU1JSBAwdiY2NT7h9Nv5VbpOPvPyZz9W4hY7u60KVJ2e9xrevUJiqe7+KMriif9Ru3knjuJIunjZPkVTQ4E3yb1Hqdnp6e7Nixw7D8zjvvkJGRYUjsli9fjouLC3v27EGv1xuSQ41GU2JitaKiIsP67du3c+jQIbZv387KlSvZtGkTCxcu5NSpU+zbt48BAwaU6rVdtGgRvXv3ZsWKFSQnJxuG1kLJxyRMTEzQarWPdK4VxfDAw373/P5monhLjuQAACAASURBVEqlKvdaPCpz8/97bZlaraawsOxe+IMHD7JkyRKio6NLHFNd5BlXIYQQQgBQXFzM/v37Wb9+veG5MRcXl0onrW/u+5VrmYW81K3xE5u0PmCiUjHa142w/s9wIvEGbyxZ/dh//AkhKhYQEEBRURGrV682rPvtM5XZ2dk0btwYExMTvvnmG3Q6HQCurq5cunSJoqIisrKyOHToEAB5eXnk5OQQGhrKvHnziIuLA+Dq1at0796dGTNm4OTkxM2bN0vEkZOTQ9OmTQHYuHFjlc7Bw8OD5ORkrl69CsB3331X5n5lxWBjY0Nubm6l69q1axeFhYVkZGQQExODj49PudcCKLf8nj17smvXLgoKCsjPz+eHH36gZ8+elY7jwoULzJw5k6+++gpn55p5xEJ6XIUQQgjBr7/+yu7du8nKysLHx4ennnqq0scWafW8fSCZ5KwixnVrTEeX+tEzaaJSERXSFRSFvTu3Mfez9bz7+hjU6ro7yZQQTzqVSsWKFSuYN28en332GU5OTlhaWjJ79mwAXnrpJSZOnEh0dDTBwcGG3khXV1ciIyMJCQnBzc3NMLQ1NzeXqKgoioqKUBSFf/zjHwDMnz+fpKQkFEUhICCAzp07ExMTY4jj1VdfZerUqfzrX/966HOjZbGwsOD9999n7NixWFpa0rNnzzKTxbJicHV1ZenSpfTr148//elPFdbVsWNHRowYQUZGBlOnTjUk22VdC4DRo0czevRomjRpQnR0tGF9ly5dGDFiBM888wwAo0aNwsvL66HDgn/r7bffJi8vjwkTJqAoCq6urqxatapSx1aWSnnCXlb2+7shdYmzs7NhpkVR/0j71l/StvWXtG3l/Pzzz/z000/Y29sbnmeqLK1eYcGB65y+lcfYri54N629nlZ7e3syMzNrvB69ovDJ94eI+elHXh09jJfDfWu8zoZOfnaNJz8/v8aHxWs0mkceXivqtqq2bVmft+bNm5dd9mNFJoQQQognllarRaPR0LZtW/Lz8/H398fU1LTSx+sVhY9jbnHqVh4jOjvVatJam0xUKv4UGYCJtQNbU21odzWbwNZVn/xECCHEo5NnXIUQQogGJjs7my1bthgmAnF0dCQwMLBKSSvAV6dT+elqNgPa2dOrZaOaCLXOUJuoeCW4E20czHn/h7N8/7+Txg5JCCEaFElchRBCiAZCp9Nx4sQJVq1axdWrV2nSpAmP+sTQ3suZbI7PoLdbI0Lb2lVzpHWTqdqE8d0bo9yIZ+GqbzkZf9nYIQkhRIMhiasQQgjRAGRkZLB27VoOHjxIy5YtGTduHH5+fhW+l7UssSn5fHr8Nu2cLBjcwfGRynhSWZmq+evYoWisbZm59L/cSEk3dkhCCNEgSOIqhBBC1GN6vR64P8slwKBBgxg8eDB2do/WS3o7p5h3f7qOo6WGsV1dUJs0nKT1AVdHG14f+xzZRVqmfrSa/IKy32kohBCi+kjiKoQQQtRDhYWF/PTTT2zcuBFFUbCysmLs2LF4eHg8cg9poVbPgp+uo1Pg5e5NsDJtuK+F6dqmKaOHDeHqrTQWRe83djhCCFHvSeIqhBBC1COFhYXExMTw5ZdfcvLkSRwcHLh37x7AYw3pVRSFz47fJjmrmNHeLrhYV20ip/ro6R4dCBs4hJO04lhyjrHDEaLeSE1NZfLkyfj7+9O/f38iIyMNk8nVpLNnzzJ37txqKWv48OH06NGjxDwCUVFRtGvXDoDbt2/zxz/+8aFlHDlyhBMnTlRLPGVxdXXl9ddfNyxrtVq6dOnC2LFjq1ROz549ycjIeOx9KiKvwxFCCCHqidu3bxMdHU1RUREeHh74+/vTuHHjail7V2ImB5KyCfewp4OLZbWUWR+MDezM0uO3WbT/Mv8IcKFLW1djhyTEE01RFKKiohgxYgRLly4F4Pr16+zevbvG6/bx8cHHx6fayrOzs+PEiRP4+fmRlZVFSkqKYVvTpk1Zvnz5Q4+PiYnB2tqaHj16VLrOB685qwwrKysuXrxIQUEBlpaWHDx4kKZNm1a6rtomPa5CCCHqNUVRKCgoIDU11fC8Z2ZmJjdu3CA1NZXMzExDj+ST6O7du/z6668AODs7065dO1588UUGDRpUbUlrQnoBy39OoYOzJf3cG8YMwpVlqjZhjI8L6WcPMOPjr7ibnWvskIR4oh06dAgzM7MSvX4tWrQgKioKgOTkZIYMGUJERAQRERGGHskjR46UOGbOnDls2LABgHfeeYegoCDCwsJ46623ANi6dSshISGEhYUxdOjQUmWcPn2ayMhIwsPDefbZZ0lMTARgw4YNTJgwgdGjR9OnTx/mz59f7rk8++yzbNmyBYCdO3cyYMAAw7bk5GRCQkIAWLZsGX/5y18AiI+PJyQkhEuXLrFmzRqWL19Ov379OHbsGFOnTmXbtm2GMh703h45coQhQ4Ywbtw4goKC0Ol0vP322zz99NOEhYWxZs2acmMMCQlh3759AHz33XcMHjzYsO3u3btERUURFhbGwIEDiYuLA+5P9jdq1CiCg4N54403SvQqR0dH88wzz9CvXz/++te/otPpyq27qiqVjp85c4aVK1ei1+sJDQ0tcUIAcXFxrF69mmvXrjF16lR69epl2Pbcc8/h5uYG3P9C/dvf/gZASkoKS5YsIScnh7Zt2/L6669X+u6AEEII8TDXr18nLi6OgoICrl69ilarBeDVV1/FysqK8+fPc/z48RLHWFpaMmHCBMzMzLh+/Tr5+fk4Oztjb2+PiUndus+r1Wq5cuUKsbGxJCUlYWdnR1RUFBqNhoiIiGqtK7dYx/v/u0kjMxNe8HbGpAHNIFxZTlamTBj6NP9a8RWzvtjEZ2+Ma1AzLYv67UHy91vt27ena9eu3Lt3j2+//bbU9s6dO+Pl5UVBQQHff/99iW3PPffcQ+u7dOkSXl5e5W53dnZm/fr1WFhYcOXKFSZPnvzQYcQZGRns3LmTgwcPolKpyMrKAmDJkiWsW7eOZs2aGdb9loeHB5s3b0aj0XDw4EEWLlxo6CGNjY1l165dmJmZERgYyPjx43F1LT3aIiAgwJC8bdmyhffff58lS5aU2m/ChAkMHz6cnTt38vHHH7Nw4UI8PT0ZM2YM1tbWvPLKKwCsX7++3PM8f/48P/74I25ubqxdu5ZGjRqxY8cOioqKGDx4MH379jXkZL81aNAgPvroI8LCwoiPj+f555/n2LFjAHzwwQd4eXnxn//8h0OHDjFlyhT27NnDRx99hJ+fH9OmTWPv3r2GuBISEvjuu+/47rvvMDU1ZdasWXz77beMGDGi3LirosJMUa/Xs2LFCt58802cnJyYNWsWvr6+tGjRwrCPs7Mzr732Glu3bi11vJmZGYsWLSq1fu3atTzzzDP06dOHZcuW8eOPPxIeHv6YpyOEEKIhKigoID4+Hk9PT2xsbMjOziYhIYG2bdvi7e2Nra0tNjY2mJmZAdClSxdatmxJcXExxcXF5Obmkp+fb9h++vRpLl26BICJiQlOTk40a9aMfv36AVUbilXdzp8/z8GDByksLMTa2pqePXvi4+NTY4nSF8dvk5Z/jz/1bIa1WcOdjKkiPTu4Edy3L/t+/JGvdh3hpf59jB2SEPXC7NmzOX78OGZmZuzYsYN79+4xZ84c4uLiMDEx4cqVKw893tbWFnNzc6ZPn05YWBhhYWEA+Pr6Mm3aNCIjI0v0hD6QnZ3N1KlTSUpKQqVSlRiZExAQgK2tLQCenp7cuHGjzMRVrVbTo0cPtmzZQmFhIS1btiwzRhMTE0Py+OKLL1ZpaPADXbt2NSSmP/30E/Hx8Wzfvh2AnJwckpKSykxcO3XqxPXr19myZYuhB/iB48ePG5L1gIAA7t69S05ODkePHuXLL78EICwsDHt7e+B+b/m5c+d4+umngftzLjg7O1f5XMpT4bduYmIiTZs2pUmTJgD07t2bEydOlEhcHwxFquyXpqIoxMbGMmXKFACCgoLYtGmTJK5CCCGqJCsri2PHjhEXF4dOp0Oj0eDt7U2HDh3o2LEjLi4upKWllTrO3t7e8EVblv79+9OjRw/S09NJT08nNTWVnJz/m3wnOjqarKwsXFxcDP+aNGmCg4NDtZ2bXq8nNTWVmzdv8uuvvxIQEICTkxPW1ta0bt2aTp060apVqxrtDT6QlMXBazkMaGdPK3vzGqunvhg/oA8XExJZFv0DvTq7075l3X1WTIjKelgPqamp6UO3W1paVtjD+nuenp7s2LHDsPzOO++QkZFhSC6XL1+Oi4sLe/bsQa/X07ZtWwA0Gk2JIatFRUWG9du3b+fQoUNs376dlStXsmnTJhYuXMipU6fYt28fAwYMKNVru2jRInr37s2KFStITk5m+PDhhm0PbnLC/aTzwaiesgwaNIiXX36Z6dOnP/S8k5KSsLa25s6dO+Xuo9FoDI+86PX6Esm0lZVViX3nz59PUFDQQ+t8IDw8nLfeeovo6Gju3r1bqWPKoigKI0eOZObMmY9cxsNUmLhmZGTg5ORkWHZyciIhIaHSFdy7d4+ZM2eiVqsZNGgQfn5+5OTkYGVlhVp9/86to6NjubNM7d27l7179wLw3nvvVWvWXt00Gk2djk88Hmnf+kva9smj1+vZtWsXp0+fBsDf35/u3bvTrFmzEvs9Ttv+vqzf6tWrF9evX+fOnTvExcWh1+tp164dzz//PAA7duzAwsICW1tbbG1tsba2xtbWlkaNGhmeuVUUBb1eT0FBAXl5edjZ2eHo6Mjdu3f57rvvuHPnjuGPkkaNGqFWq3F2dsbZ2Rk/P79HOqequJlVyBc/J9DOxZqh3VphUgff16rWqB96A8IY/v7qS/zt869ZHZvPZz5ODfI9t9VBfi8bz507d2plREl5dQQFBbFw4ULWrl3LuHHjACguLkalUqHRaMjNzaV58+aYmZmxfv16w03LVq1akZCQgE6no7CwkEOHDtGrVy+KiorIz88nIiICf39//Pz80Gg0XL16FT8/P/z8/Dhw4AApKSmo1eoS9bi6uqLRaIiOjjasV6vVmJiYGOJXqVSo1epS5/NgfZ8+fZgyZQrDhg0z7POgnAf/z87O5u9//ztbtmxh1qxZ7Ny5k8jISGxtbcnJyTEc5+bmRmxsLEOHDjX0Pj8o60F8AMHBwaxZs4a+fftiamrK5cuXadq0KdbW1mW2w+jRo7G3t6dLly4cPnzYUFavXr3YsmULf/nLXzh8+DBOTk44ODjg7+9vWL9v3z4yMzNRq9X07duXl156iUmTJuHi4sLdu3fJzc2lZcuW5V4nc3PzSv+s1/in8tNPP8XR0ZE7d+7w1ltv4ebmVuqOwMP8tksfKPPOeV3h7Oxcp+MTj0fat/6Stn0ypaam0qZNG3r27GkYsvX7dqyptm3dujWtW7cG7g8bzsjIQFEU0tLS0Ov1nD17ltzcXMOdcYDu3bsTHBxMcXExn3zySaky/f396d27NwUFBeTn59O2bVuaNWtG8+bNsbW1RaVS1drnVKdXmLvnGnq9npGdHMjOLv38V11gb29PZmamscMowc7ChAkjnmX9+TQ+P3CR57xdjB3SE0l+LxtPUVGRIamqKRqN5qG9lCtWrGDevHn8+9//xsnJCUtLS2bNmoVWq2XMmDFMnDiRDRs2EBwcjJWVFVqtliZNmjBw4EACAwNxc3PDy8sLnU5HZmYmUVFRFBUVoSgKf//739FqtcybN4+kpCQURSEgIID27dsTExODoihotVpeeeUVpk6dyocffkhoaKhhvU6nQ6/XG+JXFAWdTlfqfB6s1+l0TJw4EcCwz4NyHvz/zTff5KWXXqJVq1YsWrSIkSNH0qNHD0JCQpg0aRI7d+5k/vz5jBo1ivHjxxMUFFTi3HU6nSE+gOeff55r164RFhaGoig4Ojryn//8B3Pz0iNntFotjRs3Zvz48aXKmjZtGtOnTycoKAgLCws++ugjtFotU6dOZfLkyfzhD3/A19cXV1dXdDod7u7uzJw5k5EjR6IoChqNhgULFtCsWbNyr1NRUVGpn/XmzZuX+blQKb/tUy/DpUuX2LRpE3PmzAFg8+bNAAwZMqTUvkuXLuWpp54qMTlTWdt79uzJhAkTWLZsGWq1ulQdD3Pz5s0K9zEW+SVbv0n71l/Stk+GwsJC9u/fj5+fH05OTiiKUuEjKsZsW71eT35+PtnZ2RQUFGBra4uLiws6nY6zZ8+iUqkwMTHBwsICS0tLHBwcaNSokVFi/b2vz6ex/lwao72d6d7cxtjhlKsuJq4P/OfoVQ7u3s57UZEE+HgaO5wnjvxeNp78/PwqdTI9iooSV/HkqmrblvV5Ky9xrfDBGHd3d27dukVKSgparZYjR47g6+tbqUByc3MNw5yys7P55ZdfaNGiBSqVis6dO3P06FEADhw4UOkyhRBCNDwpKSmsW7eOixcvcuvWLaDy8yoYi4mJCTY2NjRv3hx3d3dcXO73vKnVarp37063bt3w8fGhffv2uLm51ZmkNTG9kA3n0+jezLpOJ6113XDvpmgULQtWbSY3v9DY4QghxBOvwqHCarWaqKgoFixYgF6vJzg4mJYtW7Jhwwbc3d3x9fUlMTGRxYsXk5eXx8mTJ9m4cSMffvghN27cYNmyZZiYmKDX6xk8eLBhUqfRo0ezZMkSvv76a9q0aVNqFishhBAC4PLly2zbtg0LCwtGjhxZ5syNonrc0yl8HHOTRmZqhnZyNHY4TzRbKwteGjqQT1et473/7mD+hKHGDkkIIZ5oFQ4VrmtkqLAwFmnf+kvatu66du0a3377LS4uLgwZMqTMiSUeRtq2av57LpUN59N5uXtjOjWu2aGC1aEuDxV+YOGa7zlz5jQf/uVl/uDtYexwnhjys2s8MlRYPA6jDhUWQgghjMXV1ZUePXowcuTIKietomquZBQSfSGdp5pbPxFJ65Pi1WH9sLaxZfGmvdzT6Ss+QAgje8L6tMQTriqfN0lchRBC1Dm//PILBQUFaDQaAgICSrwzT1Q/rV7h45hbWJmaMLiDDBGuTrZWlkwa8xz6dgF8fT7d2OEIUaGK3ksqRHXRarVVehd5zb+kSQghhKiCCxcusGvXLnx9fenbt6+xw2kQvo1NJymziHHdXLAyq9nXYDREPTyac7kgjehzt/G20+LTpvx3BAthbBYWFhQWFlJUVFRjk+CZm5tTVFRUI2UL46ps2yqKYphZv7IkcRVCCFFnXLt2jT179tCqVSv69Olj7HAahOSsIr6+kEbXplZ0aSLDsWvKs+0d2L/9W/5+QcumBX/GzFT+BBN1k0qlwtLSskbrkGeY66+abFsZKiyEEKJOSE9PZ+vWrTg6OhIZGYlGI3/Y1zRFUfjs+G3M1SYM7uhk7HDqNSszNYODe3PjTiqfbvnJ2OEIIcQTRxJXIYQQdcLevXvRaDQMGTIEc3NzY4fTIOxPyiY2pYBnPB1oZC5DhGtaf7/OtHV3Z+MPP5F0S3qbhBCiKiRxFUIIUScMHDiQoUOHYmtra+xQGoTsIh3/OZVCa3tz/FrYGDucBkGlUjFpxNPogbdWfieztwohRBVI4iqEEMKokpOT0ev1WFtb07hxY2OH02B8dTqFvGIdwzo7YVJDE7CI0lo1diQsMIBLKbkcScowdjhCCPHEkMRVCCGE0SQnJ7Np0yaOHz9u7FAalPiUfPZcziKwtS3NG8mrhmrbCxF98AoayH/O3qVIK+92FUKIypDEVQghhFEUFhayc+dO7O3t6datm7HDaTC0eoVPj9/GwUJNuLu9scNpkMw0aoZ2duZORjafbDtq7HCEEOKJIImrEEIIozhw4AB5eXkMGDBAJmOqRTsu3eXXrGIGd3TCXCN/BhiLu6MFjbMS2LRlO/HXbhk7HCGEqPPkG0sIIUStS0xMJDY2lh49etCsWTNjh9NgZBVqWX8ujfbOFnRuXLPvaRQVe/nZYEw0Gt5ds00mahJCiApI4iqEEKLW2djY0L59e/z9/Y0dSoOy9mwqhVo9gzo4opIJmYyuiYMt/foGEJdwhe9jzhs7HCGEqNMkcRVCCFHrmjZtysCBA1Gr5d2hteVyRiF7ErMIcLOliY1MyFRXPBfmj4OjI0s37qSgqNjY4QghRJ0liasQQohac+PGDfbu3UtRUZGxQ2lQFEVh2c93sDYzIdzDztjhiN8wN9XwQmQ4hZZOfBeXZuxwhBCizpLEVQghRK3Q6XTs2bOHK1euyDDVWva/azlcTC1gQDsHLE2ll7uu+YNPewJC+rH5Ug4ZBVpjhyOEEHWSJK5CCCFqxalTp0hPTyckJAQzMxmqWlsKtXpWnkqhha0Zfi1sjB2OKMez7R3Jv5vGe1/vMXYoQghRJ0niKoQQosbl5eVx9OhR2rZti4eHh7HDaVC+iU0no0DL4I6OmEhPd53lbG2Ku5LC/gMHORp3xdjhCCFEnSOJqxBCiBp35MgRtFotffv2NXYoDUpK7j02x2XQrZk1bRwsjB2OqMDYgUFYWFnz0fod6PV6Y4cjhBB1iiSuQgghapyvry/h4eE4OjoaO5QGZd3ZVACe8XQwciSiMuysLOkf/Acu/3qD7w6fM3Y4QghRp0jiKoQQosY5ODjQuXNnY4fRoFzOKOTA1Wz+0MoWB0uNscMRlTS0ry8OTo588e1uiu7JRE1CCPGAJK5CCCFqzNWrV9m8eTN5eXnGDqVBURSF/5xKwdrMhJC28vqbJ4mpRs3IZ8K559CKXQl3jR2OEELUGZK4CiGEqBGKonDo0CHS0tKwsJDnK2vTyZt5XLiTTz93eyxN5av+SRPk40m3Hj3ZGHuX3GKdscMRQog6Qb7NhBBC1IjExETu3LmDv78/arW8O7S26PQKq06l4Gylwb9lI2OHIx6BSqUisoMDqclJfBi939jhCCFEnSCJqxBCiGqn1+s5fPgwjo6OdOrUydjhNCj7rmSRnF3MM54OaEzk9TdPqha25rhqU9ix7wCXb6UZOxwhhDA6SVyFEEJUu4sXL5Kenk7v3r0xMZGvmtpScE/PurOptLY3p0sTK2OHIx7TmGfDUIAPv95l7FCEEMLo5K8JIYQQ1a5NmzYEBQXh6elp7FAalC3xGWQW6hjY3gGVSnpbn3QtXRzo+dRTHD9zgZMJvxo7HCGEMCpJXIUQQlQ7S0tLnnrqKUmealFmoZZv49PxbmJFGweZDKu+ePHpvpiZm/PxRul1FUI0bJK4CiGEqDZ6vZ6tW7fy66/SO1Tboi+kc0+nMKCdg7FDEdXIsZEVYcFB3DRrTnxKvrHDEUIIo5HEVQghRLVJSEjg0qVLFBYWGjuUBiU17x47EzLxdbWhsY2pscMR1WxkiB9NWrfjq7OpKIpi7HCEEMIoJHEVQghRLRRF4dixYzg6OtKuXTtjh9OgfH0+DVAId7c3diiiBphrTAhp3YjjJ35m+/F4Y4cjhBBGIYmrEEKIanHlyhVSU1Pp0aOHPNtai65nF7Hvcha9W9riYKkxdjiihvRs2QhuJ/B59A/o9HpjhyOEELVOElchhBCPTVEUjh8/jq2tLR07djR2OA3KurNpmKlVhLrbGTsUUYPMTTVEhvXldkoq6348aexwhBCi1kniKoQQolp07tyZP/zhD6jVamOH0mBczijkyK85BLa2xcZMrnt9N8DfBxcXF9Zs/ZHCYq2xwxFCiFoliasQQojHplKp8Pb2pkOHDsYOpUFZeyYVK1MT+raW3taGQG1iwrCIIO5mZrLihxhjhyOEELWqUg/DnDlzhpUrV6LX6wkNDWXw4MEltsfFxbF69WquXbvG1KlT6dWrFwBXr15l+fLlFBQUYGJiwtChQ+nduzcAS5cuJS4uDisrKwAmT55M69atq/HUhBBC1IaMjAyuXbuGl5cXpqYyo21tib2Tz6lbeQxs74ClqdyHbigCu3Zgz6kuHLyjYoJWj7lG2l4I0TBUmLjq9XpWrFjBm2++iZOTE7NmzcLX15cWLVoY9nF2dua1115j69atJY41MzPjT3/6E82aNSMjI4OZM2fi4+ODtbU1AGPGjDEkuUIIIZ5Mp06d4vz583h6ekriWksUReGrs6nYmasJcGtk7HBELTIxMWH8kAF8duIOPyRkMqijo7FDEkKIWlHhbbrExESaNm1KkyZN0Gg09O7dmxMnTpTYp3HjxrRq1arULJLNmzenWbNmADg6OmJnZ0d2dnY1hi+EEMKYCgoKiI2NpVOnToabkqLmnbyZx8XUAvq522Oqlh63hsbDyZI2NgrLvvmBrNx8Y4cjhBC1osJvu4yMDJycnAzLTk5OZGRkVLmixMREtFotTZo0Maxbv349b7zxBqtWreLevXtVLlMIIYRxnT9/Hq1WS/fu3Y0dSoOhKAprz6biZKXBr4WNscMRRtK7iZo7Cef5fOv/jB2KEELUilp54dvdu3f55JNPmDx5MiYm93PlF154AXt7e7RaLV988QVbtmxh+PDhpY7du3cve/fuBeC9997D2dm5NkJ+JBqNpk7HJx6PtG/9JW37aHQ6Hb/88gudO3eus6/AqY9te/ByOkl3i4jq6YaTo4OxwzEqtUaNvb29scMwipAe9mze245tPx1l5kuDcLKtXzcx6uPPrvg/0r71V022bYWJq6OjI+np6Ybl9PR0HB0r/zxFfn4+7733HqNGjcLT09Ow3sHh/petqakpwcHBpZ6PfSAsLIywsDDDclpaWqXrrm3Ozs51Oj7xeKR96y9p20eTnZ2NRqPB09Ozzl6/+ta2ekXhi0NXcbHS0MFORWZmprFDMip7e/sGfQ1G9AvknU9XsGDlZmaPHmDscKpVffvZFSVJ+9Zf1dG2zZs3L3N9hUOF3d3duXXrFikpKWi1Wo4cOYKvr2+lKtVqtSxevJjAwMBSkzDdvXsXuD/k6cSJE7Rs2bJSZQohhKgbbG1tGTVqFG3atDF2KA3GseRcrmYW0c/DHrWJquIDRL3mbMdwbQAAIABJREFU7dECDw93dvx0lLTsPGOHI4QQNarCHle1Wk1UVBQLFixAr9cTHBxMy5Yt2bBhA+7u7vj6+pKYmMjixYvJy8vj5MmTbPx/7N13fFT3ne//15kZ9YokVOggiqmmCBAYkCi2cYl7qpON7Xjj/JIbr3P3ZnNT9vHLb+91rnPddo3tFNshdoIdA264YZAB0zHFgEE00UGAUO/SzJzv7w/W2mCKJNDMGc28n//NOd9zznv8NTPz0fec73fhQp566inWr1/Pnj17qKurY9WqVcB/LXvzzDPPtE3U1L9/f77//e8H9I2KiEjXqaurw+12Ex8ff8HEfBIYtjG8uvMsmQkexuVoIiw55xs3FfLU4hUsKT7LA/n6/0JEwpdljDFOh+iM0tJSpyNckm57CG/q3/Clvu28ZcuWceDAAR566CE8nqBMl3BFwqlv1x2t5f+uLeXeMRmM7xVezzNeqUi/VfgLf/6sjJKKZv54Ry7JMW6n43SJcPq3KxdS/4YvR28VFhER+XvNzc3s2bOHIUOGhHTRGk78tuG1neVkJUQxVqOt8iU3Dk6lrqqcP3z4qdNRREQCRoWriIh0SnFxMT6fj2uvvdbpKBFj3bE6jte2cv3gVFy6NVu+JCcpmvSqA7z5/oecrqpzOo6ISECocBURkQ4zxrBjxw5ycnLOW5dbAsdvG/72eTnZiVFcmx3vdBwJUd+8qQCf18e8t1Y6HUVEJCBUuIqISIedPXuWqqoqxowZ43SUiLH2aC0na1u5QaOtchlD+mYzcvgwVqzfQmlFrdNxRES6nApXERHpsMzMTO6//36GDRvmdJSI8MVoa05SFKOzNNoql/fNm2bi8/l4/u1VTkcREelyKlxFRKRDvpiEvkePHkRFRTmcJjKsPlJLaZ2XG3I12irty+2dyejRo9l8qomaZp/TcUREupQKVxER6ZBPP/2Ud955B7/f73SUiOC3Da/vKqdXUhSjNNoqHfS9e24hbuC1vLOn0ukoIiJdSoWriIi0y7ZtduzYgdfrxe0Oj3UiQ93qI7WcqvPq2VbplKzEaMZkxbN43S7NMCwiYUWFq4iItOvYsWPU1dUxevRop6NEBL9tWLSr4txoa6ZGW6VzpmZanNpSxO/fWeV0FBGRLqPCVURE2vX5558TFxdHbm6u01EiwvpjdZysa2VObiqWRlulkwb37smwYUNZtm4zZ6vrnY4jItIlVLiKiMhlNTU1UVJSwvDhw/F4PE7HCXu2MSzcVU5WgmYSliv39RsLaG318od3VzsdRUSkS6hwFRGRdk2ZMkVrtwbJphP1HKtpZXZuip5tlSs2vH8OQ4cMYemaTVTWNTgdR0TkqqlwFRGRy4qLiyM/P5/09HSno4Q9Ywyvf15ORryHsdkJTseRbu5rN87Aa3lYvPmI01FERK6aClcREbmkiooK9u3bpyVwgmRraQOHq1qYPSgFt0ujrXJ1Rg3qw8w7v82qiigavfo3LCLdmwpXERG5pG3btrF06VK8Xq/TUcKeMYa/fV5OWpyHCb0SnY4jYeKGoWnUN3v526aDTkcREbkqKlxFROSivF4ve/fuZejQocTGxjodJ+xtP93IgYpmZmm0VbpQv5QYEk9sZf6Cv1HT0Ox0HBGRK6bCVURELmr//v20trYyatQop6OEvS+ebU2NdTOxt0ZbpWvdXjCRpsYmXvxgndNRRESumApXERG5qF27dpGamkqfPn2cjhL2dpU1sudsEzMHpuDRaKt0sQnDBtK/fz+WrFxPfVOL03FERK6IClcREbmA1+ulsbGRESNGYGlJloB7/fMKkmPcTO6j0VYJjLtvKKChoZEXP9zgdBQRkSuiwlVERC4QFRXFfffdx6RJk5yOEvb2nG3k8zONFA5IJsqtr2UJjEnDB9G3bx8+3LyHFp/tdBwRkU7TN6SIiJzHGIPP58OyLNxut9Nxwt7CzytIiHaR3zfJ6SgS5u772p3EjCyk6GCN01FERDpNhauIiJzn9OnT/P73v+fkyZNORwl7Byqa2HaqgYIBycR49JUsgTWiVw8GpcWycHspza1a4kpEuhd9S4qIyHmKi4vx+Xykp6c7HSXsLfy8gvgoF9f1S3Y6ikQAy7KY2tNiz7LXeWX5ZqfjiIh0igpXERFp4/f72bdvH7m5uVq7NcAOVTbz6cl6pvdPJlajrRIkYwdkkpmRxsJlq2n1+Z2OIyLSYfqmFBGRNocPH6apqYkRI0Y4HSXsLdxVTqzHYnp/PdsqweNyubht9nSqq2tYsGKL03FERDpMhauIiLQpLi4mPj6eAQMGOB0lrB2rbmHD8Xqm9U8mLkoTYElwFY4bTs+MDF5f+glejbqKSDehwlVERNpMnDiR2bNnazbhAFu0q4IYt8WM/nq2VYLP5XJx66xplFdUsWTzfqfjiIh0iMfpACIiEjpycnLIyclxOkZYK61tZc3RWgoGJpMQrT8QiDNmTxzFpxUWKytiuNMYXJbldCQRkcvSiKuIiADw6aefcvr0aadjhL03iivwuCwKB6Q4HUUimMft5pZxgzhe08r6o1rXVURCnwpXERGhqqqKNWvWcOzYMaejhLWyei8rD9UwuU8iSTEabRVnjc1JwHV8J4//4S/Ytu10HBGRy1LhKiIi7NmzB8uyuOaaa5yOEtbeLK4AoHCgRlvFeS7LYsqQbE6WlvLepj1OxxERuSwVriIiEc4Yw969e+nTpw/JyZosKFAqm3wUHawhr3ciPeI0xYSEhluvG09yUhLz31uhUVcRCWkqXEVEIlxZWRlVVVUabQ2wd/ZU4rMNswZptFVCR3SUhzkzpnD8RCkfbdMMwyISulS4iohEuJqaGhISEhgyZIjTUcJWbbOPDw9UMb5XAhnxUU7HETnP7dPzSEyIZ/57qzDGOB1HROSidK+SiEiEGzp0KEOGDMHSchgB8+6+Klp8Gm2V0BQTHcVtt97Kx6V+dpU1MjorwelIIiIX6FDhun37dubPn49t28yePZs77rjjvP3FxcW8/PLLHD16lEceeYT8/Py2fatWreLNN98E4K677qKwsBCAQ4cO8dxzz9Ha2sq4ceO4//779aNJRCTIvF4vHo9Hn78B1NDq5719VYzJiic7MdrpOCIXdVPeMLatPsnCXRUqXEUkJLV7q7Bt27z00kv84he/4Omnn2bdunWcOHHivDYZGRn88Ic/ZNq0aedtr6+vZ/HixfzmN7/hN7/5DYsXL6a+vh6AF154gYceeohnnnmG06dPs3379i58WyIi0hFr1qzh5Zdf1qQsAfTB/ioavTazczXaKqEryu1icobFyvffYe3uI07HERG5QLuFa0lJCdnZ2WRlZeHxeJg6dSqbN28+r01mZib9+/e/4C/227dvZ8yYMSQmJpKYmMiYMWPYvn07VVVVNDU1MXToUCzLYsaMGRecU0REAsu2bfbv309aWhoul6Y8CIRmn807eyoZ3jOOPskxTscRuaxpuenQUMEf31nhdBQRkQu0e6twZWUl6enpba/T09M5cOBAh07+5WPT0tKorKy86DkrKys7k1tERK7S8ePHaWhoCPpswsbnhdLjUFWBaW7EioqG1DTI6YsVFx/ULIH20YFq6lpt5ujZVukGkuLjmD45j+WfrGPT/uNMHtrX6UgiIm1CfnKmoqIiioqKAHjsscfIyMhwONGleTyekM4nV0f9G74itW83bNhAamoqEydOJCoqsDPd2k2NNK9ZTvO6j/Hu2QFeb9u+tjlMLQtP7jXETp5B7OxbcfdIv+i5OsPJvm3x2byz9yDXZCYydmC2IxnCndvjJjU11ekYYeWBu+ayZuMWXnpvNbf85p8cyxGpn8uRQv0bvgLZt+0WrmlpaVRUVLS9rqioIC0trUMnT0tLo7i4uO11ZWUlI0aM6NQ558yZw5w5c9pel5eXd+jaTsjIyAjpfHJ11L/hKxL71u/3s23bNgYOHEhNTU3ArmOamzAfvYVZ8S40NkBaBoyfitW7P6T0gJhY8PmgtgpOn8R3cC/1C/5A/WsvYk2ZiXXbN7HSel7x9Z3s2w/3V1HR6OWrI9Oorq52JEO4S01N1X/bLmYB+Xnj+GT9Joo+3cnYQb0cyRGJn8uRRP0bvrqib3v1uvjnTruFa25uLqdOnaKsrIy0tDTWr1/Pww8/3KGLjh07ltdee61tQqYdO3bwrW99i8TEROLi4ti/fz9Dhgxh9erVzJ07txNvR0RErtYNN9xAUlJSwM5vtq3HfvWPUFMJw0ZjTZkJvS+cDwGA7N4wdBTWjBsxlWcxm9dgNq7CfPoJ1u33Yl1/O5bLHbCsXc1nG94orqB/agyD02KdjiPSKV+7fiq7Krx8cKSFsYOcTiMick67havb7eaBBx7g0UcfxbZtZs6cSd++fXn99dfJzc0lLy+PkpISnnjiCRoaGti6dSsLFy7kqaeeIjExkbvvvpuf//znANxzzz0kJiYC8OCDD/L888/T2trK2LFjGTduXGDfqYiItHG73QwZMiQg5zYtLZhXf49Z/zFk98G689tYfQZ0+HgrrSfWjXdhJhdilr+NWfxnzLYNuB76GVZa97i17JPDNZxt8PG98Wlaaki6nbTkJG4uvI4Vh2o4XtNC3xRNLCYizrOMMab9ZqGjtLTU6QiXpNsewpv6N3xFWt96vV62bt3KyJEju3zE1VRXYM/733D8EFw3B2v6DVjuKx8pNcbA7m2YDxZDdAyuH/4Ca8iIDh/vRN/6bcOP3juEC4ufTM1R4RpAulU4cOpb/fzy9bVck+jjyQdvCfr1I+1zOdKof8NXIG8V1voHIiIR5tChQ6xbt46qqqouPa85cRj7N/8DTh3H+uoDuApvuqqiFcCyLKxRE7AeeASiY7Cf/lfMztBePm3dsTpO1XmZnZuiolW6rcRoNwPdNaxZt54DpSowRMR5KlxFRCLM3r17SUhIoE+fPl12TnP8MPYTvwSfD+u7P8YaOrLLzg1gZWRhfffHkJGF/dyj2BtXdun5u4ptDIt2lZOVGMXorPBa2kciz9dvnI4FPPdWaP57E5HIosJVRCSCtLS0cPjwYYYOHYrL1TVfAebkUeyn/hXcHqzv/Agru3eXnPfLrIRErG//EPoOwvzp3zHb1gfkOldj84l6jtW0MmtgCi6Ntko3l5Peg3FjRrFh6w4On+naOzRERDpLhauISAQpKSnB7/czbNiwLjmfOXsa+8lfgWVhffv/CfjkSVZMLNbXvwe9+2P/8QnMnh0BvV5nGGNYuKuC9DgP43ISnI4j0iW+dsN0jG3z/DurnI4iIhFOhauISASpra2lR48el5z4oDNMYwP2vP8F3lase39wVeutdoYVHXOueE3LwH7uUczxw0G5bnu2n26kpLKZWYNScLs02irhoV92BhMnTWJXfRwVjV6n44hIBFPhKiISQaZMmcJ999131ZMGGb8f+w+/hTMnse65Dysjq4sSdowVl4D1rYcgJvZc8VpXG9TrX8zCz8tJiXWT1zvR6SgiXeofvjKb2OwBvFVc6XQUEYlgKlxFRCKEz+cD6JJnW80bf4bi7Vg33YM1IDDrwbbHSkrBuuc+qKnE/sNvMX6/IzkAdp9ppPhsEzMHpODRaKuEmfT4KMZmxrBo+VpOV9U7HUdEIpQKVxGRCPHmm2/y4YcfXvV5zPaNmOXvQN51WOPyuyDZlbN69cO66R7Y9znmzZcdy7FwVzmJ0S4m99Voq4SnSRlQuXcLv1uy2ukoIhKhVLiKiESA+vp6Tpw4QUpKylWdx5Sfwf7Tf0BOH6w5t3dRuqtjXTsJJkzFLHsbs2tr0K9/oKKJ7acbKRiQQrRbX6sSnob0zeGaoYMpWvcp5TUNTscRkQikb1gRkQiwb98+jDFcc801V3wO4/dj//EJsP1Yd30Xy+PpwoRXx5pzG/TMwZ7/H5ja6qBee9GuCuKjXEztlxTU64oE21dvmEFLSwu/f3+d01FEJAKpcBURiQD79u2jZ8+epKWlXfE5zNI34PA+rJu/itUjvQvTXT0rKhrrznuhsf5c8WpMUK57pKqZTSfqmdYviViPvlIlvI0Y1JchgwaydPVGqhqanY4jIhFG37IiImGuurqaU6dOXd1o6/HDmHdfgxFjsUaO68J0XcfK7IU161bYtRWzdnlQrrlodwUxbovp/ZODcj0Rp91zw3RMXDLv7Cx1OoqIRBgVriIiYS42NpbZs2dfceFqfF7sPz0NsfFYc+/u4nRdbOI06JeLWfQnTHVFQC91oraFdUfrmNovifhod0CvJRIqrh06kIKb72T5CR+NXudm8haRyKPCVUQkzMXGxjJ27FiSk69sVNAsfRNOHMG65atY8QldnK5rWZYL69avgdeL/dffBfSW4Td2V+BxWRQMuLoJr0S6mzm5qdTU1rNg7R6no4hIBFHhKiISxqqrq9m5cyctLS1XdLw5exrzwUIYfi3W0FFdnC4wrLSeWAVzYcentKz7OCDXOFXXyqrDtUzpm0RSjEZbJbL0S4kh9sgG/rp4CY0tXqfjiEiEUOEqIhLGiouLKSoqwuvt/I9LYwz2q38Ay4V1fWgsfdNhkwsgpy91L/07pqmxy0+/eHcFLsuicKCebZXIdNvsaTTU1zF/2adORxGRCKHCVUQkTBlj2LdvH3369CExMbHzJ9i2AXZtxSqYi5Wc2vUBA8hyubDm3oVdXYl57/UuPfeZ+lZWHqphcp9EUmJDZ0kgkWCaMmoovXOyeGP5Gpq9PqfjiEgEUOEqIhKmysrKqKysvKJJmUxzI/bf/ghZvc9NeNQNWb37EzNpOqZoCebU8S477xu7K7EsmDVIz7ZK5LIsi9tmT6O2poZXirY4HUdEIoAKVxGRMLVv3z5cLhdDhgzp9LHm3dehpgrrpruxXN33Gc64m++G6Gjs1/7YJRM1nW3wUnSomkm9k0jVaKtEuOljR5CTlckHO47h9Qdn7WQRiVwqXEVEwlR1dTX9+/cnLi6uU8eZslOYj5fAmIlYfQYEJlyQuBKTsQpvgj07YNv6qz7fG7srwGi0VQTA5XLx4HfvxfQdw8rDNU7HEZEwp8JVRCRM3Xbbbdx2222dPs5+42Vwu88VfOFg/BTI6oW98CWMt/WKT1Pe6GX5wRom9k6kR5xGW0UARmQl0TclmgUbD+Lz207HEZEwpsJVRCQM2fa5H5AeT+cKLHOgGLatx8qfiZUUHqOKlst9blbkynLMiveu+Dxv7q7ANobZGm0VaWNZFmPj6ti9fDGvr9npdBwRCWMqXEVEwoxt27z00kts3ry5U8cZ28Ze+BIkpUB+YWDCOcQaMARyr8G8vwjTUNfp4ysavSwrqSGvVyJp8VFdH1CkG5s2Kpe01BQWfLBKo64iEjAqXEVEwszx48epra0lJaVzI4Nm8xo4cgCr8Gas6JgApXOONftWaG7EfLCo08e+tacSvzHMztVoq8iXeTwebiq8jrIzZ3hrQ7HTcUQkTKlwFREJM/v27SMqKoqBAwd2+Bjj9WLefAWye8OYCQFM5xwrsxeMmYj5+D1M+ZkOH1fV5GPpgWrG5ySQodFWkYuaO2UsqSlJvPL+SuwumMFbROTLVLiKiIQRv9/PgQMHGDx4MFFRHS+yzNplUHkWa+atWFb4fjVYBXPBsjBv/7XDx7y9pxKf3zAnNzWAyUS6t6ioKG6Yns+ZsjJW7Cl1Oo6IhKHw/XUiIhKBjhw5QnNzM8OGDevwMaalBfPe69B3EAwaGsB0zrOSU2HSDMymTzDHD7fbvrrZx4f7qxjfK4GeCRptFbmcW67LY+QNX+f9I61dsm6yiMjfU+EqIhJGMjIymDp1KgMGDOjwMeaTD6C2GqvwJizLCly4EGFNmQkxcdhLXm237Tt7Kmn1ayZhkY6IiYlmzjVZHKxsZsPhSqfjiEiYUeEqIhJGUlJSmDJlCm63u0PtTXMj5sPFMGgYVv/cAKcLDVZcPNbkGbB9E+bowUu2q2328f7+KsbmJJCVGB3EhCLdV16vRFp2r+TpV97QqKuIdCkVriIiYaK0tJSSkpK2NVw7whS9C/V15579jCSTZkBs/GVHXd/ZW0WrzzBHo60iHeZ2WUwdMZDDhw/z8c5DTscRkTCiwlVEJExs2bKFoqKiDrc3DfWYZW/B0JFYvfsHMFnosWLjsPILYedmzOEDF+yvbfbx7r5KxmTHk52k0VaRzrhz1hQS4mJ58e2PnY4iImFEhauISBhoaWnh0KFDDB06FJerYx/tZtlb0NSIVXBTgNOFqInTIC7hoqOubxZX0uoz3DhYMwmLdFZ8bCyF1+VTcugwK3dc+nZ8EZHOUOEqIhIGSkpK8Pv9HZ5N2NTXYoqWwIixWFm9ApwuNFkxsVhTCmHXVszBvW3bq5vOPds6Ts+2ilyxu2dNJj4+jvkfrnM6ioiECRWuIiJhYO/evSQnJ9OrV8eKUPPxe9DagjXt+gAnC3F50yA+8bxR1zeKK/DZhhs02ipyxeJjY7njrrup7jWBfeVNTscRkTCgwlVEpJvz+/1UV1czbNiwDi1nY5oaMR+/C8NGYWXmBCFh6LKiY7CmzoLi7ZiSPVQ0evnwQDUTeiVq3VaRq3TjmIEkxkWxYPtpp6OISBhQ4Soi0s253W4eeOABpkyZ0qH2ZtUH0NSAdV2Ej7Z+YfyUc6OuHyzijd0V+G3D9bmaSVjkasV4XEzs4eejRQso+myf03FEpJvzdKTR9u3bmT9/PrZtM3v2bO64447z9nu9Xp599lkOHTpEUlISjzzyCJmZmaxZs4YlS5a0tTt27Bi//e1vGTBgAL/+9a+pqqoiOvrc80O/+tWvSEnRDwURkc4yxmBZFlFR7Y8QmpYWzLK3z63b2qtvENKFPis6BiZN5+yG9XyUUc2k3omkx2u0VaQrzBnRm7c9Fn94s4hZ1w7p8ORxIiJf1m7hats2L730Er/61a9IT0/n5z//OXl5efTp06etzYoVK0hISGDevHmsW7eOBQsW8JOf/ITp06czffp04FzR+vjjjzNgwIC24x5++GFycyNjwXsRkUCoq6tjwYIF3HDDDQwaNKjd9mbtMqivxbrrH4KQrhvJm8bi04kY22Z2rp5tFekqcTHRXD99Km9/uIylW/Zy86QRTkcSkW6q3T97lZSUkJ2dTVZWFh6Ph6lTp7J58+bz2mzZsoXCwkIA8vPz2bVrF8aY89qsXbuWqVOndl1yERFh7969NDQ00KNHj3bbGq8Xs/QN6JeL1a/9IjeSlBHLx1kTmFO6icz6MqfjiISVOwvySExK4oW3irBt2+k4ItJNtTviWllZSXp6etvr9PR0Dhw4cMk2breb+Ph46urqSE5ObmuzYcMGfvrTn5533PPPP4/L5WLy5MncfffdF51UpKioiKKiIgAee+wxMjIyOvH2gsvj8YR0Prk66t/w1Z379uTJk+Tm5jJkyJB22zYue4e66kqSvnofUamRMarodrtJ7cB7/eOOOlwuP3efXE3ShkYa/+Gfg5BOrpbb07H+FefdOXcOf1n0Fqv3lnLPjLHttu/On8vSPvVv+Apk33boGderdeDAAaKjo+nXr1/btocffpi0tDSampp48sknWb16NQUFBRccO2fOHObMmdP2ury8PBiRr0hGRkZI55Oro/4NX921bysrKzl48CAFBQXt5jd+P/aiP0NOX+oze2NVVwcnpMNSU1Opbue9nmoyLD3uZ2qGRdTQ4URvLOLMjNvxp+pHVajrSP9KaLg+bwQbTjby1uFWZgw/i6udGdC76+eydIz6N3x1Rd9eamm/dm8VTktLo6Kiou11RUUFaWlpl2zj9/tpbGwkKSmpbf+6deu47rrrLjgGIC4ujmnTplFSUtLBtyIiIgD79u3DsiyGDRvWbluzeTWUn8Gadn2HlsyJJAuP27gtmJVlUT2+EDCkrF7S3mEi0gnRUR7unDqKI9WtrD9W63QcEemG2i1cc3NzOXXqFGVlZfh8PtavX09eXt55bSZMmMCqVasA2LhxIyNHjmz7YWTbNhs2bDivcPX7/dTWnvvQ8vl8bN26lb59NbuliEhn9O/fn+nTp5/3h8KLMcZgPlwMPXNgqCZG+XsnGg2rygxTMiySoyz8yT2ov2YCSZ9+jKu+xul4ImFlXE4CMeUlPPa7V/D59ayriHROu7cKf7E+4KOPPopt28ycOZO+ffvy+uuvk5ubS15eHrNmzeLZZ5/lxz/+MYmJiTzyyCNtx+/Zs4eMjAyysrLatnm9Xh599FH8fj+2bTN69OjzbgcWEZH29erV65K305xn9zYoPY5127ewLC1F8fcWHLWJdp0bbf1Cdd4sEou3kLz+Q6pv+IaD6UTCi8uymJ6bzsvb1vHaqm18Z3Ze+weJiPwny3x5+t8QV1pa6nSES9L9+uFN/Ru+umPfHjx4kISEBLKzs9tt63/yV3DyKNZ/+yWWOyhTG4SMyz0DeaDO8NMdfq7Psrgx5/yCPuvd+cSWHuH4z5/HRMcGI6pcAT3j2v34/H4e+e3vsP027z7+34mJuvhnUnf8XJaOU/+GL0efcRURkdBijOHjjz9m/fr17bc9dgj27sSaOD3iitb2/PWInwQ3FGRe+Mxv9YRC3E31JG5ZFfxgImHM43Zz140zqayq4uVlm5yOIyLdiApXEZFuprS0lLq6Oq655pp229rL3oboGBg/JQjJuo8d1TY7amB2lkWs+8LCtaXXQJpzBpCy5l3w+x1IKBK+CscNJycnh799uIqmVq/TcUSkm1DhKiLSzezduxe3201ubu5l25nKcti8BsZOxoqNC1K60GeM4a9HbFKjYErGpWdYrp5QSFTVWRJ2a1RIpCu5XC6+edtcYkYWsuxgndNxRKSbUOEqItKN2LbN/v37yc3NJSYm5rJtzYp3wdhYk2YEKV33sKnCcKAebsi2iHJdunBtHDSS1tSeJH+yBLrXdBAiIW/SsH6MGtibRbvKqW/xOR1HRLoBFa4iIt1IZWUlPp+v3dt5rGnUAAAgAElEQVSETVMj5pOlMPxarNS0y7aNJH5j+OtRm8wYmJDWznq2Lhc14wuIPXmI2EPFwQkoEkFuHpLKsW2reeJvHzkdRUS6ARWuIiLdSEZGBj/4wQ8YNGjQZduZtcuhuQkrvzA4wbqJVWWGE00wN8eF22qncAXqR+Thj0skZfWSIKQTiSx9U2MZkBLFR2s2cvysZocWkctT4Soi0k0YYzDGEBUVhdvtvnQ7vx9T9A70G4TVq18QE4Y2r2147ZhN3zgYndKxY4wnipqx04jf9xlRZ44HNqBIBPrOLbPw+/08vWi501FEJMSpcBUR6Sb27dvHyy+/TG1t7WXbma3roLJco61fsvSUobwFburlwurAaOsXasdMxfZEkbL63QCmE4lMA3v1ZPzYa1m7+TP2HCtzOo6IhDAVriIi3cSePXtoaWkhKSnpkm2MMZiP3oK0njBkRBDThbYmn2HRcZvBiTAksXPH2nEJ1I2cROL2NbhrKgMTUCSCffvmAtxuN08vWuZ0FBEJYSpcRUS6gcbGRo4cOcLw4cMvP1q4fzccO4g1uQDL0kf8F5aUGmp9cFNO50Zbv1AzvgBsm+T1HwQgnUhky05LZe7cmylNHc6Biian44hIiNKvGhGRbmDv3r3Yts2IEZcfRbWXvQXxCTAmL0jJQl9Vq+HNEzajU6B/QueLVgBfSjoNg8eQvHE5VnNjFycUkbunX0tKchJ/3laG0fJTInIRKlxFRLqBPXv20LNnTzIyMi7Zxpw+ATs3w4TrsKKig5gutP3tmI3Xhpt7Xd1XXs2EQlwtTSRtXtFFyUTkC7EeF9dle1j1/lss3bLX6TgiEoJUuIqIhDhjDOPHj2fKlCmXb7f8HfBEYeVdF6Rkoe9InY9lpw1TMyx6xlzZaOsXWrL70dR7EMnrPgC/v4sSisgXCgb3JMbfxPOLluK3bafjiEiIUeEqIhLiLMti+PDhDBky5JJtTG01ZsMKGD0BK+HSkzdFmj/saSDGBddnX13R+oWa8QVEVZeTsGtjl5xPRP5LTHQUX5lTyKkzZbz43hqn44hIiFHhKiISwowxfPbZZ9TX11++3aoPwOvVEjh/5/Nqmw1lXmZlWSR4uqZwbRw0gtYePc8tjaPn8ES63M1TryUzK4vf/20JDc0tTscRkRCiwlVEJISVlpayYsUKjh49esk2prUFs/J9GDISKz0ziOlCl20M8w/bpMVYTO/ZNUUrAJaLmnEFxJw8ROzhPV13XhEBwO1y8c1br6eqppbnl6x1Oo6IhBAVriIiIay4uBiPx3P524Q3rIT6Oo22/p01Zw2HGuD2/nFEubqwcAXqR+Thj0sgZfWSLj2viJwzZWQu+bNvYrMvh+pmn9NxRCREqHAVEQlRfr+f/fv3M3jwYKKjLz5LsLFtzLK3Iacv9BsU5IShqdU2/OWoTe84mNQzqsvPbzxR1Fx7HfF7txFVdqLLzy8i8IOvzMCLi1d3nHU6ioiECBWuIiIh6vDhwzQ3NzN8+PBLN9q5GcpKsfILsayuHVnsrt4rNZS3wK29XLgC9N+kdsx12J4okte+H5Dzi0S67ORYxiQ085dXXmbHoVKn44hICFDhKiISosrLy0lISGDAgAGXbGMvextSesDwMcELFsJqvYZFx22uSYYhSYEr5O34ROqH55G4bTWuuuqAXUckkt08ujcubxNPvPah01FEJASocBURCVH5+fl873vfw+W6+Ee1OXwADuzGmjQDy+UOcrrQ9Noxm2Y/3JoT+K+3mnEzsHw+kjd8FPBriUSi9JQkZk+fyp79JXywWZOhiUQ6Fa4iIiHI5zs3IUlU1KWf0TTL34aYOBg7OVixQtqRBsPSU4YpGRbZcYG/bdqblknjoBEkb/wIq1XLdogEwtfnTCUlJZnnFn6I1+d3Oo6IOEiFq4hIiDHG8Oqrr7Jy5cpLtyk/g9myFsbnY8XEBjFdaDLG8OIhP3FumJsdvGd9ayYU4m6sJ3HbJ0G7pkgkiYmO4q65szlTdpZXVmxzOo6IOEiFq4hIiDlz5gxnz54lLS3tkm3Mx++CZWFNnB7EZKFrY4VhVw3cmG0R7wle4drcayDNWf1IWfMe2HbQrisSSa6fOIqxhTexsiaFhlaNuopEKhWuIiIhZvfu3bjdboYNG3bR/aaxHrNmGYwcj5WcGuR0oafFb5h/2CY7FvIzgjyzsmVRM6GAqIrTxO/ZEtxri0QIl8vFvQXXUt9q8+pnp5yOIyIOUeEqIhJCfD4fe/fuZciQIcTGXvwWYLP6I2hpxsovCHK60LSk1FDWArf3duF2YEmghsGj8SankbL63aBfWyRS9E2JYXh0LS+9+CLbD510Oo6IOECFq4hICDl48CDNzc2MHDnyovuNz3vuNuGBQ7Gyegc5XeipaDEsPm4zKiWwy99clstNzbgZxB7dR8yx/c5kEIkAt+cNwu2yeOyVJdi6NV8k4qhwFREJIX369GHmzJn069fvovvN5rVQXanR1v/0yhEbv4Gv9HL266xu5CT8MXHnnnUVkYDISE5k7szplBw+xuK1O5yOIyJBpsJVRCSEJCQkMH78+Iuu3WqMwSx7C3pmw6BrHEgXWvbWGj45ayjItEiPcWi09T+Z6BhqR08hftcmPBWnHc0iEs7umTmZjIwM/rB4KfVNWoZKJJKocBURCRF79+5l9+7dGGMu3qB4O5w4gpVfiOXAs5yhxG8MLxzykxIFszJD479F7dhpYLlIXvuB01FEwlaUx823b7+Rmppa/lj0mdNxRCSIVLiKiIQAYwzr169n165dlyxK7Y/egqQUGDU+yOlCT9Fpw8F6uCXHIsYdGoWrPzGF+mvGk7RlJa7GeqfjiISt/JGDmX3HN1hf34PTda1OxxGRIFHhKiISAkpLS6mqqrr0pEzHDsGe7VgTp2G5PUFOF1qqWw2vHLXJTYRxPUKjaP1C9fgCXN4WkjYtczqKSFj72qRcLAt+t+ag01FEJEhUuIqIhIDPP/+cqKgohg4detH9ZtlbEB0D46cGOVnoeeWITZMP7urjCrlbpr0ZOTT2H0by+qXg8zodRyRspcR6yEts4IOFf+HdDZ87HUdEgkCFq4iIw5qamti7dy/Dhw8nOjr6gv2m8ixm8xoYl48VG+dAwtBRXGNYUWYozLTIig2tovULNRMK8dRVk7h9rdNRRMLarXlDSU3twbzX36ehWRM1iYQ7Fa4iIg6rr6+nR48ejBkz5qL7TdESAKxJM4IZK+T4bMPvD/rpEQWzs0KzaAVo6juEloxepKx+Fy410ZaIXLWYKA/fvvMmKqtqeHJhkdNxRCTAOvSg1Pbt25k/fz62bTN79mzuuOOO8/Z7vV6effZZDh06RFJSEo888giZmZmUlZXxk5/8hF69egEwZMgQvv/97wNw6NAhnnvuOVpbWxk3bhz3339/yN3yJSISDD179uQf/uEfLvoZaBrrMas/ghFjsVJ6OJAudLx/ynCsEe4b6AqZCZkuyrKomVBA5kevEbd/O03DxjmdSCRsXTdqMCtGjuD9VRu4q2A8o/rnOB1JRAKk3RFX27Z56aWX+MUvfsHTTz/NunXrOHHixHltVqxYQUJCAvPmzeOWW25hwYIFbfuys7N5/PHHefzxx9uKVoAXXniBhx56iGeeeYbTp0+zffv2LnxbIiLdQ319PS0tLZf8w51Z/RG0NGPlzwxystBS3mJ47ajNiGQYmex0mvbVDx2LLzGFlDXvOR1FJOz94O65RMXE8FzRZZYTE5Fur93CtaSkhOzsbLKysvB4PEydOpXNmzef12bLli0UFhYCkJ+fz65duy77wVFVVUVTUxNDhw7FsixmzJhxwTlFRCLB2rVr2+5o+TLj82KK3oWBQ7GyezuQLnT86ZCN38DtvUNvQqaLcnuoGTuNuJLPiT552Ok0ImGtZ2oSD37ve5TG9mHl4Vqn44hIgLR7q3BlZSXp6eltr9PT0zlw4MAl27jdbuLj46mrqwOgrKyMf/mXfyEuLo5vfOMbDB8+/KLnrKysvOj1i4qKKCo699zCY489RkZGRiffYvB4PJ6QzidXR/0bvpzq26amJo4fP864cePIzMy8cP+KD6itqSTxa/cRnZoa9HyhYt3pFtZX1HF7/1gGpMd26li320VSYlKAkl2eyb8e+9OPydj0EY0P/E9HMoQ7t8dNagT/2whnne3bm8em8HlFCc999BnTf3QDOWnd4NaMCKbfVOErkH0b0MUAe/TowfPPP09SUhKHDh3i8ccf58knn+zUOebMmcOcOXPaXpeXl3d1zC6TkZER0vnk6qh/w5dTfbtt2zZqa2sZMGDABdc3xmC/8Qpk5tCQ2ZvG6uqg5wsFjT7Dkzv95MTC1NRW6uo7t8RMUmISdfV1AUrXPs/ISaRs+YQzs+7Bn6ofaV0tNTWV6gj9txHurqRv5+TA//u3Zfx372nmPfyNACWTrqDfVOGrK/r2i/mRvqzdW4XT0tKoqKhoe11RUUFaWtol2/j9fhobG0lKSiIqKoqkpHN/6R40aBBZWVmcOnWqQ+cUEQlnxhh27tzZ9ijGBXZtg9JjWPmF3ePW2AB55YhNVSt8ta8Ldzf871A7bjoYQ/L6D52OIhL2hvTuyXWTJ7Fh6w6Wb93rdBwR6WLtFq65ubmcOnWKsrIyfD4f69evJy8v77w2EyZMYNWqVQBs3LiRkSNHYlkWtbW1bc9tnTlzhlOnTpGVlUWPHj2Ii4tj//79GGNYvXr1BecUEQlnZWVlVFRUXHIJHHvZW5CUAiMjd0ba4hrD0tOG6T0t+iV0v6IVwJecRsOQMSRvKsJqbnQ6jkjYu/8rM0lJTeWJv7xNfZPWdhUJJ+3eKux2u3nggQd49NFHsW2bmTNn0rdvX15//XVyc3PJy8tj1qxZPPvss/z4xz8mMTGRRx55BIDi4mIWLlyI2+3G5XLxj//4jyQmJgLw4IMP8vzzz9Pa2srYsWMZNy5yf5yJSOTJysriO9/5zkWf4TJHDsDenVizb8VyB/SJjpDVahueK/GTFg03ZnfPovUL1RMKSdy/naTNK6idfqvTcUTCWnxsDN+96xae+dMCHnt1Kf/7e7c7HUlEuohlutm84aWlpU5HuCTdrx/e1L/hK9T61v+7/wPF27F+/K9YMZ2bjChcvHrUz8Ljhn8c5GJY8pUXrk4/4/qFnEXPE1VfzfF/mQcR+seIQNAzruHravv2ydfe52BzLM/dP4sh6XFdmEy6Qqh970rXcfQZVxER6VqbNm1i6dKlF18Cp/QYbNsAedMitmg90mB444RhQg/rqorWUFIzoRBPTQUJO9c7HUUkIvzwnpvIHjiUZzacwuu/8LNWRLofFa4iIkHk9/v57LPPaGxsxOW68CPYLH0DoqKxJk13IJ3z/Mbw7AE/sW64rXd4FK0AjQOvoTU9m9SVb8FF/mAhIl0rLsrF3cPTKP58J48vXuF0HBHpAipcRUSCaP/+/TQ0NFz0uX5Tfgaz6RMYPwUrPtGBdM5784ShpB7u7G2R4AmfwhXLRdXE2USXnSS+eIvTaUQiwsisBPpY1byzdCXbSk44HUdErpIKVxGRIPrss8/o0aMHAwYMuGCf+egtsCysyQXBDxYCDtcb/nbM5tpUi7E9wu/rqWHotXhTMkhd+SZ0r+klRLqtH33tFmJiY/n1C4to8fqcjiMiVyH8fhmIiISoU6dOcerUKcaOHXvB2qympgqzdjmMmYiVfOFMw+HOaxv+fb+feDfc1SeMRlr/nstN9cRZxJw8RNyBHU6nEYkIaSlJfOO2myg9XcYTry9zOo6IXAUVriIiQZKYmEheXh4jR468YJ9Z/g74fVhTZjqQzHkLj9scbYR7+rrC6xbhL6kbPgFfUiopK950OopIxLhh0ihGjhjOOx+vZefRs07HEZErpMJVRCRIkpKSKCgoICYm5rztpqEes+oDGD4WK62nQ+mcs7/O8MZxQ14Pi5Ep4Vu0AuD2UD2hkLgje4k9VOx0GpGI8eOv30r//Ln8YWcdLT5NkCbSHalwFREJgn379nH06NGL7jMr34OWZqzrZgc5lfNa/Ib/2O8nOQpuD6NZhC+nblQ+vvhEUlZq1FUkWFKT4rmvcDQnalt5Yd0hp+OIyBVQ4SoiEmA+n48VK1awdevWC/aZ5kZM0bswZARW1sUX3A5nC47anGyCr/Z1ERfGtwj/PeOJomZcAfEHdhJ9vMTpOCIRY1hGHCM9lbw8/098tHWf03FEpJNUuIqIBFhxcTGNjY3k5eVdsM+s/AAa6rCmXe9AMmd9VmWzpNQwNcNiWHJkFK1fqB0zFX9s/Ll1XUUkaL4+bQSpKSn89s9vUlHb4HQcEekEFa4iIgFkjGHr1q1kZmbSt2/f8/c1N55bAid3OFbv/g4ldEZ1q+Hf99tkx8JXekVW0QpgYmKpvXYaCcWbiTp9zOk4IhEjITaWh755B/X1dfz8j29gtDSVSLehwlVEJIAOHjxIZWUleXl5Fy6B88Vo64wbHErnDGMM8w74afDBvf1dRLkir3AFqBk3HTs6lh5Fi5yOIhJRxg7pz5yC6Wz7fA8vL9vodBwR6SAVriIiAeT1esnJyWHo0KHnbT832vomDI680db3Txm2VsGtvSxy4iKzaAWwY+OpGTeDhF2biC494nQckYjynZtm0Ld/f17ffprS2lan44hIB6hwFREJoOHDh/PNb34Tt9t93naz4n1oqMeaHlmjrUcaDH8+bDMiGa7LiNyi9Qs142bgj4kjtWih01FEIkqUx83P//FeUgeO5PG1J/H6tUSOSKhT4SoiEiDHjx/Htu0LbxFubsQsewsGj4io0dYWv+GJvX7i3fC1fq4L/rtEIjs2jprxBSQUbyH6hJboEAmmtPgovjE6g+J9Jfy/r3zodBwRaYcKVxGRAKisrGTRokVs3rz5gn1to60R9mzri4dsTjTBN/q7SIyQpW86ombs9HMzDGvUVSToRmbGM9hTzbJVa3lz3edOxxGRy1DhKiISAJs2bcLtdjN69OjztpvGhv8abe3Vz6F0wffxGZvlZwyzMi2GJqlo/XsmJpbqCYUk7N1GzLEDTscRiTg/vGcumZkZPPXKmxw8VeF0HBG5BBWuIiJdrLq6mr179zJmzBji4+PP22eWv31utLXgRofSBd+RBsPvS2wGJ8KNOSpaL6b22mn44xI06irigJjoKH7y3a9i2zY/nbeA5lav05FE5CJUuIqIdLEvbg/Oy8s7b7uprcIsextGjMXK6XuxQ8NOg8/w2B4/ce5zS9+49VzrRZnoGKrzZhK/fwcxR/Y6HUck4gzMyeRbd97KsZOneHKJlsgRCUUqXEVEupBt25w+fZpRo0aRlJR03j7zwWLwebEK5jqULrjOrddqU9YM3x7gIilKRevl1I65Dl9CMmkf/hWMcTqOSMSZO3k0N9/5VTY1Z7D6SK3TcUTkS1S4ioh0IZfLxbe//W0KCgrO227Kz2BWfQDXTsJKz3QoXXAtKTVsrDDc3MtiUKKK1vaYqGiq8m8k9uh+4ou3OB1HJCLde91wctNiebpoD5+VnHA6joj8HRWuIiJdpKmpiebmZizLIjo6+rx9ZsmrYLkiZt3W3TWGlw/bjEqBgp4qWjuqbuREWntk0uOjV8HvdzqOSMRxuyzuHZNO9c5V/Muzr1JZ2+h0JBH5TypcRUS6yIYNG5g/fz6tra3nbTcnj2I2roKJ07CSU50JF0RlzYbf7vGTHgNf13qtneNyU3ndzUSXnSRx6yqn04hEpNS4aB762leoqanmn+b9Db9tOx1JRFDhKiLSJWpra9mxYwe5ubkXjLbab/0FYmKxps5yKF3wNPsN/2ePn1Yb7h/oIs6torWzGnNH0Zzdnx7LF2K1tjgdRyQi5Q0fxG03zmbP/gP8r7984HQcEUGFq4hIl9i48dwslPn5+edtN/t3wY5PsfJnYsUlOBEtaIwxPHPA5kgD3DvARWasitYrYllUTr8VT10Vyes/dDqNSMT62uwpjLt2DO+tWMcb63Y6HUck4qlwFRG5SlVVVezevZtrr72W5OTktu3GtrFffxGSU2HyDAcTBsfiE4b15YabcyyGJ6tovRrNvQfRMGgkKavextVQ53QckYjkcrn4p2/eyrAJU3jtqIcDFU1ORxKJaCpcRUSuUklJCS6Xi8mTJ5+33WxYCccOYc26FSsq+hJHh4eNFTYLjtqM72FRmKmitStUXnczrpYmUj9e5HQUkYgVGx3F/7hnNkmx0fx/y0o4cqbK6UgiEUuFq4jIVZo4cSL3338/CQn/dSuwaW7CvPky9O4PI8c5mC7w9tcZntpn0y8evtrX0mRMXcSbnk3t6Ckkb1xG1OljTscRiVhJMW4eGN+TYxs/4sdP/pnqBo28ijhBhauIyFWor68HOO8WYQDz4RtQW411wx1hXcidaTY8WuwnyXNuMqYoV/i+VydUTbkROyqGtPf+DMY4HUckYuUkxXDfV2Zzpuws//Qfr+H1abkqkWBT4SoicoWOHDnCCy+8wPHjx8/bbirKMMveglHjsXr3dyhd4NX7DP+2+9wMwt8b5CIpSkVrV7PjEqmaciPxJbuIL97sdByRiDZj3HBumzub3fsO8MsX38boj0kiQaXCVUTkCti2zSeffEJSUhI5OTnn7TNvvAwWWDNvcShd4Hltw2PFfs40w30DNYNwINWOmUprejZp772C5W1t/wARCZhvXH8dUyZNZMWGLTzz9hqn44hEFBWuIiJXoLi4mPLycqZPn47H42nbbvbvxmxeA/kzsVJ6OJgwcPzG8Mx+m1218LV+FrmJKloDyuWmvOB2oqrKSF77vtNpRCLef/vaXMZNnMzy6iRWHqpxOo5IxFDhKiLSSa2traxdu5acnByGDh3att34fNgLnoeUNKypsxxMGDjGGF44aLOm3HBLL4vxPfQ1EgzN/YbSkDuK1JVv4q6pdDqOSERzu1z889dv5JqcHjyz/gRFO484HUkkIugXh4hIJ508eZKmpiYKCgrOm3jJfLwESo9j3Xhn2C5/89oxm6WnDYWZFjMz9RUSTBUzbgO/n7T3X3E6ikjE87gs7huXiev4Dv71mT+xZtdBpyOJhD396hAR6aSBAwfy4IMP0rt377ZtpvIsZslrMHQU1tCRDqYLnHdP2iw8bpiUZnFLjm4PDjZfSjrVk+aQuHM9cfs+czqOSMSL9bj456/OIT4hnp/P+wvbDpxwOpJIWPO03wS2b9/O/PnzsW2b2bNnc8cdd5y33+v18uyzz3Lo0CGSkpJ45JFHyMzMZOfOnSxYsACfz4fH4+E73/kOo0aNAuDXv/41VVVVREefG5X41a9+RUpKShe/PRGRrmOMoaysjKysLJKSks7bZ//tBTA21g13XOLo7m3FGZuXDtuMToG7tVarY6onzCRx3zbS336Rkz95ChMd43QkkYiW2SOZXz70Hf7tuZf5yVPzee5nDzJqQE77B4pIp7U74mrbNi+99BK/+MUvePrpp1m3bh0nTpz/F6UVK1aQkJDAvHnzuOWWW1iwYAEASUlJ/OxnP+PJJ5/kRz/6EfPmzTvvuIcffpjHH3+cxx9/XEWriIS8gwcP8te//pWSkpLztpsdm+GzjVjTbsBKTXMoXeCsKrOZd8BmSCLc29+FW0WrczweymfdQ1TVWVI/Xux0GhEB+mWl8z8fuheD4Z/+/a+U1rY4HUkkLLVbuJaUlJCdnU1WVhYej4epU6eyefP5a8lt2bKFwsJCAPLz89m1axfGGAYOHEha2rkfcX379qW1tRWv19v170JEJMC8Xi8rV64kPT2dgQMHtm03jQ3Yf30OemZDfoGDCQPjkzKb/9hvk5sI9w9y4XGpaHVac59cakdOImXNe0SdPuZ0HBEBBvfO4qf/+G1SR07jXz8+Tlm9fu+KdLV2C9fKykrS09PbXqenp1NZWXnJNm63m/j4eOrq6s5rs2nTJgYNGkRUVFTbtueff56f/vSnLF68WIs4i0hI27RpE7W1tcyZMwe329223SyeDzVVWF/5Bpa7Q09fdBtrz9r8+36bQYlw/0AX0SpaQ0bltFuxY2LJWPw78PudjiMiwPD+vfjxnBE0tNr8txeXcrC03OlIImElKL+yjh8/zoIFC/jlL3/Ztu3hhx8mLS2NpqYmnnzySVavXk1BwYWjFUVFRRQVFQHw2GOPkZGREYzIV8Tj8YR0Prk66t/w1V7fnj59mt27d5Ofn8/YsWPbtrfs2Ez1mmXEFs4lfsSYYEQNmlWlLTy1v47ByW5+PDKRGHf3LFrdbhdJiUntN+xuEpNomPstkt96gazNy2mZ+w2nEznC7XGTmprqdAwJgO7at6mp8EOXi1/85jV+8NvdLPjNPzOiv555/TL9pgpfgezbdgvXtLQ0Kioq2l5XVFS03f775Tbp6en4/X4aGxvbJi6pqKjgiSee4Ec/+hHZ2dnnHQMQFxfHtGnTKCkpuWjhOmfOHObMmdP2urw8dP96lZGREdL55Oqof8NXe3175MgRoqOjGT9+fFs709yEPe9RSOtJy+RCWqurgxU34D4+Y/PsAZsBCfDd/obWpnpanQ51hZISk6irr2u/YTdU128YriFjSHj3FSoHjMCb3c/pSEGXmppKdRj925P/0p37tldiFP/9e/fy1EsL+PrP/i/P/vRBRvbPcjpWSNFvqvDVFX3bq1evi25v91bh3NxcTp06RVlZGT6fj/Xr15OXl3demwkTJrBq1SoANm7cyMiRI7Esi4aGBh577DG+9a1vcc0117S19/v91NbWAuDz+di6dSt9+/a90vcmIhJQgwcP5rvf/S5xcXFt28ybr0DlWaxbvx5Wa7a+V/qfEzElwYODXMR205HWiGBZlM+8Gzs6lp6Lnge/z+lEIvKfxuT24Wff/zY+n48f/vaPWipHpAu4f/3rX//6cg1cLhfZ2dnMmzePpUuXMn36dPLz83n99ddpbm6mV69e9OvXjxJa7H0AACAASURBVLVr1/Lqq69y5MgRvv/975OYmMi7777Lli1bOHbsGMuXL2f58uVMnjwZy7L4t3/7Nz766COWLVvGgAEDuOuuuzq0vMKXn50NJfHx8TQ2NjodQwJE/Ru+LtW3FRUVlJSUkJmZicv1X3/nM7u2Yv7/9u47Po7yTvz4Z2Zni1YrrbTqki3bkmxwbzLuMcY2/WhHCST0HCEGEuCSX8Il4cgRAlwCJMEhkMMQjiOJA9i0O2xMjDFu4G5j4yp3NUurvto6z++PlWU1Wza2tZL8fb9eemlnnmdmv7uPZjTfmWee+dt/wQVT0cdM6spQzxqlFG8dUry2z2SYG27vr2PrBUmr3WYnGOyp14s7p6w2wm4P7nVLwWLgzxsS65C6lMPhwO/3xzoMcRb0hrZNS0pk6KCBrFi/hdVeC2ML+pDitHa+4DlAjql6rzPRtm0fOXiUpnrYqEjFxcWxDuG4pNtD7ybt23t11LamaTJv3jy8Xi933XVX89VWVVeD+dgDYI9Du/tBNKPnH4SYSvHf+0zeOawYk6xxU67Wax5505u7CreU/n+vE797C8Xfe5xg34JYh9NlenJ3UnFivaltS2saeGVjFfXBCD8Y52FyflqsQ4o5OabqvWLaVVgIIc5Fn3/+OcXFxVx00UXHklalMF97Hhrq0a75Vq9IWkOm4rkd0aR1UqrGN3tR0nouqbjon4nEJ5D+19+hBRpjHY4QooVMdzz3T8gkPuDl4V/9nlc/+jzWIQnRI0niKoQQbRw+fJhVq1YxePBgBg8e3DxffbYINn2BdtEVaBkdnw3sSerDise+jPBZheLyLI1rczR0SVp7JNPhpPzSb2FUlZPyztxYhyOEaCPRbvDARefTJyuTP7zxDk/9dRGmacY6LCF6FElchRCihUgkwocffkhiYiIzZsxonq8O7UP97WUYMAgumBrDCM+MI37FTzZF2F4Ht/TTuChDP6lxBkT35c/Jo2r8LBI2LCN+w7JYhyOEaCPRGcd/3Hcbw4YM5q2FS3lwzjz8QRlUTYiTJYmrEEK0YLFYmDlzJpdffjl2ux0A1ejD/OOTYHegXX0Lmtazd53baxU/3BShMgj/kqczJrlnfx5xTPUFM2nMySN1wcsYR7rvmBBCnKvsNiuP3HU9F06ZxMp1m/n+nz+mNhCJdVhC9AhytCKEEE0aGhoA6N+/f/PAAEopzD//HirK0K67Dc2VGMsQT9vHpSY/2xLBqsF9A3UKEuQqa6+iWyi/9BaUxULG67+W+12F6IYsus69183i7m9/kwpXf364cB97vTLCrhCdkcRVCCGAkpISXn75ZXbs2NFqvvr4PVi/Em36FWi5eTGK7vRFlOLloghzdpsMiIfvD9LJdEjS2htFEpIpu/xWrEeKSXvrj9CzHh4gxDlj1pjzuG98JrV1ddz26O94+7ONsQ5JiG5NElchxDmvoaGB9957D6fTSW5ubvN8tWML6q0/w3nDYMKFMYvvdFUHo4MwfVCsmJqm8Z18HachSWtv5u87EO+UK4jfshr3svdiHY4Q4jj6JTn4XmEmic44nnx5Hv/+2v8SCkvXYSE6IomrEOKcZpom8+fPp7GxkauuuurYo2/Ki6P3tXpS0f7p5h47cNGXNYoHN0TYXgs39dW4OkeXx92cI2rGXEj9wJEkL/wLcTs3xTocIcRx5KQl8fSDdzF65HD+d8lybntiLsXe2liHJUS3I4mrEOKcpZRi6dKl7Nu3j5kzZ5KRkRGd31CP+fvHwVRoN96N5oiLcaSnzlSKNw+a/HxLBKsODwzSGZciu/xziqZxZNZNBFMySX/jWaylB2IdkRDiOBx2Gz++/Tpu/KdL2XvgAPfMeYdt5XLfqxAtyVGMEOKclpiYyIQJExg2bBgAKhzGfOlpqChFu+FONE9qjCM8dd6g4vGtEd7YbzIySeMHg3Sy4+Qq67lI2eyUXn03pmEj89UnsdR6Yx2SEOIErps+nn+bfSep54/l3xYf4JXV+6XrsBBNJHEVQpyTQqEQmqZRWFjIrFmzgOgVWPU/f4CvNqFdfkOPHIxpVYXJD9ZH+LIGruuj8a1+Gg6LJK3nskhCMqVX34XuqyPjz0+jBfyxDkkIcQJD+ufww2kDGJMZx8tv/J2b/v2P7Dh0JNZhCRFzkrgKIc45xcXFzJ07l8OHDzfPU0qh3vozasU/YOrFaCMviGGEp64xrHh+Z4Snt5u4rfDQeTqTUvUee2+uOLOC6X0ou/xWbCX7SP/LsxAOxzokIcQJOAydW0ZlcNOsKZSVH+HOXzzPnz5YjmmasQ5NiJiRxFUIcU4pLy9nwYIFGIZBcnJy83y1cD7qowVQOBntG5fEMMJTt6HK5AcbIiwpV8zM0HhgkE66POpGtNE4YAgVM67HuWMj6fN+DxHpfihEd3fF5NE88a/3kpGZxZ/e/F++9cv/Yt+RmliHJURMSOIqhDhnVFZW8vbbb2MYBtdffz1OpxMA30fvoOa/BkNHo11ybY+5SlkTUjy3I8IvtkbPwN83UOfSLBk1WBxf3bAJVHzjauK3rCb17T+CXL0Rotvrk5bMU9+/nWsvv5jD3np+9NFh3t/uJWLKM5rFucWIdQBCCNEVamtrefPNN9E0jRtuuIGkpCQAzKUfUvfGHyF/MNpVN6Np3f98nlKKpUcUrxSZ+CIwK0NjRoaGoUvCKjpXO+Yb6KEAnlULUVY7lVffDXr3/7sX4lxm0XVumjmRWZMKeXublz99foj/+ft8/t/NlzB2YN9YhydEl5DEVQhxTnC5XBQUFDB69Gg8Hg8A5pIPUH/9E9bBIwlffQuapfvvEvc3KOYWRdhcA/2d8N2+OpkyYrA4RdUXzEQLBUn+fDFaOETFP98ryasQPYDHaeU7Y9NZsrme11cV870n/sj0SYX8+OZL8SQ4Yx2eEGdV9z9KE0KI01BWVobT6SQhIYGZM2c2zzcXv4v6+1w4bxiu2++jpr4+hlF2rjak+NsBk4UlCrsFrs3RmJiqoUu3YPF1aBpVky9HGQae1R+hBf0cuen7YMhhgRDdnaZpzBiZz6gB9/HKu4tZsnItqzdu49YrZ3DnpeOxyEko0UvJfyghRK914MAB3nnnHfr06cN1110HgDJN1PzXUIsWwOCRaNd8G60bH6yHTcXCUsVfD5j4wjAxVeOSTI14QxJWcZo0jeoJl6BsDlKWvYceDFD+rYdRNnusIxNCnISUxHh+dOs1bJo4mtfeWcSfP17LV1oWd41JZ2iGXH0VvU/3PVoTQojT8NVXX7Fo0SKSkpKOPac1FEK9+lvUms+iowdffC1aNz0zbSrFygrFXw6YFDfCQBdclaOTJd2CxRlWM2YaptVG6pK3yXzp3ym//cdEEpM7X1AI0S2MLOjHbx7+Dp/v97J4XyP/7/1tpJRv5gc3XsLYgpxYhyfEGSOJqxCiV1FK8dlnn7FmzRr69OnDVVddRVxcHKquFvOPT8KurWgzroQJ07vl6MFKKb7wKv6y32S/DzIdcOcAnSGJdMt4Re9QN3wikfhE0j/8H7L/8Ahlt/+EYHb/WIclhDhJuq4zcUAqhbkm76yq4oMv9nPvL//AmOGD+cH1MxnSLyvWIQpx2iRxFUL0KqFQiL179zJ8+HBmzJiBxWJB7d+N+cKTUFsV7Ro8bEysw2xHKcW6quh9rLvrIc0O3+qnMTJJ7mMVXcOXN5TiG+4n871XyHrxUY7ceD++YRfEOiwhxCmwWnRumDKCmSPy+NtHy1m5Zj13/PtXFI4axo/vuI5+SY5YhyjE16YppXrUQ6CKi4tjHcJxpaamUlFREeswxFki7du9lZWVkZSUhN1uJxAIYLPZ0DQNc8U/UP/zAjhdaNffgZbd/rEBSUlJVFdXxyDq6D2snx1RLDhscsAHybbo423GejR5HusZkOBKoK6+LtZh9CiW+hoy3n8VR9lBaqZcifeyW6Cbjrgdy21XnF3StmdGeXUtf1u0gm3l9cQPHMf4Pi5mZML48/rENC45puq9zkTbZmdndzi/e/4nEkKIk6SUYv369SxbtowRI0YwY8YM7HY7yu/D/Nt/oVb8A/oPRLv2VrR4V6zDbeYLKz4uU7x32KQiCFkOuDlXY1SyJKwitiIuN8U33E/K8vdxL/8A+4EdlN/yEJGk1FiHJoQ4RelJiXz/psuoD0ZYvr+WJZuLmP/KB+T168NNsyZz1cThWA1LrMMU4qRI4iqE6LFqa2tZvHgx+/bto6CggEmTJgGgdm/DnPscVJbDlFlo37gYTe8e/5j3NSgWlph8Uq4ImJAfD1f30Tk/Qe5hFd2IYVB54bX4sweQ9vHfyfntD/H+0x3Uj5kG8ncqRI/jslm4dGAyk3KG8E68j89WreHJl+fx4luLuHTqOO64ZCIpCXGxDlOIE5LEVQjRIxUVFfHBBx8AMGPGDEaOHAnBIObbr6EWzQe3B+22+9H6DohxpBA0FasqFAtLTb6qBasGI5M1JqVq5DolCRDdV8OgUQTS+5C2eB5pb76A88vPqbz2Hhl1WIgeKtEZx22XTuXmmZNYsvZLFq/4gnkLP2VFIIML8zxMzjIYnZuC3k1H3BfnNklchRA9SiQSwWKxkJaWRv/+/Zk2bRputxu1ZS3mGy9Gr7KOGo8262o0e+wGoVBK8VUtfFJusqJC4YtAqg2uzNYY55HnsIqeI5yUSsn138O94TOSV35IzrMPUXXxN6kbPwss3aMngxDi1FgNC5dMGMklE0ayt6yKtUcifLavhr++/g6JDgvTLxjJjReOJS8zJdahCtFMBmc6g+RG895N2je26urqWL58OXV1ddxwww3N3WpVeTHm26/B+lWQmoF22fVo/fJPad1nahAQpRQHfbCiwmTpEUWZH2w6DHdHB1sqcCEjBHcxGZzpzLJ6y0n5ZD7Og7sIZPWj8uq7CfQ/P2bxyAA+vZe0bddrDIZZtGoDy9du5NDhYjQgr19fLps+mWsmDiXJceaud8kxVe8lgzMJIc5ZPp+PtWvXsmHDBkzTZOzYsZimid5Qh/rgb6hPF4FhQbvwcph4IVoXj36qlGJPPayqNFlVoSj2gwYUuGB6rsZwt4bdIsmq6B1CnnRKr/su8bs3k7LsPbJffJSGIYVUXfxNQpm5sQ5PCHEa4mwG10wbxzXTxrG/tIKPVm9i/Zat/GV9Me+U2Clwhhmge7lqwnD6psvtAqLryRXXM0jOHvVu0r5d7+DBg8yfP59IJMJ5553H5MmTcWOiPn4XtfRDCAVh1ITo4EuuxK/9Pqd6Zj8QUWypUWysUnzuVRwJgA7ku2B4ksYwt0aiVZLV7kCuuJ49WiiAe/0yktYtRQsGaBg1meqL/plQek6XxSBX5XovadvuwTRNSutDbC7zseyL9RStW44G9MnOpHD4+cwsHEJhQQ6WU7wnVo6peq+zecVVEtczSDbC3k3at2uUlJQQDAbp168fwWCQTz/9lLFjx5Ic8kcT1uUfQyQMg0ehTbsELSX9tN+zswMkpRT7fbCxSrG+SrGtVhFW0UGWChKiXYGHuuW+1e5IEtezT/f7SFq7hMSNy9HDIRoGj6XmG1dFuxCf5a7xktz0XtK23Y9Sil0Hy1i26Su+3L6TspJS0DQGXnIL43I95DsaGdXXQ+5JXI2VY6reS7oKCyF6hZ8u3g/AE7P6tZrf2NjI9u3b2fLcE5TV+8i5/tv069cPq2Gh9LCf8k+fwl22PVp59AS0idPRPGkAmL/+t1br0n/0K366ORx9nxEd7+I6Kv/p5jCUFfN4yfuoW2eztx621UaT1G3lPmot0ccE9Al4ubTyK0Yf2cbgmr2orNafxXbkMMG0HEpumN08L+vNF3AUF0Xjtdqby7PefAHH4T2gafiz81qtQw8FmuvajhxuLjtg85BbX4KO4tFR30VpOj8t/4hfpc1EAf6c6P2931v1InneItA0fnT50zy+8Gf8avidmOh8Z9c7zetL91fhiAQ4EJ/Jz0ffS//6Yn6x8U8ABA07+5yZDKrdT9Bi5VfD7wTgsY1/Qkfh123YzFD0teHAFg6ga3Dr5MfwWRzk1R/m1+t+D5qGabWjB/0UuXII6lbsKkSGr5Knht+OX7dx1+73+fno75JbX8q/7HqX3IZSnGaAR0feA8B/bHypOea98dk8Ovq7zbE+OvpegroVRyTAY9teY789hT8OvA6A2TveAjSy/BXYIiEeG3VsfSYa357yHwQtVjQUT6+bQ7bvSLSdNI2fj7qX+3a8RbbvCCVxqfRrKIm+vyubH499AFskxE+3vMLQmr08Ouq7bE3K5/zqvQR1K6amYVEmDjPIHlcf5q58HIcZ5NYpv8BnaRo0TNP4jw0vAvDoqOjnAdiXEL1i6dAhYMIvNrzIKwVXUerOJjsOihrAHg4AEDDs0Pb8s6bRvy66rlJ3NgETFNC/vpg4HfbbU/DptmhSebxz100Jp0Z0WXs4QFZjJfsSso8tc7SOshIXdxHGBVO57PBKLtu1kuyv1lHkymZx1ng+Sx9FozWOE50l15QJgNK6biRTR9OJJn+4dWRHP3Nbx5t/PBqQl2ynqCraVqd6leBU3+94yx597TA0npjZj39duA+AZy7tD8APm6aPxmpv8b20/AwKyE+2M3t8Vrtljjr6ng5Da/e9OgyNnARb8/ThumDze7Rc9mjMDqvO0WsrLdfVMoa238/Rz3i0LD/Z3vxeOQm25lh/0/TZX/i8pFU8LeM7Og8gJ8HG7PFZAK2+vxc+L2lVp6Nl25bNHp/V/L4tX7d19P1aalv3aJ2OPkdHy58uTdMYlJvJxyWK+Lg8xhWGGOz0U+lws+ZwPW+t+j8C3jIcCUm407O4ccL5TBmeR37WqT0Duu2xwvGOHcS5RxJXIURMrVq1itWrV2OaJilKMb1/NsPGjsSc9zLqi2XcUltNjcMNbg+4EtCvuPGMvr9SivLGCLUhaDQS+EXm5exYFcEfPY4mxQYjGw5wfmMJGdMvYvC7b0YTy6AfAP8ZjUaInqvW5mLegItZkHsh00vXcXHxar67awF37PmAVWnD+SRzLFvdeZjd5JnKQojT43DGc9n4AgBMpdia90+s/2oPK7bsouLAbp7f8xX/lZLNeVMvY2BKHM7q/Qzpm8bYgTmknlouKwQgiasQogsFGxuoLT3A/PnrmDVrFgkJCaSmpjJq5EiGuuyklG+Fumr4zSMoiwEDB/OX865iR8ZgfrH55dN+/4hSlDZCbQj8EXh8a4Td9YqaUFW0guHCaQYZnayR54K8eA23TSPrzaUAlNhmnHYMQvR2QYuNRTkTWZQ9gbz6w8ws+YKpZRu5sGw9tYaTL1KHsjptGJuTBxLW5TBEiN5A1zSGD8hh+IAcGlIGYkYiTM+CvV4/VRYb20tq2PrhByil0DSNtNQU+mWlM2bEMAqHDaJvopUUp7X5iQFCdET+Ywghzqra2lpWr17NoUOH2Lp2LwCVo3OpKdpNfIOX/K0byNu6HhrqowvY49CuvAnOH4HmiOOrpm69p8JUEDRhZYXJQR8c9CkO+hSHG6Fl7zVDVwxM0ChIcrCsuJEEbwmPHn6fkvGzj79yIcTJ0TSKEvrwp4Q+vJZ/JaO8O5lwZAsTj2xhZuka/LqVrUl5bE4eyJakfA7EZ3RpV2EhxNmjWywMHZDF0AFH56RTPvYHfLX3MLsOlnCw9Ajbig6yM+jig/I4wr5aatd9SEZqKo1WFy53Mh8k1jM0r29zsiuEJK5CiNOmlMLn81FeXk5ZWRmlpaXk5eUxYsQIdF1nx7at5DhsXOK2MUrVMvmrZWhfvB+9Pyk+AfIHoxUMRn2xDHQL2qjxJ3y/sKmoCkJFEEpTR1Dq8FBmT6bM4aH0izBVTbcW/ef2aH/fFBukO2BKmkaGA5aXK6w63D8o2mUxwWXni7JG5JBZiLMjaLHxedowPk8bhmGGGVG1izHeHYyo2s2dez4AoNZwsiuxL7sSc9mV0JddiX2pt8bHOHIhxJmSnuwmPdnNtDFDmgffqguEKasPsbfUYGP9IEqPVFB5aC9lRdt5bMMqPCOnEZfRH1VXwbdXLyIzzUNOajI5aR5yMzzk5WSQ6rLLM9LPEZK4CiFOmt/vp76+nurqagzDoH///iileOmll2hoaECFgtBQj8fQ6HtwJ5EV/0fcob3cW17SnBRWOZPR8grQ+g6A3DxITUfTooNwNKxbQ53FQX2dojakqA5CeQDCJvwq41IqDSfez8PUhloM5DHwejRl4gnWkeH3UuDS2F2nMHT4Vj+dNDvtnqO6pjLSlV+bEKKFsG6wPmUw61MGA5Dir2Zk1U4G1+yjoO4Qo/d9jN60hZfbkzkYn8GB+AwOOTM4GJ9BcVwqPmtcLD+CEOIMSbAbJNgNClL6MGtoHwD+sLqYgL+RKWkQtCfwyaEAtTVhvHU+9h4sJhA4NiBY2vgriEtKxVp9kEjxTjxuFylJiaQmJZCa5Ob8QfmkJjhxWkySHAZxdtvxQhE9wEklrhs3buTVV1/FNE1mzJjBNddc06o8FAoxZ84cioqKSEhI4MEHHyQ9PfqIigULFrBkyRJ0XefOO+9k1KhRJ7VOIUTXaWxspKGhgcbGRnw+H42NjVitVoYOHQrA/PnzObx/H4HaGvD7IdBI34R4cocPgiovw7ZsxFFXjduM4IpPwLTF4auJ48ukDHxZo2kcfDGNCSks9LnxGXEUejRqw1BXCrUHTWpDJrVhiAyIjlrLptaJpQ6UGy6Sww2kJ2q4reC2QqJVY/hbz5EWqMGmol2K9097ghd2RZfv45QzsEJ0d5WOJD7JLOSTzEKUpuMI+8mvO8ygugP0qy8ht6GMEVW7sKpj+4UGi4MjjmTKm36OOJKpsiVQY3NRY43+rrM6MaXrsRA9jqZpOOKcXDA0OjLy5soSUvr3Y/ZNE1BKUVXfyOEjXkorq7B5sqgPW9hXVM6+Yo2iQ6Vs3rGHcCgEQOa0G7HY46jdvYG6os3Y7TaccXEkxDuJd8YxfOolJMbZaKw4TKShmkSng4Q4O26XE3e8nX59++AwdAxl4rRbsRmyT4mlThNX0zSZO3cuP/vZz0hJSeGRRx6hsLCQPn36NNdZsmQJ8fHxPP/886xYsYI33niDhx56iEOHDrFy5UqeffZZqqqqePzxx/nd734H0Ok6hTiXmKZJJBJBKYVpmiilUEphdzhQaDQ2+mj0BwiFTcLhEKFwhFAkQnZ2DmbEpLSkGG9FJeFQkFA4RCgQImyaDD5/GJFIiC2bNlF8+BDBQAC/P0Ao4Ee3WJg68ULC4QifrVhKWXkpZsRERSJETBO7I56xoycRNGH3vt34Iyaaw40W1weSXOx1eVhSHU/QsBEcP4OAbiWknWC00KZbWAnCZ0cU8QY4DUiwQqZDI96AzK0rSYj4CU2cgcuAeAP+ut9E1+DxXW8BUDK59f2nOf7Ks9QqQohY8BsOtibnszU5v3mebkbI8HvJbSgjw19Jur+KNH8V6X4vw6t3ExcJtltPBI06azy11ngaDTs+iwOf0fRjcRC0OWg0HNRpNoK6laBuENINQro1uj/TjaZ5VoK6lYiuE9F0IpoFE42IpkcTY+miKESX0TQNT4ITT4KT4Xkt8obzLoDLLgCity/VNfqpqKpDOVw0huFg2lCKc93UNfioa/BR72vEW9vA9soA/rCf4s2bqD+4q/V7WQyyZ3wLgKotn+ErKcJiGFgNA6vVIC7exaBpV+EwdCp3bSBcW4ndZsNqNbBZrThdLgqGjcHQNY4c2E0k0IjNsGA1LNgMC/Hx8eT2G4DVolFVXgpmGMNiwWq1YFgsxNkduJPc6LpGY30dKIVh0bHoGtamek5HHLoOkVAIi65h0XUMi46u673yvuBOE9fdu3eTmZlJRkYGAJMmTWLNmjWtksy1a9dyww03ADBhwgReeeUVlFKsWbOGSZMmYbVaSU9PJzMzk927dwN0us6eZsPqTfzuw8+pOVLcar6mW8gbMw2AsqKt1HlLW5VbDBv9R00FoGTXJhqqWz+w1+pwkjt8IgDFO9bjq61qVW53JtBnaPR+wENbv8Dvq231XD5HQjI5g8cBcGDzCkJ+H6A4+uQ0Z1IqWYNGA7Bvw6eEQ4FW63d5ssgoGAHA3rX/wIyEWz03LTEth9S8YQAUfb6I1k+Og8TMfqT2Ox8zEmbvmo/bfW/unAI8fQcSDvo5sP6TduWe3PNIys4j2Ojj4KZPozObPp8CUgcMIzGzH4H6Gg5vWQ5tttG0glG40vrQWFtJyZer2q0/4/xCnJ5MfN5yyrZ/0apMAZlDJhKXlEr9kcMc2bm+3TPjskdMwZbgoa50PxV7NtJ2BTljLsLmTKD60E4q9mwBokPGK6I71pwp16I74qkq2kJ10eZou2ig0FBoZF10M7pho2bHGur3b2sXf/as29A0jeptK2k41NEOty8ANTuqCFTWoFmsaEYcuiUR3e5kvTc6Hr0/eyoqLYhus6Nb7ehWBzarwScaWDGxec7DqimsuobVomE07XgNw8Cqg1UHmw4Oi4ZDB4el6UePdtM9+vq1vRE04L5BHe96MldsBqDUPat5nkUzW3+uto3chtaixvHqti23WTR0LVqga2DXWy93dNKuR+tpHZS1XHfbdbSs07K85fy2dVrW7exfz/HWc6J5Xa2rQugGH/WsOu227KFfkKlbKHGmUeJMa1+oFPHhRpKCdbhD9SQF63EH60kKRX8nhHzERfzEhxtJC1TjDDcSFw4QZ7ZPdr+OCDqmph1LajUNU4smubrFQshs2qc37dubw26aB1r0fwIaaC1eA6ppurlO04Foczk0LdN5w+oa5Oyw8XR103NH90S7TD7VNG03NAJh1fw3Zjb9wzs6/+jrzLXWdsu0pWnHlm/5/jbLsatVwYjZrk7b+ke1rOcwdDLXWnmyOtDhMjk77M1ljqarY8GIidWiN8eaUxR9vuv9LZ61GowobG1uKwlGovVtFo2sdU3fV9OzYHOKwHBuWQAADlJJREFU7NxfF2xVp6Nl25ZlrbNF31fTyFpv5f66UIefP2u9td28tnWP1mk5PxgxsVn0Dpc/U+6vC53wfY7Gc6IYLBaD+EjHgy+2Xf5k1nfKrE0/26LP0jbtivp8nSplUGtaqFU6dVhwlS8hoBscsFdTkWbHr8BvQoAwZsSH59BO/LpBw8ED1NR4CZqKiKkImRH0uAS2aX0JaxYq1qwmUFXeKgSbO5W08VcAUL7qfUJ13lbldk8mqYWXAFC2/G3CvvpW5XHpffGMugiA0qXziASPPaBPQ+HKGkDKsMnoSnH407+DaaJpMHXKxTx557Qz9112oU4TV6/XS0pKSvN0SkoKu3btOm4di8WC0+mkrq4Or9fLwIEDm+t5PB68Xm/zek60zqM+/vhjPv44mvA89dRTpHbTBz+53YnEGxZCttZfqa7pJOjRA+8Gq46ytd7oDMOK+2i5zcDSpu+9zWaQbDlWbnO0LrfbDVIs0e5TDQ4rAWVHa5FeORwGaZYwGgpfnJWQbm+1fLzdQqYlukNojLcTDrXuApHg0Mkyojt2n8uBMlt34XTHWchoKvcnHLvn6OhZnqQ4C2lGEFOLEHK3H2QjxamTagQIqxC4XceWb/qdGqfhsQYImkEsSa52y2c4NZKsfvyOMA5PYrvyzDhFouGj0R7B1bZcg+w4E5e1kXpnhKRUd7vl+8SFcFobqHWalKa5mxZT0YU1jdy4MA5bA9UuRXl6cqtldSDXEcBu06lOsuHNzULXNHSi/1x1TSPXVo3d0kBtqkGNkRst10An+jtXP4ih6dRmGzQk5WHoOoZhYFgMLIaFDFcpukUnMKYPalQWhsWInsmzWrHbDOzWeqyGjnHeGKy28VhtVgy7HcNmw2q3Y7XbsDocWO0jMRx2rBY9mpjq0TN6Z9qszipMmwBA/xazJjTvW6e0K4sus6TVZHarZTrWttwwDF4bdfQf6I3ktl1gajSu3BavO9NqHR0scyrraqnlI+X/1Pzq1havm9w9pfnlawB3L21Vv61sYHHz1O3Nr/p3+H63cyKLW019s9VUZpu6x9Z5M581v77phO+ZSetYWz8o6XaygIlNU4ZxK+HwsYOjtutr/dfT+g/jjRax9G8xPxtY3jw1o816p9Chf43WW9yuYEqLtRzPlE7Kz12GYbRq31OhImFUoy/6EwxG79EPBlHBACoUQAWDEAygQtF5hMMo04RIBMwIquk3kaOvzY7nK9V0wlUdO7Hc4rdqLqNNXVov13TSs/U6FO3OqHZiaNqJp7/OOlpRpxjQcWia1vTddGzYCQ4H25d1vJ6BKR3OPqFhLZb5OssDDPQci2dg8gkqtl3uOHVPZR1nQmfvdzLxnKh92y7fVZ8vkej+vbWmZDId2qdNEQitib7sC/RtWx6GmgUowJ+v8JsuAkDQ1PADESLYaxcT0nTK+0NjJImwig5AGVZgsUZIbVhFBI2DfWz4I25MFe1NElFgd4RJ820komkU5TgJRRxE0DCVIgI4403SA3tQwM50F6apMFGMj2s8q/mUYRhnbf3dfnCmmTNnMnPmzObpioqKE9SOnbzBA5j3ux93Et+FnazlbJdPj3H5zE7KL+mk/PJOyq/spPyqTspPLDU1tdv+/Z02M0DEFyAC+Dut3Pv06rY9x0nb9m5npH01A+wG2J1nJihxRsi227uda+1ra/o5noKzHsHVrabO5nd/Jto2O7v9KQSg86c/eDweKiuP3UNWWVmJx+M5bp1IJILP5yMhIaHdsl6vF4/Hc1LrFEIIIYQQQggh4CQS1/z8fEpKSigvLyccDrNy5UoKCwtb1Rk7dixLly4FYPXq1QwdOhRN0ygsLGTlypWEQiHKy8spKSmhoKDgpNYphBBCCCGEEELASXQVtlgs3HXXXTzxxBOYpsn06dPp27cv8+bNIz8/n8LCQi666CLmzJnDAw88gMvl4sEHHwSgb9++TJw4kYcffhhd17n77rvR9Wiu3NE6hRBCCCGEEEKItjR1ojvfu6Hi4uLOK8XIudZf/1wj7dt7Sdv2XtK2vZu0b+8lbdu7Sfv2XjG9x1UIIYQQQgghhIglSVyFEEIIIYQQQnRrkrgKIYQQQgghhOjWJHEVQgghhBBCCNGtSeIqhBBCCCGEEKJbk8RVCCGEEEIIIUS3JomrEEIIIYQQQohuTRJXIYQQQgghhBDdmiSuQgghhBBCCCG6NU0ppWIdhBBCCCGEEEIIcTxyxfUM+slPfhLrEMRZJO3be0nb9l7Str2btG/vJW3bu0n79l5ns20lcRVCCCGEEEII0a1J4iqEEEIIIYQQoluzPPbYY4/FOojeJC8vL9YhiLNI2rf3krbtvaRtezdp395L2rZ3k/btvc5W28rgTEIIIYQQQgghujXpKiyEEEIIIYQQolszYh1AT/b3v/+df/zjHyQmJgJw8803M2bMmHb1Nm7cyKuvvoppmsyYMYNrrrmmq0MVX8Prr7/OunXrMAyDjIwMZs+eTXx8fLt69913Hw6HA13XsVgsPPXUUzGIVpyMzrbFUCjEnDlzKCoqIiEhgQcffJD09PQYRStOVkVFBX/4wx+orq5G0zRmzpzJ5Zdf3qrO1q1b+c///M/m9hw/fjzXX399LMIVX0Nn+1mlFK+++iobNmzAbrcze/Zs6YbYAxQXF/Pcc881T5eXl3PjjTdyxRVXNM+TbbdneeGFF1i/fj1ut5tnnnkGgPr6ep577jmOHDlCWloaDz30EC6Xq92yS5cuZf78+QBcd911XHjhhV0ZuuhER23b5cfKSnxt8+bNU+++++4J60QiEXX//fer0tJSFQqF1A9/+EN18ODBLopQnI6NGzeqcDislFLq9ddfV6+//nqH9WbPnq1qamq6MjTxNZzMtrhw4UL10ksvKaWUWr58uXr22WdjEao4RV6vV+3Zs0cppZTP51Pf//7327Xtl19+qZ588slYhCfOgM72s+vWrVNPPPGEMk1T7dixQz3yyCNdGJ04EyKRiPrOd76jysvLW82Xbbdn2bp1q9qzZ496+OGHm+e9/vrrasGCBUoppRYsWNDh8VRdXZ267777VF1dXavXovvoqG27+lhZugqfZbt37yYzM5OMjAwMw2DSpEmsWbMm1mGJkzBy5EgsFgsAgwYNwuv1xjgicTpOZltcu3Zt8xneCRMm8OWXX6JkGIBuLzk5ufnqWlxcHDk5ObK9nmPWrl3LN77xDTRNY9CgQTQ0NFBVVRXrsMQp2LJlC5mZmaSlpcU6FHEahgwZ0u5q6po1a5g2bRoA06ZN6/A4eOPGjYwYMQKXy4XL5WLEiBFs3LixS2IWJ6ejtu3qY2XpKnyaFi1axLJly8jLy+O2225r16Ber5eUlJTm6ZSUFHbt2tXVYYrTtGTJEiZNmnTc8ieeeAKAWbNmMXPmzK4KS5yCk9kWW9axWCw4nU7q6uqabwcQ3V95eTl79+6loKCgXdnOnTv50Y9+RHJyMrfeeit9+/aNQYTi6zrRftbr9ZKamto8nZKSgtfrJTk5uUtjFF/fihUrmDx5codlsu32bDU1Nc3bYlJSEjU1Ne3qtP0f7fF45ARkD9MVx8qSuHbi8ccfp7q6ut38b37zm1x88cXN91nMmzeP//7v/2b27NldHaI4DSdq33HjxgEwf/58LBYLU6dOPe46PB4PNTU1/PKXvyQ7O5shQ4ac1biFEO35/X6eeeYZ7rjjDpxOZ6uyAQMG8MILL+BwOFi/fj2//vWv+f3vfx+jSMWpkv1s7xYOh1m3bh233HJLuzLZdnsXTdPQNC3WYYgzrKuOlSVx7cTPf/7zk6o3Y8YMnn766XbzPR4PlZWVzdOVlZV4PJ4zFp84PZ2179KlS1m3bh2PPvrocXe0R9vT7XYzbtw4du/eLQdU3dDJbItH66SkpBCJRPD5fCQkJHR1qOJrCIfDPPPMM0ydOpXx48e3K2+ZyI4ZM4a5c+dSW1srV9N7iM72sx6Ph4qKiuZp+V/bs2zYsIEBAwaQlJTUrky23Z7P7XZTVVVFcnIyVVVVHbadx+Nh27ZtzdNer1eOpXqIrjxWlntcT0PL+2e++OKLDruu5OfnU1JSQnl5OeFwmJUrV1JYWNiVYYqvaePGjbz77rv8+Mc/xm63d1jH7/fT2NjY/Hrz5s3k5uZ2ZZjiJJ3Mtjh27FiWLl0KwOrVqxk6dKicGe4BlFK8+OKL5OTkcOWVV3ZYp7q6uvl+5d27d2OappyU6CFOZj9bWFjIsmXLUEqxc+dOnE6ndBPuQU7UTVi23Z6vsLCQTz/9FIBPP/20uUdbS6NGjWLTpk3U19dTX1/Ppk2bGDVqVFeHKk5RVx8ra0pGHvnann/+efbt24emaaSlpXHPPfeQnJyM1+vlpZde4pFHHgFg/fr1vPbaa5imyfTp07nuuutiHLk4GQ888ADhcLj5vuWBAwdyzz33tGrfsrIyfvOb3wAQiUSYMmWKtG831tG2OG/ePPLz8yksLCQYDDJnzhz27t2Ly+XiwQcfJCMjI9Zhi05s376dRx99lNzc3OYTDTfffHPzFbiLL76YhQsX8tFHH2GxWLDZbNx2222cd955sQxbnKTj7Wc/+ugjINq+Sinmzp3Lpk2bsNlszJ49m/z8/FiGLU6S3+9n9uzZzJkzp/nqasu2lW23Z/ntb3/Ltm3bqKurw+12c+ONNzJu3Diee+45KioqWj0OZ8+ePSxevJh7770XiN4juWDBAiD6OJzp06fH8qOINjpq2wULFnTpsbIkrkIIIYQQQgghujXpKiyEEEIIIYQQoluTxFUIIYQQQgghRLcmiasQQgghhBBCiG5NElchhBBCCCGEEN2aJK5CCCGEEEIIIbo1SVyFEEIIIYQQQnRrkrgKIYQQQgghhOjWJHEVQgghhBBCCNGt/X+T6IPXiKwhxQAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 1152x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ps8koqaeKZgN"
      },
      "source": [
        "## [EM-algorithm](https://en.wikipedia.org/wiki/EM_algorithm_and_GMM_model)\n",
        "Evaluate the parameters for the Gaussian distributions with observed data $\\textbf X$ by integrating out latend variables $\\textbf Z$. \n",
        "We want to maximize $ l(\\boldsymbol{\\theta}) =\\sum_{i=0}^{N-1} log p(x_i;\\boldsymbol{\\theta})$, where  $\\boldsymbol{\\theta} = (\\pi,\\mu,\\sigma)$.\n",
        "<br/>\n",
        "\n",
        "_Here is the derivation of EM-algorithms for general case:_\n",
        "\n",
        "\n",
        "\\begin{align}\n",
        "l(\\boldsymbol{\\theta})&=\\sum_{i=0}^{N-1} log p(x_i;\\boldsymbol{\\theta})\\\\\n",
        "&=\\sum_{i=0}^{N-1}log \\sum_{j=0}^{M-1}p(x_i,z_j;\\boldsymbol{\\theta})\\\\\n",
        "&=\\sum_{i=0}^{N-1}log \\sum_{j=0}^{M-1} q(z_j)\\frac{p(x_i,z_j;\\boldsymbol{\\theta})}{q(z_j)}\\\\\n",
        "&=\\sum_{i=0}^{N-1}log \\mathbb{E}_{z_j \\sim q(z)}\\frac{p(x_i,z_j;\\boldsymbol{\\theta})}{q(z_j)}\\\\\n",
        "&\\leq ^{Jensen's inequality} \\mathbb{E}_{z_j \\sim q(z)} log \\frac{p(x_i,z_j;\\boldsymbol{\\theta})}{q(z_j)} \\\\\n",
        "&=\\sum_{i=0}^{N-1}\\sum_{j=0}^{M-1} q(z_j)log \\frac{p(x_i,z_j;\\boldsymbol{\\theta})}{q(z_j)} \n",
        "\\end{align}\n",
        "\n",
        "Let $q(z_j)$ be the posterior $p(z_j|x_i;\\boldsymbol{\\theta})$:\n",
        "\n",
        "\n",
        "\\begin{align}\n",
        "l(\\boldsymbol{\\theta})&= \\sum_{i=0}^{N-1}\\sum_{j=0}^{M-1} p(z_j|x_i;\\boldsymbol{\\theta})log \\frac{p(x_i,z_j;\\boldsymbol{\\theta})}{p(z_j|x_i;\\boldsymbol{\\theta})}\\\\\n",
        "&=\\sum_{i=0}^{N-1}\\sum_{j=0}^{M-1} p(z_j|x_i;\\boldsymbol{\\theta}) \\left ( log p(x_i,z_j;\\boldsymbol{\\theta}) -logp(z_j|x_i;\\boldsymbol{\\theta}) \\right )\\\\\n",
        "\\end{align}\n",
        "\n",
        "So we use current parameters $\\boldsymbol{\\theta}^k$ to calculate $p(z_j|x_i;\\boldsymbol{\\theta}^k)$, which is also called *responsibility*, then to update $\\boldsymbol{\\theta}^{k+1}$ by maximizing $l(\\boldsymbol{\\theta}^{k+1})$:\n",
        "\n",
        "\n",
        "\\begin{align}\n",
        "\\boldsymbol{\\theta}^{k+1} &= \\underset{\\boldsymbol{\\theta}^{k+1}} {\\mathrm{argmax}}\\,l(\\boldsymbol{\\theta}^{k+1})\\\\\n",
        "&=\\underset{\\boldsymbol{\\theta}^{k+1}} {\\mathrm{argmax}}\\, \\sum_{i=0}^{N-1}\\sum_{j=0}^{M-1} p(z_j|x_i;\\boldsymbol{\\theta}^k) \\left ( log p(x_i,z_j;\\boldsymbol{\\theta}^{k+1}) -logp(z_j|x_i;\\boldsymbol{\\theta}^k) \\right )\\\\\n",
        "\\end{align}\n",
        "\n",
        "The last term on the RHS above is constant w.r.t. $\\boldsymbol{\\theta}^{k+1}$\n",
        "\n",
        "_In summary:_\n",
        "1. Initialize parameters $\\boldsymbol{\\theta}^0 $ randomly<br/>\n",
        "2. until convergence:<br/>\n",
        "  1. E-step:<br/>\n",
        "      calculate responsibilites $p(z_j|x_i;\\boldsymbol{\\theta}^k)$ with current parameters <br/>\n",
        "  2. M-step:<br/>\n",
        "      update parameters with responsibilites:\n",
        "\\begin{align}\n",
        "\\boldsymbol{\\theta}^{k+1}=\\underset{\\boldsymbol{\\theta}^{k+1}} {\\mathrm{argmax}}\\, \\sum_{i=0}^{N-1}\\sum_{j=0}^{M-1} p(z_j|x_i;\\boldsymbol{\\theta}^k) \\left ( log p(x_i,z_j;\\boldsymbol{\\theta}^{k+1}) -logp(z_j|x_i;\\boldsymbol{\\theta}^k) \\right )\\\\\n",
        "\\end{align}\n",
        "\n",
        "_Here is the EM-algorithm in GMM, derivation is omitted:_\n",
        "1. Initialize parameters $\\boldsymbol{\\theta}^0 $ randomly<br/>\n",
        "2. until convergence:<br/>\n",
        "  1. E-step:<br/>\n",
        "      calculate responsibilites $p(z_j|x_i;\\boldsymbol{\\theta}^k)$ with current parameters: <br/>\n",
        "      \\begin{align}\n",
        "      \\widetilde{\\rho} ^{i}_{j} &= \\pi_j \\frac{1}{\\sigma_j \\sqrt{2\\pi}}exp(-\\frac{1}{2\\sigma^2}(x_j^i - \\mu_j)^2)\\\\\n",
        "      \\rho ^{i}_{j} &= \\frac{\\widetilde{\\rho} ^{i}_{j} }{\\sum_{c=0}^{M-1} \\widetilde{\\rho} ^{i}_{c}}\\\\\n",
        "      \\gamma_{j} &= \\sum_{i=0}^{N-1}\\rho ^{i}_{j}\\\\\n",
        "      \\end{align}\n",
        "  2. M-step:<br/>\n",
        "      update parameters with responsibilites:\n",
        "      \\begin{align}\n",
        "      \\pi_j &:= \\frac{\\gamma_j}{N}\\\\\n",
        "      \\mu_j &:=\\sum_{i=0}^{N-1}\\frac{\\rho_j^i x_j^i}{\\gamma_j}\\\\\n",
        "      \\sigma_j &:= \\sqrt{\\frac{1}{\\gamma_j}\\sum_{i=0}^{N-1}\\rho_j^i(x^i-\\mu_j)^2}\n",
        "      \\end{align}\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7K5FLxrvHr-g"
      },
      "source": [
        "def EM(n_iter, dataset):\n",
        "    \"\"\"EM algorithms.\"\"\"\n",
        "    # random initilization\n",
        "    pis = np.array([0.2, 0.8])\n",
        "    mus = np.array([1.5, 3.0])\n",
        "    sigmas = np.array([2.0, 2.0])\n",
        "\n",
        "\n",
        "    pis_set = np.zeros((n_iter,2))\n",
        "    mus_set = np.zeros((n_iter,2))\n",
        "    sigmas_set = np.zeros((n_iter,2))\n",
        "    pis_set[0] = pis\n",
        "    mus_set[0] = mus\n",
        "    sigmas_set[0] = sigmas\n",
        "\n",
        "    for i in tqdm(range(n_iter)):\n",
        "        # calculate responsibilites with current parameters\n",
        "        # E-step\n",
        "        responsibilities = (\n",
        "            tfd.Normal(loc=mus, scale=sigmas).prob(dataset.reshape(-1, 1)).numpy() * pis\n",
        "        )  # n_samples * c_classes matrix\n",
        "        responsibilities /= np.linalg.norm(\n",
        "            responsibilities, axis=1, ord=1, keepdims=True\n",
        "        )  # normalzied over rows(n samples)\n",
        "\n",
        "        # responsibilities for each class\n",
        "        class_responsibilities = np.sum(responsibilities, axis=0)\n",
        "\n",
        "        # M-step update parameters\n",
        "        for j in range(len(pis)):\n",
        "            pis[j] = class_responsibilities[j] / sample_size\n",
        "            mus[j] = (\n",
        "                np.sum(responsibilities[:, j] * dataset) / class_responsibilities[j]\n",
        "            )\n",
        "            sigmas[j] = np.sqrt(\n",
        "                np.sum(responsibilities[:, j] * (dataset - mus[j]) ** 2)\n",
        "                / class_responsibilities[j]\n",
        "            )\n",
        "        \n",
        "\n",
        "        pis_set[i] = pis\n",
        "        mus_set[i] = mus\n",
        "        sigmas_set[i] = sigmas\n",
        "\n",
        "    param_dict = {\"pis\":pis_set, \"mus\":mus_set,\"sigmas\":sigmas_set}\n",
        "\n",
        "    return pis, mus, sigmas, param_dict"
      ],
      "execution_count": 123,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Cgyg8ILeH4-T",
        "outputId": "199c4c18-8ccc-4941-8092-f64cb56887db"
      },
      "source": [
        "n_iter = 100\n",
        "pis, mus, sigmas, param_dict = EM(n_iter , x)\n",
        "print()\n",
        "print(pis)\n",
        "#==> [0.25, 0.75]\n",
        "print(mus)\n",
        "#==> [2.0, 5.0]\n",
        "print(sigmas)\n",
        "#==>[1.0, 1.8]"
      ],
      "execution_count": 135,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "100%|██████████| 100/100 [00:00<00:00, 874.19it/s]"
          ],
          "name": "stderr"
        },
        {
          "output_type": "stream",
          "text": [
            "\n",
            "[0.26413991 0.73586009]\n",
            "[2.07717848 4.88057437]\n",
            "[1.02806032 1.86795329]\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "stream",
          "text": [
            "\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8BleZ6oYpQtD"
      },
      "source": [
        "## Visualization of EM process, a .gif<br/>\n",
        "![SegmentLocal](EM-process.gif \"EM-process\")<br/>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6rRiIjp4YsGi"
      },
      "source": [
        "from celluloid import Camera\n",
        "from matplotlib.animation import writers"
      ],
      "execution_count": 136,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "me75Yj8MH7U0",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 391
        },
        "outputId": "9e49c876-6e98-451a-997a-578e19ef099d"
      },
      "source": [
        "xx = np.linspace(-5, 12, 1000)\n",
        "fig = plt.figure(figsize=(16, 6))\n",
        "\n",
        "camera = Camera(fig)\n",
        "for i in range(n_iter):\n",
        "    #ground truth plot\n",
        "    plt.title(\"Gaussian Mixture Model\")\n",
        "    plt.plot(xx, Gaussian_models[0].prob(xx) * true_pis[0], label=\"Gaussian distribution 1\",color='skyblue')\n",
        "    plt.plot(xx, Gaussian_models[1].prob(xx) * true_pis[1], label=\"Gaussian distribution 2\",color='salmon')\n",
        "    plt.plot(\n",
        "        xx,\n",
        "        (\n",
        "            Gaussian_models[0].prob(xx) * true_pis[0]\n",
        "            + Gaussian_models[1].prob(xx) * true_pis[1]\n",
        "        ),\n",
        "        label=\"Gaussian Mixture Model\",\n",
        "        ls=\"dashed\",\n",
        "        color=\"black\",\n",
        "        alpha=0.5,\n",
        "    )\n",
        "    plt.fill_betweenx(Gaussian_models[0].prob(xx) * true_pis[0], xx, color='skyblue', alpha=0.2)\n",
        "    plt.fill_betweenx(Gaussian_models[1].prob(xx) * true_pis[1], xx, color='salmon', alpha=0.2)\n",
        "    plt.vlines(x, 0, 0.005, color=np.array([\"skyblue\", \"salmon\"])[states])\n",
        "    # update evaluation\n",
        "    univariate_gmm = tfd.MixtureSameFamily(\n",
        "    mixture_distribution=tfd.Categorical(probs=param_dict['pis'][i]),\n",
        "    components_distribution=tfd.Normal(loc=param_dict['mus'][i], scale=param_dict['sigmas'][i]),\n",
        ")\n",
        "    plt.plot(\n",
        "    xx, univariate_gmm.prob(xx), color='yellowgreen', label='EM evaluation', alpha=0.5\n",
        "    )\n",
        "    camera.snap()\n",
        "\n",
        "#Creating the animation from captured frames\n",
        "animation = camera.animate()\n",
        "\n",
        "#Saving the animation\n",
        "animation.save('EM-process.gif', writer = 'pillow')\n",
        "\n"
      ],
      "execution_count": 137,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA64AAAF2CAYAAABu9MUZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de1hU9b7H8c8wkxcuITMjEIlao9ZW87KbzMhUhOq0q5PHrG0X22lt6+C2sMsuyjqe2uwobx1vZW4fMvUU5s7HLrsbGbm3aKGJZVhJ2hWKYCjBSwKzzh89zWkCYwyVn6v363l6HtZav/Vb38XXofmw1hoclmVZAgAAAADAUFHtXQAAAAAAAD+H4AoAAAAAMBrBFQAAAABgNIIrAAAAAMBoBFcAAAAAgNEIrgAAAAAAoxFcAQC/StOnT1evXr3au4wWmVybKXr27Km//OUvh7SPw+HQ8uXLj1BFAIAjieAKADhsAoGAcnJy1LdvX0VHRyshIUGDBg3S3Xffrc8++6y9ywtz2223aePGjUflWD179pTD4dDcuXObbZs6daocDocyMzPbVJvL5dLjjz/e1lLb5OOPP5bD4ZDD4dC7777bbPvgwYPlcDgOOXACAEBwBQAcFp999pkGDx6slStXKicnRxs3blRpaakefvhh1dTUaObMme1dYpjY2Fh5vd6jdrzu3bvrb3/7W9i6/fv364knnlCPHj3atbafamhoaNP+3bt31+LFi8PWvfXWW9qxY4c8Hk+b5gYA/DoRXAEAh0VWVpYOHDigLVu2aPz48RowYIB69OihkSNH6tFHH9XDDz8cGvvqq69q5MiRcrvdio+P14gRI/TWW2+FzdfSbZ2ZmZm69tprQ8tr1qzR4MGDFR0drS5dumjIkCHasmWLpO/D1y233KJu3bqpY8eOOuGEEzRu3LjQvj+9HXfXrl0aM2aMUlJSFB0drdNOO03Lli0LO/7IkSN1/fXX6/7771dycrLcbreuueYa1dfXt/r9GTdunHbu3Kk333wztG7VqlVKSEjQiBEjwsb+uDbLsnThhRfqjDPOCAXKYDCozMxMDR8+XE1NTerZs6eampo0YcKE0BVPSXr88cflcrnC5v7888/lcDhUVFQkSSoqKpLD4dALL7ygYcOGqVOnTqGAPW/ePJ166qnq1KmTevfurdzcXDU2NrZ6rtddd52WL1+u/fv3h9Y99thj+v3vf6/Y2NiwsXV1dbrhhhvUtWtXdezYUX6/X6+88krYmK1btyotLU0dO3ZU7969tXLlymbHrK+v180336wTTzxR0dHRGjx4sJ555plWawUAHBsIrgCANgsEAvrHP/6hKVOm6Pjjj29xzA9hSvo+ZGRlZWnDhg0qLi5W79699W//9m+qqamJ+JhffvmlLrvsMl1xxRV67733tGHDBmVnZ4eC2rx587Ry5UotX75cO3bs0LPPPquhQ4cedL76+nqNGjVKL774ot59911NmjRJEyZM0Ouvvx42btWqVQoEAioqKtJTTz2l559/Xg8++GCr9cbFxWncuHFhVyIfe+wxXX/99WHfm59yOBx6/PHH9cUXXygnJ0eS9MADD2jLli1asWKFnE6nSkpK5HQ69fDDD6uyslKVlZWt1vNTt956q+644w5t375dF198saZPn66ZM2fqgQce0Pbt2/U///M/WrRokf77v/+71bnS09Pl9Xq1atUqSd+H06eeekp//OMfm42dOHGiXn75ZS1fvlylpaU6++yzddFFF+n999+XJO3bt0+/+93v1KVLF7311lt64oknNGPGDFVVVYXmsCxLF198sbZu3aqCggJt27ZN//mf/6lx48bptddeO+TvBQDAQBYAAG305ptvWpKsZ555Jmz9WWedZcXExFgxMTFW3759D7p/U1OT1aVLF2v58uWhdZKsZcuWhY3LyMiw/vCHP1iWZVlvv/22JcnatWtXi3PedNNNVnp6uhUMBlvc/l//9V+Wz+f72fP693//d+v6668PLY8YMcIaMGBA2Jgbb7zRGjp06M/O06NHD+v++++33nzzTSsmJsbavXu3tX37duu4446zvvzyS+sPf/iDlZGR8bO1rV271nI6ndb06dMtl8tlrV69Omy70+m08vPzw9bl5+dbTqczbN1nn31mSbJef/11y7Is6/XXX7ckWU888URozJ49e6zOnTtbL774Yti+S5cuteLj4w96nrt27bIkWf/85z+tBx980Bo+fLhlWZb1yCOPWKeddlrY98KyLGvHjh2WJOuFF14Im2fw4MHWhAkTLMuyrMWLF1sxMTFWIBAIbX/33XctSaF5Xn/9datjx47WN998EzbPhAkTrEsuuSS03NK/KQDAscF10EQLAMAhsiwrbLmgoEDfffedFi5cGHbb5q5du3Tvvfdqw4YNqqqqUjAY1N69e/XJJ59EfKwBAwbo/PPPV//+/XXuuedq5MiRGjNmjFJTUyVJEyZM0LnnnqtevXrp3HPP1bnnnquLL75YHTp0aHG+vXv36r777tNzzz2nyspKHThwQN99953S09PDxg0cODBsOSUlRS+//HJENQ8ZMkS9evXSk08+qffff18XX3yxkpKSIto3PT1dt956q6ZPn64bb7xRo0ePjmi/SA0ZMiT09Xvvvad9+/bp0ksvDbsa3NTUpP379+vrr79W165df3a+a6+9Vvfcc48++OADLV68uMWrrWVlZZKk4cOHh60fPny4NmzYEBrzm9/8RgkJCaHt/fv3V3x8fGi5pKREBw4c0Iknnhg2z4EDB9S7d+/WTh0AcAwguAIA2qxXr16KiorS9u3bw9b/ECLdbnfY+osuukher1cLFixQamqqOnTooGHDhunAgQOhMQ6Ho1kQ/vGHBjmdTr344osqKSlRYWGh/v73v+vOO+/U008/rYsuukiDBg3Srl279Oqrr+r111/XzTffrHvuuUcbN25s8Xbm22+/XWvWrNHs2bN1yimnKCYmRrfeequ+/fbbsHE/Db4Oh0PBYDDi79WkSZP0yCOP6LPPPtOKFSsi3q+pqUnr16+X0+nURx99JMuyfvYWY0mKimr+RNDBPngpJiYm9PUP5/P000+rT58+zcb+tJ8tSUxM1CWXXKLJkydr+/btGj9+fKv7/FLBYFDx8fEqKSlptu1gv6gAABxbeMYVANBmbrdbF1xwgebNm9cs6P1UTU2NysrKdOedd+r8889X37591alTp7BnFqXvg09FRUVo+bvvvgtdofuBw+HQkCFDdNddd2ndunUaMWKE8vPzQ9tjY2P1H//xH5o7d642bdqk7du364033mixrnXr1umqq67S5ZdfroEDB+rkk0/Whx9+eKjfilZdffXV2rFjh+Li4nTuuedGvN/06dNVXl6u9evX66233tJDDz0Utr1Dhw5qamoKW5eYmKimpiZ99dVXoXVvv/12q8fq16+fOnXqpJ07d6pXr17N/nM6nRHVfMMNN+i1117T2LFj1aVLlxaPI33/vf+xdevWqX///pKkvn37avv27frmm29C2997772wf2d+v1/ffPON9u/f36zW7t27R1QrAMBsXHEFABwWCxcu1Nlnn63Bgwdr+vTpGjRokGJjY/XBBx/o+eefD4WdhIQEde3aVYsXL5bP51NNTY3+/Oc/q3PnzmHzZWZm6tFHH9Xw4cMVFxen3NzcsCuyxcXFeu2113TeeefphBNO0I4dO/TOO+/ouuuukyTNmDFDKSkpGjRokKKjo/Xkk0/K6XS2eAVRkk455RStWbNGl156qWJjYzV79mxVVFREfCtvpI4//nh98cUXioqKavGKaEveeOMN5eXl6fnnn9eZZ56pxx57TFdffbXS09NDt/iedNJJev3113XBBReoQ4cO8nq9GjJkiOLi4nTnnXfqrrvu0kcffaT77ruv1ePFxsbqrrvu0l133RX6G7ONjY169913tWXLlog+jEqSMjIy9PXXXzf7JOEf+Hw+XXbZZcrKytKiRYvUo0cPPfLII9q2bZv+93//V5J05ZVX6p577tHVV1+t3Nxc7du3TzfffHPYv5dRo0YpMzNTY8aM0UMPPaQBAwaotrZWxcXF6tSpU4u3KQMAji1ccQUAHBbdu3fXli1bdNlll+mBBx7QmWeeqX79+unWW2/VWWedFfp016ioKD399NP66KOPNGDAAF177bXKzs7WCSecEDbfzJkz1b9/f51//vm64IILNHz4cJ1xxhmh7fHx8dqwYYMuueQS9e7dWxMnTtRVV12le+65R9L3AXH27Nk666yzdNppp2n16tX6+9//rlNOOaXF+ufMmaMePXooPT1dGRkZOvHEEzV27Ngj8r2Kj49XXFxcRGMDgYCuvvpq3XzzzTr//PMlSZdffrmuvfZaXXHFFaqrq5MkzZo1S5s3b1bPnj1Dz5+63W49+eST2rhxowYMGKD777+/2ZXag7nnnns0e/ZsLV68WAMHDtSwYcM0Z84c9ezZ85DO1ev1qlOnTgfd/re//U3nn3++rr76ag0cOFDr16/X888/r1NPPVWSFB0drX/84x+qqanRkCFDdNVVV2nq1KlKTEwMzeFwOPTss89qzJgxmjp1qk499VRdeOGFeuGFF+Tz+Q6pXgCAmRzWTx8gAgAAAADAIFxxBQAAAAAYjeAKAAAAADAawRUAAAAAYDSCKwAAAADAaARXAAAAAIDRCK4AAAAAAKO52ruAQ1VRUdHeJRyU1+tVdXV1e5eBI4T+2he9tS96a2/0177orb3RX/s6HL1NSUlpcT1XXAEAAAAARiO4AgAAAACMRnAFAAAAABiN4AoAAAAAMBrBFQAAAABgNIIrAAAAAMBoBFcAAAAAgNEIrgAAAAAAoxFcAQAAAABGI7gCAAAAAIxGcAUAAAAAGI3gCgAAAAAwGsEVAAAAAGA0gisAAAAAwGgEVwAAAACA0QiuAAAAAACjEVwBAAAAAEYjuAIAAAAAjEZwBQAAAAAYjeAKAAAAADAawRUAAAAAYDSCKwAAAADAaARXAAAAAIDRCK4AAAAAAKMRXAEAAAAARiO4AgAAAACMRnAFAAAAABiN4AoAAAAAMBrBFQAAAABgNIIrAAAAAMBoBFcAAAAAgNEIrgAAAAAAo7kiGVRaWqr8/HwFg0FlZGRo9OjRYdvLysq0dOlSffLJJ8rOztbQoUMlSdu2bdPSpUtD4yoqKnTzzTdryJAhWrBggcrKyhQdHS1Jmjx5snr27HmYTgsAAAAAYBetBtdgMKglS5Zo2rRp8ng8ysnJkd/vV7du3UJjvF6vsrKy9Nxzz4Xt279/f82YMUOSVF9frylTpmjgwIGh7ePHjw+FXAAAAAAAWtJqcC0vL1dycrKSkpIkSWlpaSopKQkLromJiZIkh8Nx0Hk2btyowYMHq2PHjm2tGQAAAADwK9LqM66BQEAejye07PF4FAgEDvlA69ev19lnnx227sknn9Rtt92mxx9/XA0NDYc8JwAAAADA/iJ6xrWtamtr9emnn4bdJnzllVeqS5cuamxs1KJFi7RmzRqNHTu22b6FhYUqLCyUJOXl5cnr9R6Nkn8Rl8tldH1oG/prX/TWvuitvdFf+6K39kZ/7etI9rbV4Op2u1VTUxNarqmpkdvtPqSDbNiwQUOGDJHL9f+HS0hIkCQdd9xxSk9Pb/Z87A8yMzOVmZkZWq6urj6kYx9NXq/X6PrQNvTXvuitfdFbe6O/9kVv7Y3+2tfh6G1KSkqL61u9Vdjn86myslJVVVVqbGxUcXGx/H7/IR28pduEa2trJUmWZamkpESpqamHNCcAAAAA4Neh1SuuTqdTEydOVG5uroLBoNLT05WamqqCggL5fD75/X6Vl5dr5syZ2rNnjzZv3qyVK1dq9uzZkqSqqipVV1erb9++YfPOnTtXu3fvliT16NFDkyZNOgKnBwAAAAA41jksy7Lau4hDUVFR0d4lHBS3Pdgb/bUvemtf9Nbe6K990Vt7o7/21a63CgMAAAAA0J4IrgAAAAAAoxFcAQAAAABGI7gCAAAAAIxGcAUAAAAAGI3gCgAAAAAwGsEVAAAAAGA0gisAAAAAwGgEVwAAAACA0QiuAAAAAACjEVwBAAAAAEYjuAIAAAAAjEZwBQAAAAAYjeAKAAAAADAawRUAAAAAYDSCKwAAAADAaARXAAAAAIDRCK4AAAAAAKMRXAEAAAAARiO4AgAAAACMRnAFAAAAABiN4AoAAAAAMBrBFQAAAABgNIIrAAAAAMBoBFcAAAAAgNEIrgAAAAAAoxFcAQAAAABGI7gCAAAAAIxGcAUAAAAAGI3gCgAAAAAwmiuSQaWlpcrPz1cwGFRGRoZGjx4dtr2srExLly7VJ598ouzsbA0dOjS07fe//726d+8uSfJ6vbrjjjskSVVVVXr44YdVV1enk08+WVOmTJHLFVE5AAAAAIBfkVaTYjAY1JIlSzRt2jR5PB7l5OTI7/erW7duoTFer1dZWVl67rnnmu3foUMHzZgxo9n65cuX68ILL9TZZ5+txx57TGvXrtV5553XxtMBAAAAANhNq7cKl5eXKzk5WUlJSXK5XEpLS1NJSUnYmMTERPXo0UMOhyOig1qWpffeey90ZXbkyJHN5gQAAAAAQIrgimsgEJDH4wktezwe7dixI+IDNDQ06M4775TT6dQll1yiIUOGqK6uTtHR0XI6nZIkt9utQCDwC8oHAAAAANjdEX+odOHChXK73frqq6903333qXv37oqOjo54/8LCQhUWFkqS8vLy5PV6j1SpbeZyuYyuD21Df+2L3toXvbU3+mtf9Nbe6K99Hcnethpc3W63ampqQss1NTVyu90RH+CHsUlJSerbt68+/vhjnXnmmdq7d6+amprkdDoVCAQOOmdmZqYyMzNDy9XV1REf+2jzer1G14e2ob/2RW/ti97aG/21L3prb/TXvg5Hb1NSUlpc3+ozrj6fT5WVlaqqqlJjY6OKi4vl9/sjOmh9fb0aGhokSbt379YHH3ygbt26yeFwqF+/ftq4caMkqaioKOI5AQAAAAC/Lq1ecXU6nZo4caJyc3MVDAaVnp6u1NRUFRQUyOfzye/3q7y8XDNnztSePXu0efNmrVy5UrNnz9YXX3yhxx57TFFRUQoGgxo9enTo04ivuuoqPfzww3rqqad00kknadSoUUf8ZAEAAAAAxx6HZVlWexdxKCoqKtq7hIPitgd7o7/2RW/ti97aG/21L3prb/TXvtr1VmEAAAAAANoTwRUAAAAAYDSCKwAAAADAaARXAAAAAIDRCK4AAAAAAKMRXAEAAAAARiO4AgAAAACMRnAFAAAAABiN4AoAAAAAMBrBFQAAAABgNIIrAAAAAMBoBFcAAAAAgNEIrgAAAAAAoxFcAQAAAABGI7gCAAAAAIxGcAUAAAAAGI3gCgAAAAAwGsEVAAAAAGA0gisAAAAAwGgEVwAAAACA0QiuAAAAAACjEVwBAAAAAEYjuAIAAAAAjEZwBQAAAAAYjeAKAAAAADAawRUAAAAAYDSCKwAAAADAaARXAAAAAIDRCK4AAAAAAKMRXAEAAAAARnNFMqi0tFT5+fkKBoPKyMjQ6NGjw7aXlZVp6dKl+uSTT5Sdna2hQ4dKkj7++GMtXrxY+/btU1RUlMaMGaO0tDRJ0oIFC1RWVqbo6GhJ0uTJk9WzZ8/DeGoAAAAAADtoNbgGg0EtWbJE06ZNk8fjUU5Ojvx+v7p16xYa4/V6lZWVpeeeey5s3w4dOuhPf/qTTjjhBAUCAd15550aOHCgYmJiJEnjx48PhVwAAAAAAFrSanAtLy9XcnKykpKSJElpaWkqKSkJC66JiYmSJIfDEbZvSkpK6Gu32634+Hjt3r07FFwBAAAAAGhNq8+4BgIBeTye0LLH41EgEDjkA5WXl6uxsTEUgCXpySef1G233abHH39cDQ0NhzwnAAAAAMD+InrGta1qa2s1b948TZ48WVFR32flK6+8Ul26dFFjY6MWLVqkNWvWaOzYsc32LSwsVGFhoSQpLy9PXq/3aJT8i7hcLqPrQ9vQX/uit/ZFb+2N/toXvbU3+mtfR7K3rQZXt9utmpqa0HJNTY3cbnfEB9i7d6/y8vJ0xRVXqE+fPqH1CQkJkqTjjjtO6enpzZ6P/UFmZqYyMzNDy9XV1REf+2jzer1G14e2ob/2RW/ti97aG/21L3prb/TXvg5Hb3/8uOmPtXqrsM/nU2VlpaqqqtTY2Kji4mL5/f6IDtrY2KiZM2dq+PDhzT6Eqba2VpJkWZZKSkqUmpoa0ZwAAAAAgF+XVq+4Op1OTZw4Ubm5uQoGg0pPT1dqaqoKCgrk8/nk9/tVXl6umTNnas+ePdq8ebNWrlyp2bNnq7i4WNu3b1ddXZ2Kiook/f+fvZk7d652794tSerRo4cmTZp0RE8UAAAAAHBscliWZbV3EYeioqKivUs4KG57sDf6a1/01r7orb3RX/uit/ZGf+2rXW8VBgAAAACgPRFcAQAAAABGI7gCAAAAAIxGcAUAAAAAGI3gCgAAAAAwGsEVAAAAAGA0gisAAAAAwGgEVwAAAACA0QiuAAAAAACjEVwBAAAAAEYjuAIAAAAAjEZwBQAAAAAYjeAKAAAAADAawRUAAAAAYDSCKwAAAADAaARXAAAAAIDRCK4AAAAAAKMRXAEAAAAARiO4AgAAAACMRnAFAAAAABiN4AoAAAAAMBrBFQAAAABgNIIrAAAAAMBoBFcAAAAAgNEIrgAAAAAAoxFcAQAAAABGI7gCAAAAAIxGcAUAAAAAGI3gCgAAAAAwGsEVAAAAAGA0VySDSktLlZ+fr2AwqIyMDI0ePTpse1lZmZYuXapPPvlE2dnZGjp0aGhbUVGRnnnmGUnSmDFjNHLkSEnSzp07tWDBAh04cECDBw/WhAkT5HA4DtNpAQAAAADsotUrrsFgUEuWLNFdd92lOXPmaP369fr888/Dxni9XmVlZWnYsGFh6+vr67Vq1Sr99a9/1V//+letWrVK9fX1kqTFixfrhhtu0Ny5c/Xll1+qtLT0MJ4WAAAAAMAuWg2u5eXlSk5OVlJSklwul9LS0lRSUhI2JjExUT169Gh2xbS0tFQDBgxQbGysYmNjNWDAAJWWlqq2tlb79u1Tnz595HA4NHz48GZzAgAAAAAgRRBcA4GAPB5PaNnj8SgQCEQ0+U/3dbvdCgQCbZoTAAAAAPDrEtEzru2psLBQhYWFkqS8vDx5vd52rujgXC6X0fWhbeivfdFb+6K39kZ/7Yve2hv9ta8j2dtWg6vb7VZNTU1ouaamRm63O6LJ3W63ysrKQsuBQEB9+/Y9pDkzMzOVmZkZWq6uro7o2O3B6/UaXR/ahv7aF721L3prb/TXvuitvdFf+zocvU1JSWlxfau3Cvt8PlVWVqqqqkqNjY0qLi6W3++P6KCDBg3S1q1bVV9fr/r6em3dulWDBg1SQkKCOnfurA8//FCWZWndunURzwkAAAAA+HVp9Yqr0+nUxIkTlZubq2AwqPT0dKWmpqqgoEA+n09+v1/l5eWaOXOm9uzZo82bN2vlypWaPXu2YmNjdemllyonJ0eSNHbsWMXGxkqSrr/+ei1cuFAHDhzQoEGDNHjw4CN7pgAAAACAY5LDsiyrvYs4FBUVFe1dwkFx24O90V/7orf2RW/tjf7aF721N/prX+16qzAAAAAAAO2J4AoAAAAAMBrBFQAAAABgNIIrAAAAAMBoBFcAAAAAgNEIrgAAAAAAoxFcAQAAAABGI7gCAAAAAIxGcAUAAAAAGI3gCgAAAAAwGsEVAAAAAGA0gisAAAAAwGgEVwAAAACA0QiuAAAAAACjEVwBAAAAAEYjuAIAAAAAjEZwBQAAAAAYjeAKAAAAADAawRUAAAAAYDSCKwAAAADAaARXAAAAAIDRCK4AAAAAAKMRXAEAAAAARiO4AgAAAACMRnAFAAAAABiN4AoAAAAAMBrBFQAAAABgNIIrAAAAAMBoBFcAAAAAgNEIrgAAAAAAo7kiGVRaWqr8/HwFg0FlZGRo9OjRYdsbGho0f/587dy5U3FxccrOzlZiYqL++c9/6tlnnw2N+/TTT/Xggw+qZ8+emj59umpra9WhQwdJ0rRp0xQfH38YTw0AAAAAYAetBtdgMKglS5Zo2rRp8ng8ysnJkd/vV7du3UJj1q5dq5iYGM2bN0/r16/XihUrNHXqVJ1zzjk655xzJH0fWmfMmKGePXuG9rvpppvk8/kO/1kBAAAAAGyj1VuFy8vLlZycrKSkJLlcLqWlpamkpCRszKZNmzRy5EhJ0tChQ7Vt2zZZlhU25l//+pfS0tIOX+UAAAAAgF+FVoNrIBCQx+MJLXs8HgUCgYOOcTqdio6OVl1dXdiYDRs26Oyzzw5bt3DhQt1+++1atWpVs6ALAAAAAIAU4TOubbVjxw516NBB3bt3D6276aab5Ha7tW/fPs2aNUvr1q3TiBEjmu1bWFiowsJCSVJeXp68Xu/RKPkXcblcRteHtqG/9kVv7Yve2hv9tS96a2/0176OZG9bDa5ut1s1NTWh5ZqaGrnd7hbHeDweNTU1ae/evYqLiwttX79+fbOrrT/M0blzZw0bNkzl5eUtBtfMzExlZmaGlqurqyM8taPP6/UaXR/ahv7aF721L3prb/TXvuitvdFf+zocvU1JSWlxfau3Cvt8PlVWVqqqqkqNjY0qLi6W3+8PG3P66aerqKhIkrRx40b169dPDodD0vcf7vTT24Sbmpq0e/duSVJjY6M2b96s1NTUX3RiAAAAAAB7a/WKq9Pp1MSJE5Wbm6tgMKj09HSlpqaqoKBAPp9Pfr9fo0aN0vz58zVlyhTFxsYqOzs7tP/27dvl9XqVlJQUWtfQ0KDc3Fw1NTUpGAzqtNNOC7uqCgAAAADADxzWMfapSBUVFe1dwkFx24O90V/7orf2RW/tjf7aF721N/prX+16qzAAAAAAAO2J4AoAAAAAMBrBFQAAAABgNIIrAAAAAMBoBFcAAAAAgNEIrgAAAAAAoxFcAQAAAABGI7gCAAAAAIxGcAUAAAAAGI3gCgAAAAAwGsEVAAAAAGA0gisAAAAAwGgEVwAAAACA0QiuAAAAAACjEVwBAAAAAEYjuAIAAGI4P/8AABLGSURBVAAAjEZwBQAAAAAYjeAKAAAAADAawRUAAAAAYDSCKwAAAADAaARXAAAAAIDRCK4AAAAAAKMRXAEAAAAARiO4AgAAAACMRnAFAAAAABiN4AoAAAAAMBrBFQAAAABgNIIrAAAAAMBoBFcAAAAAgNEIrgAAAAAAo7kiGVRaWqr8/HwFg0FlZGRo9OjRYdsbGho0f/587dy5U3FxccrOzlZiYqKqqqo0depUpaSkSJJ69+6tSZMmSZJ27typBQsW6MCBAxo8eLAmTJggh8NxmE8PAAAAAHCsazW4BoNBLVmyRNOmTZPH41FOTo78fr+6desWGrN27VrFxMRo3rx5Wr9+vVasWKGpU6dKkpKTkzVjxoxm8y5evFg33HCDevfurQceeEClpaUaPHjwYTw1AAAAAIAdtHqrcHl5uZKTk5WUlCSXy6W0tDSVlJSEjdm0aZNGjhwpSRo6dKi2bdsmy7IOOmdtba327dunPn36yOFwaPjw4c3mBAAAAABAiuCKayAQkMfjCS17PB7t2LHjoGOcTqeio6NVV1cnSaqqqtKf//xnde7cWePGjdNvfvObFucMBAKH5YQAAAAAAPYS0TOuv1RCQoIWLlyouLg47dy5UzNmzNCsWbMOaY7CwkIVFhZKkvLy8uT1eo9EqYeFy+Uyuj60Df21L3prX/TW3uivfdFbe6O/9nUke9tqcHW73aqpqQkt19TUyO12tzjG4/GoqalJe/fuVVxcnBwOh4477jhJ0sknn6ykpCRVVlZGNOcPMjMzlZmZGVqurq4+tDM8irxer9H1oW3or33RW/uit/ZGf+2L3tob/bWvw9HbHz7Y96dafcbV5/OpsrJSVVVVamxsVHFxsfx+f9iY008/XUVFRZKkjRs3ql+/fnI4HNq9e7eCwaAk6auvvlJlZaWSkpKUkJCgzp0768MPP5RlWVq3bl2zOQEAAAAAkCK44up0OjVx4kTl5uYqGAwqPT1dqampKigokM/nk9/v16hRozR//nxNmTJFsbGxys7OliSVlZVp5cqVcjqdioqK0h//+EfFxsZKkq6//notXLhQBw4c0KBBg/hEYQAAAABAixzWz338r4EqKirau4SD4rYHe6O/9kVv7Yve2hv9tS96a2/0177a9VZhAAAAAADaE8EVAAAAAGA0gisAAAAAwGgEVwAAAACA0QiuAAAAAACjEVwBAAAAAEYjuAIAAAAAjEZwBQAAAAAYjeAKAAAAADAawRUAAAAAYDSCKwAAAADAaARXAAAAAIDRCK4AAAAAAKMRXAEAAAAARiO4AgAAAACMRnAFAAAAABiN4AoAAAAAMBrBFQAAAABgNIIrAAAAAMBoBFcAAAAAgNEIrgAAAAAAoxFcAQAAAABGI7gCAAAAAIxGcAUAAAAAGI3gCgAAAAAwGsEVAAAAAGA0gisAAAAAwGgEVwAAAACA0QiuAAAAAACjEVwBAAAAAEZzRTKotLRU+fn5CgaDysjI0OjRo8O2NzQ0aP78+dq5c6fi4uKUnZ2txMREvfPOO1qxYoUaGxvlcrk0fvx49e/fX5I0ffp01dbWqkOHDpKkadOmKT4+/jCfHgAAAADgWNdqcA0Gg1qyZImmTZsmj8ejnJwc+f1+devWLTRm7dq1iomJ0bx587R+/XqtWLFCU6dOVVxcnO644w653W59+umnys3N1aJFi0L73XTTTfL5fEfmzAAAAAAAttDqrcLl5eVKTk5WUlKSXC6X0tLSVFJSEjZm06ZNGjlypCRp6NCh2rZtmyzL0kknnSS32y1JSk1N1YEDB9TQ0HD4zwIAAAAAYFutXnENBALyeDyhZY/Hox07dhx0jNPpVHR0tOrq6nT88ceHxrz55ps6+eSTddxxx4XWLVy4UFFRUTrzzDN16aWXyuFwtPmEAAAAAAD2EtEzrm312WefacWKFbr77rtD62666Sa53W7t27dPs2bN0rp16zRixIhm+xYWFqqwsFCSlJeXJ6/XezRK/kVcLpfR9aFt6K990Vv7orf2Rn/ti97aG/21ryPZ21aDq9vtVk1NTWi5pqYmdPvvT8d4PB41NTVp7969iouLC42fOXOmJk+erOTk5LB9JKlz584aNmyYysvLWwyumZmZyszMDC1XV1cf4ikePV6v1+j60Db0177orX3RW3ujv/ZFb+2N/trX4ehtSkpKi+tbfcbV5/OpsrJSVVVVamxsVHFxsfx+f9iY008/XUVFRZKkjRs3ql+/fnI4HNqzZ4/y8vJ05ZVX6tRTTw2Nb2pq0u7duyVJjY2N2rx5s1JTU3/puQEAAAAAbKzVK65Op1MTJ05Ubm6ugsGg0tPTlZqaqoKCAvl8Pvn9fo0aNUrz58/XlClTFBsbq+zsbEnSSy+9pC+//FKrVq3SqlWrJH3/Z286duyo3NxcNTU1KRgM6rTTTgu7qgoAAAAAwA8clmVZ7V3EoaioqGjvEg6K2x7sjf7aF721L3prb/TXvuitvdFf+2rXW4UBAAAAAGhPBFcAAAAAgNEIrgAAAAAAoxFcAQAAAABGI7gCAAAAAIxGcAUAAAAAGI3gCgAAAAAwGsEVAAAAAGA0gisAAAAAwGgEVwAAAACA0QiuAAAAAACjEVwBAAAAAEYjuAIAAAAAjEZwBQAAAAAYjeAKAAAAADAawRUAAAAAYDSCKwAAAADAaARXAAAAAIDRCK4AAAAAAKMRXAEAAAAARiO4AgAAAACMRnAFAAAAABiN4AoAAAAAMBrBFQAAAABgNIIrAAAAAMBoBFcAAAAAgNEIrgAAAAAAoxFcAQAAAABGI7gCAAAAAIxGcAUAAAAAGI3gCgAAAAAwmiuSQaWlpcrPz1cwGFRGRoZGjx4dtr2hoUHz58/Xzp07FRcXp+zsbCUmJkqSVq9erbVr1yoqKkoTJkzQoEGDIpoTAAAAAAApgiuuwWBQS5Ys0V133aU5c+Zo/fr1+vzzz8PGrF27VjExMZo3b54uvPBCrVixQpL0+eefq7i4WLNnz9bdd9+tJUuWKBgMRjQnAAAAAABSBMG1vLxcycnJSkpKksvlUlpamkpKSsLGbNq0SSNHjpQkDR06VNu2bZNlWSopKVFaWpqOO+44JSYmKjk5WeXl5RHNCQAAAACAFMGtwoFAQB6PJ7Ts8Xi0Y8eOg45xOp2Kjo5WXV2dAoGAevfuHRrndrsVCARC8/zcnD8oLCxUYWGhJCkvL09erzfSczvqXC6X0fWhbeivfdFb+6K39kZ/7Yve2hv9ta8j2duInnFtT5mZmcrMzAwtV1dXt2M1P8/r9RpdH9qG/toXvbUvemtv9Ne+6K290V/7Ohy9TUlJaXF9q7cKu91u1dTUhJZramrkdrsPOqapqUl79+5VXFxcs30DgYDcbndEcwIAAAAAIEUQXH0+nyorK1VVVaXGxkYVFxfL7/eHjTn99NNVVFQkSdq4caP69esnh8Mhv9+v4uJiNTQ0qKqqSpWVlerVq1dEcwIAAAAAIEVwq7DT6dTEiROVm5urYDCo9PR0paamqqCgQD6fT36/X6NGjdL8+fM1ZcoUxcbGKjs7W5KUmpqqs846S7fccouioqJ03XXXKSrq+6zc0pwAAAAAAPyUw7Isq72LOBQVFRXtXcJBcb++vdFf+6K39kVv7Y3+2he9tTf6a1/t+owrAAAAAADtieAKAAAAADAawRUAAAAAYDSCKwAAAADAaARXAAAAAIDRCK4AAAAAAKMRXAEAAAAARiO4AgAAAACMRnAFAAAAABjNYVmW1d5FAAAAAABwMFxxPYzuvPPO9i4BRxD9tS96a1/01t7or33RW3ujv/Z1JHtLcAUAAAAAGI3gCgAAAAAwmnP69OnT27sIOzn55JPbuwQcQfTXvuitfdFbe6O/9kVv7Y3+2teR6i0fzgQAAAAAMBq3CgMAAAAAjOZq7wKOZStXrtRrr72m448/XpJ0xRVX6Le//W2zcaWlpcrPz1cwGFRGRoZGjx59tEvFL7Bs2TJt3rxZLpdLSUlJysrKUkxMTLNxkydPVqdOnRQVFSWn06m8vLx2qBaRaO212NDQoPnz52vnzp2Ki4tTdna2EhMT26laRKq6uloLFizQN998I4fDoczMTP3ud78LG/Pee+/poYceCvXzzDPP1NixY9ujXPwCrf2ctSxL+fn52rJlizp27KisrCxuQzwGVFRUaM6cOaHlqqoqXX755brwwgtD63jtHlsWLlyot99+W/Hx8Zo1a5Ykqb6+XnPmzNHXX3+trl27aurUqYqNjW22b1FRkZ555hlJ0pgxYzRy5MijWTpa0VJvj/p7ZQu/WEFBgbVmzZqfHdPU1GT96U9/sr788kuroaHBuu2226zPPvvsKFWItigtLbUaGxsty7KsZcuWWcuWLWtxXFZWlvXtt98ezdLwC0TyWnzppZesRYsWWZZlWf/617+s2bNnt0epOESBQMD66KOPLMuyrL1791o33XRTs95u27bNeuCBB9qjPBwGrf2c3bx5s5Wbm2sFg0Hrgw8+sHJyco5idTgcmpqarOuvv96qqqoKW89r99jy3nvvWR999JF1yy23hNYtW7bMWr16tWVZlrV69eoW30/V1dVZkydPturq6sK+hjla6u3Rfq/MrcJHWHl5uZKTk5WUlCSXy6W0tDSVlJS0d1mIwMCBA+V0OiVJffr0USAQaOeK0BaRvBY3bdoU+g3v0KFDtW3bNll8DIDxEhISQlfXOnfurBNPPJHX66/Mpk2bNHz4cDkcDvXp00d79uxRbW1te5eFQ/Duu+8qOTlZXbt2be9S0AZ9+/ZtdjW1pKREI0aMkCSNGDGixffBpaWlGjBggGJjYxUbG6sBAwaotLT0qNSMyLTU26P9Xplbhdvo5Zdf1rp163TyySfrmmuuadbQQCAgj8cTWvZ4PNqxY8fRLhNttHbtWqWlpR10e25uriTp3HPPVWZm5tEqC4cgktfij8c4nU5FR0errq4u9DgAzFdVVaVdu3apV69ezbZ9+OGHuv3225WQkKDx48crNTW1HSrEL/VzP2cDgYC8Xm9o2ePxKBAIKCEh4ajWiF9u/fr1Ovvss1vcxmv32Pbtt9+GXotdunTRt99+22zMT/8f7Xa7+QXkMeZovFcmuLbi/vvv1zfffNNs/bhx43TeeeeFnrMoKCjQE088oaysrKNdItrg5/p7xhlnSJKeeeYZOZ1OnXPOOQedw+1269tvv9Vf/vIXpaSkqG/fvke0bgDN7d+/X7NmzdK1116r6OjosG0nnXSSFi5cqE6dOuntt9/WjBkzNHfu3HaqFIeKn7P21tjYqM2bN+vKK69sto3Xrr04HA45HI72LgOH2dF6r0xwbcU999wT0biMjAw9+OCDzda73W7V1NSElmtqauR2uw9bfWib1vpbVFSkzZs369577z3oD9of+hkfH68zzjhD5eXlvKEyUCSvxR/GeDweNTU1ae/evYqLizvapeIXaGxs1KxZs3TOOefozDPPbLb9x0H2t7/9rZYsWaLdu3dzNf0Y0drPWbfbrerq6tAy/689tmzZskUnnXSSunTp0mwbr91jX3x8vGpra5WQkKDa2toWe+d2u1VWVhZaDgQCvJc6RhzN98o849oGP35+5q233mrx1hWfz6fKykpVVVWpsbFRxcXF8vv9R7NM/EKlpaVas2aN7rjjDnXs2LHFMfv379e+fftCX7/zzjvq3r370SwTEYrktXj66aerqKhIkrRx40b169eP3wwfAyzL0qOPPqoTTzxRF110UYtjvvnmm9DzyuXl5QoGg/xS4hgRyc9Zv9+vdevWybIsffjhh4qOjuY24WPIz90mzGv32Of3+/XGG29Ikt54443QHW0/NmjQIG3dulX19fWqr6/X1q1bNWjQoKNdKg7R0X6v7LD45JFfbN68efr444/lcDjUtWtXTZo0SQkJCQoEAlq0aJFycnIkSW+//baWLl2qYDCo9PR0jRkzpp0rRySmTJmixsbG0HPLvXv31qRJk8L6+9VXX2nmzJmSpKamJg0bNoz+Gqyl12JBQYF8Pp/8fr8OHDig+fPna9euXYqNjVV2draSkpLau2y04v3339e9996r7t27h37RcMUVV4SuwJ133nl66aWX9Morr8jpdKpDhw665pprdMopp7Rn2YjQwX7OvvLKK5K+769lWVqyZIm2bt2qDh06KCsrSz6frz3LRoT279+vrKwszZ8/P3R19ce95bV7bHn44YdVVlamuro6xcfH6/LLL9cZZ5yhOXPmqLq6OuzP4Xz00Ud69dVXdeONN0r6/hnJ1atXS/r+z+Gkp6e356ngJ1rq7erVq4/qe2WCKwAAAADAaNwqDAAAAAAwGsEVAAAAAGA0gisAAAAAwGgEVwAAAACA0QiuAAAAAACjEVwBAAAAAEYjuAIAAAAAjEZwBQAAAAAY7f8A8th2umns1WUAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 1152x432 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2N6Pv_29pS3O"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}