{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 4-6: Mixture Models+Model orden selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this lab session is to study mixture models. In the first part you will code the EM algorithm to estimate the parameters of a GMM given the number of mixed distributions and in the second part you will try different model order selection methods. You will send only one notebook for both parts.\n",
    "\n",
    "You have to send the filled notebook named **\"L4_6_familyname1_familyname2.ipynb\"** (groups of 2) by email to aml.centralesupelec.2019@gmail.com before November 28 at 23:59 and put **\"AML-L4-6\"** in the subject. \n",
    "\n",
    "We begin with the standard imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline\n",
    "sns.set_context('poster')\n",
    "sns.set_color_codes()\n",
    "plot_kwds = {'alpha' : 0.25, 's' : 80, 'linewidths':0}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will have two toy datasets to try the different methods:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GMM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Gaussian mixture model is a probabilistic model that assumes all the data points are generated from a mixture of a finite number of Gaussian distributions with unknown parameters. After estimation of those parameters we get an estimation of the distribution of our data. For the clustering task, one can think of mixture models as generalizing k-means clustering to incorporate information about the covariance structure of the data as well as the centers of the latent Gaussians. \n",
    "\n",
    "### First part\n",
    "\n",
    "Fill in the following class to implement a multivariate GMM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class my_GMM():\n",
    "    \n",
    "    def __init__(self, k):\n",
    "        '''\n",
    "        Parameters:\n",
    "        k: integer\n",
    "            number of components\n",
    "        \n",
    "        Attributes:\n",
    "        \n",
    "        mu_: np.array\n",
    "            array containing means\n",
    "        Sigma_: np.array\n",
    "            array cointaining covariance matrix\n",
    "        cond_prob_: (n, K) np.array\n",
    "            conditional probabilities for all data points \n",
    "        labels_: (n, ) np.array\n",
    "            labels for data points\n",
    "        '''\n",
    "        self.mu_ = None\n",
    "        self.Sigma_ = None\n",
    "        self.cond_prob_ = None\n",
    "        self.labels_ = None\n",
    "        \n",
    "    def fit(self, X):\n",
    "        \"\"\" Find the parameters mu_ and Sigma_\n",
    "        that better fit the data\n",
    "        \n",
    "        Parameters:\n",
    "        -----------\n",
    "        X: (n, p) np.array\n",
    "            Data matrix\n",
    "        \n",
    "        Returns:\n",
    "        -----\n",
    "        self\n",
    "        \"\"\"\n",
    "        def compute_condition_prob_matrix(X, mu, Sigma):\n",
    "            '''Compute the conditional probability matrix \n",
    "            shape: (n, K)\n",
    "            '''\n",
    "        \n",
    "        # TODO:\n",
    "        # initialize the parameters\n",
    "        # apply sklearn kmeans or randomly initialize them\n",
    "        \n",
    "        # While not(convergence)\n",
    "        #     Compute conditional probability matrix\n",
    "        #     Update parameters\n",
    "        \n",
    "        # Update labels_\n",
    "        \n",
    "        # Return self\n",
    "        \n",
    "    def predict(self, X):\n",
    "        \"\"\" Predict labels for X\n",
    "        \n",
    "        Parameters:\n",
    "        -----------\n",
    "        X: (n, p) np.array\n",
    "            New data matrix\n",
    "        \n",
    "        Returns:\n",
    "        -----\n",
    "        label assigment        \n",
    "        \"\"\"\n",
    "        # TODO\n",
    "        \n",
    "    def compute_proba(self, X):\n",
    "        \"\"\" Compute probability vector for X\n",
    "        \n",
    "        Parameters:\n",
    "        -----------\n",
    "        X: (n, p) np.array\n",
    "            New data matrix\n",
    "        \n",
    "        Returns:\n",
    "        -----\n",
    "        proba: (n, k) np.array        \n",
    "        \"\"\"\n",
    "        # TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate your own mixture of Gaussian distributions to test the model, choose parameters so that GMM performs better than K-Means on it. Use `np.random.multivariate_normal`. \n",
    "\n",
    "Plot data with colors representing predicted labels and shapes representing real labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bonus (not graded): Implement a mixture of asymmetric generalized Gaussians (AGGD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Second Part\n",
    " \n",
    "- Implement the information criterions from the lecture (AIC, BIC, etc.) to select the number of clusters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Implement the merge criterions \n",
    "    - Correlation coefficients\n",
    "    - Measuring Error \n",
    "    - Comparing the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Implement cross-validation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the model selection criterions to choose the number of clusters for the two given datasets (data-MM-i.csv). Compare the results and the computational time. Try to visually validate your results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You are going to work with the following data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape: (1797, 64)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAHZCAYAAAC1llPzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3cFrVdma9/Enb6TKogKnQgVKMGDLDZRgwMgVSugCdVaDBu/9C2znF2xnd2Y7qHE68wKvM2f3Cs69Ql+wwItKWU0aUthCBBsUK9URYojsd2DyVtzr9+Q8e5+1snLe/n4mUqtOdp69ztp75Zz17GdNNE1jAACgnv9TOwAAAP63YzIGAKAyJmMAACpjMgYAoDImYwAAKmMyBgCgMiZjAAAqYzIGAKAyJmMAACpjMgYAoDImYwAAKmMyBgCgMiZjAAAqYzIGAKAyJmMAACpjMgYAoDImYwAAKjuU60ATExMPzey4ma2b2Uqu4xYwZ2ZTZvZ0+7+JuYxxi/n/xds0zekxGc/j1sdmxLxfiLm8D+4Zox5somma0UMys4mJiZ/NbJDlYPtjbftfYi5r3GJea5rmszEbz+PWx2bEvF+Iuby1pmk+G/Ug2T4Z2/u/YAaDwcAWFhb2fOHLly9l+3/9138lbZ99lp7jP/zDPyRthw7FTuXRo0e2trZm9j5es2DMex2vbWtrK2lTMc/MzIR/R86Yf/7556TtyZMnSdvU1FTSFv2dfWNeXV2V7WpsHD58OBRfZGyIeMPj2aPGwfLyctI2Pz/f6/h9+1iNWTPdnydOnOgV216/u8b1d+bMmV7H3/kdOceyik/dE9+8eZO0TU5OymOePXv2g//+4Ycf7Jdffukc88qK/hCq4jty5EjSNjs7m7SVvjere5eZ7ue+40wR8Y4k52S8YmZHFxYW7K9//eueL/zLX/4i2//1X/81aTt//nzodWrSVs6fP2/37t0z+/Wrj1DMex2vTU12169fT9ouXrwY/h05Y1Y/97vf/S5pUwM3+jv7xvxv//Zvsl295+oPnDt37iRtkbEh4g2PZ48aB//8z/+ctHnXwzB9+1iNWTPdn3/60596xbbX765x/fU9/s7vyDmWVXxqDDx+/DhpU38gm6Xj/p/+6Z/sb3/7m1nHmP/lX/5Ftqv41FhWP1/63qzuXWb5x0GbiHckJHABAFAZkzEAAJXl/Jo6TH3laKbXBdVXDV2+TvO+wshFfQWz/dXFB+7evZu0Rb+m7stbG7xw4ULSNhikuRLq/chJjQPvK1v1lZ/6Skyds/e1bGlqTOZcs+rLe1/VuL1582bSduzYsfAxc7l9+7ZsVzFfu3ataCyjUvcMNb5Vm7ofqmNG12nbvHuGosa3+ho451fDapx5Y0OZmJhI2k6dOpW0demHXPhkDABAZUzGAABUxmQMAEBlTMYAAFRWPIFLLYR7yR7quTqVLKKSsrwF91wJXN7xo8kJNRJ3vGQolbCg+kk9G53T5cuXkzbvOUeVhKUS+Woka3lJNSrBRZ1fl+Qndc5dec99Pnv2LGlTiX3RZ3v3+l1ddUnKKp20GeWNZUUlM6pxkTMZSvHuU9GkWfV+ezH3uVa9caacO3cuaVPnUbpPo/hkDABAZUzGAABUxmQMAEBlTMYAAFTGZAwAQGXFs6lV9puXsacyp5XSmcmqDJ1XwnN7C62hamT5etmcKqNQvbZ0uU71fquMXjOdWRrN6s2V0evxSrGqmEfd6cYbh114GdnqaQY1vtX1V7qPvSxa9WRAjScXRi0D6e3w1OY9IaHGVR/ecU6fPp20jbLlbV9djqX6SmXad8nQLolPxgAAVMZkDABAZUzGAABUxmQMAEBlVRK4Rk1mKp2koxJqvMSG6enp0DFLJwmo43tJIV4SSJuXmFSSl8SnkkVUMoZq8863z5hRe6devXpVvvbSpUuhYy4tLSVtN27c6BZYB15/qIQjVQbWO1+lS0nIvXjXj0roUeNejYvSiUWjltBV71PpRNAu9ym1l/TTp0+Ttpz9rK5ZlcRnpu/NV65cSdq6lGzOeS5tfDIGAKAyJmMAACpjMgYAoDImYwAAKiuewKUW3L3EBkUlFKifPyh7mHpUzDkrBanKTCoxyKOSRUpXVepCxaJiVglDXiJbn2pWan9f1WZmdvPmzaQtOvZrjOdRkoO67Mnch5c4o5KI1D1DJZ09fPhQHrPPdani8xLlJiYmQq8tnaylxuKFCxfka9V+0tGkSq8fciVDedfUKPdcL/EwmvzaB5+MAQCojMkYAIDKmIwBAKiMyRgAgMqKJ3AdP348afMW3NXieHTBPFeln3GlKoR5lX7UVnkq8UJtoXj58mV5zFzbLXpJVdHtEtU550yGisZhpse5+nlVqatk8pyqImamE9GiSW6lE868CngqMUslBqlkI+/ekiux0rsnqX4+d+5clt/ZheonLxlRnYvqU7XVolfJL8d2oHtR76M6DxVfyUQtD5+MAQCojMkYAIDKmIwBAKiMyRgAgMqYjAEAqKx4NrXan9bLolPtKiMuuh9oTl52q8oiVtmqKmYvQ7QP1U9dysSpvlfn4ZWwy5VN7fVzNFteZfV65TBLU+eytraWtOUcBxF3796V7dHyqSr7u3TpRq+PVEavyo5V8ZXOAPfuU6pMao3Ss+p3eu+j2htYZV6r+0DpJ12846v7XPQJjJyliqP4ZAwAQGVMxgAAVMZkDABAZRNN0+Q50MTEqpkdHQwGQ79vf/HihWxfXV1N2qamppK2EydO9IrR7P06wva63fPtplDMnidPniRtr169Stq++OKLpC16HrljXl9fT9rU2ps6j6NHj8pjzs3NffDffWNWY8DMHzNtah2sHZuyO96maWa7jGfPxsZG0vb9998nbadOnUraImuIfft4ZWVFtj9//ly2t9UYy1tbW7JdjVs1VlR/evkP7XtO35gfPHgg29XvnZmZ2fNYXfWNWd3PzPS9YHJyMmlT/eyNjUOHPkxZyj2e1X1OjSMVn5p32tr3jKE/METOyfhnM9O11A6mnUwaYi5r3GJea5rmszEbz+PWx2bEvF+Iuby1pmlGzsDLmU390fa/78ws/ZPk4Jgys0n7NV4zYi5h3GJuxzsO43nc+tiMmPcLMZen4u0t52S8aWaf2PvgxuEvms3tf4m5rHGLeXPXv+MW87jEa0bM+4WYy9sc/pLhck7G62Y2qLX+Mz8/Hwpy1/f8O39thWLu4v79+0lbe33EzH+WbY+1lE4xv3z5UrardVnVfyrmqEjMak21y5qxik+tvR05ckQec/e6kIg3PJ49ajyr8zt79mzSFun7SB+rMeD1sbou37x5MzQOM7OvvvpKth8+fLhzzF2U7mOz/jF79zkVn2pTY7nH2nynmJeXl2W7WkON3psjORtm/WP28knU2FD91/cZbxHvSHJOxitmdnRhYWFoUQ5vyzlVeEI9wP/1118nbdEtr86fP2/37t0zex+vWTDmLlSChnrD79y5I3++/dq+MXtb5S0uLiZtqv9GKUQQiVldLF6BDjUOVHyqkINXMGL3RS7iDY9njxrP6vzUOIj0faSP1RhQ77+Zvi7VdpvKrVu3ZHv7Wsh9/ZXuY7P+MXv3ORWfavvmm2+SNm87wra+MUeulb1iUUVDokV3+sbcZYtGFUvfgjUi3pHwaBMAAJUVL4epeH99qb/ir127lrSpv4S8v45KlxpUMT979izU5v3lnKs0nipb6B1f9V/pMnbqk7H3V7CKRfWfKuno9WeupQnvfVR96j1OEz1m17Fx48aNpG37r/mEKm+orj/1SSJ6Xrmp8aL6qHS5SVV6sUsJTxXfQSr7q85v1PtIrjHjfSuq7rnRT/M18MkYAIDKmIwBAKiMyRgAgMqYjAEAqKx4ApdKVvAeuVEJRyo9XSW3eHv3lnblypXQ686dO5e0lU568Y6vEkPUI0GlE7hU4oT3PqrECzU2VBJS6X1rvX5S41Qlm6j3yUsqiT7Ct6PLPtfqtercauy968WsktG8R7dKUslCXoJgNNlL3TtL864V9UhQ9BHO0ve5Lv2s9pJW95EaCYl8MgYAoDImYwAAKmMyBgCgMiZjAAAqYzIGAKCy4tnUXTIvo6UrS2dzqixYL2NWZVHWoDIvvSxD1X81Mje7iGYRqwzKnJmRKqtUZWia6axeFcv2zi8fyLmLWJs3ZqMZwTXGSpenJUpnzysXL15M2o4dOyZfq54mUeNbnYfX97nGuDfuVMzq6ZfoRhY5efdm9dSI6if1812fWsiBT8YAAFTGZAwAQGVMxgAAVMZkDABAZcUTuGqVqRyFSpLwEidUkkaX0ni5qMQEVebNE91vuUYpRDOdOBUt35gzGaNL8pJKZlHnoZw+fTr8e/ai+qNLss/ly5ezxDEqb39n5fjx40nbqVOnkrbr16/Ln1fJWH2M+h6qxEBv/OXa+9hLflP9pxJua9wfvN8Z7RN1zl1KxubCJ2MAACpjMgYAoDImYwAAKmMyBgCgsuIJXF0WvFUloujexV2SlYZRMXvJAKoyjUoIGCWZpy8vUUfFovYBrpWspahzUeOgy3vn7Rm8l+j+2mY6cUyNcZUEmCuJSL2HXqU7b5/xtmi/59Tl+o7uMe69rk/fqzFw7do1+Vo1HlVilnqfalQXM4vvgZ4rkWw/qSRFb7yVrMzFJ2MAACpjMgYAoDImYwAAKmMyBgCgsipbKJ47d06+Vm059+c//zl0zNIJJB6V+KTUSIbythZbWlpK2tR5qJ/3zqOdbLKxsTE0PpX0cu/ePfna169fJ20qAU4lSOXc8k+dv7dtnDq/6enppK1PIllUlz5WW+Kpyks1rjUvMSia0KTGsroOzNLxEhnLalx4SUAqMUu9TzmTUhX1O72kT/Xag7LtqpdAGa3++PTp06TNS2bcfc6RcdEFn4wBAKiMyRgAgMqYjAEAqIzJGACAypiMAQCorHg2teKVFFMZjyojzsterUFllqoM1MePHydtXhZgrsxrr+yhyoJU56HeJy+2dkZw32xqlVHfhSpl6PVDaWo8q6z1kvGp60dlTZvpTPSS5f+68DK4VR+rLGSVOe2VvWxnFB8+fHh4gB2ocV8yo96jrmXvnqTi87LR95s3H1y9ejX08+p+7Y2N3X126FDe6ZNPxgAAVMZkDABAZUzGAABUNtE0TZ4DTUysmtnRwWAwtELP1taWbF9ZWUna1tfXk7YTJ04kbVNTU6E4Hz16tLM29ny7KRSzR52LWqd78+ZN0vaP//iP8pjttYi+Mau+M9Nrxqr/Xr58OTS2He11tpWVlZ1zdmNW68rLy8vy+Irqe7W+Nz8/P/RYu/u4aZrZLuPZo85F9amKL5I3EBkXag3wyZMn8njv3r1L2r766qukbZQ11NzX3+rqatKmxrc6t88//1wes/1+5I5Z9b/q07m5uV7HN+sfszc21L1ExTczM9Mx0l/1jVmNATOzn376KfR7P/3006TNG+O7554ffvjBfvnlF7Pte0bol+0h52T8s5nFakMeDDvZKsRc1rjFvNY0zWdjNp7HrY/NiHm/EHN5a03TjJx1mzMd7KPtf9+Zmf5IdjBMmdmk/RqvGTGXMG4xt+Mdh/E8bn1sRsz7hZjLU/H2lnMy3jSzT+x9cOPwF83m9r/EXNa4xby5699xi3lc4jUj5v1CzOVtDn/JcDkn43UzG4yylqnW2NR392o9bXY29pX9rnWJnSBCMXvUuuf3338f+lm1HmeWnnPfmL1dVZ49e5a0nTx5MmnLtP7jxqzWfL31H7X2qcaRWtNWOQZmH44jEW94PHeh8glUfJF12b7jwtvNJvp8u+rj0tefl2eizkW9Vq3L98gzyXLPUNflixcvQj/r/c5c9wwvZ0P1qXrPR6mR0DdmlWtkpu8ZR44cSdqiY7dNxDuSnJPxipkdXVhYcLc72+HdDFTxA7Wll3oA3dsuUP3s9hZyO+9gKGaPurCOHz8e+tlbt27J9vY5943Z24Lt+vXrSdu3336btHkPvkdEYlYXi9oW0UxvoafGkboZeMfcPY5EvOHx3IUau6pogbeVXftYfcaFV2AiWnhC9XHp669LMQr12u+++y5pi06mue8Z6rqMFjIqfc/wCtCoPlXv+SjFS/rG7I099XPq/KJjt03EOxIebQIAoDImYwAAKqtSm9r7+lTVb1Ztt2/fTtp+97vfyWNGvu4bhbcuu9/U10heXWH19bPqv1yPvXlU33lLGOrrL9WmvpryxlvOr58V9dWjOudctcgV1Z/bX62F2tVYqVFH2VtqUPcHVWu4ZB93FV16U9ev95Vq6Rri6lqJ3vu86yzXe+LdM9TYUPWq1b2v9Lyh8MkYAIDKmIwBAKiMyRgAgMqYjAEAqKx4ApdavFcJWGZmV65cSdpU8k3OIgyKSobykgS85KC2c+fOJW05kwRUMoTXTyqxSL1WnXPOvlfH6pKIooqXqJ/3np3MxUtQuXz5ctK2uLiYtKnkpOi4GkaNi2PHjsnXRp/bLk0lBqln4z19n9veL2o8RmsslH4/vH6KJpipn++SlJmTSuRTSV3R55FL45MxAACVMRkDAFAZkzEAAJUxGQMAUBmTMQAAlVUph+nxSt61qSzanFQ2piqjdtB5pfNUtrjKMjxIGajqPVfZkipDu3RmpNfP6ukA9dqJiYmkzev7ruei3ldPtBxmad4OTUrppxSiVAa4txOTyi5W56zGfOlz87L41XWlMrujpT5z8q4J9TSDos6ZbGoAAP4XYjIGAKAyJmMAACpjMgYAoLLiCVxdFu9VEoNKElBJG16yRJ+ygirJxjsPlXR28+bNpK3Gvsde6Up1fqpPD9IesKqEo0oWUeeW8zzUOPOSpFT/e/tut+VKIFFJbl6ClIpNJaFFEy376pKgqa6rGoloalx0KeEZlXMsq3HgjbtoIqBX+rIkL2bVrsbL8ePHkzbvPEqWYuaTMQAAlTEZAwBQGZMxAACVMRkDAFBZlQpcg8FAtqtkK5UsohIPSlem8Rbuo7+3RlUgLwlBJZuopLODTiXVqTEUTZqKUEkh09PT8rV//vOfk7b9TuRT485L4Iomt6gkuZzj29tvWVHJXtH3+8aNG7K9T/KcGnddkkdVzDX22fX25n748GHSpu4jKj4vubYGNU67JCmWPBc+GQMAUBmTMQAAlTEZAwBQGZMxAACVVUng8pIQ/vKXv4R+XiWg5EzS6SKauKKqAnnJPLmSYbr0iUrQUG3eMXNVpvESJ9R7rsZLl+33cvGqO6l2lQAS3eotFy9Jx2tvU+M2ZwKXOpaX1DXKdqpeMk6NLTdv376dtC0uLiZtOStwqWN5x1fJoKr/atyHvURVNZ7VPUONZ29c7U7K29jYCMUXxSdjAAAqYzIGAKAyJmMAACpjMgYAoDImYwAAKquSTa2yCc10VpzKiFNZfLX23lUlGdXewOrcSmdTe/2sMgqjmbRexnv757e2tkLHa+uSTa36fmlpqdfv3S9q7F67dm1fY/D6WI1RVSqwyx7lfahr2Rt3aoyr81DHrPUEhopP7TtdOqtb8X5ntNxnjZi9e1f06Rx1v42UPz58+HDo+FF8MgYAoDImYwAAKmMyBgCgsommafIcaGJi1cyODgaDodWYvMoly8vLSZtae5ybm0vaomvGjx49srW1NTOz59tNoZg96lzUeayvrydt8/Pz8pjtc+kb8+rqqmx/+fLlnj+3F28tuB3LDz/8YL/88otZx5jv378f/r3qPVdjI7K2s7uPm6aZ7TKeu4iuZ0byBvqOiydPnsh2NUZnZmaSNtXHUX1jVrGZma2srIRee+hQmh4zOzsrj9luz33PUGNAjW91fHUe3u/oE7OXx/LixYukTfXfkSNHkrbSMZe4z3kx775nt+8ZvX/ZtpyT8c9mpjcqPpjWtv8l5rLGLea1pmk+G7PxPG59bEbM+4WYy1trmmbkDOKc2dQfbf/7zsz0n7EHw5SZTdqv8ZoRcwnjFnM73nEYz+PWx2bEvF+IuTwVb285J+NNM/vE3gc3Dn/RbG7/S8xljVvMm7v+HbeYxyVeM2LeL8Rc3ubwlwyXczJeN7NB5Dt+taZqpp8lnZqaStrUepp6nbLre/6dv7ZCMXvr3Gr9Ta2zjfLscCRmFZ+3NvjmzZvesXz++eeyvb3+HYlZrUN56z+q/7x1xDZvbXD3upCINzyePWrNSq1xquN3XOd2+1j1kbfLjddPkdjUWqHS9/pTY8VMr3FG1+Cjz4n2jfnBgweyXd2r1Dp8dK1V6Ruzd/0p6j1R9xb1DLXZnrkxbsxqbd1b537+/HnS9umnnyZtauxGrgUR70hyTsYrZnZ0YWFhaAEJ78Fw9XNq8KiHzaM3zPPnz+9sZ7hzVwzF7L3hqnCAalMxR0ViVvF5RQ0eP37cO5avv/5atrcfsI/ErApgeAUprl+/nrSpLR4Vr/jJ7puBiDc8nj1qWzxVROPWrVtJW+SPt0gfq4nXK9oRLdigYov+bN/rz9vuUF1X6vzU66J/IPeN2bsnqXY17kcpZNQ3Zu/6U9R7ou4t3jHb71MkZvWBzbu3qgJA6o8eNXa9e8aQeEfCo00AAFTGZAwAQGXFa1Orr0Vu3rwpX6vWFqJfA3vrYLlqVntf76ivZVTbKF/tRUTjMDO7dOlS0vb73/8+aRsM0ryJnM/cqq/WvZhHqSPsfR2Zq//VV2dmuua06r9ctcgVFdv2OldCLQUo6jr1vvouXWddUfcC1e/eNd1njKslCW8sq/dEfeXb5Zxz8ZYDFBWf+nnv3tynxrk6vleD+u7du71/vkbf88kYAIDKmIwBAKiMyRgAgMqYjAEAqKx4AlcXaiFdJYCo13mJB7kW4qenp2W7SnKKxpwzgev169fh16oElWPHjoVel1OX5B71fG40GapPokgX288aJlQCT3TD81y6JDCqPlbXT8mEMzOd2OclnalkRHUvUGPFu2d0edZ2h7oPeNR1r35njSQi75pXMav+U+Mt531OHd9LoFTXmkoevnjx4uiBZcAnYwAAKmMyBgCgMiZjAAAqYzIGAKAyJmMAACqrUg7TE83SVBl1x48fD/+ePryMO7VjyNWrV5M2b9enXLySc4qKT7lx40bSljMzsgu1A4vKYPVKrZbkZbKr+KKlUnPpkk2t+lhlzEZ3V+urS8yqjGv0mBcuXAj/nmHUe6ieUDCLlx1VfV/6+vPGouorlcle+mkBdf7evU/dcxcXF5O2PtnzJfDJGACAypiMAQCojMkYAIDKmIwBAKjsQJXDjFIJJLUSi1TJumgJQS/xoE8yjDp/LxFGlYeMluPL2c/qWF4CiSqHqBJIVIJGznJ3qvSeV7YwGrN6P7xSjV3351b9qfrIO7aKQ7XlTILpco5eklSbet+6lJAdpsv1rdpV4pPqU29f7z77tqs4uuxn3OW1JY069lQSrpdwWzLZkk/GAABUxmQMAEBlTMYAAFTGZAwAQGUHPoFLJV48e/YsaSu99+6oVJKPShww61fFRp2/1yeqT709Qfebt/ewSjZRCVKlx4HqJy+pQ71WnYdKyvHGRo5EKS/hTMWr4ii9n7GKw9svWN0LotXovGSoXLokUKo2NZa9e0OfxMoue6yrPr19+3bSdlD2Bu5C9b1XObJkojCfjAEAqIzJGACAypiMAQCojMkYAIDKiidwqcVxbwsxlbihFsxVMkefCjQ5qMQGdR4qGcPbyrBd/WVjY6NXbF4SgkrgUVv7qS0Uc1L95CUoqWQTNQ5qbDHnxay2nVNVi/Y7eU4l3piZXblyJWmLJnXlpK5llaxnpseouibVfaj0PcPr54cPHyZtp0+fTtpUzF7f9xn3XSq/qXGvzqN0Apd6b71zVwl6qp/VMS9fvtw9uBHxyRgAgMqYjAEAqIzJGACAypiMAQCojMkYAIDKqmRTexl309PTSdu5c+eSNi9LuCRvf0t1firzUGWlnjp1asSofqXiU9m8ZjoL+dq1a0lb6cxklcXoZdqrmNU4qFEW1SvhqTLAVYZnyaxwNe7UnspmOrtYjasaTy54ezCrJxLU/SXnfstR6pryqPhUBnnOe4a6VrzM/ps3byZtpZ+2UFTM3vWnMsPVNal+Xs07pfHJGACAypiMAQCojMkYAIDKJpqmyXOgiYlVMzs6GAyGrts9efJEtr969SppU+tpc3NzSdvU1FQozkePHu2sxTzfbgrF7FXBevDgQdJ2+PDhpG1raytpO3RIL9nPz89/8N//8R//Yf/zP/+zZ8wqvu+//14ef3JyMmmbnZ1N2kbZKi/Sz2p9SlUCM9Mxq/csOg7adsfbNM1sl/HsefnyZdL2448/Jm3RczP78PwifazG3f379+Wx1Xg8c+ZM6HVRfa+/1dVV2f7TTz8lbZ9//nnSduLEiaQteh59Y1b3Bo+6ft+9e5e0ffrpp/Ln2+9T37GxvLwsj6/uzV9++WXSduTIEfnzEX37eWVlRbar60+95+qeoeaY9s+37xl7BhmQczL+2cz0pqMH0052BDGXNW4xrzVN89mYjedx62MzYt4vxFzeWtM0I2c15sym/mj733dmtp7xuLlNmdmk/RqvGTGXMG4xt+Mdh/E8bn1sRsz7hZjLU/H2lnMy3jSzT+x9cOPwF83m9r/EXNa4xby5699xi3lc4jUj5v1CzOVtDn/JcDkn43UzG4yyZqy+z1drKWotM/rs467v+Xf+2grF7FHrFS9evEjazp49m7T1WLPKErP3O9rUOvLMzEz4eDX6Wa3/RH6niDc8nr01q+fPn8v2iJMnT8r23f3ft4+7xKvWKNW4iK4VRmLuss7dzq8wy/8cdO6xHF2rVbkn3lpmW9+YR82NUWvzPfJ5OsWs+tNMjxkVszp+5N4s4h1Jzsl4xcyOLiwsDC3KoQofmMWLDqhtxLwHv9Xr7t27Z/Y+XrNgzB61HaF62PzOnTtJW/SmkTtm73e0RQsqeMer0c/qwor8ThFveDyr2MzMlpaWhv5ez7fffivbd/d/3z7uEq+6+aufjxYpicSsEvu8ZEJVLCN6L4jKPZajW8V22a6zrW/MXnEjdV2p+L777rvQzyp9Y/YKlaj4VFvfe7OIdyQ82gQAQGXFy2Eq3l8y6q8y9VeuKvP4+vVrecxcX1l5f52pTxOqlFqNEoIe1c/BxPRLAAAgAElEQVTbf+ENVXrzcK+f1V+06huS0mUP1dhVJfbMdNlJdR6qBKjauN0sT/97n1Si5TvVxus5S6eqPlalIc38kq9tx44dS9q6jLWc1Dc6t2/fTtpylr6MUteUR/WTul979/tcVH+a6TGjYlH3jC79kAufjAEAqIzJGACAypiMAQCojMkYAIDKqiRweY82qYV4L9W+rXSClBezSgxR56F+3ksSyPX8sJc4EU22qZF05j2WEn1cRfW9eobarF8/Rx+/80QTzLzxloP3/qvxqGrDq71tc+ryyJBKaIu+r6UTi7zz6PJ4UknquvDeW7V3sRpHpfdAV2PU6+crV64kbeqxPHWteedRMrmPT8YAAFTGZAwAQGVMxgAAVMZkDABAZUzGAABUViWb2stUUxsTqIzHu3fv5g7pAyo7zyvHp85FZdeqcndetrJX3m0v6nd6/RwtfVk6m1r1s1dacpQs5pwZkCoD1YtZvTaaKexl2nu/K4fTp08nbarf1RMEOU1PT4/08+o8SpdxVdeal5msYnn27FnSVvr665JNHr3WVLayN5b7lJxUfeJlz6vjq59XMXvXaclscT4ZAwBQGZMxAACVMRkDAFAZkzEAAJVVSeBSC+YeteBeukxclyQZldgQPb+ce2aqhAMvAUSdn9p7t3Q/K2p/aDNdltFLqmvz3s8+yRjqWGo/Yo86D5XQU2M/VRWHSpZU48IrOdoneU7F4SVgqd+ryiCW3mNcJV16iZgqZpXgWXpfZXXP8JLzouVZS5cYPX78eNLmJZdF31+VANYniXZUfDIGAKAyJmMAACpjMgYAoDImYwAAKquSwOVVN1lcXEzaVMKM2g+0S1LYMNH9Rs38hKM2lRiRM0Gjy96iKslCJUN5STm5qESgpmnCP6/6T42XXPtDm+nEKq+fVbKJ+vmcYzc3dS10qfbWZW/ivXhJeCrxSSUbqYScg9zvZt32ye5D9ZN374veq1TSVM5EOZXIp5JPzfT7G71nlL73KXwyBgCgMiZjAAAqYzIGAKAyJmMAACornsDVJYFDJdqoxfUayRheYoNKTlAVmWpUVPJEq+SUTiDpIpqMkTNZK6pLtZ6SW7BFedekao9uA1m6370+vnz5cujna4wLj9ouUSl9/ak+8ZKt1LhVCWDqPhdNcu3Lu6bU+al7xuPHj5O2GzdujBpWZ3wyBgCgMiZjAAAqYzIGAKAyJmMAACpjMgYAoLLi2dQqo83LbFRZcSr7rUu5ytKimck19gb2RGM+SBmoKrM0usdqaV6pRlWmL2dpwL6860dlTqvrT2W2l35awDu+KjOrMq8P0vV36tSppE2dR42Yvaz1aBlXNba8vahz8cazyvpX975r164lbTWeeuCTMQAAlTEZAwBQGZMxAACVTXTZsm7PA01MrJrZ0cFg8MFa49bWVvLa5eVleYz19fWk7fDhw0mbWseKrsU9evRoZ7vA59tNScxdrKysJG3Pnz9P2r766qukTZ2bkjtm1c9///vfk7ZRtn3MHfOTJ0+SNvWez87O9jr+7nibppn1xrPy4MED2T41NZW0nThxold8St8+Vn1pFr/+VL+XHhf3798PHd9M9/Eoa/W5x/LGxob8HW2qT48cORL6HX1j9vJJVM6GGi9zc3NJW+mYPeq6VPORii8yntv3jB4hfiDnZPyzmQ2yHGx/7GzgS8xljVvMa03TfDZm43nc+tiMmPcLMZe31jTNyJmZObOpP9r+952ZpX8yHRxTZjZpv8ZrRswljFvM7XjHYTyPWx+bEfN+IebyVLy95ZyMN83sE3sf3Dj8RbO5/S8xlzVuMW/u+nfcYh6XeM2Ieb8Qc3mbw18yXM7JeN3MBqN8x//y5cukTa3JqnUsby2u/dpd3/Pv/LU1Usxq3USt/8zMzCRtan3FzOzQoQ/flr4xq9jMdJ9uH3+oL7/8Ura31136xuytWal2tfvNyZMnkzbV920i3pHHs1qf8taX2+bn52X77nXovn2s1ix3jtemzkFda5E+3vkdOa8/73e0qfPwfmeu68+jYlFrsup98sZFW+6YVZ6Pur+o47f70xOJWfWdup+Z6flE5Q6o8RyJWcQ7kpyT8YqZHV1YWOi0beJut2/fTtquXLmStKnFde9h9fZrz58/b/fu3TN7H6/ZiDGrC189rP/NN98kbdFtGfvGrGIz0w/rbx9/qD/+8Y+yvf2QfN+Yu2zvp7aq/Pbbb5O2SNEBEe/I41n9ARG9GX733XeyfffP9+1jb2s+NW7VOfTt453fkfP6835HmzqPO3fuyJ/Pdf15VCyqqIZ6n7wCM225Y1ZFMNT9RfVpNHkuErPqO2/7XNVXX3/9ddKm5o5IzCLekfBoEwAAlTEZAwBQWfHa1Ir39Wy0/qn6WsH76i36/GNfKmZVM/nmzZtJm1f/NFdN2ujX4GZmi4uLSdvVq1eTNu9rsly1XL3jLy0tJW2qpuyNGzeSttK1cT1q7Kqvqbt8HZljbdWr86zW4BVVc9u7/krX4lbnopZcBoM0B8jLTygds7q+Hz9+XPR3Kur8ves4+trSfafGmXePV9ef+mpd3a9zLptE8ckYAIDKmIwBAKiMyRgAgMqYjAEAqKxKApe3yK8W11XigEpuyVk0oAsVn0owUzF7SS+5eM9eKyoWlQBWOkHjwoULsj36bGY0QcqsfHKfiiWa1JWLuqZUMqGZTsxS1M97z7TnSkb0Eou8c2lTfVz6/fdEEyi9BMxc1DWlaj2YmZ07dy5p857vLUm9j971o85F3RNVP9RI7uOTMQAAlTEZAwBQGZMxAACVMRkDAFAZkzEAAJVVyabuUnJNlSqrkeXrUbF4maVtObM5VeagKg1pprMPo1sUls6g9EpXqtKXqmRd6Qz1LlRfqfGiYs7Vz9ESl2bx7PvS/a6OH82a9njZsTWo81P3gtLXWpf3rNb9dRTq/qfGgbpf1zhfPhkDAFAZkzEAAJUxGQMAUBmTMQAAlVVJ4PKSKVRilyqnl6vEXg5q8T+auJPzPO7evZu0eaXtvPYIL8mndOk+1VeqdKZK9MqZKKeSb7y9T1W7Gi+qVGqNBBL13qprUvVn6QQpVS7STCcjqnFRo3Sjd52pc6mRwNWFGsuqHO1Bul+r/quxT3EUn4wBAKiMyRgAgMqYjAEAqIzJGACAyooncKlF/uvXr8vXnjp1KmlTyS2lqUV+VQnMzGxtbS1pu3LlStLmVR3LRfWzF7Pq06WlpaRNVbCpcR5mOrno2LFjSdvp06czR/QhVbXIG8+K6tOSCS5qH9rBYCBfq5LwoslaORPOuiQBRRPHSldmU3139erV8M971fIOCnWfU/cRdS1451b6XqLGjEr4U/ecLnvB58InYwAAKmMyBgCgMiZjAAAqYzIGAKCy4glcKqFGJd6YmT1+/Dhpi26hmLPKklr4946vYo4u/nsJDCrJYBiVQOMlvaj+U+9J6QQLxUuGUskiKlHOS07KJVolznutGhtqbHnjres47zIuVMUoVc1K9XGNsXKQqGtWjU8znSB6+fLlpE2Nea8qV58kQPUzXqUzlYwWHd9eUmauMeNtWauS+9T1o7bn9CoKlqyMxydjAAAqYzIGAKAyJmMAACpjMgYAoDImYwAAKiueTX3x4sVQm5nOxFNtKsPay2Duk5mseBl7KjNSZUGqNi8zNlfM3t6dKmu2RtlRxSudp95zVeoxV9914b2Pqk+j+zJ770eOpwa67D2tykiW3ru6C5XdqsZF6X1sR93HN7r3bul9rrvsoexlSbepbOWcvKxs9aSLehLg0qVLSVuN/cT5ZAwAQGVMxgAAVMZkDABAZRNN0+Q50MTEqpkdHQwGvdftXrx4EWrb2tpK2k6cOCGPOTU19cF/P3r0aGc7sOfbTSPFrKq8vHz5MtQ2Ozsrj9lu7xuzt72cWks5efJk0jYzM7Pn8ffSN2b1fpuZ/ed//mfSptZ/5ufnk7ZDh4anRuyOt2ma2RzjWY1TlXvw5s2bpE29H2Yfvid9+3hjY0O2r6yshF47NzeXtEXX2HJff97viIj+zv2IWfW9un5r3edWV1eTtui2lF7M7ftL35gfPHgg29V1NTk5OTSOvWLerX3PGPoDQ+ScjH82s7K1CPPa2aCTmMsat5jXmqb5bMzG87j1sRkx7xdiLm+taZqRM75yZlN/tP3vOzNbz3jc3KbMbNJ+jdeMmEsYt5jb8Y7DeB63PjYj5v1CzOWpeHvLORlvmtkn9j64cfiLZnP7X2Iua9xi3tz177jFPC7xmhHzfiHm8jaHv2S4nJPxupkNIt/xLy8vy3a1rqrW+9T3+T3WrHb+2grF7K1lqrUUtVbx5ZdfJm1HjhwZFq6Z9Y/5/v37sv3t27dJ28cff5y0qWdac8a8vp7+wfv3v/89dHwzHbNah/di3j22RLzh8exRY0atfZ86dSppi4znvuPCo9Ytu+RnROSOWd1LVMwqlyCqb8zqfmam+/ns2bO941P2Y2woKp8gKve9OZpv1OWesZuIdyQ5J+MVMzu6sLAw9AF77yFtVehA3ZRU0YHoA/bnz5+3e/fumb2P1ywYs1dURMWiEqT++Mc/Jm3RLcT6xuwViHj27FnSpgZf6ZhVoo3actOjYv7DH/6QtHkx7x5bIt7wePaoMaO2yus7nvuOC48q+KCSiKJbhCq5Y1bvrYp5lKI2fWNWxXXM9NaKuYuS7MfYUEYpCJP73qza1dhQ5xa5z4l4R8KjTQAAVMZkDABAZcVrUyveVwDqqzn1tYeq5fv06VN5zD61fNXXp+rrRTOzY8eOJW3q+Vf186reslm+uqjeV0vq/FT9WBWzt37TZy1KvTfqKzyP+hrq6tWrSZsXW5fawXvxnudWY1eNjf2up+19lbi0tJS0LS4ulg4nxBvLatx69c33m9d33nV/EHjPaKux0eVazSW69GOm44ve2717Q47a8B4+GQMAUBmTMQAAlTEZAwBQGZMxAACVVUng8hbH1eK6emZXLcznXFhXCTXec4oXL15M2tSm29evX0/avMSf0glc6vnH6AbgOfs5+gy51779wP0H1MbyOROk1HvWZXPz6LP0JXnPZariI9Hnyksbx5i9ZKhRnsXNScXnJZddunQpaVPnoTaPyHnP8O6ZUeqc1Rgqmajl4ZMxAACVMRkDAFAZkzEAAJUxGQMAUBmTMQAAlVXJpvZEM9j2u3ygmc6aHtWomYF9qexB5dq1a0nbfmf+7lBZmorKus0Zs8og9XbnURmoahyp7F/viYMcmcJeX6qnAGq9321eHDXuBYq6llW2v9nBiVmNZbWjm5nZ73//+9DPqzHkjbc+Y0s9IeLt6qTuBeo9OSjZ7XwyBgCgMiZjAAAqYzIGAKAyJmMAACo7UAlcakFfJcGohfmcJfBUyTQvocZL0ojwSvyVTihQiXKqjKQq3eiV2Cyd6KPeX9V/XWLuw0sWUVQ/q1hUKVIv6aXrOFfxemNWJd+o/lQlE3P2sUqG8pIdVT+p81D9VqPkoZmOWfWzSvTKtQe3me4nr5/VfViNo2hyaE5eqWLVV+o9z9mno+CTMQAAlTEZAwBQGZMxAACVMRkDAFDZgUrgUlQy0/T0dNLmJdb0WZxXVWi8ZA+V8KB+Plp5qRaVzHHhwoWkzUs6y5nAo6hkFtV/6jxyxqbO36uopPawjvL2le2qS5KSOg/Vn+qa9Pbu9cbLXlQyoNfH9+7dC7Wp90IlJZml57e1tSVft1u0QpyZ2fHjx8OvbVN7uZv1S/pUY8NLhlLHv3r1atJWugKe4o09NQ7UeD4o+GQMAEBlTMYAAFTGZAwAQGVMxgAAVFY8gUslOHmL6Cph6PXr16Hf4y3i90ngUslW3haK6vjqnGskNnhb+929ezdp8xI3SlL91CXZSlWuUryx0WcrO5X04lUtUkmFaoyrsZUr6UzFq7bGNItvOdcl8ScXL6FNJemoanKKN37aiYHr6+tDj6XG0mAwCMVhpvtPvR9LS0vy50tX7VP3KnV+NbaH7JKkWKvqWgSfjAEAqIzJGACAypiMAQCojMkYAIDKmIwBAKiseDa1ysLzMi9V1qtqU9mnucoHerxymyqbc3FxMWkrnTmteFmzjx8/Dv28KheYs4SnykL2Mp9VzCprVpXoq5HhaeafS9t+l+jzfp/KNFVjWZV+LH39ednlKpZoydLoPrh9r13v+OoaUln2anyr92M/qLGhMu3VExzekyi5eP2sfi/Z1AAAwMVkDABAZUzGAABUNtE0TZ4DTUysmtnRwWAwdI3u/v37sv3w4cNJm6p+o9Zw5ubmQsd89OjRzlrH8+2mUMxelSW1lvmb3/wmaZudnd3z+HvpG/ODBw9k+5s3b0K/94svvkjavH4+dOjD9INIzBsbG8lxnjx5Eo5ZVQBS/TwzMyOP6cXbNM1sl/HsWV1dTdp++umnpO23v/1t0jY1NTX0+H3HhefFixehNvW+eX3cHi+5Y15ZWUnaVMyqP71rsn0uue8Zy8vLSdvbt2+TNjW+vX5un0vuflbnou59J0+eTNoi159Z/5jV+21m9vLly6Rtfn4+FEtE+54x6vFyTsY/m1m8/lt9O9kHxFzWuMW81jTNZ2M2nsetj82Ieb8Qc3lrTdOMnKGbM5v6o+1/35nZ8GKu9UyZ2aT9Gq8ZMZcwbjG34x2H8TxufWxGzPuFmMtT8faWczLeNLNP7H1w4/AXzeb2v8Rc1rjFvLnr33GLeVziNSPm/ULM5W0Of8lwOSfjdTMbjLIuoaj1ZbXWE12T3fU9/85fW0nMW1tb8ucUtX4WWe8z67Vm5casqPPYOV6bivnEiRN7Hn8vffvZWzNWr1Xr132fCRXxhseztzavqPWzI0eOJG0qf6Kt77joQq3JqvXDM2fOhI4XiVnliajrzIuvnb9gpq811e99Yx6VGveqn701z/a47xuzWs820/2vxmjpe4bSZW1e3ef6riOLeEeSczJeMbOjCwsLboGMPtRD2n/4wx+StuiWc+fPn98p1LFzFScxqzfX24pRFR2IXqSqQIVZ+rB6JGbFG6TqXFTMqnhCVN9+9opHqNeqbeP6bJm583OteMPjuctNWZ2fKgIRKU7Qd1x0oa4rdfzo74zErP5YfPbsmTzelStXkjb1B5k6j2gBm/3oZzUu1PG9rRLb475vzF6fqPucGqOl7xmK9//VuXQp/jKMiHckPNoEAEBlxcthdqH+alF/Eau/0nJSf5l7JSS9kpNt6q/InGXi1F/M3idvRZ2fej9ybXpvpvtZlRf1qBKCDx8+TNpyfp2oyv1Fy4t6rx3lk2ZO6rryNrNv876F6bNsoD5deXGokpHRT2w5S7t2oca9GleKV8q0z3hRcdy8eTP88+paVX3a99sqRcXcpRSr6md1T/O+gSiJT8YAAFTGZAwAQGVMxgAAVMZkDABAZVUSuLzEiWgiTOk9KVWSwKlTp+Rr93svWo+XQKOoc1FJTurcciZwlaASL0Z53KJNJQypGsJmuk9rJGZFqetPPTqk+jPnft1dEu7U+63GaI39xD3qUZpjx44lbd7jXLmoseiNZfVadR6q76P7ekeo3+nNB+r3qjaVYObd10uOIz4ZAwBQGZMxAACVMRkDAFAZkzEAAJUxGQMAUFnxbGqV7eiVXFOlJa9fv560dSl/1ofKTPZKcKrMT/Xag5SZHC2Mrl7nZUb2KTmpshhv3LghX3v58uWkTWWgqrGVM5taZVN65x7d3EI9XVAj+zdanrVvYf0odX175QlVuyrT6I2rGk6fPp20qXuGGss5S3iq8emNZdUe3VwmZ6nUUY8TLfHrjfGSJVT5ZAwAQGVMxgAAVMZkDABAZUzGAABUVjyBSy14Hz9+XL729evXoWOq0mw5F9bVgr5Xcm16ejppU3vqRkvH9aWSQrzSdoo6vy7lHHPtGey9j9H3d2JiImnzYs61z6p3HJV8qPaYVmOjRtlMlZx36dKlpK10CdguSXIqyUnFXHrvYpWk5CUOqvdW7bOrytbW2oM5qsvewAellLDiJeyWxCdjAAAqYzIGAKAyJmMAACpjMgYAoLLiCVwqGcOr9HNQ9npVMXdJnFhcXAwdMye1z65H9bNKQlKVjLxkjINCJb14VcNyJXB5iSjq+CpZa2lpKWnLVelMvV8qicyjrtXSleNUMlTpql+jUmNAva8eVU2u9P1QJX2O2s+l9/DuErMaRwdpX+s2PhkDAFAZkzEAAJUxGQMAUBmTMQAAlRVP4OrCq3J1EHgVWdR2bypxwEvIyUUlJngJSl0qjLXlqrRlpvvJ225NtavEkBp978W8traWtKlkE1UpLde1oN5rL6FG9ZOqDKWSwryx1me8qLGs+tJMJy6W3kZT6bINporl2bNnSZtKCsuZQKmS865cuSJfq2JRbercct7XVYJml+1tVSxq3NdIVOWTMQAAlTEZAwBQGZMxAACVMRkDAFAZkzEAAJUd+GxqlS2pMlJz7vOpsutU1rSZziJVmXg1yrB5GaSqn1VG8LVr1zJHNPx3eu+jKs2pSgiqcytdvtE7vsrqVW7cuJG0lRwvXUoeqnNTWbTe+5YrK1XtEW4WL2mqssK9srx9qMxdL5Nc9anqvy5lUnOVnPTGshoz6vxGLSU8jLq+vess+qRL9P0ojU/GAABUxmQMAEBlTMYAAFQ20TRNngNNTKya2dHBYJC1SpNaIzl0KF3qnp+fDx9vu5rP8+2mJOb19fXk5548eSKPp9ZI5ubmQjFHRWJWtra2ZPv9+/dD8R05ciRpi1bTicS8sbGR/Nzy8rI8nqrA9PHHHydthw8fTtrU+2FmNjU1JeNtmma2y3j2Yv7v//7vPX9ux5dffpm0qb5v6zsuulhZWUnaXrx4kbR58bb7vm/M6prcOV7bu3fvkraTJ08mbTMzM3v+zt2/I2c/q+tSjaFXr14lbapam1m6fts35tXVVdn+8uXLpE2dh7qPePeM9r2zb8wqNjOzH3/8MWmbnJxM2mZnZ5O2yH2ufc8Y+gND5JyMfzYzPVIOpp27OzGXNW4xrzVN89mYjedx62MzYt4vxFzeWtM0I2dc5sym/mj733dmpv+MPRimzGzSfo3XjJhLGLeY2/GOw3getz42I+b9QszlqXh7yzkZb5rZJ/Y+uHH4i2Zz+19iLmvcYt7c9e+4xTwu8ZoR834h5vI2h79kuJyT8bqZDXKvpTx48CBpU+sS3u9sv3bX9/w7f22FYvbWrNRa8tu3b5M2tb559uzZPX/njr4xe9T6lFp3Ucffvc66l0jMas3YW5tXr1VrfmrtMvLMrohXjmc1Rr1nP1XMZ86cSdrUOndE7nGhzu1vf/tb0pZp/bVTzKovzfS4Vbv4qD4ufc/wqHX458+fJ21qh6Lo8+f7cc9Q4yWau6P0jVnNEWZmb968SdrUmrGKuec9YyQ5J+MVMzu6sLDQ+wF09UB29MHyO3fuyGO2X3v+/PmdAhI7V0QoZu+Gqx4sV9uhqUki2k99Y/aoB9rVQ/3fffdd0ha9mCMxq5umV1xFvfabb75J2tS5RQpDiHjleFZj1Du+ivnWrVtJW98t5nKPC3Vu09PTSdu3336btEULaPSN2dsmT41bVZRE9XHpe4ZHFZlQBT66bFWpXlf6nqHGS5eCMm19Y/buSY8fP07a1IeJvv0s4h0JjzYBAFAZkzEAAJUdqNrUivrKV7V5X2PleubSq1WqvjJXbeorEy/mvl9btnl1gVUtV/U1Y87nVRV1ntevX5evVc9Xqq/71Nfc3hJDn35WX8N576Oi3pNc9ZtHpep/K6qWds46z4r3taF6D6PLMF6/q6+5h1Ff2XrHUV9Jq/XhXPeBLrxrRd0zLl26VDqcEG9sqHuBGgcXLlxI2rxa6CXviXwyBgCgMiZjAAAqYzIGAKAyJmMAACornsClNvVWCSBmo22QXTrZyItNJWupxCKV5FM6QeNPf/qTbFfJUN5r91uXRCDVf6UT5VSCi7chu/q9uTaBL2FxcbF2CC6v37psXFKSuv+oRFMzs2PHjiVtKtmoRgJXl+eED8o9o0sCZPT8vPFGAhcAAP8fYzIGAKAyJmMAACpjMgYAoDImYwAAKiueTa0yd71Mte3tqHodszSvtJ3K0lRt6ty8Epu5shS9LGJVPi66NdtBojL1z507l7RFd7qJUKXzVHasmc6y3u9+VpmmXkZptBxmzv6MGjWzWJWrHOWYW1tbHxxTHd+7T6lymCqbWh2zdOlU797sjfGDTPWfetpC3TO8JyRK4pMxAACVMRkDAFAZkzEAAJUxGQMAUFnxBC6V7KEW1s10kpTa37bPfqNdqCQGlXTRhSrzWDpJwEuIU4ldqk8PQgLSXtQesKXjU++jV2pRJUR5pWBLUeX7vOtPUedQuvSsSsxTCU6jmp6e7v2zhw4d+mCsqf1+vZij51LjWvOSPlVpTxWfSjDzElVLU/GppLoaZUcVPhkDAFAZkzEAAJUxGQMAUBmTMQAAlRVP4OrCSx5oK10BSB3/7t278rWqmpFqU4kNpRMHvApAqgqNalMxexWAaiRpqESYLvux9qGSiy5duiRfq/q/RCLSXtRY9q4flbCnErhKj1uVmOe9r0+fPk3aVAU7lbSmfk9fKrFPxWamx4C6/mpUgfLeW5XApRL5usRc456h4itd1SyKT8YAAFTGZAwAQGVMxgAAVMZkDABAZQcqgSuaGKIqHpWuCuTFppJFVGJCjSovXrKQqhaktkhTSRte9bN2UtDGxsbQ+FRSjZeoo5L7VN+X7ufXr18nbV6ls1ESuLzqSyUT1LpU5ipJvYfe+6piVmNUjfnSY8U7fvT31qjA5Y3PaCJf9P0wy5fA5VXAU+0qWUtdv15SV8mkOj4ZAwBQGZMxAACVMRkDAFAZkzEAAJUxGQMAUNmByqZW2XUqe1RltHkl/nJlTHqlOtXxS++3HOVlBKosTZWZrHhZ6+1+OHz48NBjqcxLLw7V/yqW6Hn01aUEp4pZZUTCmSUAAAjbSURBVKUq586d6xRXDuq6UhmpqnRjrT1h1T1DtalylUh52cJqLKt939U4KF0C1rv3e085tKmxUaNUJ5+MAQCojMkYAIDKmIwBAKhsommaPAeamFg1s6ODwaB3NSxVtenJkyeh1505c0Yes712+ejRo521hOfbTaGYvepEKysrSZs61qFD/Zfn+8a8tbUl29X6z4sXL5K2d+/eJW2ff/65POb8/PwH/x2JWb2Py8vL8vjqtVNTU0nbiRMnkrZI3++Ot2maWW88qz7tEvObN2+GxmLmb3+5O5a+48Kzvr6etKlzU+uCMzMzod+RO2Z1f1D5CnNzc72Ob7Y/MXe5p0Xkjlnd59Q9UfW9NzaOHDnywX/3jfnf//3fZbu6fynqnqbuI2Yf3kva94zQL9tDzsn4ZzPTd5CDaWd1n5jLGreY15qm+WzMxvO49bEZMe8XYi5vrWmakWuX5sym/mj733dmlv6JfXBMmdmk/RqvGTGXMG4xt+Mdh/E8bn1sRsz7hZjLU/H2lnMy3jSzT+x9cOPwF83m9r/EXNa4xby5699xi3lc4jUj5v1CzOVtDn/JcDkn43UzG4yyLqHW49T6inpOdnZWf2XfXi/c9T3/zl9bScxqLcRbF3z79q1sb1NrgO111lFiVtQ6j5k+P9X3ap2tx9pgp5g9ahy8evUqaTt58mTSFolZxBsez2qt1Uz3v3r28YsvvkjavDWr3Wr1sdrhq8uOa8NiVmPRe87/+fPnSdvk5GTSpsZApI/N8vezGi+q79WacTT3JBKzWqd+8OCBPF50/fXTTz9N2rz73B75PJ362cuNuX///tDfadZ/bV7EO5Kck/GKmR1dWFiQD4NHqElCPTCuHvL2HlZvT9znz5/fKbywc6dMYlbxew+Bq20GFTWgvGIRfWJWvD5RP6f6/ttvv03aosUT+sbsUePg9u3bSVvfmEW84fHsbeGm+l8V/fjmm2+Stkjxklp9rK6FaKGbSMxdtuFbWlpK2lRiX98+Nsvfz2q8qL6/c+dO0hbdVjESc7SQjlm8gIb6A/7WrVvyte0/4Pr2s5dcG92Ks+/7KOIdCY82AQBQGZMxAACVHaja1OrrG/W1nmrzvpLtI/r1hvd7VU1oFbP3NXWfuqjqqxb1FZ6Z2aVLl5K26FdnuR6F6yr69Vxp6qs9r/auGjOq5rS3HlqK9/Ws+kr61KlTSVvpWsPquvDG8rVr15I2NVauXr2atHn3jL5rwVHq+q5R21v1U5exrMaRqlvujbfSNfzVV+sqPjXua9Qy55MxAACVMRkDAFAZkzEAAJUxGQMAUFmVBC7vubBoEpZaXM+Z4NPlWTS1+K+ePVZJUzk3sFbPXqvkGzOdIBPtP++9y9X/Xj/fvHkzy/FHpRJAvM3NVeKKes+9Pi2lS8KYSuipkeDkJTuqPlbXn7pnlE6aUomcZnoMqfPr8hxwH+qajT57bWZ2/fr1pE0VN8qZXKt0iVmpkayl8MkYAIDKmIwBAKiMyRgAgMqYjAEAqIzJGACAyopnU6tMUS+LWGUZRo9Zi8oeVEpnoCrebkKqXWUxqxKCXiZun/OL7hLUhdreLyeVeellY6pxqvpZ9YPXzzkygL3MZCVnxn+U6jevP1Tm9I0bN5K20uehrimv3GM0e/f06dNJ29OnT+Vra5TTVFTmdOlStl2yqb0nTA4CPhkDAFAZkzEAAJUxGQMAUBmTMQAAlVVJ4PISsNTepKrk2kFJVujioOzHa6aTraJJcTkT0S5fvhw+frTsYY1+9hJ11NhV1F693vsxauk/s27lMFVSZY0ykl1irpEsqcai2k/Xo5LOFK9cbI1EO5UsqcZnzn2L1diIJv4edHwyBgCgMiZjAAAqYzIGAKAyJmMAACornsDVZW9gRS3YH6QELpUsoqpyqapHNZIuPNGkF6+qV5+kmS7VrHIkLpXSZTyqCkBqHOTaA1YlPnVJLFLVwdRYzrknrKo85SXmqfhUwlCXqmN9qPP3krLUWL53717SpsZKjeQ0j4pFVdDLuQe6ula67DV/kPHJGACAypiMAQCojMkYAIDKmIwBAKiseAJXF2ohXiWbnD9/vmgcKuFAJZKZxRNkSleGUjF7iUU3b95M2h4+fBj6PbUqialkEVX1qAYvEU9VGFNVuXImP7Wp9+vSpUvytWpcKFeuXEnacp6DOpYXsxoXautPlXhYOhnKGxcq6Uxdq+rnc8as7l3evVWNo2jlq5zbrqo4vHuzSuA6yNW6+GQMAEBlTMYAAFTGZAwAQGVMxgAAVMZkDABAZQcqm1plyqmScKWp0nTRvWnNdDlML+MvF9V33j6iKjNVZYCrrNlapUijpfdK76eqeGX31NgtmTmtqHHhlRZVWa/qWojufZ2Tl+WrMqfVvugHqYyk6r8uJUpL8t5b9eSC2s9Y7c1duu+9JzxUfOo8amTaK3wyBgCgMiZjAAAqYzIGAKCyiaZp8hxoYmLVzI4OBoPe37evr68nbcvLy0nb3Nxc0hatDPXo0aOd9Znn201JzC9fvkx+7scffwwd38xscnIyaTt79mzSduhQbMk+ErOyuroq29Xa4Lt375K2o0ePJm2q75W+MXtUzNF1rMg69+54m6aZ7TKe1XjxYj5z5szQWKJy97FaO1NrmWp8f/311+Hf0SfmJ0+eyPZXr14lbX3HgCd3P29sbCRt33//fdL2m9/8JmmbnZ0N/Y5IzCqOBw8eyOOp+8PHH3+ctKn7w8zMTLaYla2tLdmuzuXt27dJ229/+9ukbWpqali4yT1j6A8MkXMy/tnM0sylg2vnLkPMZY1bzGtN03w2ZuN53PrYjJj3CzGXt9Y0zch1gnNmUz81s+Nmtm5mKxmPm9ucmU3Z+3jNiLmUcYu5He84jOdx62MzYt4vxFxeO96RZPtkDAAA+iGBCwCAypiMAQCojMkYAIDKmIwBAKiMyRgAgMqYjAEAqIzJGACAypiMAQCojMkYAIDKmIwBAKiMyRgAgMqYjAEAqIzJGACAypiMAQCojMkYAIDKmIwBAKjs/wIOF+DAb83uvAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 100 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.datasets import load_digits\n",
    "digits = load_digits()\n",
    "print(\"shape:\", digits.data.shape)\n",
    "\n",
    "def plot_digits(data):\n",
    "    fig, ax = plt.subplots(10, 10, figsize=(8, 8),\n",
    "                           subplot_kw=dict(xticks=[], yticks=[]))\n",
    "    fig.subplots_adjust(hspace=0.05, wspace=0.05)\n",
    "    for i, axi in enumerate(ax.flat):\n",
    "        im = axi.imshow(data[i].reshape(8, 8), cmap='binary')\n",
    "        im.set_clim(0, 16)\n",
    "plot_digits(digits.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model your data with your GMM class using a model order selection method to produce new synthetic handwritten numbers. Explain why you used that model selection method in this case. You should use PCA to reduce the dimension as GMM doesn't perform well in high-dimensional contexts. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# TODO"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}