gaussian mixture model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian Mixture Model\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Gaussian Mixture Model is another form of unsupervised learning. Mathematically, it is somewhat similar to Gaussian Discriminant Analysis. In GDA, we can predict a new data point based upon p(y) and p(x|y=label) where p(y) is the prior probability of a data point belongs to one specific label and p(x|y=label) is the probability density function of Gaussian distribution under the label. The concept holds in GMM. \n",
    "\n",
    "From the perspective of GMM, the metadata consists of data points belong to different labels. Each label follows a multivariate Gaussian distribution with unique mean and covariance matrix. However, GMM is an unsupervised learning which indicates we do not have training dataset to compute p(y). p(y) is crucial to a follow-up computation of μ and Σ in Gaussian distribution, yet it remains a latent variable we do not observe. In order to estimate p(y), we need to take an initial guess and gradually converge to the optimal number throughout iterations. This is when the algorithm called Expectation-Maximization kicks in and saves the day.\n",
    "\n",
    "The material from Stanford university is the Bible of machine learning. It did an excellent job on the derivation of EM algorithm, but it didn’t cover much about the mathematics of Gaussian Mixture Model. I will recommend the following materials as a supplementary to the Bible.\n",
    "\n",
    "Detailed explanation on Gaussian Mixture Model\n",
    "\n",
    "https://towardsdatascience.com/gaussian-mixture-models-explained-6986aaf5a95\n",
    "\n",
    "Detailed notebook on EM algorithm and Jensen’s inequality\n",
    "\n",
    "http://krasserm.github.io/2019/11/21/latent-variable-models-part-1/\n",
    "\n",
    "Reference to Gaussian Discriminant Analysis\n",
    "\n",
    "https://github.com/je-suis-tm/machine-learning/blob/master/gaussian%20discriminant%20analysis.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import sklearn.datasets\n",
    "import sklearn.mixture\n",
    "import sklearn.decomposition\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#official sklearn gmm\n",
    "def skl_gmm(df,xcols,num_of_clusters,diagnosis=False):\n",
    "    \n",
    "    data=df.copy()\n",
    "    \n",
    "    clf=sklearn.mixture.GaussianMixture(num_of_clusters,\n",
    "                                        init_params='random')\n",
    "    clf.fit(data[xcols])\n",
    "    \n",
    "    if diagnosis:\n",
    "        print(f'{clf.n_iter_} iterations to reach convergence\\n')\n",
    "    \n",
    "    data['label']=np.argmax(clf.predict_proba(data[xcols]),axis=1)\n",
    "    \n",
    "    #compute accuracy\n",
    "    erreur=0\n",
    "    checked=[]\n",
    "    for i in data['y'].unique():\n",
    "        erreur+=get_accuracy(data,i,checked)\n",
    "        checked.append(i)\n",
    "    accuracy=1-erreur/len(df)\n",
    "    if diagnosis:\n",
    "        print('\\naccuracy: %s'%(accuracy))\n",
    "\n",
    "    return clf.means_,clf.covariances_, \\\n",
    "clf.weights_,clf.predict_proba(data[xcols])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for unsupervised learning, clf.score doesnt return the accuracy\n",
    "#there is no cross validation, no known labels\n",
    "#the only way to detect the accuracy is vote of the majority\n",
    "#for each label given\n",
    "#we check which iris type is the majority\n",
    "#we consider the majority as the correct classification\n",
    "#all we need to do is to count the minority\n",
    "def get_accuracy(data,class_,checked):\n",
    "    \n",
    "    df=data.copy()\n",
    "    \n",
    "    #use dictionary to keep track of everything\n",
    "    d={}\n",
    "    \n",
    "    #counting\n",
    "    for i in df['label'][df['y']==class_].unique():\n",
    "        if i not in checked and i!=-1:\n",
    "            d[i]=df['label'][df['y']==class_].tolist().count(i)\n",
    "\n",
    "    #comparison\n",
    "    large=-1\n",
    "    lbl=None\n",
    "    for i in d:\n",
    "        if d[i]>large:\n",
    "            lbl=i\n",
    "            large=d[i]\n",
    "\n",
    "    return len(df['label'][df['y']==class_][df['label']!=lbl])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize random miu, sigma and weight for each gaussian\n",
    "#you can use random number like i do\n",
    "#or use kmeans centroids which is more computation expensive\n",
    "def init_params(df,num_of_clusters):\n",
    "    \n",
    "    mean_vec=[];cov_mat=[];w_vec=[]\n",
    "\n",
    "    for _ in range(num_of_clusters):\n",
    "        mean_vec.append(\n",
    "            df.loc[np.random.choice(len(df),10)].mean().tolist())\n",
    "        cov_mat.append(np.diag([1]*len(df.columns)))\n",
    "        w_vec.append(1/num_of_clusters)\n",
    "    \n",
    "    return (mean_vec,cov_mat,w_vec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pdf of multivariate gaussian\n",
    "#equivalent to scipy.stats.multivariate_normal(miu,sigma).pdf(X)\n",
    "def multivariate_gaussian(data,miu,sigma):\n",
    "    \n",
    "    numerator=np.exp(-0.5*(\n",
    "        np.mat(data)-miu)*np.linalg.inv(sigma)*(np.mat(data)-miu).T)\n",
    "\n",
    "    denominator=(2*np.pi)**(data.shape[1]/2)*np.linalg.det(sigma)**0.5\n",
    "\n",
    "    return np.diagonal(numerator/denominator)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#e step\n",
    "#compute the posterior probability with given params\n",
    "def e_step(data,num_of_clusters,params):\n",
    "    \n",
    "    #unpack params\n",
    "    mean_vec,cov_mat,w_vec=params\n",
    "    \n",
    "    collections={}\n",
    "    posterior={}\n",
    "    \n",
    "    #Bayes rule\n",
    "    for i in range(num_of_clusters):        \n",
    "        \n",
    "        collections[i]=multivariate_gaussian(data,mean_vec[i],cov_mat[i])*w_vec[i]\n",
    "    \n",
    "    #compute posterior\n",
    "    for i in range(num_of_clusters):\n",
    "        \n",
    "        posterior[i]=collections[i]/np.sum([collections[i] for i in collections],axis=0)\n",
    "        \n",
    "    return posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#m step\n",
    "#use posterior probability to update params\n",
    "def m_step(data,num_of_clusters,params,posterior):\n",
    "    \n",
    "    #unpack params\n",
    "    mean_vec,cov_mat,w_vec=params\n",
    "    \n",
    "    #solving lagrangian by partial derivatives\n",
    "    for i in range(num_of_clusters):\n",
    "\n",
    "        mean_vec[i]=np.sum(\n",
    "            np.multiply(\n",
    "                np.mat(data),\n",
    "                posterior[i].reshape(-1,1),),axis=0)/sum(posterior[i]).T\n",
    "        \n",
    "        cov_mat[i]=np.multiply(\n",
    "            posterior[i].reshape(-1,1),\n",
    "            (np.mat(data)-mean_vec[i])).T*(\n",
    "            np.mat(data)-mean_vec[i])/sum(posterior[i])\n",
    "        \n",
    "        w_vec[i]=posterior[i].mean()\n",
    "        \n",
    "    return (mean_vec,cov_mat,w_vec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#using jensens inequality to compute the lower bound\n",
    "#which is the function of the expectation\n",
    "def get_lower_bound(data,num_of_clusters,params,posterior):\n",
    "    \n",
    "    #unpack params\n",
    "    mean_vec,cov_mat,w_vec=params\n",
    "    \n",
    "    lowerbound=0\n",
    "\n",
    "    for i in range(num_of_clusters):\n",
    "\n",
    "        #log pdf of multivariate gaussian\n",
    "        logpdf=np.log(multivariate_gaussian(data,mean_vec[i],cov_mat[i]))\n",
    "        \n",
    "        #use bayes rules to decompose the numerator of the lower bound\n",
    "        lowerbound+=np.sum(\n",
    "            posterior[i]*(logpdf+np.log(w_vec[i])-np.log(posterior[i])))\n",
    "\n",
    "    return lowerbound"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#using mle for training\n",
    "def training(data,num_of_clusters,\n",
    "             tolerance=0.001,num_of_itr=50,diagnosis=False):\n",
    "    \n",
    "    #initialize\n",
    "    lower_bound_old=None\n",
    "    lower_bound=None\n",
    "    counter=0    \n",
    "    \n",
    "    #initial parameters\n",
    "    params=init_params(data,num_of_clusters)\n",
    "\n",
    "    #cap the maximum number of iterations\n",
    "    while counter<num_of_itr:\n",
    "        \n",
    "        try:\n",
    "            \n",
    "            #e step\n",
    "            posterior=e_step(data,num_of_clusters,params)\n",
    "            \n",
    "            #m step\n",
    "            params=m_step(data,num_of_clusters,params,posterior)\n",
    "            \n",
    "            #use lower bound to determine if converged\n",
    "            lower_bound_old=lower_bound\n",
    "            lower_bound=get_lower_bound(data,num_of_clusters,params,posterior)\n",
    "\n",
    "            if lower_bound_old and np.abs(lower_bound/lower_bound_old-1)<tolerance:\n",
    "                \n",
    "                if diagnosis:\n",
    "                    print(f'{counter} iterations to reach convergence\\n')\n",
    "                return params\n",
    "            \n",
    "            counter+=1\n",
    "        \n",
    "        #reboot everything when singular matrix is encountered\n",
    "        except np.linalg.LinAlgError:\n",
    "            \n",
    "            params=init_params(data,num_of_clusters)\n",
    "            lower_bound_old=None\n",
    "            lower_bound=None\n",
    "            counter=0\n",
    "    \n",
    "    if diagnosis:\n",
    "        print(f'{counter} iterations to reach convergence\\n')\n",
    "        \n",
    "    return params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#the main function\n",
    "def gmm(df,xcols,num_of_clusters,\n",
    "        tolerance=0.001,num_of_itr=50,diagnosis=False):\n",
    "    \n",
    "    data=df.copy()\n",
    "    \n",
    "    #obtain parameters of gmm\n",
    "    params=training(data[xcols],num_of_clusters,\n",
    "                    tolerance,num_of_itr,diagnosis)\n",
    "    \n",
    "    #compute the prediction probability\n",
    "    posterior=e_step(data[xcols],num_of_clusters,params)\n",
    "    \n",
    "    #make forecast based upon probability\n",
    "    data['label']=np.argmax(\n",
    "          np.concatenate(\n",
    "              [posterior[i].reshape(-1,1) for i in posterior],axis=1),1)\n",
    "    \n",
    "    #compute accuracy\n",
    "    erreur=0\n",
    "    checked=[]\n",
    "    for i in data['y'].unique():\n",
    "        erreur+=get_accuracy(data,i,checked)\n",
    "        checked.append(i)\n",
    "    accuracy=1-erreur/len(df)\n",
    "    if diagnosis:\n",
    "        print('\\naccuracy: %s'%(accuracy))\n",
    "    \n",
    "    return params,posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize number of clusters\n",
    "num_of_clusters=3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load data\n",
    "iris=sklearn.datasets.load_iris()\n",
    "\n",
    "#create dataframe\n",
    "df=pd.DataFrame(iris.data,columns=iris.feature_names)\n",
    "df['y']=iris.target\n",
    "\n",
    "xcols=iris.feature_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15 iterations to reach convergence\n",
      "\n",
      "\n",
      "accuracy: 0.84\n"
     ]
    }
   ],
   "source": [
    "#the accuracy entirely depends on the initial parameters\n",
    "#to obtain a better result\n",
    "#the initial parameters should be generated by kmeans centroids\n",
    "params,posterior=gmm(df,xcols,num_of_clusters,diagnosis=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14 iterations to reach convergence\n",
      "\n",
      "\n",
      "accuracy: 0.78\n"
     ]
    }
   ],
   "source": [
    "#to level the playing field\n",
    "#we are using random initial parameters in sklearn\n",
    "miu,sigma,weight_gaussian,prob_gaussian= \\\n",
    "skl_gmm(df,xcols,num_of_clusters,diagnosis=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[matrix([[5.00600031, 3.42800669, 1.46199645, 0.245995  ]]), matrix([[6.33380974, 2.68772305, 4.62903768, 1.41648825]]), matrix([[6.24236692, 2.92232959, 4.98157584, 1.74684726]])] \n",
      "\n",
      " [[5.00623485 3.42852345 1.46207491 0.24597442]\n",
      " [6.3016641  2.75354375 4.77784322 1.46295965]\n",
      " [6.23091878 2.96243017 5.00260651 1.83849055]]\n"
     ]
    }
   ],
   "source": [
    "#mean comparison\n",
    "print(params[0],'\\n\\n',miu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[matrix([[0.12176663, 0.09723391, 0.01602788, 0.01012392],\n",
      "        [0.09723391, 0.14081754, 0.01146455, 0.00911303],\n",
      "        [0.01602788, 0.01146455, 0.02955576, 0.00594712],\n",
      "        [0.01012392, 0.00911303, 0.00594712, 0.01088312]]), matrix([[0.56529097, 0.14708111, 0.65033649, 0.22170147],\n",
      "        [0.14708111, 0.09357417, 0.07091504, 0.02652935],\n",
      "        [0.65033649, 0.07091504, 0.94237864, 0.3174626 ],\n",
      "        [0.22170147, 0.02652935, 0.3174626 , 0.1111077 ]]), matrix([[0.39759079, 0.11838905, 0.40078107, 0.15664252],\n",
      "        [0.11838905, 0.10219061, 0.14284535, 0.07698576],\n",
      "        [0.40078107, 0.14284535, 0.57520702, 0.25229624],\n",
      "        [0.15664252, 0.07698576, 0.25229624, 0.17367612]])] \n",
      "\n",
      " [[[0.1217028  0.09701231 0.01599739 0.01014138]\n",
      "  [0.09701231 0.14029193 0.01138451 0.0091445 ]\n",
      "  [0.01599739 0.01138451 0.02955853 0.00595474]\n",
      "  [0.01014138 0.0091445  0.00595474 0.01088854]]\n",
      "\n",
      " [[0.63630405 0.19544266 0.75715975 0.24144962]\n",
      "  [0.19544266 0.14277862 0.20658867 0.07444944]\n",
      "  [0.75715975 0.20658867 1.08097842 0.33723951]\n",
      "  [0.24144962 0.07444944 0.33723951 0.11844495]]\n",
      "\n",
      " [[0.27977572 0.07067894 0.22215519 0.1197936 ]\n",
      "  [0.07067894 0.06540801 0.07187768 0.04916496]\n",
      "  [0.22215519 0.07187768 0.34701922 0.21183261]\n",
      "  [0.1197936  0.04916496 0.21183261 0.16424845]]]\n"
     ]
    }
   ],
   "source": [
    "#covariance matrix comparison\n",
    "print(params[1],'\\n\\n',sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.33332396280474647, 0.14300775583353068, 0.5236682813617229] \n",
      "\n",
      " [0.3331781  0.28909928 0.37772262]\n"
     ]
    }
   ],
   "source": [
    "#prior probability comparison\n",
    "print(params[2],'\\n\\n',weight_gaussian)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[9.99999988e-001 1.20274201e-008 4.70792632e-015]\n",
      " [9.99991268e-001 8.73143775e-006 1.73316927e-010]\n",
      " [9.99998946e-001 1.05387904e-006 6.06216278e-012]\n",
      " [9.99999987e-001 1.19230384e-008 7.55110179e-010]\n",
      " [1.00000000e+000 6.58718716e-013 1.37859263e-015]\n",
      " [1.00000000e+000 2.07968759e-015 5.27580548e-017]\n",
      " [1.00000000e+000 6.06936521e-015 3.50570830e-012]\n",
      " [9.99999997e-001 3.13047774e-009 2.91760507e-013]\n",
      " [9.99993973e-001 6.01165782e-006 1.51959895e-008]\n",
      " [9.99999941e-001 5.89408539e-008 2.65089953e-010]\n",
      " [9.99999999e-001 7.11578311e-010 5.36101035e-017]\n",
      " [1.00000000e+000 2.75470166e-015 7.09083463e-012]\n",
      " [9.99999408e-001 5.91338630e-007 4.89288429e-010]\n",
      " [9.99997245e-001 2.75356546e-006 1.12655397e-009]\n",
      " [9.99999982e-001 1.83813966e-008 5.26903077e-023]\n",
      " [1.00000000e+000 2.51143585e-023 3.37683608e-025]\n",
      " [9.99999999e-001 7.00732678e-010 7.97072229e-020]\n",
      " [9.99999760e-001 2.40292809e-007 5.65485655e-015]\n",
      " [9.99999968e-001 3.20118021e-008 6.08416726e-017]\n",
      " [1.00000000e+000 4.04454947e-018 6.06457269e-017]\n",
      " [9.99999877e-001 1.23022844e-007 3.10505303e-012]\n",
      " [1.00000000e+000 5.60186091e-013 1.29140622e-015]\n",
      " [1.00000000e+000 3.15039604e-014 5.42771586e-015]\n",
      " [9.99296727e-001 7.03272817e-004 5.84622844e-010]\n",
      " [9.99999993e-001 8.85065349e-023 7.25794890e-009]\n",
      " [9.99995907e-001 4.08995150e-006 3.15717953e-009]\n",
      " [9.99999800e-001 2.00216121e-007 4.88876325e-012]\n",
      " [9.99999987e-001 1.25398997e-008 1.24115192e-014]\n",
      " [9.99999210e-001 7.90041093e-007 2.10281215e-014]\n",
      " [1.00000000e+000 4.94302342e-011 4.13394711e-010]\n",
      " [9.99999640e-001 3.58719082e-007 1.15119280e-009]\n",
      " [9.99979484e-001 2.05164518e-005 8.57898024e-014]\n",
      " [1.00000000e+000 4.50257541e-034 7.07984803e-020]\n",
      " [1.00000000e+000 1.99427895e-021 3.29492224e-023]\n",
      " [9.99993502e-001 6.49810931e-006 1.25533395e-010]\n",
      " [9.99991830e-001 8.16988651e-006 3.43934954e-013]\n",
      " [9.99999971e-001 2.92919408e-008 1.91072324e-016]\n",
      " [1.00000000e+000 1.38305530e-017 5.63135936e-015]\n",
      " [9.99998670e-001 1.32852935e-006 1.47657685e-009]\n",
      " [9.99999950e-001 4.99900429e-008 1.63222562e-013]\n",
      " [9.99999863e-001 1.36747975e-007 3.97417534e-015]\n",
      " [9.99985013e-001 2.52493351e-012 1.49871361e-005]\n",
      " [1.00000000e+000 3.44996177e-012 1.02051396e-010]\n",
      " [9.99999997e-001 3.10641877e-009 8.34255943e-011]\n",
      " [1.00000000e+000 1.60774680e-026 8.03143471e-013]\n",
      " [9.99568175e-001 4.31824779e-004 3.62894083e-010]\n",
      " [1.00000000e+000 1.11052424e-022 3.24096905e-016]\n",
      " [9.99999999e-001 1.24644049e-009 4.20305428e-011]\n",
      " [1.00000000e+000 2.27935085e-011 9.41867526e-017]\n",
      " [9.99998795e-001 1.20491096e-006 4.34212961e-013]\n",
      " [1.22819335e-092 9.65557238e-001 3.44427615e-002]\n",
      " [7.88415573e-083 6.64764718e-006 9.99993352e-001]\n",
      " [9.61020142e-105 8.95108864e-001 1.04891136e-001]\n",
      " [2.85563719e-064 6.71719289e-001 3.28280711e-001]\n",
      " [1.79744556e-092 7.34681485e-001 2.65318515e-001]\n",
      " [4.56976318e-079 3.56019484e-014 1.00000000e+000]\n",
      " [2.42689241e-093 4.27934862e-018 1.00000000e+000]\n",
      " [1.46830280e-033 1.71267805e-002 9.82873219e-001]\n",
      " [1.91502500e-086 7.49377050e-001 2.50622950e-001]\n",
      " [1.67592339e-060 1.36326411e-016 1.00000000e+000]\n",
      " [3.17456025e-042 8.23701805e-001 1.76298195e-001]\n",
      " [3.37745958e-072 3.46819676e-010 1.00000000e+000]\n",
      " [2.31424494e-060 7.30528823e-011 1.00000000e+000]\n",
      " [1.80316331e-090 3.30922178e-008 9.99999967e-001]\n",
      " [2.16564440e-046 2.44975928e-005 9.99975502e-001]\n",
      " [1.83834223e-079 9.59597953e-001 4.04020472e-002]\n",
      " [1.72030889e-083 4.50989124e-030 1.00000000e+000]\n",
      " [3.11707982e-058 1.88942964e-004 9.99811057e-001]\n",
      " [1.78738554e-092 1.82496461e-007 9.99999818e-001]\n",
      " [1.26768950e-054 7.17060900e-001 2.82939100e-001]\n",
      " [2.54018876e-106 2.19704477e-040 1.00000000e+000]\n",
      " [5.70371948e-062 8.43735012e-001 1.56264988e-001]\n",
      " [2.79460091e-108 8.99939397e-001 1.00060603e-001]\n",
      " [1.42398547e-086 4.57606293e-008 9.99999954e-001]\n",
      " [1.27541760e-073 9.35599971e-001 6.44000294e-002]\n",
      " [5.60177050e-080 9.30801063e-001 6.91989366e-002]\n",
      " [9.11891994e-100 2.49005851e-001 7.50994149e-001]\n",
      " [3.35989427e-115 1.71310085e-002 9.82868992e-001]\n",
      " [7.00095058e-085 5.69075670e-008 9.99999943e-001]\n",
      " [1.35555344e-039 8.36354180e-001 1.63645820e-001]\n",
      " [3.64316300e-052 8.73276669e-001 1.26723331e-001]\n",
      " [1.29036745e-046 7.32110713e-001 2.67889287e-001]\n",
      " [7.39559493e-056 7.84063873e-001 2.15936127e-001]\n",
      " [1.50629142e-116 9.93845516e-012 1.00000000e+000]\n",
      " [1.60158615e-083 1.16860484e-042 1.00000000e+000]\n",
      " [1.04176230e-083 3.54792048e-034 1.00000000e+000]\n",
      " [4.14694266e-094 6.56651994e-001 3.43348006e-001]\n",
      " [1.99229296e-082 1.32661577e-007 9.99999867e-001]\n",
      " [6.47547482e-062 4.09655406e-017 1.00000000e+000]\n",
      " [2.78232308e-062 3.64402878e-002 9.63559712e-001]\n",
      " [1.14712900e-073 1.56633014e-011 1.00000000e+000]\n",
      " [3.45887530e-085 1.03907705e-009 9.99999999e-001]\n",
      " [2.19900416e-060 8.54281949e-001 1.45718051e-001]\n",
      " [5.80872650e-035 8.27680203e-001 1.72319797e-001]\n",
      " [5.19983745e-068 4.47215933e-007 9.99999553e-001]\n",
      " [4.95408225e-063 6.35460693e-016 1.00000000e+000]\n",
      " [1.65527612e-066 4.45321005e-011 1.00000000e+000]\n",
      " [6.87607335e-072 4.78218687e-001 5.21781313e-001]\n",
      " [3.76952730e-028 4.97033092e-002 9.50296691e-001]\n",
      " [1.03513031e-063 7.93079787e-006 9.99992069e-001]\n",
      " [1.39206900e-203 2.83842899e-086 1.00000000e+000]\n",
      " [2.28067343e-128 9.42624062e-025 1.00000000e+000]\n",
      " [5.59261525e-181 7.73146184e-007 9.99999227e-001]\n",
      " [1.85747114e-148 1.00874410e-021 1.00000000e+000]\n",
      " [9.88020068e-178 4.25499074e-030 1.00000000e+000]\n",
      " [5.87114399e-228 7.80615665e-002 9.21938434e-001]\n",
      " [4.79633889e-095 1.45387278e-041 1.00000000e+000]\n",
      " [1.00476083e-195 1.50006562e-003 9.98499934e-001]\n",
      " [5.93349412e-167 9.11726232e-001 8.82737685e-002]\n",
      " [2.28989198e-207 1.14218835e-047 1.00000000e+000]\n",
      " [7.76267459e-130 2.26229765e-022 1.00000000e+000]\n",
      " [1.11981076e-139 6.39211045e-007 9.99999361e-001]\n",
      " [1.08016048e-157 2.26827749e-012 1.00000000e+000]\n",
      " [3.60210601e-130 6.07741181e-022 1.00000000e+000]\n",
      " [3.37748833e-156 1.00226029e-057 1.00000000e+000]\n",
      " [1.26703704e-155 2.70023736e-044 1.00000000e+000]\n",
      " [6.14308898e-143 2.07683246e-015 1.00000000e+000]\n",
      " [4.49920241e-228 8.75667152e-039 1.00000000e+000]\n",
      " [7.96751246e-266 9.78456123e-001 2.15438774e-002]\n",
      " [1.54394607e-114 9.67061353e-001 3.29386468e-002]\n",
      " [5.39173328e-177 1.55990594e-025 1.00000000e+000]\n",
      " [1.33718089e-123 4.02458066e-041 1.00000000e+000]\n",
      " [4.88632048e-236 9.70560955e-001 2.94390446e-002]\n",
      " [3.76418107e-116 1.26652300e-005 9.99987335e-001]\n",
      " [1.83424561e-164 3.38937934e-032 1.00000000e+000]\n",
      " [2.03867515e-172 2.55341509e-009 9.99999997e-001]\n",
      " [1.94847808e-110 4.35817276e-009 9.99999996e-001]\n",
      " [1.10524539e-112 2.77525492e-020 1.00000000e+000]\n",
      " [1.85015889e-163 3.83100202e-018 1.00000000e+000]\n",
      " [4.52499367e-157 2.72364548e-003 9.97276355e-001]\n",
      " [1.38918116e-190 9.61628572e-001 3.83714277e-002]\n",
      " [7.04926218e-201 2.27103675e-017 1.00000000e+000]\n",
      " [9.17523022e-169 4.74287702e-022 1.00000000e+000]\n",
      " [1.39740755e-113 2.35839150e-006 9.99997642e-001]\n",
      " [1.61436013e-137 9.44229310e-021 1.00000000e+000]\n",
      " [7.02602729e-207 1.75754312e-007 9.99999824e-001]\n",
      " [9.32993653e-175 2.00750785e-077 1.00000000e+000]\n",
      " [1.68765702e-141 1.54870520e-025 1.00000000e+000]\n",
      " [3.71109079e-108 6.52585677e-023 1.00000000e+000]\n",
      " [7.91006814e-152 1.66654912e-012 1.00000000e+000]\n",
      " [3.11166083e-178 4.01704516e-034 1.00000000e+000]\n",
      " [2.57070399e-148 2.79944555e-026 1.00000000e+000]\n",
      " [2.28067343e-128 9.42624062e-025 1.00000000e+000]\n",
      " [1.34435594e-187 4.33144112e-032 1.00000000e+000]\n",
      " [1.14430499e-187 1.67195567e-051 1.00000000e+000]\n",
      " [1.58374617e-153 2.25886174e-026 1.00000000e+000]\n",
      " [4.96827858e-127 1.31518701e-006 9.99998685e-001]\n",
      " [1.03994643e-136 1.46831479e-015 1.00000000e+000]\n",
      " [1.56194051e-158 9.73622864e-074 1.00000000e+000]\n",
      " [1.65328816e-121 2.31477311e-035 1.00000000e+000]] \n",
      "\n",
      " [[1.00000000e+000 1.14273414e-010 7.77496975e-019]\n",
      " [9.99999823e-001 1.77461434e-007 1.51077192e-014]\n",
      " [9.99999956e-001 4.38418937e-008 1.68629010e-016]\n",
      " [9.99999815e-001 1.84901535e-007 2.53860195e-013]\n",
      " [1.00000000e+000 5.40506533e-011 1.88265355e-019]\n",
      " [1.00000000e+000 8.09594967e-011 9.44135727e-021]\n",
      " [9.99999890e-001 1.10268343e-007 6.21198302e-017]\n",
      " [9.99999999e-001 7.03863901e-010 1.50346743e-016]\n",
      " [9.99995535e-001 4.46512399e-006 9.88431891e-013]\n",
      " [9.99999994e-001 6.13658713e-009 4.97307119e-013]\n",
      " [1.00000000e+000 2.09246484e-012 3.96342360e-020]\n",
      " [9.99999997e-001 2.52517159e-009 9.61447916e-015]\n",
      " [9.99999971e-001 2.94639993e-008 2.26865781e-013]\n",
      " [9.99998012e-001 1.98799301e-006 7.95998772e-015]\n",
      " [1.00000000e+000 8.13099591e-015 1.24460720e-027]\n",
      " [1.00000000e+000 4.20041788e-015 1.67180874e-030]\n",
      " [1.00000000e+000 2.90119066e-011 6.21379988e-026]\n",
      " [9.99999999e-001 1.40312167e-009 1.47719310e-019]\n",
      " [1.00000000e+000 6.76861445e-012 1.23168526e-019]\n",
      " [1.00000000e+000 3.52927279e-011 4.17654533e-021]\n",
      " [9.99999999e-001 6.56799371e-010 2.31426470e-014]\n",
      " [9.99999998e-001 2.29521948e-009 1.67399013e-020]\n",
      " [9.99999975e-001 2.47068900e-008 1.50036723e-021]\n",
      " [9.99974805e-001 2.51952343e-005 2.19613373e-014]\n",
      " [9.99999972e-001 2.76550626e-008 1.77511786e-010]\n",
      " [9.99999745e-001 2.55036803e-007 3.48002083e-012]\n",
      " [9.99999769e-001 2.31470484e-007 2.79100064e-016]\n",
      " [1.00000000e+000 8.03010281e-011 7.91315575e-018]\n",
      " [1.00000000e+000 3.03589866e-010 3.91817012e-018]\n",
      " [9.99999953e-001 4.69676447e-008 5.14076489e-013]\n",
      " [9.99999893e-001 1.06800777e-007 1.36250532e-012]\n",
      " [9.99999968e-001 3.19598044e-008 1.85093389e-018]\n",
      " [1.00000000e+000 1.90722139e-015 6.46765545e-023]\n",
      " [1.00000000e+000 8.51562510e-016 2.69236897e-027]\n",
      " [9.99999947e-001 5.34597562e-008 4.76921476e-014]\n",
      " [9.99999983e-001 1.66605494e-008 2.72293593e-018]\n",
      " [1.00000000e+000 3.50449027e-011 1.49544622e-020]\n",
      " [1.00000000e+000 9.37270224e-012 3.68786215e-018]\n",
      " [9.99997848e-001 2.15205690e-006 2.90601503e-014]\n",
      " [1.00000000e+000 4.70229489e-010 9.34084908e-017]\n",
      " [9.99999997e-001 2.93492726e-009 2.21816644e-020]\n",
      " [9.78014973e-001 2.19850269e-002 2.74323369e-012]\n",
      " [9.99999667e-001 3.32799715e-007 2.55783678e-015]\n",
      " [9.99946800e-001 5.31997288e-005 9.60115840e-017]\n",
      " [9.99999986e-001 1.35397992e-008 1.02313636e-015]\n",
      " [9.99997430e-001 2.56983194e-006 4.80612730e-015]\n",
      " [1.00000000e+000 2.90004470e-012 3.68642806e-019]\n",
      " [9.99999939e-001 6.11316199e-008 4.35513317e-015]\n",
      " [1.00000000e+000 3.10011614e-012 5.88080215e-020]\n",
      " [9.99999998e-001 2.45452039e-009 6.23925129e-017]\n",
      " [8.39332151e-092 6.39117049e-001 3.60882951e-001]\n",
      " [2.99783940e-083 4.89032910e-001 5.10967090e-001]\n",
      " [3.13110253e-104 6.81485158e-001 3.18514842e-001]\n",
      " [4.98795493e-064 8.89651628e-001 1.10348372e-001]\n",
      " [2.88313153e-092 5.60378092e-001 4.39621908e-001]\n",
      " [3.80124787e-079 3.48804741e-001 6.51195259e-001]\n",
      " [1.67314993e-093 2.77477333e-001 7.22522667e-001]\n",
      " [7.18881692e-034 9.57760269e-001 4.22397312e-002]\n",
      " [3.63568221e-086 5.98936639e-001 4.01063361e-001]\n",
      " [3.95566206e-060 1.58149982e-001 8.41850018e-001]\n",
      " [2.83983281e-042 9.94314521e-001 5.68547921e-003]\n",
      " [3.38354783e-072 2.42488782e-001 7.57511218e-001]\n",
      " [1.67282112e-060 9.07029189e-001 9.29708106e-002]\n",
      " [1.26099053e-090 5.07670785e-001 4.92329215e-001]\n",
      " [1.74259685e-046 8.84488430e-001 1.15511570e-001]\n",
      " [8.78005594e-079 6.45320425e-001 3.54679575e-001]\n",
      " [1.34783520e-083 6.29881399e-002 9.37011860e-001]\n",
      " [3.45881709e-058 2.46093882e-001 7.53906118e-001]\n",
      " [7.50583072e-093 9.40229484e-001 5.97705156e-002]\n",
      " [2.15785764e-054 7.80768054e-001 2.19231946e-001]\n",
      " [2.28759407e-106 2.31214657e-003 9.97687853e-001]\n",
      " [1.18080385e-061 7.76359162e-001 2.23640838e-001]\n",
      " [1.13154776e-107 8.82129309e-001 1.17870691e-001]\n",
      " [2.92463551e-086 5.04206042e-001 4.95793958e-001]\n",
      " [5.65459172e-073 6.39719582e-001 3.60280418e-001]\n",
      " [1.95559387e-079 6.24130443e-001 3.75869557e-001]\n",
      " [5.31615083e-100 7.32122777e-001 2.67877223e-001]\n",
      " [2.15933414e-115 3.06034574e-001 6.93965426e-001]\n",
      " [4.67306533e-085 2.63535216e-001 7.36464784e-001]\n",
      " [2.09936801e-039 9.28673699e-001 7.13263006e-002]\n",
      " [1.12210274e-051 8.90719347e-001 1.09280653e-001]\n",
      " [2.25519262e-046 8.67483424e-001 1.32516576e-001]\n",
      " [1.28019249e-055 7.79981645e-001 2.20018355e-001]\n",
      " [1.79858995e-116 5.39427689e-001 4.60572311e-001]\n",
      " [1.21365812e-083 2.41887836e-002 9.75811216e-001]\n",
      " [1.14398540e-083 1.44102534e-001 8.55897466e-001]\n",
      " [3.71421015e-094 5.75782885e-001 4.24217115e-001]\n",
      " [5.05727743e-083 9.48338037e-001 5.16619628e-002]\n",
      " [4.23173928e-062 3.42521976e-001 6.57478024e-001]\n",
      " [2.23902663e-062 6.99039269e-001 3.00960731e-001]\n",
      " [1.14614405e-073 4.14543755e-001 5.85456245e-001]\n",
      " [2.05755694e-085 4.43501631e-001 5.56498369e-001]\n",
      " [6.11826471e-060 7.71953942e-001 2.28046058e-001]\n",
      " [1.00312403e-034 9.84169601e-001 1.58303991e-002]\n",
      " [3.43150190e-068 4.37526206e-001 5.62473794e-001]\n",
      " [3.09305207e-063 3.25915805e-001 6.74084195e-001]\n",
      " [9.72904023e-067 3.79806381e-001 6.20193619e-001]\n",
      " [4.63414198e-072 5.96934061e-001 4.03065939e-001]\n",
      " [1.23657363e-028 9.98808279e-001 1.19172055e-003]\n",
      " [6.07103577e-064 5.00652051e-001 4.99347949e-001]\n",
      " [1.28603290e-203 8.52254165e-010 9.99999999e-001]\n",
      " [1.72624989e-128 8.42567054e-004 9.99157433e-001]\n",
      " [1.82116656e-180 1.52524270e-001 8.47475730e-001]\n",
      " [5.10972990e-148 6.61120869e-001 3.38879131e-001]\n",
      " [1.38267021e-177 2.38675062e-004 9.99761325e-001]\n",
      " [1.24877176e-228 9.99877746e-001 1.22254419e-004]\n",
      " [3.43255357e-095 1.63351502e-004 9.99836648e-001]\n",
      " [1.50381706e-196 9.99970303e-001 2.96972562e-005]\n",
      " [1.21250829e-166 9.96710054e-001 3.28994644e-003]\n",
      " [4.95030904e-207 1.56578052e-006 9.99998434e-001]\n",
      " [5.91591599e-130 5.99061108e-004 9.99400939e-001]\n",
      " [1.47561570e-139 4.39512903e-002 9.56048710e-001]\n",
      " [9.22191630e-158 1.26138764e-003 9.98738612e-001]\n",
      " [3.53829114e-130 2.80512513e-005 9.99971949e-001]\n",
      " [5.05103221e-156 5.61110771e-013 1.00000000e+000]\n",
      " [1.01081211e-155 1.90610925e-008 9.99999981e-001]\n",
      " [1.24631289e-142 5.42390305e-001 4.57609695e-001]\n",
      " [3.45255273e-228 9.99904095e-001 9.59053349e-005]\n",
      " [4.45751152e-265 9.99966916e-001 3.30835057e-005]\n",
      " [7.35805718e-114 9.86210267e-001 1.37897329e-002]\n",
      " [4.08012330e-177 8.44117727e-006 9.99991559e-001]\n",
      " [9.51628097e-124 2.25812474e-006 9.99997742e-001]\n",
      " [5.54494581e-236 9.99998884e-001 1.11631250e-006]\n",
      " [3.93045187e-116 3.17032006e-002 9.68296799e-001]\n",
      " [4.05616109e-164 6.33538205e-003 9.93664618e-001]\n",
      " [9.62609598e-173 9.97522334e-001 2.47766619e-003]\n",
      " [1.83817150e-110 1.16377806e-002 9.88362219e-001]\n",
      " [8.66244311e-113 8.11703230e-003 9.91882968e-001]\n",
      " [2.19393120e-163 1.09772047e-003 9.98902280e-001]\n",
      " [3.00040602e-157 9.97403580e-001 2.59641987e-003]\n",
      " [7.30903873e-190 9.98356146e-001 1.64385429e-003]\n",
      " [1.10114165e-201 9.99926571e-001 7.34288123e-005]\n",
      " [7.43529023e-169 2.66653348e-005 9.99973335e-001]\n",
      " [1.16828599e-113 8.14615953e-001 1.85384047e-001]\n",
      " [3.28369573e-137 9.97595755e-001 2.40424474e-003]\n",
      " [1.30698041e-206 7.30312037e-003 9.92696880e-001]\n",
      " [8.58854561e-175 1.84481906e-009 9.99999998e-001]\n",
      " [4.66309900e-141 4.32191958e-001 5.67808042e-001]\n",
      " [2.97480670e-108 3.73767448e-003 9.96262326e-001]\n",
      " [5.39877370e-152 7.26063124e-004 9.99273937e-001]\n",
      " [1.80227038e-178 1.30468985e-008 9.99999987e-001]\n",
      " [3.14299324e-148 4.72141409e-008 9.99999953e-001]\n",
      " [1.72624989e-128 8.42567054e-004 9.99157433e-001]\n",
      " [1.88049626e-187 5.23407525e-005 9.99947659e-001]\n",
      " [7.70390871e-188 3.42762442e-010 1.00000000e+000]\n",
      " [1.51887823e-153 3.63216445e-008 9.99999964e-001]\n",
      " [8.01330639e-127 1.08284088e-002 9.89171591e-001]\n",
      " [7.62443324e-137 9.94641796e-004 9.99005358e-001]\n",
      " [1.57562234e-158 1.76283888e-008 9.99999982e-001]\n",
      " [1.54485867e-121 7.85752725e-003 9.92142473e-001]]\n"
     ]
    }
   ],
   "source": [
    "#posterior probability comparison\n",
    "print(np.concatenate(\n",
    "    [posterior[i].reshape(-1,1) for i in posterior],\n",
    "    axis=1),'\\n\\n',prob_gaussian)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use pca to reduce dimensions\n",
    "#check the link below for more details of pca\n",
    "# https://github.com/je-suis-tm/machine-learning/blob/master/principal%20component%20analysis.ipynb   \n",
    "dims=2\n",
    "x=sklearn.decomposition.PCA(n_components=dims).fit_transform(df[xcols])\n",
    "\n",
    "df['dim1']=x[:,0]\n",
    "df['dim2']=x[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#prepare labels for viz\n",
    "df['skl']=np.argmax(prob_gaussian,axis=1)\n",
    "\n",
    "df['gmm']=np.argmax(\n",
    "    np.concatenate(\n",
    "        [posterior[i].reshape(-1,1) for i in posterior],\n",
    "        axis=1),axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAFNCAYAAABST1gVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X98XHWd9/3XJ5lMMvlRaZNWaGMJFZAWQSpJabWywIXiLr3gWgSJslYu2UfdmwXXtd0fiF50ueFeV1rXRddrrSs3FFmD96oLFAX0wop1221aKSClWCyhpKW0TQvk92SS7/3HmQmTNJnMJDNzzkzez8cjj8mcOTnnM9NCPv3++HzMOYeIiIiIBFuJ3wGIiIiIyMSUtImIiIgUACVtIiIiIgVASZuIiIhIAVDSJiIiIlIAlLSJiIiIFAAlbSIiaTKzi8ys3e84RGR6UtImIgXDzJrN7L/MrNvMDse/v9HMLP76vWbmzOyKUT/39fjx6+PPr48//9qo8/5H/Pi9+XpPIiLpUtImIgXBzFYD/wTcBZwMvBP4M+CDQDjp1N8Bn076uRBwDfD7UZf8PXBt/PWElfGfFxEJHCVtIhJ4ZvYO4HbgRufcvzvnOp3naefcdc65/qTTHwE+aGYz488/CjwLHBp12UPAc8Bl8XvMAj4APJxGPF80s6Nm1mZm18WPNZnZ68lJoJl9zMx2jXONzWb2p0nPrzezLfHvzcz+MT6a+KaZPWtm750oLhEpbkraRKQQLAPKgYfSOLcPL/Fqjj9fCWwc59yN8deJn/8Q0D/OuQknA3XAPLwRvQ1m9h7nXCvQAXw46dw/Ae5PI+bRPgJcCJwJnARcG7+2iExjStpEpBDUAUedc7HEATP7TzN7w8x6zezCUedvBFbGR+j+APiPca77Y+Ci+HmpkrvRvuyc63fO/RJ4FPh4/Ph9eIlaYuTuMuDf0rxmsgGgBjgLMOfcC8651yZxHREpIkraRKQQdAB1yVOPzrkPOOdOir824v9lzrktwGzgS8Am51zvWBeNH380fl6dc+7XacRy3DnXnfT8FWBu/PvvAf/dzKrxErlfTSbZcs49CXwT+GfgdTPbYGYzMr2OiBQXJW0iUgi24k1bXpnBz3wPWM3Eo2cb4+elO40508yqkp7PBw4COOcOxGP9Y+BTE1yzG6hMen5y8ovOubudc+cDZ+NNk/5VmvGJSJFS0iYigeecewP4O+BbZna1mVWbWYmZnQdUjfNjd+OtL3tqgsv/Mn7eNzII6e/MLGxmHwJWAP9f0msbgb8GzsGbfh3PLuAqM6s0s9OBGxIvxDc1XGBmZXjJXR8wmEF8IlKEQhOfIiLiP+fcV83sAF5CtBEvmdkH/A3wn2Ocfwz4P2lc16VzXpJDwHG80bUe4M+cc3uSXv8x8L+BH4+aRh3tH4Em4HW83a0PAJfGX5sRf30BXsL2OLAugxhFpAiZ9/8rERHJFjP7PfBZ59zP/Y5FRIqHpkdFRLLIzD4GOOBJv2MRkeKi6VERkSwxs83AIuBTzrkhn8MRkSKj6VERERGRAqDpUREREZECoKRNREREpAD4uqbNzO7Bq3F02Dl3QjNkM7sIrxfgy/FDP3LO3Z7qmh/96EfdY489lu1QRURERHLB0j3R740I9+K1aklVsfxXzrkV6V7w6NGjU41JREREJHB8nR51zj0FHPMzBhEREZFCUAhr2paZ2TNm9lMzO9vvYERERET84Pf06ER+A5zqnOsysz8C/gM4Y/RJZrYKWAUwf/78/EYoIiIikgeBHmlzzr3lnOuKf/8ToMzM6sY4b4NzrtE51zh79uy8xykiIiKSa4FO2szsZDOz+PdL8OLt8DcqERERkfzzu+TH94GLgDozawduA8oAnHP/AlwN/F9mFgN6gWanFg4iIiIyDfmatDnnPjHB69/EKwkiIiIiMq0FenpURERERDxK2kREREQKgJI2ERERkQIQ9DptIiIiUkRaD7XSsqeF9s526mvqaT6rmaaTm/wOqyBopE1ERETyovVQK+ta19HR20FdpI6O3g7Wta6j9VCr36EVBCVtIiIikhcte1qIhCJUh6sxM6rD1URCEVr2tPgdWkFQ0iYiIiJ50d7ZTlVZ1YhjVWVVtHe1+xRRYdGaNpnQtn0dbNzaxv5jPcyfVcnKZQ0sXVDrd1giIlJg6mvq6ejtoDpcPXyse6Cb+up6H6MqHBppk5S27evgjk27OdoZZXZ1OUc7o9yxaTfb9qmbmIiIZKb5rGZ6Y710RbtwztEV7aI31kvzWc1+h1YQlLRJShu3tlEZDlFTEaLEjJqKEJXhEBu3tvkcmYiIFJqmk5tY07SG2kgtR/uOUhupZU3TGu0eTZOmRyWl/cd6mF1dPuJYVXkp+4/1+BSRiIgUsqaTm5SkTZJG2iSl+bMq6e4fHHGsu3+Q+bMqfYpIRERkelLSJimtXNZATzRGZ1+MIefo7IvRE42xclmD36GJiIhMK0raJKWlC2r50opF1NWEOdLVT11NmC+tWKTdoyIiInmmNW0yoaULapWkiYiI+EwjbSIiIiIFQEmbiIiISAFQ0iYiIiJSAJS0iYiIiBQAJW0iIiIiBUBJm4iIiEgBUMmPIrVtXwcbt7ax/1gP82dVsnJZg8p2iIhI0Wg91ErLnhbaO9upr6mn+azmom+PpZG2IrRtXwd3bNrN0c4os6vLOdoZ5Y5Nu9m2r8Pv0ERERKas9VAr61rX0dHbQV2kjo7eDta1rqP1UKvfoeWUkrYitHFrG5XhEDUVIUrMqKkIURkOsXFrm8+RiYiITF3LnhYioQjV4WrMjOpwNZFQhJY9LX6HllNK2orQ/mM9VJWXjjhWVV7K/mM9PkUkIiKSPe2d7VSVVY04VlVWRXtXu08R5YeStiI0f1Yl3f2DI4519w8yf1alTxGJiIhkT31NPd0D3SOOdQ90U19d71NE+aGkrQitXNZATzRGZ1+MIefo7IvRE42xclmD36GJiIhMWfNZzfTGeumKduGcoyvaRW+sl+azmv0OLaeUtBWhpQtq+dKKRdTVhDnS1U9dTZgvrVik3aMiIlIUmk5uYk3TGmojtRztO0ptpJY1TWuKfveoOef8jiGrGhsb3Y4dO/wOQ0RERCQdlu6JqtM2Dahmm4iISOHT9GiRU802ERGR4qCkrcipZpuIiEhxUNJW5FSzTUREpDgoaStyqtkmIiJSHJS0FTnVbBMRESkOStqKnGq2iYiIFAeV/JgGli6oVZImIiJS4HwdaTOze8zssJn9dpzXzczuNrOXzOxZM3t/vmMUERERCQK/p0fvBT6a4vU/BM6If60C/nceYhIREREJHF+TNufcU8CxFKdcCWx0nm3ASWZ2Sn6iExEREQkOv0faJjIPeDXpeXv8mIiIiMi0EvSkbawmqid0uDezVWa2w8x2HDlyJA9hiYiIiORX0JO2duBdSc/rgYOjT3LObXDONTrnGmfPnp234ERERETyJehJ28PAyvgu0qXAm8651/wOSkRERCTffK3TZmbfBy4C6sysHbgNKANwzv0L8BPgj4CXgB7gf/oTqYiIiIi/fE3anHOfmOB1B/x5nsIRERERCaygT4+KiIiICEraRERERAqCkjYRERGRAqCG8ZIT2/Z1sHFrG/uP9TB/ViUrlzWoab2IiMgUaKRNsm7bvg7u2LSbo51RZleXc7Qzyh2bdrNtX4ffoYmIiBQsJW2SdRu3tlEZDlFTEaLEjJqKEJXhEBu3tvkcmYiISOFS0iZZt/9YD1XlpSOOVZWXsv9Yj08RiYiIFD4lbZJ182dV0t0/OOJYd/8g82dV+hSRiIhI4VPSJlm3clkDPdEYnX0xhpyjsy9GTzTGymUNfocmIiJSsLR7tEgldm++8Npb9EYHqSgrYdHcd+RlF+fSBbV8acWiUbtHz9TuURERkSkwr1NU8WhsbHQ7duzwOwxfJXZvxoYc7cd6MDMcjvqTIoRKS/jSikVKoERERILB0j1RI21FKLF7s62jm1BpCaESIzY0xLHuARrqqti4tU1Jm4iIjKv1UCste1po72ynvqae5rOaaTq5ye+wpj2taStCid2bvQODlJZ4CXypGb0Dg9rFKSIiKbUeamVd6zo6ejuoi9TR0dvButZ1tB5q9Tu0aU9JWxFK7N6MlJUyOORNfw86R6SsVLs4RUQkpZY9LURCEarD1ZgZ1eFqIqEILXta/A5t2lPSVoQSuzdnVYWJDQ7RHxsiNuSYVVWmXZwiIpJSe2c7VWVVI45VlVXR3tXuU0SSoDVtRSh592Z/bIjeaIyKshIWzKnO+u5R9RgVESku9TX1dPR2UB2uHj7WPdBNfXW9j1EJaPeoTEFil2plOERVuTf12hONaXeqiEgBS6xpi4QiVJVV0T3QTW+slzVNa7QZITfS3j2q6VGZNPUYFREpPk0nN7GmaQ21kVqO9h2lNlKrhC0gND0qk7b/WA+zq8tHHNPuVBGRwtd0cpOStADSSJtMmnqMioiI5I9G2mTSVi5r4I5NuwFGrGlbuexMnyMTkbxp2wLbN8DxV2DmqbBkFTQs9zsqyTMV480PjbTJpCV2qdbVhDnS1U9dTVibEESmk7Yt8Pit0H0Uqud4j4/f6h2XaUPFePNHI20yJUsX1CpJE5mutm+AcBWU13jPE4/bN2i0bRpJLsYLDD+27GnRaFuWaaRNREQm5/grXtKWLFzlHZdpQ8V480cjbdOciuOKyKTNPNWbEk2MsAFEu73jMm2oGG/+aKRtGksUxz3aGWV2dTlHO6PcsWk32/Z1+B2aiBSCJau8JK2/E9yQ9xjt9o7LtNF8VjO9sV66ol045+iKdtEb66X5rGa/Qys6StqmMRXHFZEpaVgOl90JVXXQddh7vOxOrWebZlSMN380PTqNqTiuiExZw3IlaaJivHmipK2ITbRebf6sSo52RqmpePuvgYrjioiIBJOmR4tUOuvVVi5roCcao7MvxpBzdPbF4sVxG/wLXERERMakkbYilbxeDRh+TKxXS4zAVZWHAMeRrv74aNyZ2j0qIiISQEraitR469V2H3yTOzbtpjIcYnZ1ebz11KA6GYiIiASckrYiNd56tb6BIebUjD0Cp6RNRKSwqOfn9KI1bUVqvPVqkXApVeWlI87VjlERkcIzXXp+th5qZfXm1Vz7yLWs3ry66N5fJpS0FanxmrkvPGUG3f2DI87VjlERkcKT3PPTzKgOVxMJRWjZ0+J3aFkzXRLTdGl6tIiN18z9jk27AW+EzVvTFmPlsjPzHZ6IiExBe2c7dZG6EceKreenmtGP5OtIm5l91MxeNLOXzOxvx3j9ejM7Yma74l9/6kecxWDbvg5ufGAndzy6m6ryUswYMQKn9WwiIoWlvqae7oHuEceKreenmtGP5FvSZmalwD8DfwgsAj5hZovGOPVB59x58a9/zWuQRWJ0zTbnjNff7KUqHGL/sR42bm1Tv1ERyZ+2LfCDlfDtP/Ae27b4HVFBmg49P6dDYpoJP0falgAvOef2OeeiQAtwpY/xFK3RPUYHBgc53BVl7+udahQvIvnVtgUevxW6j0L1HO/x8VuVuE3CdOj5OR0S00z4uaZtHvBq0vN24IIxzvuYmV0I/A74S+fcq2OcIymMrtl24HgfZaXGwJAbbhQPKvshInmwfQOEq6C8xnueeNy+QT1MJ6HYe34mEtOWPS20d7VTXz29y5r4mbTZGMfcqOePAN93zvWb2Z8B9wGXnHAhs1XAKoD58+dnO86CN7pmW+/AICUGkbK3S3+o7IdIEWrb4iVDx1+BmafCklX+J0bHX/FG2JKFq7zjImMo9sQ0E35Oj7YD70p6Xg8cTD7BOdfhnOuPP/0OcP5YF3LObXDONTrnGmfPnp2TYAvZ6JptoVJjYNBRPzMyfI7KfogUmaBOQ848FaIj1ygR7faOT0OqQTZF02x9pJ9JWytwhpmdZmZhoBl4OPkEMzsl6ekVwAt5jK9ojK7ZduY7q3nnO8oJlZSoUbxIsUqehrQS7zFc5R3305JVXpLW3wluyHuMdnvHpxnVIJuioP7DJId8mx51zsXM7CbgcaAUuMc597yZ3Q7scM49DHzOzK4AYsAx4Hq/4i10yTXbtu3r4B9/9iIvvPYWAAtPmZHVsh/b9nUMN6T3mtA3aK2cSL4FdRqyYTlcdufIaduLv+j/tK0PVINsiqbh+khfi+s6534C/GTUsf+V9P0twC35jquYJcp/VIZDnH/qTLr7B+nuj+Xk+sk7U1ULTiTPZp7qjTwkfpFBcKYhG5YX7S/VTEyH4rg5FdR/mOSQ2lhNM6PLf9RUhKgMh9i4te2EcxMFeVd841fc+MDOtEqCZHJ9EckhTUMGnmqQTdE0XB+ppG2a2X+sJ62G8aML8qZbyy3d64tIjiWmIavqoOuw93jZnVMf4SqShd9B2ACgGmRTNA3/YaKkbZqZP6syrYbxE42YjTcKl+71RSQPGpbDxzfCZ3/pPWYjYSuChd9B2QAwHYrj5lSu/mESYObc6NJoha2xsdHt2LHD7zACK3nNWXLD+NFrzlZ841fMri6nxN4upzfkHEe6+vnS5YvGvQaQ1vVFpAD9YOWJ6+T6O71flh/f6F9cGVq9eTUdvR3DC/8BuqJd1EZqWX/Ret/iaj3U6hWR7WynvmZ6F5GdZsaqWzsmjbRNQ1XlIfYceoudrxzHzI2ZUKUaMUs1Cje6vIga0osUkeOveAu9kxXgwu+pNCHP1bSqr6N/RTLlPR0oaZtGEqNszsH5p87krJNnnJCYJYwuyJtcy22idWtLF9TyrevOZ9PNH+Jb152vhE2kWBTJwu/JbgDIZWKVXP7DzKgOVxMJRWjZ0zLla6dUJFPe08W4SZuZzTCzvzez+83sk6Ne+1buQ5Nsy2RnZ6oRM61bE5mmimTh92Q3AOQysZrK6N+UpCjCHITNGjJSqjpt/y+wF/gh8Bkz+xjwyXhbqaX5CE6ya3TjeEi9szO5IG+ylcsauGPT7uGfT6xbW7nszOwHLSLBUSSFcSfbhHwyddXSXadWX1N/wjq7dMt/TGkt3Di1zlrf2Mu61nVEQpERo4raKOGvVEnbu51zH4t//x9mdivwZLxDgQREJt0HRjeOh8mNkCVG4Ube90xNg4oERS4bxRdJYdzJNCHPNLFKTKemk/g0n9XMutZ1gJcIdg90pzX6l8k9xjROEeaWSEjdGgIoVdJWbmYlzrkhAOfcnWbWDjwFVKf4OcmT0d0H9h3pZtXGHcyqKuOdMyoAozsaG07msjlCNt4onIj4LLFGKVxFa2UVLT17aP/FTdTPeR/N59+kX7hTkGlilUmbqsmO/k25FdaSVd7fF/CmRaPdEO2mfcYM6vyYrpWUUm1EeAS4JPmAc+4+YDUQzWVQkp7kNWpv9g5w4HgvzsHx7ii7Xn2THW3HeOn1Trb9voNbfvQsgHZ2ihS7+Bql1rIS1g0docOgzkJ0dLyoZuRTlGldtUzWqY2Y4kwzYcv0HmMap9ZZfe1Z6tYQQOOOtDnn/nqc448BZ+QsIklb8hq19uO9lJYYpQbHe71eombGkIMhB6+/2c8//ux3PPjZZUrSRIpZfI1Sy+ABIhjVVgpWQnUsCvFF8xptm7xMplXTnU6dyhTnVNbCDRtjyru5onxS07WSWyr5UcCSd3H2DgxSWmIMOodzXqW+EoNBB6ESo6zUeOG1t/wNWERyL16Wo90NUJX4X/zQIJRF0lo0r92C2ZPuLtWp7ErNVSssdWsIplRr2iRgRm86aDx1Jj/6zQEAKkIl9A0MAWDmJW3OQWmio4EBxdX8QkTGEl+jVF9udNgQ1UPOS9pqT83aonlJT7rr1DLdlTp6t+iKBSvYdWRXRmvh0o1ff/bBoqStQCRvOgiZsfX3HTy55zCn1VVh5qiuKKM/1s/cd1Tw2pu9dPYPgnNEykqIDQ0xMOg4Z94Mv9+GiORafI1S89b1rOt7GUorqJp1Ot2hcNYWzedSsbVySifxyWSKc6zketO+TUqusy2XO7CnIK3pUTP7gJl90sxWJr5yHZiMlNh0MDA4yO+PdOOco7y0hNff7KO7f5C/v+ocNqxsZMGcamZWhSkPlVAZ9nJyM2NOdZi//PB7/H0TIpIfDctp+sSPWfOH36F2/gc5WjKU1UXzuRKURu75lskUp2+dE6aTAHeJmHCkzczuB94N7AISZfAdUDjdgYtAYtPB8we6KS2BUEkJDogODg13NUhuGZVJ/TYRKU65WDQ/kamMlKUz2ldsI3GQWbmPyRT4lQwld4mAtx+3b/B9tC2d6dFGYJFzTiuifJQojNs7MEi41FunNjjkiJSVjtnVYHQdtW37OrjxgZ1K4kTkbUlTQM0zZrEu7FVzmuxuwamui5soISnmdXfpJtjZSq4lhXG6RHD8FX/iSZLO9OhvgZNzHYiklmjgHio1Ys4RG3IMDjnqZ0Ym7GqQWA93tDPK7OpyjnZGuWPTbrbt68jjOxCRQBk1BdTU18+ajmPUDrlJ7xac6tTdRI3cgzg1mO8dt7naLSpJ4juwR4h2e8d9lk7SVgfsNrPHzezhxFeuA5OREq2jznxnNf2xIUoMTp9TRaikJN7VoGHcn82kUbyITBNjNApvKpvJ+m7HgyseZP1F6zMevZrquriJEpIgrLtL5scaPJXiyIMlq7wkrb8T3JD3GO32jvssnenRtbkOQtKzdEEtD372AyPWq51yUnjCvp+ZNooXkeCb8tquHEwBTXXqbqK1XUGbGvRrx61KceRYoktE8u7Ri7/o+3o2SCNpc8790szeCST+hmx3zh3ObViSSqZ9PyfbKF6bGUSCKStru8ZpFD6VKaDJNj1Pliohycb1R5tK8jvV+mrFsImiaI3RJSIIJpweNbOPA9uBa4CPA/9lZlfnOjDJnsR6uM6+GEPO0dkXm3BKVevgRIIrK2u7cjAFlOupu2xff6rTmxOtwcvmvUQgvenRW4GmxOiamc0Gfg78ey4Dk6lLHimrKi8F4EhXf3zULPWUavI6OGD4cePWNo22ifgsK2UfcjQFlOupu2xef6rTm5mM/AWleLEUtnSStpJR06EdqGdp4CV3UJhdXU53/yA90RhfWrEoraRL6+BEAmCcquxZW9sV0CmgfJlq8qv6apJv6SRtj5nZ48D348+vBX6Su5AkG6Y6UjbZdXAikiWJkhzhqpFV2S+7Mydru6ajqSS/o9enrWlMPU0btE0UUpgmHDFzzv0VsAE4F3gfsME59ze5DizIEoVqV3zjV9z4wM5ArvNKnhJNmGikLPl9dXT1c6SrL6N1cCKSRWOU5CBcBds3qOxDlky25tlk1qepvppkQ1oN451zPwR+mONYCsLoacfEAv10px3zJdORsrGmU92Qwyz9dXAikkUTlORQ2YfJu//5+7l397281f8WkVCE2ZHZ9A72ppzeTDaZ9WmZTKWKjGfcpM3MtjjnlptZJ16v0eGXAOecm5Hz6AIoqAv0R5fnaDx1Jj/6zQHAG2FLrGlbuezMMX9+7PcVobY6zIPXLcvPmxCRt+WgJId4CdvdT99NyEJEQhEGBgfY37mfzy3+HJ86+1NpXWOy69OUaMtUjTs96pxbHn+scc7NSPqqma4JG0xu2jHXxirP8aPfHOCq98+jribMka5+6mrCKUcDpzqdGtRpYpGCFeCq7IXs3t33ErIQ4VAYMyMcChOyEPfuvjfta2RS6kMkm9Kp0/ZuMyuPf3+RmX3OzE7KfWjBNH9WJd39gyOO+b1Af7w2VTteOc63rjufTTd/iG9dd37KkcBM35fquInkWKIkR1UddB32Hi+7c1rv9syGt/rfoqy0bMSxstIy3up/K+1raH2a+CWdNW0/BBrN7HTgu8DDwL8Bf5TLwIJq5bIG7ti0G0hv2jEfslGeI9P3FdRpYpGiks+SHOOUFykk6XQcmFE+g55oD+FQePjYwOAAM8rTn0DS+jTxSzpJ25BzLmZmfwx83Tn3DTN7OteBBVWicfvI9k7ZX6CfSQupbJTnyPR9qY6bSBFJUV6kUBK3dFt7Xb/oeu5++m6IeSNsA4MDxFyM6xddn9H9tD5N/JBO0jZgZp8APg389/ixshTnF71Me39mKtMdqtka/cvkfamOm0gRSS4vAm8/bt9QMElbujs6E5sNErtHZ5TP4PpF16e9CUHET+l0NvifwDLgTufcy2Z2GvC93IY1vY23Rm3j1rYxz0+MkqW76SAbJtPPVEQC6vgrXtKWLKm8SCFo72ynqmzkexhvR+dZtWdx3uzzWPCOBcyvmc+T+5/k2keuZfXm1eoFKoE24Uibc2438Lmk5y8DX8llUNPdZKYecz36N9b98jFNLCJ5UATlRdLtOJA8jRoqCfF8x/MYxunvOH3cKVWRoJgwaTOzDwJrgVPj5yfqtC3IbWjTV6FMPeY7URQpZuksos+ZJau8NWzgjbBFu72vi7+Yn/tnQbqtvZKnUZ/veJ5wSRiH42DPQc6uPXv4HCVtEkTpTI9+F/gasBxoAhrjj1NmZh81sxfN7CUz+9sxXi83swfjr/+XmTVk475Bl82pR9VSEwm+ybRFyqoiKC+Sbmuv5GnUvlgfpSWllFopfbE+QE3cJdjS2YjwpnPup9m+sZmVAv8MfBhoB1rN7OH4dGzCDcBx59zpZtYM/ANew/qilq2px0JpuSUy3U2mLVLW5bO8SI6ks6MzeRq1IlTBwOAADkdFqAJQkVwJtnSStl+Y2V3Aj4D+xEHn3G+meO8lwEvOuX0AZtYCXAkkJ21X4k3NAvw78E0zM+dcclutopSNqUfVUhMpDJNtiySZS55GnVs1l71v7MUwGmoaAl0k19fpcwmMdKZHL8CbEv1/gPXxr3VZuPc84NWk5+3xY2Oe45yLAW8CyjbSFMSWWyJyIrVFyp/kadSYi3F27dksmrWIGLFxp1T95vv0uQRGOrtHL87RvW2s203iHMxsFbAKYP78+VOPrEhke0NDJgV/RSR96S6iT4dGZCZWaIVxAzF9LoGQTu/Rd5rZd83sp/Hni8zshizcux14V9LzeuDgeOeYWQh4B3Bs9IWccxucc43OucbZs2dnIbTikO0NDeo1KpIb6S6in4hGZApD66FWVm9enXZtuExq0ElxS2d69F7gcWBu/PmWrB7sAAAgAElEQVTvgM9n4d6twBlmdpqZhYFmvL6myR7G68QAcDXw5HRYz5Yt2Sy6m2nBXxHJTNPJTay/aD0PrniQ9Retn9QISvKIjJlRHa4mEorQsqclBxHLZEwmsdb0uSSksxGhzjn3AzO7Bby1ZWY2ONUbx69zE15CWArc45x73sxuB3Y45x7GKzdyv5m9hDfCFrzVoQGXrVpq6jUqMkoAG6xrQ0PwTWaqM5vT51LY0hlp6zazWuJrycxsKd6GgClzzv3EOXemc+7dzrk748f+VzxhwznX55y7xjl3unNuSWKnqeTf/FmVdPePzNWDWPBXJC8SDda7j45ssN62xdewNCITfJOZ6szW9LkUvnRG2r6AN035bjP7NTAbb6pS8iQIGwCy1ZRepCgEtMF6NkZktJEht9JttzVaoW2ekNyYcKQtXo/tD4APAJ8FznbOPZvrwMQTlA0AfjSlFwmsgDZYn+qIjDYy5F7zWc30xnrpinbhnAt0bTgJHptoXX+8c8HlQANJI3POua/lNLJJamxsdDt27PA7jKy58YGdJ5Tt6OyLUVcT5lvXnX/C+UEYlRMpej9YeWKD9f5Or/3TxzeecHqhjF6t3rz6hFGgrmgXtZFa1l+03sfIisvw34euduqrg/v3QfJmrPJmY0pnevQRoA94DhiabEQyOZlsABjdtmrf4S5WbdxBbXWYhafMUAInki0ZNFhPjF5FQpERo1d+rEmaKHmc7EaGQklKg0JTnTJZ6WxEqHfOXeWcu80593eJr5xHJkBmGwCSy3K80RPlwBt9OAddfTHVVRPJpgwarAelDEc6U5+T2cigKVWR/ElnpO2nZvYR59wTOY9GTpjebDx1Jj/6zQFg4g0AyaNyB473UVoCpSUl9MWGxu07qulUkUlKs8F6UMpwpFNqYjIbGTItYaFROZHJS2ekbRvwYzPrNbO3zKzTzN7KdWDT0VibDn70mwNc9f55aW0ASB6V6x0YpNSMwSFHpMzrPzp6WjUomxxEillQynCkU2piMhsZMilhoVE5kalJZ6RtPbAMeE7dCHIreXoTGH7c8crxMTcdjJZclqOirIS+AW8JYv3MCHDitOp49xs9Gicik5erwqiZjlilW2pi9HqrRMul8e6TSQkL9dAUmZp0Rtr2Ar9VwpZ7+4/1UFVeOuJYOl0Htu3r4MYHdnLHo7upKi/FjOHHeTMjvCNSNmbf0cneT0TSl4vCqJMZsZpMqYmJ7tN6qJVjvcd4vuN5njnyDMf7jqe8bpB7aGbaD1TED+mMtL0GbI43jO9PHAxqyY9CNn9W5QnlPSbqOjB6x2hizdtXPvY+gFHr1c4cMYI2mfuJSOayvVswnRGrsUbi1jStyajURKr7AMO7Ys886Uz2d+3nd8d/x3vr3jtuUjrZwrK5FqQdviKppDPS9jLwf4AwUJP0JVm2clkDPdEYnX0xhpwbc3RstKk0cp/M/UTEfxONWI03QgZk1JQ+1X2SE7qZkZm8b/b7OLv2bGZVzErZQzOIhWWDssNXZCITjrSpvEf+JLoOpBodG228Om67D76ZVs22TO8nIv6baMQqW2vHUt0n3V2xo0f8VixYwa4juwJVWDbTHb7aASt+GTdpM7OvO+c+b2aPEG8Wn8w5d0VOI5umli6ozShpGm+Ks29giDk13gjcse5+DrzRB4ys2ZbYhaokTSR9QfiFPdHmhmyVGUl1n5Y9LRNOdY417bhp36a8TDtm8ueUybTt6Pf08hsv8xdP/gWzKmbxnlnvUQInOZVqevT++OM6vB2ko78kDxKbDFZ841fc+MDOE8pxjDfFGQmXDm8ySNRsC4fertmW7hSqiLwtKCUrJtrckK0yI6nuk85Up1/Tjpn+OWUybZv8nt7of4PXel7zfmagSyVMJOfGHWlzzu2MP/7SzGbHvz+Sr8DkxE0GRzuj3PKjZ5lTU053dHC4GO5YU5wbt7YNj8D1DgwSLk1ds01EJhakkhWpNjdks8zIePdJJHSpNjb4VVg40z+ndN5LQvJ7Oth9kBIrIVQaon+wf+K/D21bYPsGOP4KzDzVa4eWRoFmkYRU06MG3AbchNfMtMTMYsA3nHO35ym+aW10HbXY0BCvv9nPGz0DnPeuk4anOa96/7wTfjbTmm0iQRSEqchkQeluMJFMkpCp3icbteGybTJ/Tunu8E1+T32xPspKyhgcGqQiVJH6Pm1bvH614SqongPdR73n47Q/ExlLqunRzwMfBJqcc7XOuZnABcAHzewv8xLdNDe6jlr78V7KSo3YoBveKRobHOLrP997QlcDgC+tWERdTTitmm0iQROUqchkQelukI6mk5sy2imaC37tFs3ln1Pye6oorWBgcIBBN8i86nmp77N9g5ewldeAlXiP4SrvuEiaUiVtK4FPOOdeThxwzu0D/iT+muTY6GbxvQODYAxPcQJ0dEeJDbkxS34sXVDLt647n81/dQkbVjayYHbVhK2wRILC7zIMYxVbDWrJiqDKRWHhdOTyzyn5PVWVVVFiJcytnss7wu9IfZ/jr3hJWrJwlXdcJE2pSn6UOeeOjj7onDtiZmU5jEnikqc4q8pLKSsx+mJDLKirGD6nOzpIVXjirgbaJSqFxs+pyFTFVvMx7ViU8thTJ9fTw8lTqcNT+BPdZ+ap3pRoeVKZ02i3d1wkTTZedyoz+41z7v2Zvua3xsZGt2PHDr/DyJpt+zqGNxlUhUO8/mYvc2ZEqCovpbt/kD2H3mLuSRHmnRQZ/pnOvhh1NeG0+pWKBNXqzatPWA/VFe2iNlLL+otyu4E93XsHbc1d0CQnv8kbIqZlp4HkNW3hKi9hi3ZrTZuAt28gLammR99nZm+N8dUJnDP1GCUdiSnOTTd/iAc/u4yvXP0+6mrCw9Ocn7/0DEIlpq4GUnT8nIpMp0dmENfcBY3fU9yB0rDcS9Cq6qDrsPeohE0ylKrkR+l4r4l/kqc5E6NwXf0DHO3qJxIOsfCUGnU1kKKQ8x2QKcovpLPrMUjlP4KqUHbb5k3DciVpMiXpNIyXAEqu4XZaXfVwo/hEeyqRYpDtRuvDJii/kE6dMyUkScZJgIPaIF6kUKXTMF6ybKIuB+mYSqN4kWkvjfILlaFKfnf8dzx95GnM7IR1WIVU/iOnEglw99GRCXDbFu22FckyJW15lhghG11XLdPEbXQNN1CXA5G0pSi/kFir5nAsnrOYM086kyPdR/jWrm/lvPzHWGVGAi9FAuxXyQ+RYqWkLc+yNUI2uoYbqMuBSNpmnurt3EsWL78wevH8wNAAR/qO8NIbL43YcABkNSEp2I0NE9QfC0KRX5FioTVtebb/WA+zq8tHHJvMCNnoGm5vr2k7M2uxiuSUn30Yl6zypvBgZPmFi79I+3PfGLFW7WD3QUIlIWJDseEdkOBtOMhmEjKZjQ2BKDmSpfpjgXgvIgGnkbY8y9YI2dIFtcNtqtTlQApOinVQeZGi/MLotWp9sT5wDPeWhNxsOEinzEiywIzMLVnlJWn9neCGvMdot3c8TYF5LyIBp5G2PMvmCJm6HEjBSl4HBW8/bt+Qv9G2ccovjN45GioJ0T/Yz2nVpw2fk4sNB5nutAxMyZFEApw8anrxFzP6cwzMexEJOCVteZYYIUt0OZg/q1J11WT6Of6KN8KWLCB9GJv6+lnTFaWldy/toRCn19RzhCghC+GcG7P8RzakU2YkWaBKjkyx/lig3otIgClp84FGyGTaC2ofxvi0bVO4iqbyd3sxHe6gdcmnaXnz+Zz3G60MVfLi8RfB4D0z35NyY0Mx1UArpvcikktK2kQk/1JsBPDVONO2TXs30/TxjTm7bXKPzsVzFtM90E3PQOrNSZmOzPkh3c0FhfBeRIJAGxFEJP+C2odxgvIVuTKZHp1Br4GWyeaCoL8XkaDQSJuI+COIfRh9mrad7JqunLX5yoJMNxcE+b2IBIVG2kREErJQvmIyirElVqYlTERkYkraREQSfJq2LcYencWYiIr4TdOjIiLJ8jRtO3qR/ooFK9h1ZFfOd6jmizYXiGSfOefyf1OzWcCDQAPQBnzcOXd8jPMGgefiT/c7566Y6NqNjY1ux44d2Qs2R7bt6xhVq61BZUBEponk3aLJCU2xLb4fTkyLJBEVyRFL+0SfkravAsecc18xs78FZjrn/maM87qcc9UnXmF8hZC0bdvXwR2bdlMZDo3oiqA2VCKFYyq9MldvXn1CXbKuaBe1kVrWX7Q+VyFPC+phKgUo7aTNrzVtVwL3xb+/D/gfPsWRsW37OrjxgZ2s+MavuPGBnWzb15HxNTZubaMyHKKmIkSJGTUVISrDITZubct2uCKSA1PtlalF+rmhHqZS7PxK2t7pnHsNIP44Z5zzKsxsh5ltMzPfE7vECNnRziizq8s52hnljk27M07c9h/roaq8dMSxqvJS9h9LXUxTRIJhMnXVkmmRfm5M9c9FJOhythHBzH4OnDzGS7dmcJn5zrmDZrYAeNLMnnPO/X6Me60CVgHMnz9/UvGmI3mEDBh+3Li1LaNpzfmzKjnaGR3+eYDu/kHmz6rMZrgikiNT7ZWpRfoTm8w0p3qYSrHL2Uibc+5S59x7x/h6CHjdzE4BiD8eHucaB+OP+4DNwOJxztvgnGt0zjXOnj07J+8HsjdCtnJZAz3RGJ19MYaco7MvRk80xsplDVmMVkSmqvVQK6s3r+baR65l9ebVw9NsUx0py2cHgPHeQ5BNdppTI5hS7PyaHn0Y+HT8+08DD40+wcxmmll5/Ps64IPA7rxFOIb5syrp7h8ccWwyI2RLF9TypRWLqKsJc6Srn7qasDYhiARMqsQhG3XVmk5uYv1F63lwxYOsv2h9zhK2QlzjNdlpzmKsdyeSzK+k7SvAh81sL/Dh+HPMrNHM/jV+zkJgh5k9A/wC+IpzztekLZsjZEsX1PKt685n080f4lvXna+ETSRgUiUOhdIrs1DXeE12o0ah/LmITJYvxXWdcx3Afxvj+A7gT+Pf/ydwTp5DSykxQjayvtqZWU+4VMNNJEvatsD2DV7D95mneu2o0iycO9H6qELolRmkNV73P38/9+6+l7f632JG+QyuX3Q9nzr7U2OeW19Tf0JJlHSnOQvhz0VkstQRIUNLF9TmNIFKruGWvENV06ciGWrbAo/fCuEqqJ7jNYJ//Na021JNJXEIiqC8h/ufv5+7n76bkIWIhCL0RHu4++m7AcZM3LRRQ2Rs6j0aMKrhJpIl2zd4CVt5DViJ9xiu8o6noRjWRwXlPdy7+15CFiIcCmNmhENhQhbi3t33jnm+pjlFxqaRtoDZf6yH2dXlI46phpvIJBx/xRthSxau8o6nIZE4FHIbpmy+h6l0Gnir/y0iociIY2WlZbzV/1bK2AvpsxbJByVtAaMabiJZMvNUb0q0vObtY9Fu73iaiiFxyMZ7SO6VmrwLNd3RrxnlM+iJ9hAOhYePDQwOMKN8xpTiEpluND0aMKrhJpIlS1Z5SVp/J7gh7zHa7R2XjEx1F+r1i64n5mJEY1Gcc0RjUWIuxvWLrs9t4CJFRklbwKiGm0iWNCz3Nh1U1UHXYe8xzU0IMtJUe6V+6uxP8bnFn6MyXElvrJfKcCWfW/y5cXePisjYzDnndwxZ1djY6Hbs2OF3GCIiRWP15tUn7ELtinZRG6ll/UXrfYxMpChYuidqTZuISJZNZdF+EKkEh0gwaKRNRIJpCoVxM5HtBCt50X5yglPoJSuGP6cC3UkrEmBpj7QpaROR4EkujBuu8jYQRLuzviYtFwmWphJFJENpJ23aiCAiwTPFwrjpykVvzqku2hcRGY/WtIlI8KRbGHeKU6i56M0ZlNZRIlJ8lLSJSPCkUxh3jN6irU/8FS3zzqR9qC+t9Wm5SLC0aF9EckXToyISPOkUxh01hdpaVsK6cD8dHS+OqNrfeqh13Nvkojen+maKSK5oI4KIBNNEU5/f/gNvhM28f3uujr1Kh4tRPRiDUz8ApLcBQLsiRcRnqtMmMt0UW20wGpanXp82agq13Q1Q5xyUvd2YPJ31acXQX1REpgdNj4oUgUTpio7ejrSnBgveqCnUemd0Dw3CSW+ve9MGgDxq2wI/WOmNgP5gpfdcRLJKSZtIEchF6YrAG9VbtLl8Hr3vmEtXKJy19WmSpsSmkO6jUD2H1u52Vj/5Oa794eWs3ry6uP/xIJJHmh4VKQK5KF1REJKmUJuANVqf5o+kTSGtQ92sK+kk4oy67uN0VHqjvtqMITJ1StpEioBqg3m0Ps0nSXX1WoaOEcGoLimDWN/w38mWPS36sxGZIk2PihSBXJSuEEnbzFO99YV4G0KqKIGhweFNIdNi1FckD5S0iRQB1QYTXyVtCqm3EN1uwEva4ptCpuOor0guaHpUpEhoanAKptgOa9pLbArZvoHmN46zrqIU3jGPqshMujXqK5I1Kq4rItNbcjuscJU3YhTt9pIQJW6TooLFIhlRcV0RkbQkt8OCtx+3b1DSNkka9RXJDa1pE5Hp7fgrXtKWLFzlHRcRCRAlbSIyvSXtfBwW7faOi4gEiKZHRWR6W7LKW9MGI9e0XfxFf+MSKUADAwO0t7fT19fndyiBU1FRQX19PWVlZZO+hpI2EZneknY+Du8evfiLWs8mMgnt7e3U1NTQ0NCAWdrr64uec46Ojg7a29s57bTTJn0dJW0iIkntsERk8vr6+pSwjcHMqK2t5ciRI1O6jta0iYiISNYoYRtbNj4XJW0iIiJSNKqrqyc+KW7t2rWsW7cuK9d/7LHHeM973sPpp5/OV77ylYyumS5NjwbYtn0dbNzaxv5jPcyfVcnKZQ0sXVDrd1giIiKSZHBwkD//8z/nZz/7GfX19TQ1NXHFFVewaNGirN5HI20BtW1fB3ds2s3Rziizq8s52hnljk272bavw+/QREREsmLbvg5ufGAnK77xK258YGfOfsc98sgjXHDBBSxevJhLL72U119/ffi1Z555hksuuYQzzjiD73znO8PH77rrLpqamjj33HO57bbbUl5/+/btnH766SxYsIBwOExzczMPPfRQ1t+HkraA2ri1jcpwiJqKECVm1FSEqAyH2Li1zefIREREpi6fgxPLly9n27ZtPP300zQ3N/PVr351+LVnn32WRx99lK1bt3L77bdz8OBBnnjiCfbu3cv27dvZtWsXO3fu5Kmnnhr3+gcOHOBd73rX8PP6+noOHDiQ9feh6dGA2n+sh9nV5SOOVZWXsv9Yj08RiYiIZE/y4AQw/Lhxa1vWlwK1t7dz7bXX8tprrxGNRkeU3bjyyiuJRCJEIhEuvvhitm/fzpYtW3jiiSdYvHgxAF1dXezdu5cLL7xwzOuP1cc9FxsyNNIWUPNnVdLdPzjiWHf/IPNnVfoUkYiISPbsP9ZDVXnpiGO5Gpy4+eabuemmm3juuef49re/PaL47+jkysxwznHLLbewa9cudu3axUsvvcQNN9ww7vXr6+t59dVXh5+3t7czd+7crL8PX5I2M7vGzJ43syEza0xx3kfN7EUze8nM/jafMfpt5bIGeqIxOvtiDDlHZ1+MnmiMlcsa/A5NRERkyvI5OPHmm28yb948AO67774Rrz300EP09fXR0dHB5s2baWpq4rLLLuOee+6hq6sL8KY/Dx8+PO71m5qa2Lt3Ly+//DLRaJSWlhauuOKKrL8Pv6ZHfwtcBXx7vBPMrBT4Z+DDQDvQamYPO+d25ydEfy1dUMuXViwatXv0TO0elaxpPdRKy54W2jvbqa+pp/msZppObvI7LBGZJlYua+COTd6v9KryUrr7B+ODE2dO6bo9PT3U19cPP//CF77A2rVrueaaa5g3bx5Lly7l5ZdfHn59yZIlXH755ezfv58vf/nLzJ07l7lz5/LCCy+wbNkywCvz8b3vfY85c+aMec9QKMQ3v/lNLrvsMgYHB/nMZz7D2WefPaX3MRYbax42X8xsM7DGObdjjNeWAWudc5fFn98C4Jz7+1TXbGxsdDt2nHA5EUnSeqiVda3riIQiVJVV0d15kN7OA6zpK6XppDO8fpzqECAiGXrhhRdYuHBh2udPt9JW43w+aS9+C/JGhHnAq0nP24ELfIpFpKi07GkhEopQHa6GnmNUH38FSoyWilKauo96DdQvu1OJm4jk1NIFtUWdpGVbzta0mdnPzey3Y3xdme4lxjg25rCgma0ysx1mtmOqfb1EpoP2znaqyqq8J2+8AiWlVFkZ7S4G5TUQrvIaqIuISGDkbKTNOXfpFC/RDrwr6Xk9cHCce20ANoA3PTrF+4oUvfqaejp6O7yRtoFeCIXpdkPUW5l3QrgKjr/ib5AiIjJCkEt+tAJnmNlpZhYGmoGHfY5JpCg0n9VMb6yXrmgXLlRB19AAvTiaS2Z5J0S7Yeap/gYpIiIj+FXy44/NrB1YBjxqZo/Hj881s58AOOdiwE3A48ALwA+cc8/7Ea9IsWk6uYk1TWuojdRytGomtUOONUM1NFkE+ju9pG3JKr/DFBGRJL5sRHDO/Rj48RjHDwJ/lPT8J8BP8hiayLTRdHLT2yU+2rZ4a9iOv+KNsF38RW1CEBEJmCDvHhWRfGlYriRNRIpCdXX1cFHciaxdu5bq6mrWrFkz5et/5jOfYdOmTcyZM4ff/va3aV8vE0Fe0yYiIiJSEK6//noee+yxnN5DSZuIiIj4o20L/GAlfPsPvMe2LTm5zSOPPMIFF1zA4sWLufTSS3n99deHX3vmmWe45JJLOOOMM/jOd74zfPyuu+6iqamJc889l9tuu23Ce1x44YXMmjUrJ/EnKGkTERGR/Gvb4hXy7j4K1XO8x8dvzUnitnz5crZt28bTTz9Nc3MzX/3qV4dfe/bZZ3n00UfZunUrt99+OwcPHuSJJ55g7969bN++nV27drFz506eeuqprMeVKa1pExERkfzbvsGrCVle4z1PPG7fkPU1tu3t7Vx77bW89tprRKNRTjvttOHXrrzySiKRCJFIhIsvvpjt27ezZcsWnnjiCRYvXgxAV1cXe/fu5cILL8xqXJlS0iYiIiL5d/wVb4QtWY4Ke99888184Qtf4IorrmDz5s2sXbt2+DWzkQ2YzAznHLfccguf/exnsx7LVGh6VERERPJv5qleTchkOSrs/eabbzJv3jwA7rvvvhGvPfTQQ/T19dHR0cHmzZtpamrisssu45577hneJXrgwAEOHz6c9bgypaRNRERE8m/JKi9J6+8EN5S1wt49PT3U19cPf33ta19j7dq1XHPNNXzoQx+irq5uZBhLlnD55ZezdOlSvvzlLzN37lw+8pGP8MlPfpJly5ZxzjnncPXVV9PZ2Znyvp/4xCdYtmwZL774IvX19Xz3u9+d0vsYizlXXK06Gxsb3Y4dO/wOQ0REZNp54YUXWLhwYfo/MLqw95JVRV0zcpzPx8Y6dyxa0yYikmOth1pp2dNCe2c79TX1NJ/V/HY3CpHpTIW9M6LpURGRHGo91Mq61nV09HZQF6mjo7eDda3raD3U6ndoIlJglLSJiORQy54WIqEI1eFqzIzqcDWRUISWPS1+hyYiBUZJm4hIDrV3tlNVVjXiWFVZFe1d7T5FJCKFSkmbiEgO1dfU0z0wsqxB90A39dX1PkUkIoVKSZuISA41n9VMb6yXrmgXzjm6ol30xnppPqvZ79BEpMAoaRMRyaGmk5tY07SG2kgtR/uOUhupZU3TGu0eFcmR6urqtM9du3Yt69atm/L1X331VS6++GIWLlzI2WefzT/90z9ldM10qeSHiEzeNKuxNFlNJzcpSRMpYqFQiPXr1/P+97+fzs5Ozj//fD784Q+zaNGirN5HI20iMjltW+DxW6H7qNc/sPuo97xti9+RiUiBaD3UyurNq7n2kWtZvXl1zkrhPPLII1xwwQUsXryYSy+9lNdff334tWeeeYZLLrmEM844g+985zvDx++66y6ampo499xzue2221Je/5RTTuH9738/ADU1NSxcuJADBw5k/X0oaRORydm+wWvuXF4DVuI9hqu84yIiE8hnDcPly5ezbds2nn76aZqbm/nqV786/Nqzzz7Lo48+ytatW7n99ts5ePAgTzzxBHv37mX79u3s2rWLnTt38tRTT6V1r7a2Np5++mkuuOCCrL8PTY+KyOQcf8UbYUsWrvKOi4hMILmGITD82LKnJevLCdrb27n22mt57bXXiEajnHbaacOvXXnllUQiESKRCBdffDHbt29ny5YtPPHEEyxevBiArq4u9u7dy4UXXpjyPl1dXXzsYx/j61//OjNmzMjqewCNtInIZM081WvunCza7R0XEZlAPmsY3nzzzdx0000899xzfPvb36avr2/4NbORrT/NDOcct9xyC7t27WLXrl289NJL3HDDDSnvMTAwwMc+9jGuu+46rrrqqqy/B1DSJiKTtWSVl6T1d4Ib8h6j3d5xEZEJ5LOG4Ztvvsm8efMAuO+++0a89tBDD9HX10dHRwebN2+mqamJyy67jHvuuYeuri4ADhw4wOHDh8e9vnOOG264gYULF/KFL3wh6/EnKGkTkclpWA6X3QlVddB12Hu87E7tHhWRtOSqhmFPTw/19fXDX1/72tdYu3Yt11xzDR/60Ieoq6sbcf6SJUu4/PLLWbp0KV/+8peZO3cuH/nIR/jkJz/JsmXLOOecc7j66qvp7Owc956//vWvuf/++3nyySc577zzOO+88/jJT34ypfcxFnPOZf2ifmpsbHQ7duzwOwwREZFp54UXXmDhwoVpn996qJWWPS20d7VTX11P81nNRV0eZ5zPx8Y6dyzaiCAiIiK+UA3DzGh6VERERKQAKGkTERERKQBK2kRERCRrim2tfLZk43NR0iYiIiJZUVFRQUdHhxK3UZxzdHR0UFFRMaXraCOCiIiIZEV9fT3t7e0cOXLE71ACp6Kigvr6qdWgU9ImIiIiWVFWVjaiRZRkl6ZHRURERPiUJ8UAAAafSURBVAqAkjYRERGRAqCkTURERKQAFF0bKzM7ArzidxxpqAOO+h1EgOnzmZg+o9T0+aSmzyc1fT6p6fNJLZPP56hz7qPpnFh0SVuhMLMdzrlGv+MIKn0+E9NnlJo+n9T0+aSmzyc1fT6p5erz0fSoiIiISAFQ0iYiIiJSAJS0+WeD3wEEnD6fiekzSk2fT2r6fFLT55OaPp/UcvL5aE2biIiISAHQSJuIiIhIAVDS5iMz+7/N7Fkz22VmT5jZXL9jChIzu8vM9sQ/ox+b2Ul+xxQkZnaNmT1vZkNmpl1ccWb2UTN70cxeMrO/9TueoDGze8zssJn91u9YgsjM3mVmvzCzF+L/ff2F3zEFiZlVmNl2M3sm/vn8nd8xBZGZlZrZ02a2KZvXVdLmr7ucc+c6584DNgH/y++AAuZnwHudc+cCvwNu8TmeoPktcBXwlN+BBIWZlQL/DPwhsAj4hJkt8jeqwLkXSKsm1DQVA1Y75xYCS4E/19+hEfqBS5xz7wPOAz5qZkt9jimI/gJ4IdsXVdLmI+fcW0lPqwAtMEzinHvCOReLP90G1PsZT9A4515wzr3odxwBswR4yTm3zzkXBVqAK32OKVCcc08Bx/yOI6icc685534T/74T7xfvPH+jCg7n6Yo/LYt/6XdXEjOrBy4H/jXb11bS5jMzu9PMXgWuQyNtqXwG+KnfQUjgzQNeTXrejn7hyiSZWQOwGPgvfyMJlvjU3y7gMPAz55w+n5G+Dvw1MJTtCytpyzEz+7mZ/XaMrysBnHO3OufeBTwA3ORvtPk30ecTP+dWvCmLB/yL1B/pfD4ygo1xTKMAkjEzqwZ+CHx+1KzItOecG4wv66kHlpjZe/2OKSjMbAVw2Dm3MxfXD+XiovI259ylaZ76b8CjwG05DCdwJvp8zOzTwArgv7lpWJ8mg78/4mkH3pX0vB446FMsUqDMrAwvYXvAOfcjv+MJKufcG2a2GW+NpDa2eD4IXGFmfwRUADPM7HvOuT/JxsU10uYjMzsj6ekVwB6/YgkiM/so8DfAFc65Hr/jkYLQCpxhZqeZWRhoBh72OSYpIGZmwHeBF5xzX/M7nqAxs9mJnfxmFgEuRb+7hjnnbnHO1TvnGvD+//NkthI2UNLmt6/Ep7qeBT6Ct9tE3vZNoAb4Wbwsyr/4HVCQmNkfm1k7sAx41Mwe9zsmv8U3rtwEPI63gPwHzrnn/Y0qWMzs+8BW4D1m1m5mN/gdU8B8EPgUcEn8/zu74qMm4jkF+EX891Yr3pq2rJa1kPGpI4KIiIhIAdBIm4iIiEgBUNImIiIiUgCUtImIiIgUACVtIiIiIgVASZuIiIhIAVDSJiIFxcwG42UYnjezZ8zsC2ZWEn+t0czu9imu/8zSda6Jv7chM2vMxjVFpDio5IeIFBQz63LOVce/n4PXTeTXzrmi6CZiZgvxehZ+G1jjnNvhc0giEhAaaRORguWcOwysAm4yz0VmtgnAzNaa2X1m9oSZtZnZVWb2VTN7zswei7cqwszON7NfmtlOM3vczE6JH99sZv9gZtvN7Hdm9qH48bPjx3aZ2bOJziZm1hV/NDO7K144+zkzuzZ+/KL4Nf/dzPaY2QPx6vuj39MLzrkX8/H5iUhhUdImIgXNObcP7/9lc8Z4+d3A5cCVwPeAXzjnzgF6gcvjids3gKudc+cD9wB3Jv18yDm3BPg8b/cF/jPgn+INsxvx+p0muwo4D3gfXoufuxKJILA4fq1FwAK86vsiImlRw3gRKQYnjFjF/dQ5N2BmzwGlwGPx488BDcB7gPfitUojfs5rST+faBa+M34+eC2gbjWzeuBHzrm9o+65HPi+c24QeN3Mfgk0AW8B251z7QBmtit+zS2ZvlkRmZ400iYiBc3MFgCDwOExXu4HcM4NAQPu7UW8Q3j/aDXgeefcefGvc5xzHxn98/Hrh+LX+jfgCrzRusfN7JLRIaUItz/p++FrioikQ0mbiBQsM5sN/AvwTTe5XVUvArPNbFn8emVmdvYE91wA7HPO3Q08DJw76pSngGvNrDQe34XA9knEJiIygpI2ESk0kUTJD+DnwBPA303mQs65KHA18A9m9gywC/jABD92LfDb+PTmWcDGUa//GHgWeAZ4Evhr59yhdGMysz82s3ZgGfComT2e7s+KSHFTyQ8RERGRAqCRNhEREZECoKRNREREpAAoaRMREREpAEraRERERAqAkjYRERGRAqCkTURERKQAKGkTERERKQBK2kREREQKwP8P7mhkYTKUSAkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#self implementation viz\n",
    "fig=plt.figure(figsize=(10,5))\n",
    "ax=fig.add_subplot(111)\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "for i in set(df['gmm']):\n",
    "    plt.scatter(df['dim1'][df['gmm']==i],\n",
    "                df['dim2'][df['gmm']==i],\n",
    "                label=f'Label {int(i)}',alpha=0.7)\n",
    "\n",
    "plt.title('GMM by us')\n",
    "plt.xlabel('Dimension 1')\n",
    "plt.ylabel('Dimension 2')\n",
    "plt.legend(loc=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAFNCAYAAABST1gVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xt8XPV95//XR5fRHbAtE7AVI5xwMwnBwTJW4hCgJE4DC82FWAkbl0324WQpSVvMtoWQhfKDX7PBZNMkzRZnw4IJrcKvDQVMCqSlDnFr17KDgWBMTI1tZGNsy8bW/fr5/XFG8khIoxlpZs5c3s/HQ4+jOefMOd85AubD9/L5mLsjIiIiItmtKOwGiIiIiMjkFLSJiIiI5AAFbSIiIiI5QEGbiIiISA5Q0CYiIiKSAxS0iYiIiOQABW0iUlDM7FIzaw3rHpm4v4jkJwVtIpJRZtZkZv9uZp1mdjD6+w1mZtHjD5iZm9nVY9733ej+66Ovr4++/s6Y834vuv+BTH0mEZFMUNAmIhljZquAvwTuAU4D3gV8FfgwEIk59bfA78e8rwS4FviPMZf8D2B59PiwFdH3F5Qxz0BE8pCCNhHJCDM7GbgTuMHd/87d2z3wvLtf5+69Mac/AXzYzGZEX38CeBE4MOayB4CXgGXRe8wEPgQ8nkB7bjWzw2a228yui+5rMLO3YgMgM/uMmW2b4BqfNLPtZtZuZvvM7OYJzvt69Ly6cY7NMbO/N7NDZva6mX095thiM9toZm+b2Ztm9gMzi8QcdzP7AzPbCeyM2fdVM9tpZkfN7K+GezFFJLcpaBORTGkEyoDHEji3hyDwaoq+XgGsneDctdHjRM9/DOid4NxhpwG1wFyCHr01ZnaOu7cAbcDHYs79z8BDE1znx8BX3L0GeB/w7NgTzOybwPXAR929dcyxIoIA9YVoW34H+CMzWxY9ZRD442hbG6PHbxhzi98DLgYWxOy7CmgAPgB8jmhQKyK5TUGbiGRKLXDY3QeGd5jZv0V7kbrN7JIx568FVkR76D4K/MME130UuDR6Xrzgbqxvunuvu/8SeJIguAF4kCBQG+65Wwb8zQTX6AcWmNlJ7n7U3X8dc8yi8+2WAZe5+6Fx3t8AzHb3O929z913AT8iGqy6+1Z33+TuA+6+G7iP4FnE+gt3P+Lu3TH7vuXub7v7XuBfgAsTeSAikt0UtIlIprQBtbFDj+7+IXc/JXps1H+P3H0DMBu4DVg3JiiJPa+bIOi6Dah1939NoC1H3b0z5vUeYE70958A/8nMqgkCuV+5+5sTXOczwCeBPWb2SzNrjDl2CrCSIKg6NsH7zwDmRAPXt83sbeBWgrl+mNnZZrbOzA6Y2XHg/yUIfmO9Mc51Y4eRu4DqCe4vIjlEQZuIZMpGgmHLa5J4z0+AVUzee7Y2et5Ew5hjzTCzqpjX84D9AO6+L9rWTwFfjHdNd29x92uAUwl6Ah+JOXyUYJjy/5rZhye4xBvA6+5+SsxPjbt/Mnr8fwM7gLPc/SSCgG7s/DSf/OOKSD5Q0CYiGeHubwN/DvzQzD5rZtVmVmRmFwJVE7ztewTzy56b5PK/jJ73/SSa9OdmFjGzjxAEV/9fzLG1wJ8A7ycYfn2H6HuvM7OT3b0fOE4wB22Eu68HrgMeNbOLx7nMZuC4mf2pmVWYWbGZvc/MGqLHa6LX7TCzc4H/lsTnE5E8o6BNRDLG3b8N3EQQEB0E3iKYp/WnwL+Nc/4Rd/9nd4/bmxRdhfrP7n4kwaYcIOgJ2w88DHzV3XfEHH+UYOjy0THDqGN9EdgdHbr8KtG5cGPa9gvgvwCPm9lFY44NAv+JYM7Z68Bh4P8AJ0dPuRn4AtBOMNftpwl+PhHJQzbJfwtFRAqSmf0HwcrQfwq7LSIioJ42EZF3MLPPEMwVe0cKDxGRsCiDtohIDDNbT5Dz7IvuPhRyc0RERmh4VERERCQHaHhUREREJAcoaBMRERHJAaHOaTOz+wnyIx109/eNc/xSgjqCr0d3/czd74x3zU984hP+1FNPpbqpIiIiIukwNmH2hMJeiPAA8APiZzv/lbtflegFDx8+PN02iYiIiGSdUIdH3f05INFkmCIiIiIFKxfmtDWa2Qtm9o9mdv54J5jZSjPbYmZbDh06lOn2iYiIiKRdtgdtvwbOcPcPENQU/IfxTnL3Ne6+yN0XzZ49O6MNFBEREcmErA7a3P24u3dEf/85UGpmtSE3S0RERCTjsjpoM7PTzMyivy8maG9buK0SERERybywU378LXApUGtmrcDtQCmAu/818Fngv5nZANANNLlKOIiIiEgBCjVoc/fPT3L8BwQpQUREREQKWlYPj4qIiIhIQEGbiIiISA5Q0CYiIiKSA8IuYyUiIiKFZPcG2LwGju6BGWfA4pVQvzTsVuUE9bSJiIhIZuzeAE9/AzoPQ/WpwfbpbwT7ZVIK2kRERCQzNq+BSBWU1YAVBdtIVbBfJqXhUZnUpl1trN24m71Hupg3s5IVjfUsmT8r7GaJiEiuObon6GGLFakK9suk1NMmcW3a1cZd67ZzuL2P2dVlHG7v465129m0S4UpREQkSTPOgL7O0fv6OoP9MikFbRLX2o27qYyUUFNeQpEZNeUlVEZKWLtxd8gtExGRnLN4ZRCk9baDDwXbvs5gv0xKQZvEtfdIF1VlxaP2VZUVs/dIV0gtEhGRnFW/FJbdDVW10HEw2C67W6tHE6Q5bRLXvJmVHG7vo6b8xD8qnb2DzJtZGWKrREQkZ9UvVZA2Reppk7hWNNbT1TdAe88AQ+609wzQ1TfAisb6sJsmIiJSUBS0SVxL5s/itqsWUFsT4VBHL7U1EW67aoFWj4qIiGSYhkdlUkvmz1KQJiIiEjL1tImIiIjkAAVtIiIiIjlAQZuIiIhIDlDQJiIiIpIDFLSJiIiI5AAFbSIiIiI5QEGbiIiISA5QnrY8tWlXG2s37mbvkS7mzaxkRWO9cq2JiEj+2L0BNq+Bo3tgxhlB0fk8L4+lnrY8tGlXG3et287h9j5mV5dxuL2Pu9ZtZ9OutrCbJiIiMn27N8DT34DOw1B9arB9+hvB/jymoC0Prd24m8pICTXlJRSZUVNeQmWkhLUbd4fcMhERkRTYvAYiVVBWA1YUbCNVwf48pqAtD+090kVVWfGofVVlxew90hVSi0RERFLo6J4gSIsVqQr25zEFbXlo3sxKOnsHR+3r7B1k3szKkFokIiKSQjPOgL7O0fv6OoP9eUxBWx5a0VhPV98A7T0DDLnT3jNAV98AKxrrw26aiIjI9C1eGQRpve3gQ8G2rzPYn8cUtOWhJfNncdtVC6itiXCoo5famgi3XbVAq0dFRCQ/1C+FZXdDVS10HAy2y+7O+9Wj5u5htyGlFi1a5Fu2bAm7GSIiIiKJsERPVJ62AqCcbSIiIrlPw6N5TjnbRERE8oOCtjynnG0iIiL5QUFbnlPONhERkfygoC3PKWebiIhIflDQlueUs01ERCQ/KGjLc8rZJiIikh9CTflhZvcDVwEH3f194xw34C+BTwJdwPXu/uvMtjL3LZk/S0GaiIhIjgu7p+0B4BNxjv8ucFb0ZyXwvzPQJhEREZGsE2rQ5u7PAUfinHINsNYDm4BTzOz0zLROREREJHuE3dM2mbnAGzGvW6P7RERERApKtgdt49XjekexVDNbaWZbzGzLoUOHMtAsERERkczK9qCtFXh3zOs6YP/Yk9x9jbsvcvdFs2fPzljjRERERDIl24O2x4EVFlgCHHP3N8NulIiIiEimhZ3y42+BS4FaM2sFbgdKAdz9r4GfE6T7eI0g5cd/CaelIiIiIuEKNWhz989PctyBP8hQc0RERESyVrYPj4qIiIgICtpEREREcoKCNhEREZEcoKBNREREJAeEuhBB8temXW2s3bibvUe6mDezkhWN9SpaLyIiMg3qaZOU27SrjbvWbedwex+zq8s43N7HXeu2s2lXW9hNExERyVkK2iTl1m7cTWWkhJryEorMqCkvoTJSwtqNu0NumYiISO5S0CYpt/dIF1VlxaP2VZUVs/dIV0gtEhERyX0K2iTl5s2spLN3cNS+zt5B5s2sDKlFIiIiuU8LEfJUmAsBVjTWc9e67UDQw9bZO0hX3wArGs/OyP1FRETykQWVovLHokWLfMuWLWE3I1TDCwEqIyX0Dw6yp62Lrv5BPlB3Mn/8sXMyErxp9aiIiEhCLOETFbTlnxse3srh9j4GhobY+VYHxUWG4xSZMfeUCm67aoECKBERkeyQcNCm4dE8tPdIF7Ory3h5fyfFRUZJkeEOfYM+sopTQZuIiEyk5UALzTuaaW1vpa6mjqZzm2g4rSHsZhU8LUTIQ8MLAbr7BykuCgL4QXcqSou1ilNEROJqOdDC6pbVtHW3UVtRS1t3G6tbVtNyoCXsphU8BW15aEVjPV19A5QWGQODQwwMDTE4BHNnlGsVp4iIxNW8o5mKkgqqI9WYGdWRaipKKmje0Rx20wqegrY8tGT+LG67agFnvauG3sEhzIz3zK6itLg4uoqzPuwmiohIlmptb6WqtGrUvqrSKlo7WkNqkQzTnLY8tWT+LH76lcZRqzjn1ERY0Xh2SuezaZWoiEh+qaupo627jepI9ci+zv5O6qrrQmyVgFaPyjTEphaJzcem1akiIrlreE5bRUkFVaVVdPZ30j3Qzc0NN2sxQnokvHpUw6MyZaoxKiKSfxpOa+DmhpuZVTGLwz2HmVUxSwFbltDwqEzZcGqRWFqdKiKS+xpOa1CQloXU0yZTphqjIiIimaOeNpky1RgVEXZvgM1r4OgemHEGLF4J9UvDbpVkmJLxZoZ62mTKhlOL1NZEONTRS21NRIsQRArJ7g3w9Deg8zBUnxpsn/5GsF8KhpLxZo562mRalsyfpSBNpFBtXgORKiirCV4PbzevUW9bAYlNxguMbJt3NKu3LcXU0yYiIlNzdE8QtMWKVAX7pWAoGW/mqKetwCk5rohM2YwzgiHR4R42gL7OYL8UDCXjzRz1tBWw4eS4h9v7mF1dxuH2Pu5at51Nu9rCbpqI5ILFK4MgrbcdfCjY9nUG+6VgNJ3bRPdANx19Hbg7HX0ddA9003RuU9hNyzsK2gqYkuOKyLTUL4Vld0NVLXQcDLbL7tZ8tgKjZLyZo+HRAqbkuCIybfVLFaSJkvFmiIK2PDbZfLV5Mys53N5HTfmJfwyUHFdERCQ7aXg0TyUyX21FYz1dfQO09www5E57z0A0OW59eA0XERGRcamnLU/FzlcDRrbD89WGe+CqykoA51BHb7Q37mytHhUREclCCtry1ETz1bbvP8Zd67ZTGSlhdnVZtPTUoCoZiIjkIpURKygaHs1TExVz7+kf0opREZF8UChlxHZvgEdWwH0fDbb59vmSoKAtT000X60iUkxVWfGoc7ViVEQkB8WWEbOiYBupCvbni0IJTBOkoC1PTVTM/bzTTxq3B04rRkVEckwhlBErhMA0CZrTlscmKuZ+17rtQNDDFsxpG2BF49mZbp6IiExHIZQRO7on6GGLlW+BaRJC7Wkzs0+Y2atm9pqZ/dk4x683s0Nmti3681/DaGc+2LSrjRse3spdT26nqqwYM0b1wGkRgohIjimEMmIzzgg+U6x8C0yTEFrQZmbFwF8BvwssAD5vZgvGOfWn7n5h9Of/ZLSReWJszjZ3461j3VRFSth7pIu1G3er3qiIZI4mlqdGIZQRK4TANAlh9rQtBl5z913u3gc0A9eE2J68NbbGaP/gIAc7+tj5VrsKxYtIZmlieWrVL4XPrYWv/DLY5lPABoURmCYhzDltc4E3Yl63AhePc95nzOwS4LfAH7v7G+OcI3GMzdm272gPpcVG/5CPpP2AILjTMKmIpFXsxHI4sd28pmC/iGUSqm87IsygzcbZ52NePwH8rbv3mtlXgQeBy99xIbOVwEqAefPmpbqdOW9sjdHu/kGKDCpKT6T+UNoPkTyUjYlXNbF8tGz8G0nWCnN4tBV4d8zrOmB/7Anu3ubuvdGXPwIuGu9C7r7G3Re5+6LZs2enpbG5bGzOtpJio3/QqZtRMXKO0n6I5JlsHYbUxPITsvVvlEsKbH5kmEFbC3CWmZ1pZhGgCXg89gQzOz3m5dXAKxlsX94Ym7Pt7HdV866TyygpKlKheJF8la35rTSx/IRs/RvligIMekMbHnX3ATO7EXgaKAbud/eXzexOYIu7Pw583cyuBgaAI8D1YbU318XmbNu0q43/9YtXeeXN4wCcd/pJKU37sWlX20hB+qAIfb3myolkWrYOQw5PLI8dErzs1sIcEszWv1GuKMD5kaEm13X3nwM/H7Pvf8T8fgtwS6bblc+G039URkq46IwZdPYO0tk7kJbrx65MVS44kQzL5sSrmlgeyOa/US4owKBXFREKTGz6DyDuytGp9Jglc30RSaPFK4OhIgi+yPo6g5/Lbp3edfNk4nzLgRaadzTT2t5KXU0dTec20XBaQ2Ybka6/UaEowKBXtUcLzN4jXQkVjB+bkDfRXG6JXl9E0iwd+a3yZA5Ry4EWVrespq27jdqKWtq621jdspqWAy2ZbYhykE1PAc6PVE9bgRmb/gPGXzk6WY/ZRL1wiV5fRDIg1cOQeTKHqHlHMxUlFVRHqgFGts07mjPf26ah4qkrwPmRCtoKzIrG+oQKxo9NyDt8/t4jXXHnrSV6fRHJQXkyh6i1vZXaitpR+6pKq2jtaA2pRYHQhmxzeci7wIJeDY8WoKqyEnYcOM7WPUcx83EXCcybWUln7+CofcM9ZmPLYtWUl1AZKRnphYtNL6KC9CJ5JE9yrNXV1NHZP/pzdPZ3UlddN+l7Ww60sGr9KpY/sZxV61elbEg1tCHbPBnyLhQK2grIcA+ZO1x0xgzOPe2kdwRmw8Ym5I3N5TbZvLUl82fxw+suYt3XPsIPr7tIAZtIvsiTOURN5zbRPdBNR18H7k5HXwfdA900ndsU933pDKxih2zNjOpINRUlFTTvaJ72teOKlyuuwBLX5oIJgzYzO8nM/sLMHjKzL4w59sP0N01SLV4P2Vjxeszi9cKJSB7Lk4nzDac1cHPDzcyqmMXhnsPMqpjFzQ03TzoUmc7AqrW9larSqlH7Eh2ynVbv39E9QZAWK1IFb72sHrgsFG9O2/8FdgJ/D3zJzD4DfCFaVmpJJhonqRVvntp4YhPyxtK8NZEClidziBpOa0h6vthU5sIlOk+trqaOtu62kUURkNiQ7XDvX0VJxajev0SCUGDitBl9nVD9rpxfdJJv4g2Pvsfd/8zd/8HdrwZ+DTxrZhrryiKbdrVxw8Nbuer7v+KGh7fGTcmRqh4yzVsTyXIa1kqLZOfCJTOcOtUh22n3/k005F1aOX4PXI4tOsk38XraysysyN2HANz9bjNrBZ4DquO8TzJk7CrOXYc6Wbl2CzOrSikvLeat47109w9ySkUpKy+Zn9Iesol64UQkZMMTyyNVtFRW0dy1g9Z/uZG6Uz9A00U3Zj6lRR5pOreJ1S2rgaCHrbO/M25glUxqkeEh2+YdzbR2tFJXndjq0WmvhJ0obcbmNQWXuDYXxOtpewK4PHaHuz8IrAL60tkoSUzsHLVj3f3sO9qNOxw83sOrb3Xwdlc/pUVGZ+8g3376VbbvP6YeMpF8F51Y3lJaxOqhQ7QZ1FoJbW2vhpNANo8kOxcumXlqo4ZREwzYYHorYUfUL4XPrYWv/DLY1i/Nm0Un+WbCnjZ3/5MJ9j8FnJW2FknCYueotR7tprjIKDY40usUGViR0TvgnFReDAOw5rldbLr1CgVpIvksmkuteXAfFRjVVgxWRPVAH0SHzdTbNnXJzIVLdJ7adOalJdv7l7ACTFybC5RcN4fFVh/o7h8kUlzE4NAQDhhgBoPuAJQWG29394faXhHJgOjE8taifmqJpuYZGoTSipRNmpfEJBpQJVuhYezf6ar5V7Ht0LakhlUTkieLTvKJgrYcMrZ01KIzZvCzX+8DoLykiJ7+IQBKi4whd3AoNgOgf9A5paI0rKaLSKZEi5DXlRltNkT1kAdB26wzEpo0P+VViCmST4FjovPUkpmXNt7fad2udRn/O0k4FLTliNhFByVmbPyPNp7dcZAza6swc6rLS+kd6GXOyeX0Dw6xq60LhpyKsiJ6B4boHxpi5SUa1RbJe9FhraaN97K653UoLqdq5nvpLImkbNJ8umRL4JhKiQynJpPuIxv+TgUhS0t7JVQRwcw+ZGZfMLMVwz/pbpiMNrzooH9wkP841Im7U1ZcxFvHeujsHeQvPv1+1qxYxPxTq6koK+Gcd1VzSmUp/UNOVVkxf7LsHL60dH7YH0NEMqF+KQ2ff5Sbf/dHzJr3YQ4XDaV00nxc00g3ElpVgJAlk+4jZX8nmVgWl/aatKfNzB4C3gNsA4aTfDmwNo3tkjGGFx28vK+T4iIoKSrCgb7BoZGqBioZJSKx0jFpPq6YdCOjvuwSrJqQyDBhPg2fDksm3UdK/k4SX2xpL8iqxMKJDI8uAha4R2e0SyiGFx0ECw6CeWqDQ05FaXHcqgbDxs6HW9FYrwBPpNDFDAE1nTST1ZEgm9OUVyFO88tusoAkH4dPhyUaYKdttaicEF2BPUqWJBZOZHj0N8Bp6W6IxDdcwL2k2BhwZ2DIGRxy6mZUTFrVYHg+3OH2PmZXl3G4vY+71m2PWz1BRPLcmCGghp5ebm47wqwhT6oe5ygT1bFM8MtusmHCbBw+nVbdzymYat1UScKMM4KcdLGyJLFwIj1ttcB2M9sM9A7vjJa2kgwZLh31v37xKi+0HqOytJgzT62ipKho0qoGsUl4gZHt2o271dsmUqjG6RVrABo6HT7306ldc6I6lgl+2U02TDjt7P8pFlbP31TqpkoSoiuwgeDfkeFarJfdGm67SCxouyPdjZDELJk/i59+5UOjhjpPPyXCisaz4wZfyRaKF5HsN+25XekYAkrBl128gCTb5nNpJWeeyuLEwpMGbe7+SzN7FzD8T+Bmdz+Y3mZJPMnW/YxNwjsskULxmgcnkp1S0sMzzV6xcaX5yy4d87mmE/wm2/OXj4so8laWJhaedE6bmX0O2AxcC3wO+Hcz+2y6GyapMzwfrr1ngCF32nsGokOq9RO+R/PgRLJXSuZ2pau25Hh1LFMk1fO5hoPftu62UcFvovPSkqn7Od17iUBiw6PfABqGe9fMbDbwT8DfpbNhMn2xPWVVZUE5m0MdvdFes/hDqpoHJ5K9UjK3K4uHgOJJ5Xyu6Q5vJtPzp6FUSYVEgraiMcOhbSSYlFfCE1tBYXZ1GZ29g3T1DXDbVQsSCro0D04kC0yQlT1lc7uydAgoU6Yb/CaTXy3bFlFIbkokaHvKzJ4G/jb6ejnw8/Q1SVJhuj1lU50HJyIpEidRrXJ1pcZ0gt+x89NuXhR/mDbbFlFIbpq0x8zd/zuwBrgA+ACwxt3/NN0Ny2abdrVxw8Nbuer7v+KGh7dm5Tyv2CHRYZP1lMV+rraOXg519CQ1D05EUig2JYcVBdtIFWxeo1xdKZJM+ahYU5mfNtV7icRKqGC8u/898PdpbktOGDvsODxBP9Fhx0xJtqdsvOFUH3LMEp8HJyIpNElKDuXqmrqHXn6IB7Y/wPHe41SUVDC7Yjbdg91xhzdjTWV+WjJDqSITmTBoM7MN7r7UzNoJao2OHALc3U9Ke+uyULZO0B+bnmPRGTP42a/3AUEP2/CctomS8I7/uSqYVR3hp9c1ZuZDiMgJ6UjJITz08kN87/nvUWIlVJRU0D/Yz972vXx94df54vlfTOgaU52fpkBbpmvC4VF3Xxrd1rj7STE/NYUasMHUhh3Tbbz0HD/79T4+/cG51NZEONTRS21NJG5v4HSHU7N1mFgkZ6UrJUeBe2D7A5RYCZGSCGZGpCRCiZXwwPYHEr5GMqk+RFIpkTxt7zGzsujvl5rZ183slPQ3LTvNm1lJZ+/gqH1hT9CP7SUrMqOmvITKSAlb9hzlh9ddxLqvfYQfXndR3J7AZD+X8riJpNlwSo6qWug4GGyX3V3Qqz1T4XjvcUqLS0ftKy0u5Xjv8YSvoflpEpZE5rT9PbDIzN4L/Bh4HPgb4JPpbFi2WtFYz13rtgOJDTtmQirScyT7ubJ1mFgkr2QyJccE6UVySSIVB04qO4muvi4iJZGRff2D/ZxUlvgAkuanSVgSCdqG3H3AzD4FfNfdv29mz6e7YdlquHD76PJOqZ+gn0wJqVSk50j2cymPm0geiZNeJFcCt0RLe12/4Hq+9/z3YCDoYesf7GfAB7h+wfVJ3U/z0yQMiQRt/Wb2eeD3gf8U3Vca5/y8l2ztz2Qlu0I1Vb1/yXwu5XETySOx6UXgxHbzmpwJ2hJd0fnF87/I3uN7efS1R+kc6CRSFOFT7/1UwosQRMKUSGWD/wI0Ane7++tmdibwk/Q2q7BNNEdt7cbd454/3EuW6KKDVJhKPVMRyVJH9wRBW6yY9CK5oLW9larS0Z9hvBWdLQdaePHwi7yv9n18ZO5HeF/1u3lx52O0rGmER1YEvY4iWWrSnjZ33w58Peb168C30tmoQjeVocd09/6Nd79MDBOLSAbkQXqRRCsOjOqR6zpC9dE9UGQ0lxfTkIPDwlJYJg3azOzDwB3AGdHzh/O0zU9v0wpXrgw9ZjpQFMlrYS4EWLwyCFYg6GHr6wx+Lrs1M/dPgURLe43Ksfb2HigqpsqKafWBnBwWlsKSyPDoj4HvAEuBBmBRdDttZvYJM3vVzF4zsz8b53iZmf00evzfzaw+FffNdqkcelQuNZEcMLwQoPPw6IUAmRqqy4P0IomW9hqVY62/G4qK6WSIOotO1c6xYWEpLObu8U8w+3d3vzjlNzYrBn4LfAxoBVqAz0eHY4fPuQG4wN2/amZNwKfcfXm86y5atMi3bNmS6uZmXDKrR+NdY3hBQ+wChWwruSVS8B5Z8c7hyd72IHj63Nrw2pWHYleZVh36LZ1DvXRbETcXv4uGoio9dwmDJXpiIqtH/8XM7gF+BvQO73T3X0+hYbEWA6+5+y4AM2sGrgG2x5xzDcHQLMDfAT8wM/PJIs08kIqhR+VSE8kRk9QZldQZlWOtagZ1x97P0cWnAAAgAElEQVSkyWpoKKk4UXUiC4eFE8lBJ/kvkaBtuJdtUcw+By6f5r3nAm/EvG6Nudc7zonmijsGzAIOT/PeBUG51ERyRB4sBMglo3KsjZ1LeNmtWTcsnGgOOsl/iawevSxN9x6vO3BsD1oi52BmK4GVAPPmzZt+y/JEqhc0pGLIVkTGkcKFAOqRSVImq05MUaI56CT/JVJ79F1m9mMz+8fo6wVm9uUU3LsVeHfM6zpg/0TnmFkJcDJwZOyF3H2Nuy9y90WzZ89OQdPyQ6oXNKjWqEiapGghwHCPTFt326gemZYDLWlquExFy4EWVq1fxfInlrNq/apJ/z6J5qCT/JfI6tEHgKeBOdHXvwX+KAX3bgHOMrMzzSwCNBHUNY31OEElBoDPAs8Wwny2VEll0t1kE/6KSJLqlwaT37/yy2A7hd6f2B4ZM6M6Uk1FSQXNO5rT0GCZiqkE1qNWvEaNl4NO8l8ic9pq3f0RM7sFRuaWDU73xtHr3EgQEBYD97v7y2Z2J7DF3R8nSDfykJm9RtDD1jTxFWU8qcqlpvlxImNkYYH1UTnIotQjk12mMtSZaA46yX+J9LR1mtksonPJzGwJcCwVN3f3n7v72e7+Hne/O7rvf0QDNty9x92vdff3uvvi4ZWmknnzZlbS2Ts6Vs/GhL8iGRF2XrUJqEcm+01lqDPRHHSS/xLpabuJYJjyPWb2r8BsgqFKyZBsWACQqqL0InkhSwusp6JHRgsZ0ivRcltjjVrxKgVr0p62aD62jwIfAr4CnO/uL6a7YRLIlgUAYRSlF8laSRZYT3bi+VRNt0dGCxnSr+ncJroHuuno68Dd6ejr0FCnJCyRigjFwJVAPTE9c+7+nbS2bIrypSLCsBse3vqOtB3tPQPU1kT44XUXveP8bOiVE8l7SVQwGJWBP6b3KxuHt1atX/WOXqCOvg5mVczi3kvvDbFl+WWkN7Ojlbpq9WZKaisiPAH0AC8BQ1NtkUxNMgsAYstWza4uY9fBDlau3cKs6gjnnX6SAjiRVEkir1pW5diaZPHEVBcyaEg1ORrqlKlKZCFCnbt/2t1vd/c/H/5Je8sESG4BQGxajre7+tj3dg/u0NEzoLxqIqmURF61rMmxlcDiiaksZNCQqkjmJNLT9o9m9nF3fybtrZF3DG8uOmMGP/v1PmDyBQCxvXL7jvZQXATFRUX0DAxNWHdUw6kiU5RgJv2pTjxPuQQWT0xlIUPSPYlZmCpFJFck0tO2CXjUzLrN7LiZtZvZ8XQ3rBCNt+jgZ7/ex6c/ODehBQCxvXLd/YMUmzE45FSUFgPvHFbNlkUOIvksayaeJ7B4YioLGZLqSczSVCkiuSKRnrZ7gUbgJVUjSK/Y4U1gZLtlz9FxFx2MFZuWo7y0iJ7+YApi3YwK4J3DqhPdb2xvnIhM3XAglOqJ50nPI0uwKP3Y+VbDK18nuk9SPYlZmipFJFckErTtBH6jgC39plp1IHaIs6rsRK9a78AQc06p4OSK0pi6oyeGVVXlQCQzUj3xPHZFauw8sri9YlMoSj/ZfVoOtHCk+wgvt71MeUk582rmUVpUOnFP4tE9QQ9brDipUjJKw7aSAxIZHn0TWG9mt5jZTcM/6W5YIZpK1YGxQ5zuRmfvAN/6zAdYs2IR82dXTTisqioHIrkpoRqjuzcEqUnu+2iwhaSL0se7z3BA5zhnnxL8z+Bvj/4WM5s4eJxxRhAoxhqnty/jNGwrOSKRoO114J+BCFAT8yMptqKxnq6+Adp7Bhhyj+kdq5/wPdMp5D6V+4lI+CadRzZREAJJFaWPd5/YgG5GxQw+MPsDnD/rfGaWz4zf29fXGeS086Fg29cZ7A9T7LCtFQXbSFWwXySLTDo8qvQemTNcdWD0as6z484vm2iIc/v+YwnlbEv2fiKFLhtykk06jyxFc8fi3SfRnG7veF6Lf5+GnetPDENedmv4w5BJDttmwz8DUpgmDNrM7Lvu/kdm9gTRYvGx3P3qtLasQC2ZPyupoGnezMp3VEzo7B2kp3+IU2uCHrgjnb3se7sHGJ2zbXi4VEGaSGKmNJcsDSZNzZGiuWPx7tO8o3nSBQjjPq/2p7n5kvQ/r6QCqwQXaQxfN/Yzvf726/zhs3/IzPKZnDPzHAVwklbxetoeim5XZ6IhMr7x8rZt2XN00jxuFZHikUUJieZsE5H4sqW6waQrUpMIQqZzn8lyuoX1vJIOrqdY4eJoz1He7HoTHDr6O0IL4qVwTBi0ufvW6PaXZjY7+vuhTDVMxilLdaiTZ3ccpH5mJaefUjEqj1tsILei8WzWbtw90gPX3T9IpDh+zjYRmdxUyzylQ9wVqVNYKZrsfRJJZRLW80o6WByucBG7enSCYdvYz7S/cz9FVkRJcQm9g71KLCxpF2941IDbgRsJipkWmdkA8H13vzND7StoY/OoHenso7SoiCOd/cydUTmy/+mXDzBrzLy2ZHO2iWSlLPuSy5rqBpNJIgiZjslSmYT1vKYULE6hwkXPQA+lRaUMDg1SXlIe/z7Di0MiVaMXh0yyglckVrzVo38EfBhocPdZ7j4DuBj4sJn9cUZaV+Bi865BUOWgtNjo7j+RpqNvYJAXWo+9o6oBwG1XLaC2JkJVWTFmMHfG2Jxt9Zn+SCKJCzkNw3BS2eVPLGfV+lW0HGjJnuoGiahfmtRK0XQI63lNpYZqomI/U3lxOf2D/Qz6IHOr58a/j1aoSgrEC9pWAJ9399eHd7j7LuA/R49Jmo3No1ZRWkz/4IkhToA9R7qoLC0eN+XHkvmz+OF1F7H+v18+ac42kawT4pfcREXQgaTLPBWyqZTFSoV0Bouxn6mqtIoiK2JO9RxOjpwc/z4JlBETmUy8hQil7n547E53P2RmpWlsk0TFDnFWlRUzsyrC7p5O5laVM+ROZ+8g3f2DnHfa6LR5481X0ypRyTkhZs+PNyfq3kvvVZA2FRmsqZOu0mGx1x++1sgq1cnuk6LFIVLY4gVtfVM8JikyNo/a/NlVfG5R3ahFBx+oO4WxBcY0X03yQohfclPOQaZ0D6OEmSIl1aXDpn2fFC4OkcIVL2j7gJkdH2e/AeVpao+MMV4P2ZdipqYMrzCF0Sk/YmuMiuSkEL/kEplAny0527JZtqRIyQoZWhwi+S1eyo/iiY5JdhjO4dbR28/hjl4qIiWcd3qNqhpIfkj3l1yclamTJq9FAUkisilFSlZIcIWqyEQmLWMl2Sk2h9uZtdUxPWz1Ctgkf6TrSy6B9AuVJZW8evRVMDhnxjnv6EFTQBJjggA4Z1KkiOSIRArGS4pt2tXGDQ9v5arv/4obHt7Kpl1tSV9jOoXiRQpenJWpw8OejrPw1IWcfcrZdPW/MxF1OtNK5JQ4qVlyKkWKSA5Q0JZhwz1kY/OqJRu4jc3hBqpyIJKwOOkXYoc9zYzqSDUVA300/+ImuO+j8MiKtAUk4+WGy3pxAuCwUn6I5CsNj2bY2CoHU60DOlGheK0aFUlAnJWp7xj27DpC1dHdtBYZVL9npCepYdndKU0rkbMLGyZJzZKpVZwihUBBW4btPdLF7DElp6bSQzY2h5tWjUrOCbNEVZyVqXW7Hx09D+vtPXQWFVFXVHaiJwmCnqTPrU1ZQDKlhQ3ZUOYrRalZlD5FZHIaHs2wsVUOYGo9ZMM53GprIqpyILkn5BJVIytTq2qh42CwjS5CeMewZ38X3WY0Fc088f40JPltbW+lqnT0kG3chQ1hP8Nhi1cGQVpvO/hQsO3rDPYnaKIKFDkxPCySQeppy7BU9pCpyoHkrNh5UDCq9ypjPUUTrEx9Rzb90mqaBkppiA2o0pDkN+mVltnwDCElqVmUPkUkMQraMmxslYN5MyuVV00KT4glqia1ewMNm9fQMByAvPcL8GIz0J7WJL+J5IYbJZue4TRTsyh9ikhiFLSFQD1kUvCytQ7jePnbXmyGC5rgjY1pz2Q/WW64UbL1GcZKcM6d8rmJJEZz2kQk81IwDyotJkpf8cZG+Nxa+Movg22KA7ZEc8ONkq3PcFgSc+6Uz00kMQraRCTz4iwECFWc/G3pNG5uuJIKmnc0T/ymbH2Gw+LkbxtL+dxEEqPhUREJRzbWYQxpyHHKc7qy8RkOS3LOnfK5iUxOPW0iIsNCGnLMy5JYM84Inl2sbJtzJ5JjFLSJiAwLacgxL+d0ZfucO5EcZO4edhtSatGiRb5ly5awmyEiEtfYCgAXzr6QbYe2paQkVtbIhooNItnPEj4xjKDNzGYCPwXqgd3A59z96DjnDQIvRV/udferJ7t2rgRtm3a1jcnVVq80ICIFIrbOaGxONk2+FylICQdtYQ2P/hnwz+5+FvDP0dfj6Xb3C6M/kwZsuWLTrjbuWredw+19zK4u43B7H3et286mXW1hN01EEtRyoIVV61ex/InlrFq/KqmSS1NaLSqJ2b0BHlkB93002Ga6rJdIGoUVtF0DPBj9/UHg90JqR9I27Wrjhoe3ctX3f8UND2+dUqC1duNuKiMl1JSXUGRGTXkJlZES1m7cnermikgaTLdWZtJ1RiUx2VKPVSRNwgra3uXubwJEt6dOcF65mW0xs01mFnpgl6oesr1HuqgqKx61r6qsmL1HJkmmKSJZYbo9ZXm5WjQbJJEbTiQXpS1Pm5n9E3DaOIe+kcRl5rn7fjObDzxrZi+5+3+Mc6+VwEqAefPmTam9iYjtIQNGtms37k5qPtq8mZUcbu8beT9AZ+8g82ZWprK5IpIm062VmXSd0QI0dqFGQgszsqkeq0gapK2nzd2vcPf3jfPzGPCWmZ0OEN0enOAa+6PbXcB6YOEE561x90Xuvmj27Nlp+TyQuh6yFY31dPUN0N4zwJA77T0DdPUNsKKxPoWtFZFpm2B+1HR7yjJaASAH53hNefhZueEkz4U1PPo48PvR338feGzsCWY2w8zKor/XAh8GtmesheOYN7OSzt7BUfum0kO2ZP4sbrtqAbU1EQ519FJbE+G2qxZo9ahINokzPyoVedUaTmvg3kvv5adX/ZR7L703fQFbDs7xmvLws3LDSZ4LK+XHLOARYB6wF7jW3Y+Y2SLgq+7+X83sQ8B9wBBBcPldd//xZNdOZ8qP4TltlZESqsqK6ewdpKtvQAGXSD56ZMU7S1r1tgcJdz+39sTwXTbnVZvkM2Sr5U8sp7aiFrMTmRDcncM9h/npVT+N/2blhpPck3DKj1Bqj7p7G/A74+zfAvzX6O//Brw/w02La7iHbHR+tbNTHrAph5tIikznC3yS+VE5USszi+Z4PfTyQzyw/QGO9x7npLKTuH7B9Xzx/C+Oe25dTR1t3W1UR6pH9iU8/JzN9VhFpkkF45O0ZP6stAZQsb15sStU1ZsnkqThocFI1eihwUTLUoVUPD6lsuQzPPTyQ3zv+e9RYiVUlFTQ1dfF957/HsC4gZsWaoiMT7VHs4xyuImkyHTTP+TD/KgUfobpJBN+YPsDlFgJkZIIZkakJEKJlfDA9gfGPT+jCzVEcoh62rLM3iNdzK4uG7VPOdxEpmC6Q4PDxeNjh1cvuzW3ht5S9Bliy27FruZMNJA63nucipKKUftKi0s53nt8wvfkxPCzSIYpaMsyyuEmkiKpGBrMh/lRKfgMsas5gZFt847mhAKrk8pOoquvi0hJZGRf/2A/J5WdNK12iRQaDY9mGeVwE0mRfBjezBLTLbt1/YLrGfAB+gb6cHf6BvoY8AGuX3B9Glorkr8UtGUZ5XATSZHhocGqWug4GGwTXYQgo0w3mfAXz/8iX1/4dSojlXQPdFMZqeTrC78+4epRERlfKHna0imdedpERApR7Jy22NWcWhwgkhLZnadNRCSfTaluZhYbXs2Z9cmERfKcetpEJDtlKrN9iu+jXikRSVLCPW2a0yYi2SdTNTPTcJ8p180UEZmEgjYRyT6JJsbdvSGor3nfR4NtssHWdBPwjmO6Ky1FRCaioE1Ess/RPUHwFGtsYtxU9JIlcp8kTXelpYjIRBS0iUj2mXFGkFMt1tjEuOP0krVEiln1q1sSL7WUyH2S1HRuE90D3XT0deDudPR1qG6miKSEgjYRyT6JJMYd00vWMtTJ6qLjtPV3jCq1FDdwS0MCXtXNFJF00epREclOk63qfGTFqDJVqwbeoG2ol+ricpizEICOvg5mVczi3kvvnfp9RETSS3naRApOvgUfk9XMXLwymMMGEKmidaiH2iGHWSeGNhNaAJAP9UVFpCBoeFQkH2QqRUY2GVOmqq64is4Z9VA5c+QULQDIoOmu5BWRSamnTSQfxE7KhxPbzWvyuxcpppesKZrUlr6OUUlttQAgA4b/pyFSBdWn0tLZSvOzX6f15NOpm3WuqieIpIh62kTyQRpSV+QaLQAIUcz/NLR4N6uL2mkrMmo7jya2IEREEqKeNpF8MOOMUZPygWmnrshFDac1KEgLw9E9wbA80Dx0hAqM6qJSGOihOlId7N/RrL+NyDSpp00kH6QhdYVIwmLy3bV6P1UUwdAglFYAqgghkioK2kTywZhJ+VTVBq/zeT6bZI+Y/2mosxI6vT8I2k4Jenq1IEQkNTQ8KpIvlLpi6vItXUqmDf9Pw+Y1NL19lNXlxXDyXKoqZtCpihAiKaPkuiJS2GJXPkaqgh6jvk71VE5Dy4EWmnc009rRSl11nVaPisSn5LoiIgkp1HQpaaQFISLpoTltIlLYlC5FRHKEgjYRKWwxKx9HFGC6FBHJfhoeFZHCNqaG6cictstuDbddIjmov7+f1tZWenp6wm5K1ikvL6euro7S0tIpX0NBm4gUtpiVjyOrRy+7VfPZRKagtbWVmpoa6uvrMUt4fn3ec3fa2tpobW3lzDPPnPJ1FLSJiChdikhK9PT0KGAbh5kxa9YsDh06NK3raE6biIiIpIwCtvGl4rkoaBMREZG8UV1dnfC5d9xxB6tXr07J9Z966inOOecc3vve9/Ktb30rqWsmSsOjWWzTrjbWbtzN3iNdzJtZyYrGepbMnxV2s0RERCTG4OAgf/AHf8AvfvEL6urqaGho4Oqrr2bBggUpvY962rLUpl1t3LVuO4fb+5hdXcbh9j7uWredTbvawm6aiIhISmza1cYND2/lqu//ihse3pq277gnnniCiy++mIULF3LFFVfw1ltvjRx74YUXuPzyyznrrLP40Y9+NLL/nnvuoaGhgQsuuIDbb7897vU3b97Me9/7XubPn08kEqGpqYnHHnss5Z9DQVuWWrtxN5WREmrKSygyo6a8hMpICWs37g65ZSIiItOXyc6JpUuXsmnTJp5//nmampr49re/PXLsxRdf5Mknn2Tjxo3ceeed7N+/n2eeeYadO3eyefNmtm3bxtatW3nuuecmvP6+fft497vfPfK6rq6Offv2pfxzaHg0S+090sXs6rJR+6rKitl7pCukFomIiKRObOcEMLJdu3F3yqcCtba2snz5ct588036+vpGpd245pprqKiooKKigssuu4zNmzezYcMGnnnmGRYuXAhAR0cHO3fu5JJLLhn3+uPVcU/Hggz1tGWpeTMr6ewdHLWvs3eQeTMrQ2qRiIhI6uw90kVVWfGofenqnPja177GjTfeyEsvvcR99903Kvnv2ODKzHB3brnlFrZt28a2bdt47bXX+PKXvzzh9evq6njjjTdGXre2tjJnzpyUf45QgjYzu9bMXjazITNbFOe8T5jZq2b2mpn9WSbbGLYVjfV09Q3Q3jPAkDvtPQN09Q2worE+7KaJiIhMWyY7J44dO8bcuXMBePDBB0cde+yxx+jp6aGtrY3169fT0NDAsmXLuP/+++no6ACC4c+DBw9OeP2GhgZ27tzJ66+/Tl9fH83NzVx99dUp/xxhDY/+Bvg0cN9EJ5hZMfBXwMeAVqDFzB539+2ZaWK4lsyfxW1XLRizevRsrR6V1Nm9YXQVgMUrlWBWRDJmRWM9d60LvtKryorp7B2Mdk6cPa3rdnV1UVdXN/L6pptu4o477uDaa69l7ty5LFmyhNdff33k+OLFi7nyyivZu3cv3/zmN5kzZw5z5szhlVdeobGxEQjSfPzkJz/h1FNPHfeeJSUl/OAHP2DZsmUMDg7ypS99ifPPP39an2M8Nt44bKaY2XrgZnffMs6xRuAOd18WfX0LgLv/RbxrLlq0yLdsecflRCTW7g1Bvc1IVfBzbB8ca4Xq2fCu8xXAiciUvPLKK5x33nkJn19oqa0meD4JT37L5oUIc4E3Yl63AheH1BaR/LJ5TRCsldVAVxscewNw6DkOnYeDgG7Z3QrcRCStlsyflddBWqqlbU6bmf2Tmf1mnJ9rEr3EOPvG7RY0s5VmtsXMtky3rpdIQTi6JwjaAN7eC0XFUFIGAz1BIBepCgI7ERHJGmnraXP3K6Z5iVbg3TGv64D9E9xrDbAGguHRad5XJP/NOCPoUSurgf5uKI7A0CCUVgTHI1VBYCciIlkjm1N+tABnmdmZZhYBmoDHQ26TSH5YvBL6OqG3HUrLYbA3CNpOOSM43tcZBHYiIpI1wkr58SkzawUagSfN7Ono/jlm9nMAdx8AbgSeBl4BHnH3l8Nor0jeqV8azFmrqoVIDVgRnDwPKk4JArm+ziCwExGRrBHKQgR3fxR4dJz9+4FPxrz+OfDzDDZNpHDULz2x0GBs+o/LbtUiBBGRLJPNq0dFJFNiAzgRkRxWXV09khR3MnfccQfV1dXcfPPN077+l770JdatW8epp57Kb37zm4Svl4xsntMmIiIikhOuv/56nnrqqbTeQ0GbiIiIhGP3BnhkBdz30WC7e0NabvPEE09w8cUXs3DhQq644greeuutkWMvvPACl19+OWeddRY/+tGPRvbfc889NDQ0cMEFF3D77bdPeo9LLrmEmTNnpqX9wxS0iYiISOYNV2bpPAzVp55I7J2GwG3p0qVs2rSJ559/nqamJr797W+PHHvxxRd58skn2bhxI3feeSf79+/nmWeeYefOnWzevJlt27axdetWnnvuuZS3K1ma0yYiIiKZF1uZBU5sN69J+Rzb1tZWli9fzptvvklfXx9nnnnmyLFrrrmGiooKKioquOyyy9i8eTMbNmzgmWeeYeHChQB0dHSwc+dOLrnkkpS2K1kK2kRERCTzju4JethipSmx99e+9jVuuukmrr76atavX88dd9wxcsxsdAEmM8PdueWWW/jKV76S8rZMh4ZHRUREJPNmnBHkhIyVpsTex44dY+7cuQA8+OCDo4499thj9PT00NbWxvr162loaGDZsmXcf//9I6tE9+3bx8GDB1PermQpaBMREZHMi63M4kMpS+zd1dVFXV3dyM93vvMd7rjjDq699lo+8pGPUFtbO7oZixdz5ZVXsmTJEr75zW8yZ84cPv7xj/OFL3yBxsZG3v/+9/PZz36W9vb2uPf9/Oc/T2NjI6+++ip1dXX8+Mc/ntbnGI+551epzkWLFvmWLVvCboaIiEjBeeWVVzjvvPMSf8PYxN6LV+Z1zsgJno+Nd+54NKdNRCTNWg600Lyjmdb2Vupq6mg6t4mG0xrCbpZI+JTYOykaHhURSaOWAy2sbllNW3cbtRW1tHW3sbplNS0HWsJumojkGAVtIiJp1LyjmYqSCqoj1ZgZ1ZFqKkoqaN7RHHbTRCTHKGgTEUmj1vZWqkqrRu2rKq2itaM1pBaJSK5S0CYikkZ1NXV09o9Oa9DZ30lddV1ILRKRXKWgTUQkjZrObaJ7oJuOvg7cnY6+DroHumk6tynspolIjlHQJiKSRg2nNXBzw83MqpjF4Z7DzKqYxc0NN2v1qEiaVFdXJ3zuHXfcwerVq6d9/TfeeIPLLruM8847j/PPP5+//Mu/TOqaiVLKDxGZugLLsTRVDac1KEgTyWMlJSXce++9fPCDH6S9vZ2LLrqIj33sYyxYsCCl91FPm4hMze4N8PQ3oPNwUD+w83DweveGsFsmIjmi5UALq9avYvkTy1m1flXaUuE88cQTXHzxxSxcuJArrriCt956a+TYCy+8wOWXX85ZZ53Fj370o5H999xzDw0NDVxwwQXcfvvtca9/+umn88EPfhCAmpoazjvvPPbt25fyz6GgTUSmZvOaoLhzWQ1YUbCNVAX7RUQmkckchkuXLmXTpk08//zzNDU18e1vf3vk2IsvvsiTTz7Jxo0bufPOO9m/fz/PPPMMO3fuZPPmzWzbto2tW7fy3HPPJXSv3bt38/zzz3PxxRen/HNoeFREpubonqCHLVakKtgvIjKJ2ByGwMi2eUdzyqcTtLa2snz5ct588036+vo488wzR45dc801VFRUUFFRwWWXXcbmzZvZsGEDzzzzDAsXLgSgo6ODnTt3cskll8S9T0dHB5/5zGf47ne/y0knnZTSzwDqaRORqZpxRlDcOVZfZ7BfRGQSmcxh+LWvfY0bb7yRl156ifvuu4+enp6RY2ajS3+aGe7OLbfcwrZt29i2bRuvvfYaX/7yl+Peo7+/n8985jNcd911fPrTn075ZwAFbSIyVYtXBkFabzv4ULDt6wz2i4hMIpM5DI8dO8bcuXMBePDBB0cde+yxx+jp6aGtrY3169fT0NDAsmXLuP/+++no6ABg3759HDx4cMLruztf/vKXOe+887jppptS3v5hCtpEZGrql8Kyu6GqFjoOBttld2v1qIgkJF05DLu6uqirqxv5+c53vsMdd9zBtddey0c+8hFqa2tHnb948WKuvPJKlixZwje/+U3mzJnDxz/+cb7whS/Q2NjI+9//fj772c/S3t4+4T3/9V//lYceeohnn32WCy+8kAsvvJCf//zn0/oc4zF3T/lFw7Ro0SLfsmVL2M0QEREpOK+88grnnXdewue3HGiheUczrR2t1FXX0XRuU16nx5ng+dh4545HCxFEREQkFMphmBwNj4qIiIjkAAVtIiIiIjlAQZuIiIikTL7NlU+VVDwXBW0iIiKSEpPWqoYAAAbhSURBVOXl5bS1tSlwG8PdaWtro7y8fFrX0UIEERERSYm6ujpaW1s5dOhQ2E3JOuXl5dTVTS8HnYI2ERERSYnS0tJRJaIktTQ8KiIiIpIDFLSJiIiI5AAFbSIiIiI5IO/KWJnZIWBP2O1IQC1wOOxGZDE9n8npGcWn5xOfnk98ej7x6fnEl8zzOezun0jkxLwL2nKFmW1x90VhtyNb6flMTs8oPj2f+PR84tPziU/PJ750PR8Nj4qIiIjkAAVtIiIiIjlAQVt41oTdgCyn5zM5PaP49Hzi0/OJT88nPj2f+NLyfDSnTURERCQHqKdNREREJAcoaAuRmf0/ZvaimW0zs2fMbE7YbcomZnaPme2IPqNHzeyUsNuUTczsWjN72cyGzEyruKLM7BNm9qqZvWZmfxZ2e7KNmd1vZgfN7DdhtyUbmdm7zexfzOyV6L9ffxh2m7KJmZWb2WYzeyH6fP487DZlIzMrNrPnzWxdKq+roC1c97j7Be5+IbAO+B9hNyjL/AJ4n7tfAPwWuCXk9mSb3wCfBp4LuyHZwsyKgb8CfhdYAHzezBaE26qs8wCQUE6oAjUArHL384AlwB/on6FReoHL3f0DwIXAJ8xsSchtykZ/CLyS6osqaAuRux+PeVkFaIJhDHd/xt0Hoi83AXVhtifbuPsr7v5q2O3IMouB19x9l7v3Ac3ANSG3Kau4+3PAkbDbka3c/U13/3X093aCL9654bYqe3igI/qyNPqj764YZlYHXAn8n1RfW0FbyMzsbjN7A7gO9bTF8yXgH8NuhGS9ucAbMa9b0ReuTJGZ1QMLgX8PtyXZJTr0tw04CPzC3fV8Rvsu8CfAUKovrKAtzczsn8zsN+P8XAPg7t9w93cDDwM3htvazJvs+UTP+QbBkMXD4bU0HIk8HxnFxtmnXgBJmplVA38P/NGYUZGC5+6D0Wk9dcBiM3tf2G3KFmZ2FXDQ3bem4/ol6bionODuVyR46t8ATwK3p7E5WWey52Nmvw9cBfyOF2B+miT++ZFAK/DumNd1wP6Q2iI5ysxKCQK2h939Z2G3J1u5+9tmtp5gjqQWtgQ+DFxtZp8EyoGTzOwn7v6fU3Fx9bSFyMzOinl5NbAjrLZkIzP7BPCnwNXu3hV2eyQntABnmdmZZhYBmoDHQ26T5BAzM+DHwCvu/p2w25NtzGz28Ep+M6sArkDfXSPc/RZ3r3P3eoL//jybqoANFLSF7VvRoa4XgY8TrDaRE34A1AC/iKZF+euwG5RNzOxTZtYKNAJPmtnTYbcpbNGFKzcCTxNMIH/E3V8Ot1XZxcz+FtgInGNmrWb25bDblGU+DHwRuDz6351t0V4TCZwO/Ev0e6uFYE5bStNayMRUEUFEREQkB6inTURERCQHKGgTERERyQEK2kRERERygII2ERERkRygoE1EREQkByhoE5GcYmaD0TQML5vZC2Z2k5kVRY8tMrPvhdSuf0vRda6NfrYhM1uUimuKSH5Qyg8RySlm1uHu1dHfTyWoJvKv7p4X1UTM7DyCmoX3ATe7+5aQmyQiWUI9bSKSs9z9ILASuNECl5rZOgAzu8PMHjSzZ8xst5l92sy+bWYvmdlT0VJFmNlFZvZLM9tqZk+b2enR/evN7H+a2WYz+62ZfSS6//zovm1m9uJwZRMz64huzczuiSbOfsnMlkf3Xxq95t+Z2Q4zeziafX/sZ3rF3V/NxPMTkdyioE1Ecpq77yL4b9mp4xx+D3AlcA3wE+Bf3P39QDdwZTRw+z7wWXe/CLgfuDvm/SXuvhj4I07UBf4q8JfRgtmLCOqdxvo0cCHwAYISP/cMB4LAwui1FgDzCbLvi4gkRAXjRSQfvKPHKuof3b3fzF4CioGnovtfAuqBc4D3EZRKI3rOmzHvHy4WvjV6PgQloL5hZnXAz9x955h7LgX+1t0HgbfM7JdAA3Ac2OzurQBmti16zQ3JflgRKUzqaRORnGZm84FB4OA4h3sB3H0I6PcTk3iHCP6n1YCX3f3C6M/73f3jY98fvX5J9Fp/A1xN0Fv3tJldPrZJcZrbG/P7yDVFRBKhoE1EcpaZzQb+GviBT21V1avAbDNrjF6v1MzOn+Se84Fd7v494HHggjGnPAcsN7PiaPsuATZPoW0iIqMoaBORXFMxnPID+CfgGeDPp3Ihd+8DPgv8TzN7AdgGfGiSty0HfhMd3jwXWDvm+KPAi8ALwLPAn7j7gUTbZGafMrNWoBF40syeTvS9IpLflPJDREREJAeop01EREQkByhoExEREckBCtpEREREcoCCNhEREZEcoKBNREREJAf8/+3WAQkAAACAoP+v2xHoCqUNAGBA2gAABqQNAGAgqy0T83UiNu0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#sklearn viz\n",
    "fig=plt.figure(figsize=(10,5))\n",
    "ax=fig.add_subplot(111)\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "for i in set(df['skl']):\n",
    "    plt.scatter(df['dim1'][df['skl']==i],\n",
    "                df['dim2'][df['skl']==i],\n",
    "                label=f'Label {int(i)}',alpha=0.7)\n",
    "\n",
    "plt.title('GMM by sklearn')\n",
    "plt.xlabel('Dimension 1')\n",
    "plt.ylabel('Dimension 2')\n",
    "plt.legend(loc=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimal Number of Cluster\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Since Gaussian Mixture Model is based upon multiple Gaussian distributions. Some properties of Gaussian distributions can be borrowed as a metric to determine the optimal numbers of cluster. Both Akaike Information Criterion and Bayesian Information Criterion are up for the job, or proclaimed in the textbook, lol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute log likelihood\n",
    "#basically the expression before jensens inequality in e step\n",
    "def log_likelihood(data,params):\n",
    "    \n",
    "    mean_vec,cov_mat,w_vec=params\n",
    "    \n",
    "    likelihood=[\n",
    "        multivariate_gaussian(data,\n",
    "                              mean_vec[i],\n",
    "                              cov_mat[i])*w_vec[i] for i in range(len(w_vec))]\n",
    "    \n",
    "    return np.sum(np.log(np.sum(likelihood,axis=0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "#details of aic can be found in the following link\n",
    "# http://www.modelselection.org/aic/\n",
    "def compute_aic(data,params):\n",
    "    \n",
    "    return -2*log_likelihood(data,params)+len(data.columns)*2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "#details of bic can be found in the following link\n",
    "# http://www.modelselection.org/bic/\n",
    "def compute_bic(data,params):\n",
    "    \n",
    "    return -2*log_likelihood(data,params)+len(data.columns)*np.log(len(data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:13: RuntimeWarning: divide by zero encountered in log\n",
      "  del sys.path[0]\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:17: RuntimeWarning: divide by zero encountered in log\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:17: RuntimeWarning: invalid value encountered in subtract\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:17: RuntimeWarning: invalid value encountered in multiply\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:10: RuntimeWarning: overflow encountered in true_divide\n",
      "  # Remove the CWD from sys.path while we load stuff.\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  \n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\numpy\\linalg\\linalg.py:2022: RuntimeWarning: invalid value encountered in det\n",
      "  r = _umath_linalg.det(a, signature=signature)\n"
     ]
    }
   ],
   "source": [
    "AIC,BIC=[],[]\n",
    "\n",
    "for num_of_clusters in range(2,10):\n",
    "    params,posterior=gmm(df,xcols,num_of_clusters)\n",
    "    AIC.append(compute_aic(df[xcols],params))\n",
    "    BIC.append(compute_bic(df[xcols],params))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAFNCAYAAABFbcjcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd8leX9//HXJ4OEPQNkEAhLNgHCxr0A2ShSByoiztbOX+22rba2tWr9oraKWrciW1TcIEPAAGHvJSQBwiYJ2dfvj/vYUkRyMOckJ8n7+Xjk4Tn3Ofd1f5I+HvXtdX+u6zbnHCIiIiISGsIqugARERER+S+FMxEREZEQonAmIiIiEkIUzkRERERCiMKZiIiISAhROBMREREJIQpnIiIiIiFE4UxEQp6ZLTCzo2YWdcbxf5vZQ6e9r2FmD5rZNjPLMbPdZvaCmbUKYm0PmlmhmWX7fjaZ2djTPr/EzPadcc7VZva5mZ00sywzW2hmI4JVo4hULgpnIhLSfMHqQsABpQWY6b7v3ADUB7oDK4HLg1chAG855+o45+oAPwReNbNmZ/uimV0LvA28DCQAzYDfAsODXKOIVBIKZyIS6iYAy4B/A7d825fM7ArgSmCkc+5L51yRc+64c+4p59zzZ/n+A2Y2/Yxj/zCzJ32vbzWznb7ZrV1mdqM/xTrnPgBOAm3Ock0DHgP+6Jyb6quvxDm30Dl3hz/ji0jVp3AmIqFuAvCa7+fqb5uRAq4AVjjn9vo57hvAUDOrB2Bm4cA44HUzqw08CQxxztUFBgBppQ1onmuAGsDGs3zlAqAF3gyfiMhZKZyJSMgys0FAS2Cac24lsAPvluXZNAYy/R3bObcHWAWM8h26DMh1zi3zvS8BuphZTedcpnNuwzmGG2dmx4AcYC7wJ+fcsW+pkfOpU0SqH4UzEQlltwAfOucO+d6/zrff2jwMxJ7n+K8D3/O9vsH3HudcDnA9cBeQaWbvmlmHc4wzzTnXwDlXC+925gQzu/NbauQ71Cki1YjCmYiEJDOriXeb8WIz229m+4EfAd3NrPtZTvkY6GNmCedxmbeBS3znjMYXzsDrHXPOXYkXpDYDz/kzoHNuN/A+Z2/w3wLsBcae5TMREUDhTERC1yigGOgEJPt+OgKL8PrQ/odz7mPgI2CWmfUyswgzq2tmd5nZxLNdwDmXBSwAXgR2Oec2AZhZMzMb4es9yweyfbWUyhf0BgPfuA3qnHPAj4HfmNltZlbPzMLMbJCZPevP+CJS9SmciUiougV40Tn3lXNu/9c/wBTgRjOLOMs51wLvAW8Bx4H1QArerNq3eR1vMcHrpx0LA34CZABHgIuBe84xxvVf73MGfAksAX5/ti8656bj3TKd6Bv/APAQMOcc44tINWLef8iJiIiISCjQzJmIiIhICFE4ExEREQkhCmciIiIiIUThTERERCSEKJyJiIiIhJCzLUWvNAYPHuzmz59f0WWIiIiI+MP8+VKlnjk7dOhQ6V8SERERqUQqdTgTERERqWoUzkRERERCiMKZiIiISAhROBMREREJIQpnIiIiIiEk6OHMzMLNbLWZzfO9/7eZ7TKzNN9Psu+4mdmTZrbdzNaaWc9g1yYiIiISaspjn7P7gU1AvdOO/cw5N/2M7w0B2vl++gLP+P4pIiIiUm0EdebMzBKAa4Cpfnx9JPCy8ywDGphZbDDrExEREQk1wb6t+QTw/4CSM44/7Lt1+biZRfmOxQN7T/vOPt8xERERkWojaOHMzIYBB51zK8/46BdAB6A30Aj4+dennGUYd5ZxJ5tZqpmlZmVlBbLkbyouhKVT4Nje0r8rIiIiEgDBnDkbCIwws93Am8BlZvaqcy7Td+syH3gR6OP7/j6gxWnnJwAZZw7qnHvWOZfinEuJiYkJYvnAvlT48FfwRBd4YQikvgC5R4J7TREREanWghbOnHO/cM4lOOdaAeOBT51zN33dR2ZmBowC1vtOmQtM8K3a7Accd85lBqs+v7TsDz9Ig0t/DbmHYN6P4NH28Pp4WDcdCnIrtDwRERGpespjteaZXjOzGLzbmGnAXb7j7wFDge1ALnBbBdT2TY2S4OKfwUU/hf1rYd3bsG4GbH0fImtDx2HQdRy0vgTCK+LPKSIiIlWJOfeNtq5KIyUlxaWmppb/hUuKYc9SL6htnA15x6FWE+gyBrpeBwm9wc7WQiciIiLVmF/hQOGsrIryYfvHsHYabJ0PRXnQoKUX0rqNg5gLKrY+ERERCRUKZ+Uu7wRsnufNqO1cAK4Emnf1bnt2GQv1tTOIiIhINaZwVqFOHoANs2DdNEhfCRi0GgRdr4VOI6Fmw4quUERERMqXwlnIOLzDW925bhoc3g5hkdDuKuh2HbQfDJE1K7pCERERCT6Fs5DjHGSmwdq3Yf0MyN4PNepCx+HejFrSxVrxKSIiUnUpnAVCTn4RtaOCEJhKimH3It+Kz7mQfwJqN/Wt+BwH8T214lNERKRqUTgrqz2Hc7jy8c+57IKmjEyO49IOTYmODA/8hQrzYNuH3m3PrR9AcQE0au2t+Ox6HTRpF/hrioiISHlTOCur9GOnmLpoJ++syeRQdj51oyIY3KU5o3rE0691Y8LDgjCzdeoYbHrHm1Hb9TngIDbZC2ldxkK92MBfU0RERMqDwlmgFBWX8MXOw8xencEHG/aTnV9E07pRDO8ex8jkOLrG18eCcQvyRCZsmOntoZaZBhgkXejd9uw0AqLrB/6aIiIiEiwKZ8GQV1jMp5sPMnt1Ogu2ZFFQXELrJrUZkRzHyOR4kprUDs6FD23zZtPWToOjuyA8Ctpf5QW1dldBZHRwrisiIiKBonAWbMdzC3l/fSZz0jJYtuswzkH3hPqMTI5nWPdYmtYNQmByDtJXeUFt/QzIOQhR9aHTcO/WZ6sLISwIfXEiIiJSVgpn5Snz+Cnmrclkdlo6GzJOEGYwsG0TRnSPY3CX5tSNjgz8RYuLYPfn3tYcm96BgpNQp7nXm9b1WojroRWfIiIioUPhrKJsP3iSOWkZzEnL4KsjudSICOOKjk0ZmRzPJRfEEBURjBWfp7xne66b7q38LC6Axm3/u+KzcZvAX1NERETOh8JZRXPOsXrvMeasTmfe2kwO5xRQLzqCoV1jGZEcR7+kxoQFZcXnUW/vtHVvw+7FgIO4nt6D2DuPgbrNAn9NERERKY3CWSgpKi5h8fZDzE3zVnzmFBTTvF40w7vHMjI5ns5x9YKz4vN4utebtu5t2L8WLMx7EkHX67wnE0TXC/w1RURE5GwUzkLVqYJiPt50gDlp3orPohJHm5jajEqOZ0RyHC0bB2nFZ9YWL6StexuO7oaIaO/Znl2vg3ZXQkRUcK4rIiIioHBWORzNKeA934rPFbuOANAjsQEju8cxrHscTeoEITA5B/tSvScSrJ8JuYe8PdM6jfS25mg5EMLCAn9dERGR6k3hrLJJP3aKd9ZkMHt1Opv3nyQ8zBjYtgmjkuO4qnNz6gTjGZ/FRbBzgTebtnkeFGRD3TjoOtabUWveTSs+RUREAkPhrDLbsv8kc9LSmZOWQfqxU0RHhnFFx2aMSo7novYx1IgIwsxWQS5sfd/bmmP7R1BSBE0u8K34vBYaJQX+miIiItWHwllV4Jxj5Z6jzEnL4N11mRzJKaBBrUiGdo1lZPc4erdqFJwVn7lHYONsL6h9tdQ7ltDbu+3ZeTTUiQn8NUVERKo2hbOqprC4hMXbDjE7LZ0PNxzgVGExcfWjGZ4cx6jkeDo0rxucFZ/H9v53xeeB9WDh0OZSb0atwzUQVTfw1xQREal6FM6qstyCIj7aeIA5aRl8vtVb8dm+WR1GJsczonscLRrVCs6FD2z0rficDse/goia0PNmuOohrfYUERE5N4Wz6uJITgHvrstkzup0UvccBaBXy4aMSo5jaNdYGgdrxefe5bD6FVj9KsSnwPWvQL24wF9LRESkalA4q472Hsll7poM5qSls/VANhFhxoXtmjAyOZ4rOzWjdjBWfG6cC7Pugqg6MO5lSOwX+GuIiIhUfgpn1d2mzBPMSctgblo6GcfzqBkZzpWdmjGqRxwXtoshMjyAKz4PboI3b/D604b+FVImBm5sERGRqkHhTDwlJY7UPUeZnZbOe+syOZZbSMNakVzTLZZRyfH0TGwYmBWfp47CjDu8bTh63gJD/6Y+NBERkf9SOJNvKigq4fOtWcxZk8FHG/eTV1hCfIOajEyOY2RyPBc0L+PKy5Ji+PQhWPwYJPTxbnPWiw1M8SIiIpWbwpmcW3Z+ER9t3M/s1Rks3n6I4hJHh+Z1vRWfyXHEN6j53QffMBtm3+Nts3H9K9CiT+AKFxERqZwUzsR/h7LzeXdtJnPS0ln11TEA+rRqxMgecQztEkvD2jXOf9ADG+HN78HxdLjmUeh1a2CLFhERqVwUzuS7+epwLnPS0pmdls6OrBwiw42L28cwIjmeKzs2o2aNcP8Hyz0CM26HHZ9Cr9tgyF8h4jsEPRERkcpP4UzKxjnHxv+s+Mxg/4k8atUI5+rOzRmZHMegtk2I8GfFZ0kxfPIHWPIEtOjn9aHVbRb8X0BERCS0KJxJ4JSUOJbvOsLcNem8uzaTE3lFNK5dg2HdYrmhb0v/FhKsnwFz7oPo+nD9q5CQEvzCRUREQofCmQRHflExC7dkMSctg483HcABD43qwriUFqWfvH+9tx/ayUy45u/Qc0LQ6xUREQkRCmcSfIez8/nBm6tZsv0wt/Rvya+HdSp9c9vcIzB9Iuz8DHpPgqv/rD40ERGpDvwKZwHcIl6qo8Z1onjptj5MGpTES1/s4capyzmUnX/uk2o1ghunw4AfwJdT4eURkH2wfAoWEREJcUEPZ2YWbmarzWye732SmS03s21m9paZ1fAdj/K93+77vFWwa5PAiAgP49fDOvHE9cms2XuMEf+3mLX7jp37pPAIuOqPMPZ5yEiDf10M+1aWT8EiIiIhrDxmzu4HNp32/i/A4865dsBR4Hbf8duBo865tsDjvu9JJTKqRzwz7h6AmXHtP79gxsp9pZ/U9Vq4/UMvrL04BFa/GvxCRUREQlhQw5mZJQDXAFN97w24DJju+8pLwCjf65G+9/g+v9z3falEusTXZ+59A+mV2JCfvL2G37+zgcLiknOfFNsN7lgAiX1hzr3w3s+guLBc6hUREQk1wZ45ewL4f8DX/3ZuDBxzzhX53u8D4n2v44G9AL7Pj/u+L5VM4zpRvHJ7HyYOTOLFJbu5+fnlHC6tD612Y7hpFvS/D1Y8Cy+PhOys8ilYREQkhAQtnJnZMOCgc+70RqKzzYQ5Pz47fdzJZpZqZqlZWfqXd6iKCA/jt8M78di47qz+6hgjpixhffrxc58UHgFXPwxjpkL6Snj2EkhfVS71ioiIhIpgzpwNBEaY2W7gTbzbmU8ADcwswvedBCDD93of0ALA93l94MiZgzrnnnXOpTjnUmJiYoJYvgTCmJ4JTL9rAM45xj6zlFmr/ehD63YdTPwAzOCFwZD2RvALFRERCRFBC2fOuV845xKcc62A8cCnzrkbgc+Aa31fuwWY43s91/ce3+efusq8CZv8R9eE+sz9/iCSWzTgR2+t4Y/zNlJUWh9aXDJMXgAt+sDsu+D9B9SHJiIi1UJF7HP2c+DHZrYdr6fsed/x54HGvuM/Bh6ogNokSJrUieLVSX25dUArnl+8iwkvrOBITsG5T6rdBG6eDf3ugeXPwCujIedQ+RQsIiJSQfSEACl301fu45ez1hFTJ4pnJ/Sic1z90k9a8ya8cz/UjvGeyxmXHPxCRUREAktPCJDQdG2vBN6+sz8lvj60OWnppZ/UfTxMnA/OwQtXw5q3gl+oiIhIBVA4kwrRvUUD5t43iG7xDbj/zTQeftefPrQeXh9afArMmgzzfwnFRec+R0REpJJROJMKE1PX60Ob0L8lzy3axa0vfsnR0vrQ6sTAhNnQ9y5Y9hS8OhpyDpdPwSIiIuVA4UwqVI2IMP4wsgt/vbYbK3YdYfiUxWzMOHHuk8IjYchfYNQz8NVybz+0zDXlUq+IiEiwKZxJSBiX0oK37uxHYXEJY55ZwjtrMko/KfkGmPg+uGJ4/mpYN730c0REREKcwpmEjB6JDXnn+4PoElef77+xmj+/v4niklJWE8f38vrQ4nrAjNvhg1+pD01ERCo1hTMJKU3rRvP6Hf24qV8i/1q4k1tfXMGx3NL60JrCLXOhz2T4Ygq8NhZyv/FwCRERkUpB4UxCTo2IMB4a1ZVHxnRl+c4jjJiyhM37/ehDG/o3GPkU7FkKz14M+9eVT8EiIiIBpHAmIWt8n0TevLMfeYXFjH5qKe+uzSz9pB43wW3zvVubU6+E9TOCX6iIiEgAKZxJSOuZ2JB53x9Ep7h63Pv6Kv4yf3PpfWgJX/ehJcP0ifDRb6GkuDzKFRERKTOFMwl5TetF88Yd/bihbyLPLNjBxH9/yfHcUh6CXrcZTJgLKbfDkn/Aa9eqD01ERCoFhTOpFGpEhPGn0V350+iuLN1xiBFPLWbL/pPnPimiBgx7DIY/CbsXw3OXwoEN5VOwiIjId6RwJpXKDX0TeXNyP3ILihn99BLeX+dHH1qvW+DW96AwD6ZeARtmBb9QERGR70jhTCqdXi0bMe/7g7igeV3ufm0Vj36wpfQ+tBa94c6F0LwrvH0rfPyg+tBERCQkKZxJpdSsXjRvTu7H+N4tmPLZdia99CXHT5XWh9YcbpkHvW6DxY/D6+Pg1NHyKVhERMRPCmdSaUVFhPPnMV15aFQXFm07xKinlrDtgB99aMOfgGFPwM6F8OylcGBj+RQsIiLiB4UzqdTMjJv6teSNyf04mVfEqKeWMH/9/tJPTLkNbn0XCnO9PrSNc4JfrIiIiB8UzqRK6N3K60Nr26wud726ksc+3EJJaX1oiX1h8kJo1gmmTYBP/qg+NBERqXAKZ1JlNK8fzVuT+zEuJYEnP93OHS+nciKvlD60erHeDFrPCbDoUXhjPJw6Vj4Fi4iInIXCmVQp0ZHh/GVsN/4wsjMLt2YxasoSth8srQ8tytsL7ZrHYMdn8NxlcHBz+RQsIiJyBoUzqXLMjAn9W/HapL6cyCtk1FNL+WjjgdJOgt63w63zIP8kTL0cNr1TPgWLiIicRuFMqqy+rRsz975BtI6pzR0vp/L4R1v96EPr5+2HFnMBvHUTfPowlJSUT8EiIiIonEkVF9egJtPu7M/Yngn845NtTH5lJSdL7UOL854o0OMm+Pyv8Ob3IO94+RQsIiLVnsKZVHnRkeE8el03Hhzeic+2HGTUU0vYkZV97pMio2HEFBj6KGz/2OtDy9pSPgWLiEi1pnAm1YKZcevAJF6b1JdjuYWMmrKEj/3pQ+tzB9zyjjdz9tzlsPm98ilYRESqLYUzqVb6tW7M3O8PolWT2kx6OZUnP9lWeh9aywEweQE0aevd4vzsz+pDExGRoFE4k2onvkFN3r6rP2N6xPPYR1u569WVZOcXnfuk+glw23xIvhEWPgJv3Qh5J8qnYBERqVYUzqRaio4M5+/juvPbYZ34ZLPXh7bTnz60kU/BkL/B1g+87TYObSufgkVEpNpQOJNqy8yYOCiJV27vw5GcAkY+tYRPN/vRh9Z3MtwyF3KPeAsFtrxfPgWLiEi1oHAm1d6ANk2Ye99AEhvV4vaXUpny6TacK6UPrdUgrw+tUWvvkU8L/qI+NBERCQiFMxEgoWEtpt81gJHd43j0w63c/eqq0vvQGrSAifOh+/dgwZ9g2s3qQxMRkTJTOBPxqVkjnMevT+bX13Tko00HGP3UEnYfyjn3SZE1YdQzMPgR7/bm1Cvg0PbyKVhERKokhTOR05gZky5szcsT+3AoO58RUxazYMvB0k6CfnfDhNmQe8jrQ9v6QfkULCIiVY7CmchZDGzbhLn3DSK+YS1u+/eXPPXZ9tL70JIu8vrQGraE16+Hz/8GpZ0jIiJyhqCFMzOLNrMVZrbGzDaY2e99x/9tZrvMLM33k+w7bmb2pJltN7O1ZtYzWLWJ+KNFo1rMvHsAw7vF8bcPtnDv66vIKbUPLRFu/xC6jYNPH/L60PJPlk/BIiJSJQRz5iwfuMw51x1IBgabWT/fZz9zziX7ftJ8x4YA7Xw/k4FnglibiF9q1gjnH+OT+dXQjsxfv58xTy9lz2E/+tBG/wuu/pP3uKepV3i3ObWaU0RE/BC0cOY8X+/qGen7Odc9npHAy77zlgENzCw2WPWJ+MvMuOOi1rw0sQ8HTuYxYsoSFm7NKu0k6H8v3DzLmzl7fRw83RdSX4TCU+VTuIiIVEpB7Tkzs3AzSwMOAh8555b7PnrYd+vycTOL8h2LB/aedvo+3zGRkHBhuxjm3juI2PrR3PbiCv65cEfpfWitL4b718CYqd6M2rwfwuOd4dOH4WQpG96KiEi1FNRw5pwrds4lAwlAHzPrAvwC6AD0BhoBP/d93c42xJkHzGyymaWaWWpWVimzFyIBlti4FjPvGcDQrrE88v5m7ntjNbkFpfShhUdCt+tg8kK49V1o0ddbLPBEF5h9LxzYWD7Fi4hIpWCl/pd/oC5k9jsgxzn36GnHLgF+6pwbZmb/AhY4597wfbYFuMQ5l/ltY6akpLjU1NQgVy7yTc45/vX5Tv46fzPtm9Xl2ZtTSGxcy/8BDm2H5c/A6teg6BS0ucy7Ddrmcu+WqIiIVEV+/R98MFdrxphZA9/rmsAVwOav+8jMzIBRwHrfKXOBCb5Vm/2A4+cKZiIVycy46+I2vHhbHzKOnWLEU4tZvO2Q/wM0aQvX/B1+vBEu+w0c2ACvjoWn+8OqV6AwL3jFi4hISAvazJmZdQNeAsLxQuA059wfzOxTIAYvPaYBdznnsn1hbQowGMgFbnPOnXNaTDNnEgr2HM5h8ssr2XbwJA8M6cAdF7bGznf2qygf1s+EL6bAgfVQOwZ63wG9b4faTYJTuIiIlDe//uVQbrc1g0HhTEJFTn4RP5u+hvfW7WdE9zj+MrYbNWuEn/9AzsGuz72Qtu1DiIiGbtd7tzxjLgh84SIiUp4UzkTKk3OOZxbu4G8fbKFj83r86+ZetGh0Hn1oZ8raAsuehjVvQlEetLvKC2lJF6svTUSkclI4E6kIn205yP1vrCY8zJhyQ08Gti3jbcmcQ/Dl8/Dlc5CTBc26eCGty7UQUSMwRYuISHlQOBOpKLsP5TD5lVS2H8zml0M7MnFgEmFhZZztKsyDdW/DF09B1iao0xz63AEpE6FWo8AULiIiwaRwJlKRsvOL+Om0NczfsJ+kJrWZOCiJsT3jqVUjomwDOwc7PvFC2o5PIaImJN8A/e7xVoGKiEioUjgTqWjOOeatzWTqop2s2XecBrUiubFvIhP6t6JZveiyX+DARlj2FKydBsWFcMEQ75Zny4HqSxMRCT0KZyKhwjnHyj1HeW7RTj7ceICIMGN49zgmDWpNp7h6Zb9A9kH4cqr3k3sYYrtD//ug82jvCQUiIhIKFM5EQtGewzm8uGQ301L3kltQzMC2jZk0qDUXt48JQF/aKW9157Kn4dBWqBsHfe+EXrdAzYaB+QVEROS7UjgTCWXHcwt5fcVX/HvpLg6cyKdNTG1uH9SaMT3jiY78Dnukna6kBLZ/7O2XtmshRNaGHjdBv7uhUVJgfgERETlfCmcilUFBUQnvrcvkuUU72ZBxgka1a3BTv5bc3K8lMXWjyn6BzLXeTNq66VBSBB2Hebc8W/RVX5qISPlSOBOpTJxzLN91hKmLdvLxpoPUCA9jVI84bh/Umgua1y37BU5kenulffk85B2D+F7e4oGOIyG8jCtIRUTEHwpnIpXVzqxsXliyi+kr95FXWMJF7WOYNCiJC9s1Of/ndp6pIAfSXvdm047shPotoO9d0PNmiK4fmF9ARETORuFMpLI7mlPg60vbTdbJfNo3q8OkQa0ZkRwXmL60rfO9/dL2LIYadaHnBOh3FzRIDMwvICIip1M4E6kq8ouKeWeNt1/a5v0naVKnBhP6t+LGvok0rhOAvrSM1V5I2zALXAl0Gun1pSWklH1sERH5msKZSFXjnGPpjsM8t2gnC7ZkERURxpieCdw+KIm2TeuU/QLH98GKZyH135B/3Fs00P9e6DAMwso4UyciIgpnIlXZtgMneWHJLmasSqegqITLOjRl0qAk+rdpXPa+tPxsWP2q15d2bA80aOk9HqrHjRAVgMUJIiLVk8KZSHVwKDuf15Z9xSvLdnMou4COsfWYNCiJ4d3jqBERVrbBS4ph87veLc+9yyCqvrehbd87oX5CYH4BEZHqo2zhzMyigbrOuawzjjcFTjjn8spcYhkpnIn8V15hMXPS0pm6aBfbDmbTtG4Utwzw+tIa1KpR9gvsS/U2td04BzDv0VAD7oO4HmUfW0SkeihzOHsWmO+cm3nG8RuBQc65u8tcYhkpnIl8k3OOz7cdYuqinSzadoiakeFc2yuBiYOSSGpSu+wXOLrH60tb+RIUnPQest7/Xmg/BMLKOFMnIlK1lTmcbXTOdfqWzzY45zqXobiAUDgTObfN+0/wwuJdzF6dQWFJCVd0bMakQUn0SWpU9r60vBOw+hVY9k84/hU0au31pSXfADUCEAJFRKqeMoezTc65juf7WXlSOBPxz8GTebz6xR5eWbaHo7mFdI2vz6QLkxjaNZbI8DLOdhUXwaa53i3P9JUQ3QBSJkKfyVAvNjC/gIhI1VDmcLYQ+JlzbsUZx3sDf3fOXVTmEstI4Uzk/JwqKGbm6n08v3gXO7NyiK0fza0DWjG+TyL1a0aWbXDnYO8KL6RtngcWDl3Gerc8Y7sF5hcQEancyhzO+gDTgH8DK32HU4AJwHjn3PKy11g2Cmci301JiWPB1oNMXbSLpTsOU6tGOONSWjBxYBKJjWuV/QJHdsHyf8KqV6AwB5Iu8ja1bXul+tJEpDor+1YaZtYMuAfo4ju0AZjinDtY5vICQOFMpOw2ZBzn+cW7eGdNBsUljqs6NeeOi5Lomdiw7H1pp47Bqpdg+b/gRDo0ae/1pXUfD5E1A/MLiIhUHtrnTET8t/94Hi9/sZvXln/F8VOFJLdowKQLkxjcuTkRZe5LK4QNs+GL/4MVDVx6AAAgAElEQVTMNVCrMaTcDr0nQd1mAalfRKQSKPNtzXXA2T40wDnnKryJROFMJPByC4qYsdLrS9t9OJf4BjW5bWArru/dgrrRAehL27PU29R2y3sQHgldx0GfOyC2O5R1pk5EJLSVOZy1PNeJzrk936GogFI4Ewme4hLHJ5sOMHXxLlbsOkKdqAjG927BrQNbkdAwAH1ph3d4j4da/RoUnYKYjt7tzm7joF5c2ccXEQk9gb+taWZNgMMuRO6FKpyJlI+1+47x/OJdzFubCcDgLs2ZNCiJHokNyz74qaOwfiaseRP2rQAMWl8M3cZDx+EQFYAHuouIhIYyz5z1Ax4BjgB/BF4BmgBhwATn3PzA1PndKZyJlK+MY6d4aeluXl/xFSfzikhp2ZBJFyZxZafmhIcF4Jbk4R2w9i0vqB3bA5G1vIDW7XpofQmEhZf9GiIiFafM4SwV+CVQH3gWGOKcW2ZmHYA3nHMV/kA9hTORipGdX8TbqXt5Ycku9h45RWKjWkwc2IrrUlpQOyqi7BdwDvYuhzVvwIZZkHcc6sZC12u9GbXmXUofQ0Qk9JQ5nKU555J9r//niQBmtlrhTESKSxwfbtjP1MW7WLnnKHWjI7ihbyK3DmhFbP0AbZVRmAfbPvBm07Z9CCVF0KwrdL8eul4HdZsH5joiIsFX5nC2yjnX88zXZ3tfURTORELHqq+O8vziXby/LpMwM67pFsukQa3pmlA/cBfJOQwbZnozaukrwcKg9aXeQoIO1+iZniIS6soczoqBHN9ANYHc0waOds6VcU192SmciYSevUdyeWnpbt78ci/Z+UX0TWrEpAtbc3mHpoQFoi/ta4e2ebNpa6d5D16vUQc6jvCCWqsL9SQCEQlF2oRWRCrOibxCpn25lxeX7Cb92CmSmtRm4sBWjO2VQK0aAehL+1pJCXy11AtqG+dA/gmoF+/d8uw+Hpp2LH0MEZHyEZAFAUuA94EFzrm8wNUWGApnIqGvqLiE99fvZ+qinazZd5wGtSK5sW8iE/q3olm96MBerPCUt7ntmrdg+8fgir3NbbuN9xYT1Gka2OuJiJyfMoezCGAQMBi4FDgMfAC875zbWurAZtHA50AUEAFMd879zsySgDeBRsAq4GbnXIGZRQEvA71817reObf7XNdQOBOpPJxzrNxzlKmLdvHBxv1EhBnDu8cxaVBrOsXVC/wFs7Ng/XRvRi0zDSwc2l7ubcvR4Ro921NEKkJgb2uaWSwwBC+stQO+cM7dc47vG1DbOZdtZpHAYuB+4MfATOfcm2b2T2CNc+4ZM7sH6Oacu8vMxgOjnXPXn6smhTORymnP4RxeXLKbaal7yS0opnerhoztmcCQrrHUrxmEdtaDm2Gtrz/tRDpE1YNOI6D79yBxgPrTRKS8lHnmLBqo65zLOuN4M+Ak0MM5t8Svi5jVwgtndwPvAs2dc0Vm1h940Dl3tZl94Hv9hW/Wbj8Qc66nESiciVRux3MLefPLr5iWupcdWTnUiAjjyo7NGNMznovaxxBZ1geun6mkBHYv8ja63TgHCrKhfqL3yKju46FJu8BeT0Tkf5U5nD0LzHfOzTzj+I3AIOfc3aUObhYOrATaAk8BfwOWOefa+j5vgXebtIuZrQcGO+f2+T7bAfR1zh36tvEVzkSqBucca/cdZ9bqdOauyeBITgGNa9dgePc4xvSMp2t8fSzQD0UvyIXN73ozajs+BVcCcT292bQuY6F248BeT0QkAOFso3Ou07d8tsE519nvSswaALOA3wIvnhHO3nPOdTWzDcDVZ4SzPs65w2eMNRmYDJCYmNhrz54Kf/66iARQYXEJC7dkMWt1Oh9tOkBBUQltYmozpmcCo3rEE98gCL1iJ/fDuuleUNu/DsIioO2V3mxa+8EQGeCFCyJSXZU5nP3PUwH8/exbL2T2O7y90n6ObmuKiB+OnyrkvXWZzFqVzordRwDo17oRY3omMKRLc+pGB6E/7cCG/+6flr0foupD51G+/rR+EOgZPBGpTsoczhYCP3POrTjjeG/g7865i845sFkMUOicO2ZmNYEPgb8AtwAzTlsQsNY597SZ3Qt0PW1BwBjn3LhzXUPhTKT62Hskl1mr05m5ah+7D+cSHRnGVZ2aM7pnPBe2bUJEwPvTimHXQm9bjk1zoTAXGrT0ZtO6XQ+N2wT2eiJSHZQ5nPUBpgH/xusbA0gBJgDjnXPLzzmwWTfgJSAcCAOmOef+YGat+e9WGquBm5xz+b4FCK8APYAjvmvsPNc1FM5Eqh/nHKv3HmPWqnTeWZvBsdxCmtSJYmRyHKN7xNM5rl7g+9Pys2HzPO+xUTsXAg4S+njP9+w8Bmo1Cuz1RKSqKvtWGr6VmfcAXXyHNgBTnHMHy1xeACiciVRvBUUlfLblILNWpfPJ5gMUFjvaN6vj9aclx9O8fhB6xU5kwLq3vVufBzdCWCS0v9qbUWt3FUREBf6aIlJV6PFNIlJ9HMstYN7aTGau2seqr45hBgPaNGZMjwQGd2lO7agAPjIKwDlv8cDat7z+tJyDEN3AW+nZfTwk9FZ/moicqcy3NdcBZ/vQgBLnXPfvXltgKJyJyNnsOpTDrNXpzFq9j71HTlEzMpzBXZozukc8A9s2ITyQD2AHKC6CnQu8256b34WiU9CotffYqG7joFFSYK8nIpVVmcNZy28ZNAH4pXNu6HevLTAUzkTkXL5+ZNSMVem8uzaDE3lFNKsXxcjkeEb3iKdjbBAeG5V3Aja94wW13YsBB4n9vUUEnUdBzYaBv6aIVBaBu61pZsnADcA4YBfeasspZSovABTORMRfeYXFfLb5IDNWpbNgy0GKShwdY+sxpkc8I5PjaBroh7ADHNsL66Z5Kz4PbYHwKLhgsDej1vYKiKgR+GuKSCgr88xZe2A88D28B5G/BfzUOXe2GbUKoXAmIt/FkZwC5q3NYMaqdNbsPUaYwaB2MYzpEc9VnZtRq0YQ+tMy07xFBOumQ+4hqNXY60/rNh7ie6o/TaR6KHM4KwEWAbc757b7ju10zrUOWIllpHAmImW1IyubWavSmbU6nfRjp6hdI5zBXWIZ0zOefq0bB6E/rRC2f+I9jWDze1CcD43bedtydLseGiQG9noiEkrKHM5G482cDQDm4+1NNtU5FzKdrQpnIhIoJSWOFbuPMGtVOu+ty+RkfhGx9aMZmRzPmJ7xtG9WN/AXPXXMewD72rdgzxLvWMtBXlDrNBKi6wf+miJSkQLTc2ZmtYFReLc3L8PbWHaWc+7DslZYVgpnIhIMeYXFfLTxALNWp7NwaxbFJY4u8fUY3SOBEd3jiKkbhL3Mju7xtuRY+yYc3g4R0XDBUO+xUW0ug/AA32oVkYoQ+H3OzKwRcB1wvXPusu9YWMAonIlIsB3KzmduWgazVqezLv044WHGRe2aMLpnAld1akZ0ZHhgL+gcpK/0+tPWT4dTR6F2DNz4NsT1COy1RKS8aRNaEZFA2nbgJDNXpzN7dTqZx/OoGxXBkK7NGd0jgb5JjQgLdH9aUQFs/wjWz4AR/wc1agd2fBEpbwpnIiLBUFLiWLbzMDNXp/P+ukxyCoqJb1CTUT3iGN0jgbZN61R0iSISmhTORESC7VRBMR9u3M/MVeks2pZFiYPuCfUZ0zOB4d3jaFRbe5mJyH8onImIlKeDJ/KYuyaDmavS2Zh5gogw45ILYhjTM4HLOjQNfH+aiFQ2CmciIhVl8/4T/9k/7eDJfOpGRzCsWyxjeiaQ0rIhpk1nRaojhTMRkYpWXOJYuuMQs1al8/76/ZwqLKZFo5qMTo5ndM8EkpqoyV+kGlE4ExEJJTn5RXywwetPW7LjEM5Bj8QGXn9at1ga1FJ/mkgVp3AmIhKq9h/PY05aOjNXpbPlwEkiw43LOjRldI8ELu0QQ1SE+tNEqiCFMxGRUOecY2Om1582Oy2DQ9n51K8ZyfDusYzukUDPxAbqTxOpOhTOREQqk6LiEhZvP8TMVel8uHE/eYUltGpci9E9Erihb2JwHhslIuXJr3Cmh7WJiISIiPAwLrmgKZdc0JSTeYXMX+/1pz3xyVaGdY9VOBOpJjRzJiIS4g6ezKNp3eiKLkNEys6vmbOwYFchIiJlo2AmUr0onImIiIiEEIUzERERkRCicCYiIiISQhTOREREREKIwpmIiIhICFE4ExEREQkhCmciIiIiIUThTERERCSEKJyJiIiIhBCFMxEREZEQonAmIiIiEkIUzkRERERCiMKZiIiISAgJWjgzsxZm9pmZbTKzDWZ2v+/4g2aWbmZpvp+hp53zCzPbbmZbzOzqYNUmIiIiEqoigjh2EfAT59wqM6sLrDSzj3yfPe6ce/T0L5tZJ2A80BmIAz42s/bOueIg1igiIiISUoI2c+acy3TOrfK9PglsAuLPccpI4E3nXL5zbhewHegTrPpEREREQlG59JyZWSugB7Dcd+g+M1trZi+YWUPfsXhg72mn7ePcYU5ERESkygl6ODOzOsAM4IfOuRPAM0AbIBnIBP7+9VfPcro7y3iTzSzVzFKzsrKCVLWIiIhIxQhqODOzSLxg9ppzbiaAc+6Ac67YOVcCPMd/b13uA1qcdnoCkHHmmM65Z51zKc65lJiYmGCWLyIiIlLugrla04DngU3OucdOOx572tdGA+t9r+cC480sysySgHbAimDVJyIiIhKKgrlacyBwM7DOzNJ8x34JfM/MkvFuWe4G7gRwzm0ws2nARryVnvdqpaaIiIhUN+bcN9q6Ko2UlBSXmppa0WWIiIiI+ONs/fXfoCcEiIiIiIQQhTMRERGREKJwJiIiIhJCFM5EREREQojCmYiIiEgIUTgTERERCSEKZyIiIiIhROFMREREJIQonImIiIiEEIUzERERkRCicCYiIiISQhTOREREqrKSEji6u6KrkPOgcCYiIlKVrX4ZpvSBzx+F4sKKrkb8oHAmIiJSlbUfAhcMgU//CP+6GPatrOiKpBQKZyIiIlVZ3WYw7iUY/zqcOgrPXwHzfwH52RVdmXwLhTMREZHqoMM1cO9ySJkIy56Gp/vDto8ruio5C4UzERGR6iK6Hlzzd5j4AUTWhNfGwow7IOdQRVcmp1E4ExERqW4S+8Fdi+DiB2DDLJjSG9a8Cc5VdGWCwpmIiEj1FBEFl/7CC2mN28KsO+HVMdp2IwQonImIiFRnTTt6tzmHPgp7V3i9aEunQHFRRVdWbSmciYiIVHdhYdDnDm/BQNJF8OGvvFWd+9dVdGXVksKZiIiIeOonwPfehGtfhOP7vH3RPn4QCk9VdGXVisKZiIhIFbb9YDb3vr6K/cfz/DvBDLqMgXtXQPL3YPHj8MwA2PV5cAuV/1A4ExERqcI2ZZ7g440HuPKxhby6bA8lJX6uyKzVCEY+BRPmeqs4XxoOc+7zNrKVoDJXiZfNpqSkuNTU1IouQ0REJKTtPpTDL2etY+mOw/Ru1ZA/j+lG26Z1/B+g8BQseASW/h/UagxD/wqdRnmzbHI+/PqDKZyJiIhUA845pq/cx0PvbuJUQTH3XdaWuy5uQ42I87iJlrkG5v4AMtPggqHeCs/68cEruupROBMREZH/lXUyn9+/s4F5azNp36wOj4ztRs/Ehv4PUFwEy5+BTx+GsAi44neQcru34lNKo3AmIiIiZ/fJpgP8evZ69p/I45b+rfjp1RdQJyrC/wGO7IJ5P4Kdn0GLvjD8SWjaIXgFVw0KZyIiIvLtsvOL+Nv8zby8bA+x9aJ5aHQXLuvQzP8BnPMe+/TBLyA/Gy76KQz6kff0ATkbhTMREREp3co9R3lgxlq2HcxmRPc4fju8E03qnEfAys7yAtq6tyGmgzeLltg3eAVXXgpnIiIi4p+CohKeWbCDpz7bTq2ocH59TSfG9ozHzmdF5raPvFudx/dB70lw+W8hul7wiq58FM5ERETk/Gw/eJIHZqwjdc9RBrVtwp9GdyWxcS3/B8jPhk8fguX/hLqxMOwxuGBI8AquXBTORERE5PyVlDheX/EVj7y/maKSEn58ZXsmDkwiIvw8VmTuS4W534eDG6HzaBj8F6h7Hv1sVZPCmYiIiHx3+4/n8Zs56/lo4wG6xNfjkTHd6BJf3/8Bigpg6T9g4d8gMhquehh63FSdN6/16xcP2qYkZtbCzD4zs01mtsHM7vcdb2RmH5nZNt8/G/qOm5k9aWbbzWytmfUMVm0iIiJSuub1o3n25l48c2NPDpzIZ+RTS/jz+94mtn6JqAEX/QzuXgLNusDc++DlEXB4R3ALr+SCuWNcEfAT51xHoB9wr5l1Ah4APnHOtQM+8b0HGAK08/1MBp4JYm0iIiLiBzNjSNdYPv7RxYxLSeBfC3dy9ROfs2T7If8HadIObpkHw56AjDXeg9QXPw7FhcErvBILWjhzzmU651b5Xp8ENgHxwEjgJd/XXgJG+V6PBF52nmVAAzOLDVZ9IiIi4r/6tSL585huvHFHP8LDjBunLudnb6/hWG6BfwOEhUHKbXDvcmh3JXz8IDx3KWSsDmrdlVG5PGvBzFoBPYDlQDPnXCZ4AQ5o6vtaPLD3tNP2+Y6JiIhIiOjfpjHv338h917ahlmr07nisYW8syYDv3vY68XC9a96P9lZ8Nxl8MGvoCAnuIVXIkEPZ2ZWB5gB/NA5d+JcXz3LsW/8L21mk80s1cxSs7KyAlWmiIiI+Ck6MpyfXd2BufcNIr5BTb7/xmpufymV9GOn/B+k43BvFq3nLfDFFHi6P2z/JHhFVyJBDWdmFokXzF5zzs30HT7w9e1K3z8P+o7vA1qcdnoCkHHmmM65Z51zKc65lJiYmOAVLyIiIufUKa4eM+8ZyG+GdeKLHYe56rGFvLR0N8Ulfs6i1WwAw5+AW9+D8Brw6hiYdRfkHglu4SEumKs1DXge2OSce+y0j+YCt/he3wLMOe34BN+qzX7A8a9vf4qIiEhoCg8zbh+UxIc/uoherRrxu7kbuPafS9l64KT/g7QaCHct9lZ2rnsbpqTA2re9Z3dWQ0Hb58zMBgGLgHVAie/wL/H6zqYBicBXwHXOuSO+MDcFGAzkArc55865iZn2ORMREQkdzjnmpGXw+3c2kJ1fxN2XtOXeS9sQFRHu/yAHNsDcH0B6KrS90nvCQIPE4BVdvrQJrYiIiJS/w9n5PPTuJmatTqdNTG0eGduN3q0a+T9ASTGseA4++YP3/vLfQJ/JEHYeIS80KZyJiIhIxVm4NYtfzlxH+rFT3NQvkf83uAP1oiP9H+DYV/DuT2DbhxDfC0b8HzTrHLyCg0/hTERERCpWTn4Rj320lReX7CKmbhR/HNmFqzo3938A52D9DHj/55B3DAb+0OtNi4wOXtHBo3AmIiIioWHN3mP8fMZaNu8/ydCuzXlweGea1juPgJV7xNsPbc3r0LgtDP8HtBoUvIKDQ+FMREREQkdhcQnPLdrJEx9vIyoijF8N7cj1vVtg5/Mg9B2fwjs/hGN7oNetcMXvvS05KgeFMxEREQk9uw7l8IuZa1m28wh9kxrx5zFdaR1Tx/8BCnJgwZ/hi6egdlMY+jfoNCJ4BQeOwpmIiIiEJucc01L38vC7m8grKuH+y9sx+aLWRIafxxasGath7vdh/zroMAyGPuo9Hip0KZyJiIhIaDt4Mo/fz93Iu+sy6dC8Ln8Z243uLc7jNmVxoTeDtuDP3lMGrvw99LzVe9B66FE4ExERkcrhww37+e2cDRw8mcetA5L4yVXtqR0V4f8Ah3fAvB/Crs8hcYC3YCCmffAK/m4UzkRERKTyOJlXyF/nb+GVZXuIb1CTh0d34ZILmvo/gHOQ9pq3qrMwFy76fzDwfoioEbyiz4/CmYiIiFQ+X+4+wgMz1rIjK4dRyXH8ZlgnGteJ8n+A7IPevmgbZkLTTjD8SWjRO3gF+0/hTERERCqn/KJinv5sB08v2E6dqAh+O7wTo5Ljz2/bjS3z4d0fw4kM6HsnXPZriKobvKJLp3AmIiIildvWAyd5YMZaVn11jIvax/DwqC60aFTL/wHyT3rP6FzxHNSL9x6k3v7q4BV8bgpnIiIiUvmVlDheXb6Hv7y/mRIHP7mqPbcNTCI87Dxm0fau8LbdyNoMXa6FwY9AnZjgFX12CmciIiJSdWQcO8VvZq/nk80H6ZZQn0fGdKNTXD3/BygqgMWPw6JHoUZtuPpP0P17cD63SstG4UxERESqFucc767L5MG5GziaW8idF7XmB5e3Izoy3P9BsrbA3B/A3mXQ+hIY9gQ0SgpWyadTOBMREZGq6VhuAX96bxPTUvfRqnEt/jSmKwPaNPF/gJISWPkCfPQglBTBVX+EPncErV4fv8JZSG6fKyIiInIuDWrV4K/Xdue1SX1xwA3PLeeBGWs5nlvo3wBhYdB7Ety7HNpcVp63NkulmTMRERGp1E4VFPOPT7bx3KKdNKxVgz+M7MyQLs3933bj6ywU/ICmmTMRERGp+mrWCOeBIR2Yc+9AmteP4p7XVnHHyyvJPH7KvwHMQmrmTOFMREREqoQu8fWZfc9AfjW0I4u3Z3HlY5/zyrI9lJRUrruECmciIiJSZUSEh3HHRa358IcX0yOxAb+ZvZ5x//qC7QdPVnRpflM4ExERkSonsXEtXp7Yh79f153tWdkM/cdinvh4K/lFxRVdWqkUzkRERKRKMjPG9krg4x9fzJCuzXni420Me3IxK/ccrejSzknhTERERKq0JnWi+Mf4Hrx4a29yC4q59p9L+e2c9ZzM83PbjXKmcCYiIiLVwqUdmvLhjy7i1gGteGXZHq56/HM+2XSgosv6BoUzERERqTZqR0Xwu+GdmXn3AOpFR3L7S6nc+/oqsk7mV3Rp/6FwJiIiItVOj8SGvPP9Qfz0qvZ8tOEAf3pvU0WX9B8RFV2AiIiISEWoERHGfZe1Y0jXWGrXCJ1IFDqViIiIiFSANjF1KrqE/6HbmiIiIiIhROFMREREJIQonImIiIiEEIUzERERkRCicCYiIiISQoIWzszsBTM7aGbrTzv2oJmlm1ma72foaZ/9wsy2m9kWM7s6WHWJiIiIhLJgzpz9Gxh8luOPO+eSfT/vAZhZJ2A80Nl3ztNmFh7E2kRERERCUtDCmXPuc+CIn18fCbzpnMt3zu0CtgN9glWbiIiISKiqiJ6z+8xsre+2Z0PfsXhg72nf2ec7JiIiIlKtlHc4ewZoAyQDmcDffcftLN91ZxvAzCabWaqZpWZlZQWnShEREZEKUq6Pb3LOHfj6tZk9B8zzvd0HtDjtqwlAxreM8SzwrG+MLDPbE5xq/0cT4FA5XKey0t/n3PT3KZ3+Ruemv0/p9Dc6N/19Slcef6P5zrmz9eP/j3INZ2YW65zL9L0dDXy9knMu8LqZPQbEAe2AFaWN55yLCUqhZzCzVOdcSnlcqzLS3+fc9Pcpnf5G56a/T+n0Nzo3/X1KF0p/o6CFMzN7A7gEaGJm+4DfAZeYWTLeLcvdwJ0AzrkNZjYN+P/t3X2sn+Mdx/H3h9aq9bRgwthajFXIqIeZWj1OWAWRydbNbItkmxktEZtlEZPIZtj6z2KxNmOph5US2XRW2Y6GBaX1cEqpoTF0YUE9l+pnf9zX0V+P89DiuK+Tfl7Jye/+3U/ne9/n5He+57qu+/o+AqwCTrf97lDFFhEREVGrIUvObE/pY/XMAfa/CLhoqOKJiIiIGA5SIWDdXNF2AJXL/RlY7s/gco8GlvszuNyjgeX+DK6aeyS7z4ciIyIiIqIFaTmLiIiIqEiSs35I2klSl6Qlkh6WNLXtmGojaZSkBZIeLPfoF23HVCNJG0u6X9JfB997wyNpmaTuUm/3vrbjqY2krSTdIOnR8nn0pbZjqomk3TvqNT8g6RVJ09qOqyaSziqf0YslXStpVNsx1UTS1HJvHq7ldyfdmv2QtD2wve1FkjYHFgIn2H6k5dCqIUnAGNuvSRoJ3AlMtX13y6FVRdLZwH7AFraPbTue2khaBuxnO3Mw9UHSVcAdtmdI2gQYbfvltuOqUanJ/CzwRdsfxxyY1ZP0aZrP5j1sv1lmRphr+8p2I6uDpD2B62hKRr4N3AqcZvvxNuNKy1k/bC+3vagsvwosISWl1uLGa+XtyPKVbL+DpB2BycCMtmOJ4UfSFsAkypPutt9OYjagI4Ankpi9zwhgU0kjgNH0M8n7Bmo8cLftN2yvAubTzMPaqiRn60DSWGAf4J52I6lP6bJ7AHgeuM127tHapgPnAqvbDqRiBuZJWijp+20HU5mdgReAP5au8RmSxrQdVMW+AVzbdhA1sf0scCnwNE3ZxBW257UbVVUWA5MkbS1pNPBV1q5Y1IokZ4OQtBkwB5hm+5W246mN7Xdt701TcuuA0kQcgKRjgedtL2w7lspNtD0BOAY4XdKktgOqyAhgAnC57X2A14GfthtSnUqX73HA9W3HUhNJnwSOB8bRVOAZI+nkdqOqh+0lwMXAbTRdmg/STIbfqiRnAyjjqOYAV9u+se14ala6Wm4HBq0ZtgGZCBxXxlRdBxwuaVa7IdXH9nPl9XngJpqxH9F4Bnimo0X6BppkLd7vGGBRZw3nAOBI4CnbL9h+B7gROKjlmKpie6btCbYnAS8CrY43gyRn/SqD3WcCS2z/pu14aiRpW0lbleVNaT4EHm03qnrYPs/2jrbH0nS3/NN2/mPtIGlMeeCG0l13FGtq7m7wbP8X+I+k3cuqI2jK3MX7TSFdmn15GjhQ0ujyd+0ImjHUUUj6VHn9DHAiFfwefayFz4eZicC3ge4ypgrgZ7bnthhTbbYHripPSG0EzLad6SJifWwH3NT8zWAEcI3tW9sNqTpnAFeXbrsnge+1HE91ylihr1DqNccatu+RdAOwiKa77n4qmgm/EnMkbQ28Q1Pb+6W2A8pUGhEREREVSbdmREREREWSnEVERERUJMlZREREREWSnEVERERUJMlZREREREWSnEXEkL/h9hsAAAPsSURBVJBkSZd1vD9H0gUf0bmvlPS1j+Jcg3yfkyQtkdTVx7bdJM2V9O+yz2xJ20k6VNIHmlJG0rQyLUREbMCSnEXEUFkJnChpm7YD6VTm5VtXpwI/sn1Yr3OMAm6hKau0q+3xwOXAth8yvGk0hanX2XpeT0QMA0nOImKorKKZ7PKs3ht6t3xJeq28HippfmmFWirpV5K+JWmBpG5Ju3Sc5khJd5T9ji3HbyzpEkn3SnpI0g86ztsl6Rqgu494ppTzL5Z0cVl3PnAw8HtJl/Q65JvAXbb/0rPCdpfttaobSLpA0jkd7xdLGlsqI9wi6cGy7uuSzqSpfdjV01In6ShJd0laJOn6UusXScsknS/pTuAkSWdKeqRc83WD/FwionKpEBARQ+l3wEOSfr0ex3wBGE9T4+5JYIbtAyRNpZktf1rZbyxwCLALTUKzK3AKsML2/pI+AfxL0ryy/wHAnraf6vxmknagKXy8L/ASME/SCbYvlHQ4cI7t+3rFuCfwYQraHw08Z3tyiWFL2ysknQ0cZvt/pcXx58CRtl+X9BPgbODCco63bB9cjn8OGGd7ZU9JtYgYvtJyFhFDxvYrwJ+AM9fjsHttL7e9EngC6EmuumkSsh6zba+2/ThNEvd5mtqcp5SSa/cAWwOfK/sv6J2YFfsDt5fC0KuAq4FJ6xHvB9FN0/J3saQv217Rxz4HAnvQJJgPAN8BPtux/c8dyw/RlHg6mabFMiKGsSRnETHUptOM3RrTsW4V5fOnFGPepGPbyo7l1R3vV7N2a3/v2nMGBJxhe+/yNc52T3L3ej/xaV0vpMPDNC1tg3nvOotRALaXluO7gV+WLtS+4rqt41r2sH1qx/bO65lM00q5L7BQUnpFIoaxJGcRMaRsvwjMpknQeixjTXJzPDDyA5z6JEkblXFoOwOPAX8HTpM0Et57onLMQCehaWE7RNI2ZXD9FGD+IMdcAxwkaXLPCklHS9qr137LgAll+wRgXFneAXjD9izg0p59gFeBzcvy3cDE0l2LpNGSdusdiKSNgJ1sdwHnAlsBmw0Sf0RULP9dRcTH4TLgxx3v/wDcLGkB8A/6b9UayGM0SdR2wA9tvyVpBk3X56LSIvcCcMJAJ7G9XNJ5QBdNa9Vc2zcPcsyb5SGE6ZKmA+/QdC1OpelK7TGHNd2s9wJLy/q9gEskrS7HnlbWXwH8TdJy24dJ+i5wbRk/B80YtKWsbWNglqQtS/y/tf3yQPFHRN1k9+4ZiIiIiIi2pFszIiIioiJJziIiIiIqkuQsIiIioiJJziIiIiIqkuQsIiIioiJJziIiIiIqkuQsIiIioiJJziIiIiIq8n9wLgFpJyBb4AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#mmm,neither aic nor bic works in practice\n",
    "#apparently domain knowledge is the key to success\n",
    "#otherwise machine learning is merely a state of art\n",
    "#you can try gap statistic or silhouette score\n",
    "#but there is no guarantee they give you a clear signal\n",
    "#details of gap statistic or silhouette score can be found in the link below\n",
    "# https://github.com/je-suis-tm/machine-learning/blob/master/k%20means.ipynb\n",
    "fig=plt.figure(figsize=(10,5))\n",
    "ax=fig.add_subplot(111)\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "plt.plot(range(2,10),AIC)\n",
    "plt.plot(range(2,10),BIC)\n",
    "plt.xticks(range(2,10))\n",
    "plt.title('AIC vs BIC')\n",
    "plt.xlabel('Number of Clusters')\n",
    "plt.ylabel('AIC/BIC')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}