factor analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Factor Analysis\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Factor Analysis is another model solved by Expectation Maximization algorithm (we already have Gaussian Mixture Model and Naive Bayes Mixture Model). It is related to Principal Component Analysis as it is a generalized version of Probabilistic PCA. Although both PCA and FA work on dimension reduction, FA has the capability to detect the latent variable structure from the perspective of distribution.\n",
    "\n",
    "Unfortunately, there aren't many towardsdatascience articles on FA due to its complexity in derivation. You can read the paper by Ghahramani and Hinton in 1996 as an alternative material for Stanford University (though they offer a slightly different version of derivation).\n",
    "\n",
    "FA in sklearn uses a different approach called singular value decomposition. EM is slower than SVD due to the increasing computation in log likelihood (Barber 2012). For those of you who are interested in SVD implementation of FA, you can refer to the free textbook hosted on the server of UCL.\n",
    "\n",
    "http://web4.cs.ucl.ac.uk/staff/D.Barber/pmwiki/pmwiki.php?n=Brml.Online\n",
    "\n",
    "Alternative FA material\n",
    "\n",
    "http://www.cs.utoronto.ca/~hinton/absps/tr-96-1.pdf\n",
    "\n",
    "Reference to PCA\n",
    "\n",
    "https://github.com/je-suis-tm/machine-learning/blob/master/principal%20component%20analysis.ipynb\n",
    "    "
   ]
  },
  {
   "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.decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#official sklearn factor analysis\n",
    "def skl_fa(data,num_of_latents):\n",
    "    \n",
    "    clf=sklearn.decomposition.FactorAnalysis(num_of_latents)    \n",
    "    clf.fit(data)\n",
    "       \n",
    "    return clf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pdf of multivariate gaussian\n",
    "#equivalent to scipy.stats.multivariate_normal(miu,sigma).pdf(data)\n",
    "def multivariate_gaussian(data,miu,sigma):\n",
    "    \n",
    "    numerator=np.exp(-0.5*(data-miu)@np.linalg.inv(sigma)@(data-miu).T)\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": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#e step\n",
    "#use lambda_ and psi to update conditional miu and sigma\n",
    "def e_step(data,miu,lambda_,psi):\n",
    "    \n",
    "    #denote beta\n",
    "    beta=lambda_.T@np.linalg.inv(lambda_@lambda_.T+psi)\n",
    "\n",
    "    #compute conditional params\n",
    "    miu_z_given_x=beta@(data-miu).T\n",
    "    sigma_z_given_x=np.identity(lambda_.shape[1])-beta@lambda_\n",
    "    \n",
    "    return miu_z_given_x,sigma_z_given_x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#m step\n",
    "#use conditional miu and sigma to update lambda_,psi\n",
    "def m_step(data,miu,miu_z_given_x,sigma_z_given_x):\n",
    "\n",
    "    #equation 8 in stanford material\n",
    "    lambda_=(data-miu).T@miu_z_given_x.T@np.linalg.inv(miu_z_given_x@miu_z_given_x.T+sigma_z_given_x)\n",
    "\n",
    "    #last equation\n",
    "    phi=(data.T@data-2*data.T@miu_z_given_x.T@lambda_.T+ \\\n",
    "         lambda_@(miu_z_given_x@miu_z_given_x.T+sigma_z_given_x)@lambda_.T)/data.shape[0]\n",
    "    \n",
    "    #take diagonal\n",
    "    psi=np.diag(np.diag(phi))\n",
    "    \n",
    "    return lambda_,psi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get latent variables z\n",
    "#andrew ng didnt teach us how to get the latent\n",
    "#we have to borrow ideas from sklearn even though we dont know why\n",
    "#alternatively\n",
    "#x=miu+lambda_@z+epsilon\n",
    "#z=inverse(lambda_)@(x-miu-epsilon)\n",
    "#but this is incorrect\n",
    "#as we dont know epsilon\n",
    "#we can only assume its zero\n",
    "#z=np.linalg.pinv(lambda_)@(data-miu).T\n",
    "def get_latent(data,miu,lambda_,psi):\n",
    "    \n",
    "    weighted_psi=lambda_.T@np.linalg.inv(psi)\n",
    "    cov_z=np.linalg.inv(np.identity(lambda_.shape[1])+weighted_psi@lambda_)\n",
    "    z=(data-miu)@(weighted_psi).T@cov_z\n",
    "    \n",
    "    return z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#unlike gmm, only p(x|z) depends on the params\n",
    "#thus the computation can be simplified\n",
    "#although we have one more step of computing z\n",
    "def get_lower_bound(data,miu,lambda_,psi):\n",
    "    \n",
    "    #transform x\n",
    "    z=get_latent(data,miu,lambda_,psi)\n",
    "\n",
    "    #logarithm multivariate gaussian\n",
    "    logpdf=np.log(multivariate_gaussian(data-(lambda_@z.T).T,miu,psi))    \n",
    "    lower_bound=sum(logpdf)\n",
    "    \n",
    "    return lower_bound"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#using mle for training\n",
    "def training(data,num_of_latents,\n",
    "             tolerance=0.001,num_of_itr=50,diagnosis=False):\n",
    "        \n",
    "    #get dimensions\n",
    "    num_of_features=data.shape[1]\n",
    "    \n",
    "    #initialize\n",
    "    lower_bound_old=None\n",
    "    lower_bound=None\n",
    "    counter=0\n",
    "    \n",
    "    #miu is the unconditional mean\n",
    "    miu=data.mean(axis=0)    \n",
    "    \n",
    "    #lambda_ is the coefficient of the latent variables\n",
    "    #use principal components to initialize\n",
    "    #this approach improves the performance (Barber 2012)\n",
    "    _,eigvecs=np.linalg.eig((data-miu).T@(data-miu)/data.shape[0])\n",
    "    lambda_=eigvecs[:,:num_of_latents]\n",
    "\n",
    "    #psi is the covariance matrix of noise, a diagonal matrix\n",
    "    #use the diagonal of covariance matrix to initialize\n",
    "    psi=np.diag(np.diag((data-miu).T@(data-miu)/data.shape[0]))  \n",
    "\n",
    "    #cap the maximum number of iterations\n",
    "    while counter<num_of_itr:\n",
    "        \n",
    "        #e step\n",
    "        miu_z_given_x,sigma_z_given_x=e_step(data,miu,lambda_,psi)\n",
    "            \n",
    "        #m step\n",
    "        lambda_,psi=m_step(data,miu,\n",
    "                           miu_z_given_x,\n",
    "                           sigma_z_given_x)\n",
    "            \n",
    "        #use lower bound to determine if converged\n",
    "        lower_bound_old=lower_bound\n",
    "        lower_bound=get_lower_bound(data,miu,lambda_,psi)\n",
    "            \n",
    "        if lower_bound_old and np.abs(lower_bound/lower_bound_old-1)<tolerance:\n",
    "            if diagnosis:\n",
    "                print(f'{counter} iterations to reach convergence\\n')\n",
    "            return miu_z_given_x,sigma_z_given_x,lambda_,psi\n",
    "            \n",
    "        counter+=1\n",
    "    \n",
    "    if diagnosis:\n",
    "        print(f'{counter} iterations to reach convergence\\n')\n",
    "        \n",
    "    return miu_z_given_x,sigma_z_given_x,lambda_,psi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize number of latent variables\n",
    "num_of_latents=2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load data\n",
    "iris=sklearn.datasets.load_iris()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#raise error when zero is encountered in logarithm\n",
    "np.seterr(divide='raise')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3 iterations to reach convergence\n",
      "\n"
     ]
    }
   ],
   "source": [
    "miu_z_given_x,sigma_z_given_x,lambda_,psi=training(iris.data,\n",
    "                                                   num_of_latents,\n",
    "                                                   diagnosis=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf=skl_fa(iris.data,num_of_latents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.51538355 -1.81399954]\n",
      " [-0.59921545 -3.44747338]\n",
      " [ 6.20068532  0.83567308]\n",
      " [ 2.70348799 -0.22529923]]\n",
      "[[ 0.72746132 -0.36874572]\n",
      " [-0.18085204 -0.24059   ]\n",
      " [ 1.75790224  0.03551754]\n",
      " [ 0.73100458  0.04526567]]\n",
      "[[34.3143796   0.          0.          0.        ]\n",
      " [ 0.          9.41287335  0.          0.        ]\n",
      " [ 0.          0.         14.24159276  0.        ]\n",
      " [ 0.          0.          0.          1.45144496]]\n",
      "[[0.01609026 0.         0.         0.        ]\n",
      " [0.         0.09807659 0.         0.        ]\n",
      " [0.         0.         0.00405033 0.        ]\n",
      " [0.         0.         0.         0.04071555]]\n"
     ]
    }
   ],
   "source": [
    "#both lambda and psi seem to be different from the result of sklearn\n",
    "print(lambda_)\n",
    "print(clf.components_.T)\n",
    "print(psi)\n",
    "print(np.diag(clf.noise_variance_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEWCAYAAABBvWFzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXucVWW9/98fYJAJlVFEHUBCk8hMFB3vaV7BDimUeCnr0NF+nn7liaO/fIXlsUntHJRTRr86v5MnLay8H0M6Y4rXTMsUBPEWQWg6AyqpkMIgt+/vj7X2sGfPWnuvPXvty8z+vl+v9dp7Pet5nvXdm2F99/N8bzIzHMdxHKdUBlRbAMdxHKd/4ArFcRzHSQVXKI7jOE4quEJxHMdxUsEViuM4jpMKrlAcx3GcVHCF4jg1gKR2SSeUMH6gpHcljSly3BckPdLb+zpONq5QnKoj6WVJneEDMXOMLGG+/SWVJcAqM7ek75dj/t5iZtvMbGcze6VaMkh6TNKmnH/HX4bXTgm/t9tzxrSE7Q9UR2onTVyhOLXC6eEDMXOsrpYgkgbluTwDeAv4tKSGConUl/hizr/jJ7OuvQ58TFJTVtvfA3+qrIhOuXCF4tQskgZIulPSa5LWSXpE0gFZ198n6TpJr0haL+lRSTsBj4bXM7+SDw/nukLSXyS9IemnknYN+2VWHf8g6RVgYYw8Aj4HXAYImJJ1bVA4xz9KWinp7exVjKRxkh6W9Kakv0r6maRhEfcYJWlj9kNX0pHhdzBI0gfDz7k+nOfmnPuPDc8/IelFSe+E22kX5/mqB0j6j3DOFyWdGM7xaUl/yJHva5LuzDNXPjYBvwLOycgMTAdu7uV8To3hCsWpdf4HGAfsDTwH/Czr2nXABOBIYHfg68B24HiArF/JTwFfAD4LnAB8ANgNmJtzr+OBD5GlKHI4AdgLuA24g+DXdS5/BxwGTAQ+K+mUsF3A1UAz8GFgP+BfcgebWQfwGHBWVvNngVvMbCvwbaAtlH808MMYWX8CXGBmuxB8R7+J6QdwDPBHYA/gKuCXoUKbD4yXNC5Hlp/1nCIxN7Hje/s7YAnBysXpD5iZH35U9QBeBt4F1oXH/Jh+ewAGDAUGAu8BB0b02z/40+7W9hvgwqzzA8PxAzL9gTEF5PwpcGf4/rhw/PDwfFA4x1FZ/e8Cvhoz13TgqazzduCE8P15wG+y5n0DODQ8vxn4f8ConPky9x8bnq8mUKK7FPhMXwBeBZTV9jTw6fD9fwHfCt8fAvwVaIiZ6zFgY9a/4zrgm+G1U8J/ZwGrCJT6nQSrlS8CD1T779CP0g9foTi1wjQzawqPadDluXStpFWS/gasDPvuQbBSGAz8OeH8I4G/ZJ3/JRw/Iqvt1bjBkoYCZwK/CJseA9YAn87p+lrW+43AzuH4vSXdLqkj/Cw/DT9HFL8EDg49tk4D1prZ0+G1/wM0AIskPStpRswcnwTOAF4JtwqPjPtsQLuFT/2QvxB8XwDzCBQcBKuT28xsS565vpT179hkZt/Kvhje5+fATAKlfHeeuZw+hisUp5b5e4JtkZOAYQQrCQh+5b4ObCb4pZtLlIfXauD9WedjwvFruwZ1f6jmciaBcrhe0msEymRvore9oriGYEVzkJntCnw+/Bw9hTfbCPw3wYP8c2RtMZnZGjP7gpk1A18O5dk3Yo4/mNkZwJ4E24a35pFtdM75GILvCzN7DEDSsQTKs5Ttrgw3Eci+wMw2pTCfUyO4QnFqmV0IHsJvAu8jsB8AgZsswa/874W//gdKOjb0vHoDMEn7Zc11C3CJpLGSdgnnusXMtieUZQbB9s9BBFs/hxDYXFqyHQUKfJYNwHpJ+wBfLdD/JuB8AnvOzzONks6WNCo8XUegPLdlD5TUKOkzknYNVxPv5PbJoVnSRaFh/1wCJX1v1vWfEWyzbTCzJwp90EKY2UoCe9QVpc7l1BauUJxa5icEv5RXA88Dv8u5fjHwIrCYwJX3XwlsAe8A/wb8IfQOayFQBrcBvyXYw3+HYNulIOHW0wnA98zstazjSeABAmVTiG8CRwDrgQUEK5B8PEpgJ/qDmbVntR8JPCVpA4GN5ssWHXsyA/hLuL12AcFKJ47fEdiU3gJagTPN7O2s6zcBHyHZ6uQ/c+JQnozqZGa/NbM1CeZz+hDKv8p3HKdaSHoUuNHMflplOYYSrPo+YmYvVVMWp7bxFYrj1CCSjiJYFdxRbVkI7B2PuzJxCpEvIthxnCog6RcEtpN/MrMNVZalHdgCTK2mHE7fwLe8HMdxnFTwLS/HcRwnFepqy2uPPfawsWPHVlsMx3GcPsXixYv/amYjCvWrK4UyduxYFi1aVG0xHMdx+hSS/lK4V5W3vCSdJml5mJ11VsT1SyS9IGmZpAclvT/r2gxJK8IjSRyA4ziOU0aqplAkDSTIlPpxguyrn5b04ZxuS4AWM5tAkEju2nDs7gSBYkcSBIt9U9JulZLdcRzH6Uk1VyhHACvNbJWZbSbINdTNNdHMHg7zGgE8wY6cQ5OB+83srTCi936CJHqO4zhOlaimDWUU3bO7thOsOOK4APh1nrGjeowAJF0IXAgwZkxR5bYdx+lnbNmyhfb2djZt8pyUUQwZMoTRo0fT0NC7YqTVVChRmVYjg2IkfRZoAT5W7Fgzux64HqClpcWDbhynjmlvb2eXXXZh7NixBAU4nQxmxptvvkl7ezv77tsjgXUiqrnl1Q7sk3U+mjBldjZhxbtvAGeY2XvFjHUcx8lm06ZNDB8+3JVJBJIYPnx4Sau3aiqUp4BxkvaVNBg4lyALaxeSJgI/IlAmb2Rdug+YJGm30Bg/KWxzHMfJiyuTeEr9bqq25WVmWyVdRKAIBhJkVX1e0pXAIjNbAMwhKGp0R/hBXzGzM8zsLUlXESglgCvN7K0qfAzHcRwnpKpxKGZ2j5l90Mw+YGbfDtuuCJUJZnaKme1lZoeExxlZY280s/3D4yfV+gxOP2XZ7XDdR6C1KXhddnu1JXL6CTvvvHPivq2trfz7v/97KvPfe++9jB8/nv3335/Zs2cXNWdSPJeX4+Sy7Hb41Vdg/auABa+/+oorFafPsm3bNr785S/z61//mhdeeIFbbrmFF154IfX7uEJxnFwevBK2dHZv29IZtDt1xfwlHRw7+yH2ndXGsbMfYv6SjrLc51e/+hVHHnkkEydO5JRTTuH111/vuvbMM89w0kknMW7cOP7rv/6rq33OnDkcfvjhTJgwgW9+85t553/yySfZf//92W+//Rg8eDDnnnsud999d+qfwxWK4+Syvr24dqdfMn9JB5fd9Swd6zoxoGNdJ5fd9WxZlMpHP/pRnnjiCZYsWcK5557Ltdde23Vt2bJltLW18fvf/54rr7yS1atXs3DhQlasWMGTTz7J0qVLWbx4MY8++mjs/B0dHeyzzw7H2NGjR9PRkf7nqKvkkI6TiGGjw+2uiHanbphz33I6t2zr1ta5ZRtz7lvOtImRcdS9pr29nXPOOYc1a9awefPmbnEgU6dOpbGxkcbGRk488USefPJJHnvsMRYuXMjEiRMBePfdd1mxYgXHH3985PxRda/K4e3mKxTHyeXkK6ChsXtbQ2PQ7tQNq9d1FtVeCv/0T//ERRddxLPPPsuPfvSjbrEguQ9+SZgZl112GUuXLmXp0qWsXLmSCy64IHb+0aNH8+qrO34ktbe3M3LkyNQ/hysUx8llwtlw+vdh2D6AgtfTvx+0O3XDyKbGotpLYf369YwaFax65s2b1+3a3XffzaZNm3jzzTd55JFHOPzww5k8eTI33ngj7777LhBsab3xxhs95s1w+OGHs2LFCl566SU2b97MrbfeyhlnnBHbv7f4lpfjRDHhbFcgdc6lk8dz2V3Pdtv2amwYyKWTx5c078aNGxk9esf26SWXXEJraytnnXUWo0aN4qijjuKll17qun7EEUcwZcoUXnnlFf7lX/6FkSNHMnLkSF588UWOPvpoIHAV/vnPf86ee+4Zec9Bgwbxgx/8gMmTJ7Nt2zbOP/98DjzwwJI+RxR1VVO+paXFvMCW49QvL774IgcccEDi/vOXdDDnvuWsXtfJyKZGLp08PnX7Sa0R9R1JWmxmLYXG+grFcRwnhmkTR/V7BZImbkNxHMdxUsEViuM4jpMKrlAcx3GcVHCF4jiO46SCKxTHcRwnFVyhOI7jVJBqpa8///zz2XPPPfnIRz5S1HzF4ArFcRynDvj85z/PvffeW9Z7VFWhSDpN0nJJKyXNirh+vKSnJW2VND3n2jZJS8NjQe5Yx3GckqlQobVyp68HOP7449l9993LIn+GqgU2ShoI/BA4FWgHnpK0wMyyq768Anwe+GrEFJ1mdkjZBXUcpz7JFFrL1MbJFFqD1NPyZNLXS+LHP/4x1157Ld/5zncCMZYt44knnmDDhg1MnDiRKVOm8Nxzz3WlrzczzjjjDB599NHYbMOVopqR8kcAK81sFYCkW4GpQJdCMbOXw2vbqyFgJanHFA+OU9PkK7SWskIpd/r6SlHNLa9RQHbRifawLSlDJC2S9ISkaXGdJF0Y9lu0du3a3spaVipZyMcpE16Dvv9RwUJr5U5fXymqqVCiqrsUk6lyTJis7DPA9yR9IKqTmV1vZi1m1jJixIjeyFl28hXycfoAXoO+fxJXUK0MhdbKnb6+UlRTobQD+2SdjwZWJx1sZqvD11XAI8DENIWrJJUs5OOUAa9B3z8pU6G1TPr6zPHd7363K339cccdxx577NGtfyZ9/VFHHdWVvn7SpEl85jOf4eijj+aggw5i+vTpvPPOO3nv++lPf5qjjz6a5cuXM3r0aG644YaSPkcU1bShPAWMk7Qv0AGcS7DaKIik3YCNZvaepD2AY4FrCwyrWUY2NdIRoTzKUcjHKQNeg75/krGTPHhl8G85bHSgTEq0n2zfHm0Snjp1ao+21tbW2HlmzpzJzJkze7RnVi253HLLLckELIGqKRQz2yrpIuA+YCBwo5k9L+lKYJGZLZB0OPBLYDfgdEnfMrMDgQOAH4XG+gHA7BzvsD5FuQr5OEWw7PbePzjqrAZ926o25j49l9c2vMbeQ/dm5qEzmbLflGqLVR680FpRVLUeipndA9yT03ZF1vunCLbCcsf9Djio7AJWiIw3l3t5VYlS3UNPvqL7eOi3NejbVrXR+rtWNm0LjMZrNqyh9XetAP1XqTiJ8QJbNYIX8qkipbqHlmlrpBaZ+/TcLmWSYdO2Tcx9eq4rFMcViuOkYgOpk62R1za8VlS7U194Li/HqaB7aF9n76F7F9Xu1BeuUBynTO6h/ZGZh85kyMAh3dqGDBzCzEN7ehs59YcrFMeZcDac/n0Ytg+g4PX079fFFlaxTNlvCq3HtNI8tBkhmoc203pMq9tPiqAa6etfffVVTjzxRA444AAOPPBA5s6dW9ScSXEbiuOU4jJch0zZb4orkD7GoEGD+M53vsOhhx7KO++8w2GHHcapp57Khz/84VTv4ysUp77pC2lTKp0nzPOSddG2qo1Jd05iwrwJTLpzEm2r2spyn3Knr29ububQQw8FYJddduGAAw6goyP9XIGuUJz6ppi0KdV40FZa4fUFBVshMjE3azaswbCumJtyKJVM+volS5Zw7rnncu21OxJ/LFu2jLa2Nn7/+99z5ZVXsnr1ahYuXNiVvn7p0qUsXryYRx99NNG9Xn75ZZYsWcKRRx6Z+ufwLS+nvknqMlzB2hjdqGAK9arcr4apZMxNpdLXv/vuu5x55pl873vfY9ddd031M4CvUJx6J6nLcNKVTNqrmGJiZNK4t+cl66KSMTeVSF+/ZcsWzjzzTM477zw+9alPpf4ZwBWKU+8kdRlO8qAtx3ZRUoWX1r09JqeLSsbclDt9vZlxwQUXcMABB3DJJZekLn8GVyhOfZPUZTjJg7YcaeyTKry07u0xOV2UK+amGunrH3/8cX72s5/x0EMPccghh3DIIYdwzz33xPbvLTIrpqZV36alpcUWLVpUbTGcvkiuDQWCB2228mltIrpGnKB1XWn3LuTWnOa9+7Eb9YsvvsgBBxyQuH9dZVYOifqOJC0OCxrmxY3yjpOEJAkgy5XGPkmesDTvXSd5yZLgMTfF4QrFcZJS6EFbzTT2dZRC36ldqmpDkXSapOWSVkqaFXH9eElPS9oqaXrOtRmSVoTHjMpJ7TgxVDOFi6ePSUw9bfMXS6nfTdVsKJIGAn8CTiWoL/8U8OnsyouSxgK7Al8FFpjZnWH77sAioIVg43gxcJiZvZ3vnm5DcSpGP7ZD9GVeeukldtllF4YPH97DHbfeMTPefPNN3nnnnW5xMNA3bChHACvNbBWApFuBqUCXQjGzl8NruUWYJwP3m9lb4fX7gdOA8hdNdpxCVCsI0inI6NGjaW9vZ+3atdUWpSYZMmQIo0f33uZXTYUyCsi2IrYDSXMBRI2NLHco6ULgQoAxY8YUL6XjFItHm9csDQ0NPX59O+lRTRtK1Hoz6f5b4rFmdr2ZtZhZy4gRIxIL5zi9po6izSuVPNHpG1RzhdIO7JN1PhpYXcTYE3LGPpKKVI5TKmm58C67HX79Neh8Kzhv3B0+fk3NrHIyyRMz+a4yyRMBd7WtU6q5QnkKGCdpX0mDgXOBBQnH3gdMkrSbpN2ASWGb45RGGvmw0og2X3Y73P3lHcoEgvfzv1QzmX/zJU906pOqKRQz2wpcRKAIXgRuN7PnJV0p6QwASYdLagfOAn4k6flw7FvAVQRK6SngyoyB3nF6TVr5sNJw4X3wSti2uWf79i2lpXJJkTUb1hTV7vR/qhrYaGb3APfktF2R9f4pgu2sqLE3AjeWVUCnvkjTmF5qtHk+e0uN2GIGaADbLdcBM2h36hP/l3ecDLVkTM9nb6mRzL9RyiRfu9P/cYXiOBlqKXX7yVfAwME92wc0VDedSpaNqXlbtFNm89DmCgvl1AquUBwnQ5QxHQW2lErXVp9wNkz9YeDZlaFxd5j2H9Xz8sqxMc18802GbO+uVNJI7+70XTw5pONk6JZR+FWCcKfwgVko2r0cqVZqLetvjo1pyoaNAMwdPpzXBg6om/TuTjxeD8VxorjuIzGxJPvAxc91b0tSKyUffSXvV7nqvTg1T9JcXr7l5ThRFGOgL6VaYjnKBpeLWrIxFcAj+KuDKxTHiaKYh2cp3mHlKBtcLvpIeeBMBP+aDWswrCuC35VK+XGF4jhRFPPwLOWXey25Kheij9Rc8Qj+6uFG+Qoxf0kHc+5bzup1nYxsauTSyeOZNjEyQbJTCyQp+ZuhlGqJ5SobXC5qzVEggtc2vFZUu5MevVIokp41s4PSFqa/Mn9JB5fd9SydW7YB0LGuk8vuehbAlUotk/ThWYzyycVL96bO3kP3jkz/svfQvasgTX0Rq1AkfSruEuD/MkUw577lXcokQ+eWbcy5b7krlP5Cb3+5F6uM+opHWJVoW9VG59bOHu0eH1MZ8q1QbgN+QbSf4JDyiNM/Wb2u5x94vnanzkiqjLwSZF5y0+lnGDZ4GJcdeZnHx1SAfAplGfDvZvZc7gVJp5RPpP7HyKZGOiKUx8im3Khsp64odrUR5xH266+5QiHaGA/wvob3uTKpEPm8vP4Z+FvMtU+WQZZ+y6WTx9PYMLBbW2PDQC6dPL5KEjlVpzfxJ3GeX51v1V7cShp1ZfIQFWfixvjqE6tQzOy3ZvZKzDUPNy+CaRNH8W+fOohRTY0IGNXUyL996iC3n9QzvYk/yef59euvpSNXGpQ5WDMuzmTXwbtG9ndjfOVwt+Ei6a3777SJo1yBODvoTfzJyVfAXf8r+lrnW9A6LIgNqbahPs26MhHExZkMGTSEIQOHdLvmxvjKUtXARkmnSVouaaWkWRHXd5J0W3j9D5LGhu1jJXVKWhoe/1kJeTPuvx3rOjF2uP/OX9JRids7/YneBENOOLt79uEoaiF1S5mDNeO2sNa/t57WY1ppHtqMEM1Dm2k9ptXtJxWkaisUSQOBHwKnAu3AU5IWmNkLWd0uAN42s/0lnQtcA5wTXvuzmR1SSZkLuf9GrV4y4zyg0elGb+NPPn5N/ColQ4qrgV5R5mDNfHEmU/ab4gqkihRcoUjaS9INkn4dnn9Y0gUp3PsIYKWZrTKzzcCtwNScPlOBeeH7O4GTJSmFe/eKfO6/UauXS+94hkvvfMZXNE5PepvGJMkqBaqbuqXMOb9mHjqTIQO7Ry741lZtkGSF8lPgJ8A3wvM/EcSo3FDivUcB2T9j2oEj4/qY2VZJ64Hh4bV9JS0h8ES73Mx+G3UTSRcCFwKMGTOmJIHzuf9GrV62bO8ZwuMBjU4XxQRDZrsYN+4WVG7cviW+fzVTt5SSOSABmRXI3Kfn8tqG17wOSw2RRKHsYWa3S7oMuh7s2woNSkDUSiP3CRzXZw0wxszelHQYMF/SgWbWw83ZzK4HroegHkopAl86eXy3FCqww/334tuWJp4ns9Lx/F5OInIDGjvfCsoDN+4evM8uBAa1kbqlzDm/fGurNklilN8gaTjhX6yko4D1Kdy7Hdgn63w0sDquj6RBwDDgLTN7z8zeBDCzxcCfgQ+mIFNe8rn/FhOkaMAh31ro22FOMqK8prZthsFDoXU9fOr6qmQArmbNkaufuJqDbzqYg+YdxME3HczVT1xdsXs78SRZoVwCLAA+IOlxYAQwPYV7PwWMk7Qv0AGcC3wmp88CYAbw+/CeD5mZSRpBoFi2SdoPGAesSkGmgsS5/0atXhoGCARbtvVcGK3r7Lld4dthTiSFvKaqkAE4N81JJhYEKPvK4eonrua25bd1nW+37V3nlx91eVnv7eSn4ArFzJ4GPgYcA/wjcKCZLSv1xma2FbgIuA94EbjdzJ6XdKWkM8JuNwDDJa0kUGwZ1+LjgWWSniEw1n/RzN4qVaZSiFq9zDnrYOZMP5hRRaxeOtZ1su+sNo6d/ZCvVpyAGqyUWM2aI3f86Y6i2p3KEVtTPk+2YQDM7K6ySFRGKlVTPso2cvFtSyOzbOajsWGgR9Q7pdesLwMT5k3AIv6ihVg2o+Tfm3k5aF585YxnZzxb1nvXK0lryufb8jo9zzUD+pxCqQRxtU+GNTZEbnPlw7fAnC4GNe5QKI27B/EoVYyGr2bNkQEawHbb3qNdiEl3TnLPryqSL5fXP+Q5zq+kkH2JuOBHKdplrRCe4r7OyaxOOrN2dDvfDoIb0066WERCx3LFgiQx9J/1wbMixwp5HfmQajlMFDTKhx5e3wQ+SrAyeQy4MuNl5XQnTgGs27il6C0v8BT3dU+Uh1fmL6mUeii5qfPHTYJnbk5ca6UcsSBJDf0Zw/sdf7qD7badARrATgN2onNb9+8pY9Opt1VKNR0mYm0oXR2k+4FHgZ+HTecBJ5hZn6uJUgkbyrGzH4oMfswY5qOuxeE2FIfWJqJr3GUxbB+4uEfZoniibDK5sSy9nbsEJt05KXIbrXloMwunL8w7tpo2nVqjlO8xjqQ2lCRxKLub2VVm9lJ4XA009UqqOiBf7ZNLJ4/Pu+3V1NjgKe6d7iTx5Co2zUq+VU+PuV8tSz2TKEqpZxJnu6nH1PXVrAuTRKE8LOlcSQPC42ygPjcmE5Av+HHaxFF5f2u2nnEgj886iZdmT+HxWSe5Mqk3omwYUXmxemDFPfSLVUAVymBcilLw/F47qKZyjVUokt6R9DeC2JObgc3hcStwcdkl68NMmziKx2edxHXnBMmQL75taVdcSVxMytDBA12B1CNdSmQY3HVhz6JUkJVEEmJdO3If+vkM7LGrnjzr50LFv1KgFKUwZb8pnro+pJrKtaANpT9RqTgU6Ok+DLG71EAQVT/nrINdqdQTkbaMHHJtGF3G9Ij08LDDpThf3ErUfQc0wKDBsHlDHoEFresSfbTe0raqrShDf9uqNv7tD//G+s1BNqimnZqYdcSsulQk2RT7PRYiqQ0lkUKRtBtBepMutWdmj/ZauipRSYUSZ5zPx6imRh6fdVKZJHJqjus+Eq8YulCQqys3c+9dFxL786QraWQO2copWzFpAETEdeQdnyK9ffi1rWrj8scuZ6tt7dbeMKCBq469qu6VSpqkZpSX9AUCL6/7gG+Fr62lCtjf6U38iMec1BlJbBmNu0XXZ2/cLX5MlDLJvd+Es3fYZ5IokzJlMI6rD58kbmLu03N7KBOALdu3lJQCphwxHNVMpFlJkiSHnAkcDjxhZidK+hCBYnHyEFc7pdAYp46Iq2yYIWOMj6rPPqgXfyu5tpNIb69clHo9k2zy5QQrtMLI57VUjEdT9gpp18G7snHrRraEtWbSiOGoZlxIpUni5bXJzDZBUOPdzP4IjC+vWH2fSyePp2Fg8tj4jGuxU0dEenCFfzOZNPSdb0eP7Xw7vnJj4+7xFROzjfWFttuG7RPYTC5+rmxpXsrhKlzoWja5K6T1m9d3KZMMm7Zt4uuPfb3Xq4tqJtKsNEkUSrukJmA+cL+ku+lZt8TJYdrEUQwdHL8AbBgomhobPOaknokqA/yp64MaJ5mHeL5Mwx+/JlpxfPya6PLC0H37LC+qSJGuUl2FB6nn/7GGAQ2JPZqiHvZRbLftvU7pUs24kEpTcMvLzD4Zvm2V9DBBkat7yypVP2F9nmSQc6a7R5dD4VomJ18R7bGVvQUVV2o3d97rPpJgiyuDVST55MxDZ3bbDoLiXIWBkry8evNQLzalSzUTaVaaWIUiaVcz+5uk7HV1Jjf0zkBV64/0BeIyDI9qanRl4iQjidJI+uAvJqBx2D6F+6RAqTnBSi0FHPewL0S2IirkpVaK0uxr5Fuh3Ax8AlhMsD5Wzut+pd5c0mnAXGAg8GMzm51zfSfgJuAw4E3gHDN7Obx2GXABsA34ipndV6o8aTJ/SQcbNvf0QGkYILeVOMWRVkXGWCeA6takj1IKacdRxBH1sB+kQew8eGfWv7ceSZGp8jOriyQG93Ik0qxVYhWKmX1CkoCPmdkrad9Y0kDgh8CpBLXjn5K0wMxeyOp2AfC2me0v6VzgGuAcSR8mKBl8IDASeEDSB82se974KjLnvuWRpX93HjLIVydOdYjbPjv4M7BiYfQKqArklvgtp1dUoYd9rsKA7quLpF470tI4AAAeSklEQVRqpa6k+gp5bShh/fZfEqwQ0uYIYKWZrQKQdCswFchWKFPZEfNyJ/CDUMlNBW41s/eAl8ISwUcQ1J6vCfKlsXecqjDhbHjlCVj8U7BtoIGBMvnEd6stWRdtq9q6KZMM5UxFn+9hX0jh1JPBPQlJ4lCekHS4mT2V8r1HAdnr73bgyLg+ZrZV0npgeNj+RM7YmvrZHxeH4rEmTtVYdntQ8ySzkLdtwfmYo6q6Iskmnytt7kO6Utti+RROPRnck5DEbfhE4PeS/ixpmaRnJaVRYCAqSCN3jyiuT5KxwQTShZIWSVq0du3aIkXsPVFp7BsGiI2bt7LvrLauZJHZzF/SwbGzH4q97ji9Ztnt8MsvRgdJljnpYzHk+2Wf/ZAuJcI+TTzLcXeSrFA+XqZ7twPZriSj6RnfkunTLmkQgcvyWwnHAmBm1wPXQ5DLKxXJE5Cxk8y5bzmr13UyrLGBDZu38na45ZWpNZ/pG1eLPnsux+kVmWSQcSbGYtPZl5F8XlfZD+lSIuzTJG5LDKjL+vZJ4lD+AiBpT7KSQ6bAU8A4SfsCHQRG9s/k9FkAzCCwjUwHHgrtOguAmyV9l8AoPw54MkXZUiFTAwWCZJG5LsSdW7Yx577lTJs4KrYWfea64/SaQilWkhTxKhOZbas1G9YwQAMiPaoAzhl/TrcHci3ZLnK3xOop1UouSZJDniFpBfAS8BvgZeDXpd7YzLYCFxEkm3wRuN3Mnpd0paQzwm43AMNDo/slwKxw7PPA7QQG/HuBL9eSh1cUcUb6jnWdeTMTe8JIp2TyrUAq7CKcTfa2FRCpTJqHNjP7uNlddeQz1HKFxnpKtZJLki2vq4CjgAfMbKKkE4FPp3FzM7sHuCen7Yqs95uAs2LGfhv4dhpyVIJ8ySI71nXG1kpxI75TMnHxJxq4o0ZKpVl2O3MXXcmmPPnu8tVAr+VgwVpaPVWaJEb5LWb2JjBA0gAzexg4pMxy9TuijPTZRHkaeMJIJxWiklA2NMIn/7NqyoRffYXXCjx9Cj2AhwzasQM/bPCwmqnQWMurp3KTRKGsk7QzQU2UX0iaC/QMAXfykqk1v9v7GmL7GDBQgVrxhJFOakQloazWygS6bDp7b82/Sx33AM5sla17b0f1yPe2vZeqiKVQz55fSba8pgKbCOrIn0fgaVU7foZ9jE1b8hcz2mbWtTJJS5nMX9LR5W02sqkx1bmdPkJa6VvSILTpzHx7Ha177M6mAT1/1+Z7ANeKh1cc9ZRqJZd8ySF/ANxsZr/Lap5XfpH6L1GeXFGk6d3l7sgOkFXyt8j0Krnjxk0qPU1LaNOZsmEjAHN3a2LNoIEMALZLNA9tzvsArnUbRdKAy0oFZlaSfCuUFcB3JDUDtwG3mNnSyojVPynGY6vYao9xuDuy0xWHknEdzpQRhvzKIGrcoht2XE86Ty5ZOcWmbNgYKJaGxsTbcMVGp1fywZ3UZbi/uhbH2lDMbK6ZHQ18jCCY8CeSXpR0haQPVkzCfsSwxnj7SS4ZW0qpxCkxd0euI6LiUJJEyCcpEdybSPsSbTrF2CjKHVGfWyt+9pOzE7kM91fX4qSBjdcA10iaCNwIfJMg5byTkLh09nFss3SC+j2nmBMbh1IoQj5pBH1vIu1LsOkUY6Mop70lapURR+52XK1v2/WWggpFUgNwGkEk+8kEwY3fKrNc/Y64dPZxjCrwwE9qaL908vhuNhRwd+S6Iy4OpVCEfGz9lCLnKQNJ08GX+uDOt12WtHww9NyO669JJWO3vCSdKulGgrxZFxIEIH7AzM4xs/mVErC/kG+LqWFA9+2tQg/8jKG9Y10nxg5De1QyyYy78qimRq9fX6/ExaEUipCPGpdLFSPtk1BKTEih7bKkSilqO66/uhbni0P5OkEOrQPM7HQz+4WZbaiQXP2OuC2mUU2NzDnr4KIe+PkM7VFMmziKx2edxEuzp/D4rJNcmdQbSW0Wy24P6s63NgWv0HNcywW1E8+SgFIe3IXsHHFKadjgYTQPbUYEHmtRAZdT9ptC6zGtBfv1NfJVbDyxkoL0d/JtPWUnkUyCG9qdoomzWXS5Bb9Kt1LAGQ+u078PFz9XHpl668pcBJkH9OwnZ3cFQu40cKdEYwttl8Wlf7nsyMsSKYb+WMUxSaS8kwLTJo7izMNGdfPeyqwqiq17ErfacUO7UxQZt+AuO0mOja+ctVK63dt2KLBlt5fldpu27njor9+8PpGnV6Htsv66yigFVygVYv6SDv57cUcP76189o84Lp08vofdpWGA3NDuFEcSt+By1UrprStzL+iti26S7bIp+01h4fSFLJuxjIXTF9a1MoFk6euvSdJW7xSqtpgvSj6f/SOW3DCVdMJWnHoiibIolwdXb12Ze0FvPb18BVI8SXJ5nQp8Laft4xFtdUuS9CaF7Bsd6zqZeOVCzGB955a8rsBRLshbtplHvzvFUcgtuJweXL11Ze4Fpbjo9kc7RznJ5zb8vyU9C4wPa8lnjpeANGrK9xuSeF0lsW+8vXEL6zq3FHQFdqO8kwqRbsHhUrfcHly9dWXuBf3VRbcWybfldTNwOkEZ3tOzjsPM7LOl3FTS7pLul7QifN0tpt+MsM8KSTOy2h+RtFzS0vDYsxR5SiXJA/7ED40oelcqbivMjfJOKkS5E3/qemhdH3h2ldMduIIp9XuzdZWbUqWQAb/Y/oVIe75KIUuQ4kPSQGAvsrbIzOyVXt9UuhZ4y8xmS5oF7GZmX8vpszuwCGghcD9ZTKDM3pb0CPBVM1tUzH1bWlps0aKihiQiroTvqKZGHp91EvOXdHDpnc8UFSmfQcBLs7v/4edusUHgguwBi07FqIDLb7XITakCwYomTgkV2z/t+1cCSYvNrKVQvyRG+YuA14H7gbbw+J8S5ZvKjlT484BpEX0mA/eb2Vtm9nZ4/9NKvG9ZiFp9ZEe7f+tXz0cqkwEJlixRqw6PfneqSoVdfitNsV5haSd67MuJI5MY5f8ZGB+WAU6LvcxsDYCZrYnZshoFZFvt2sO2DD+RtA34b+Bqi1lqSbqQIHUMY8aMSUP2bmTcgbNvLuDMw3YEK769cUvk2O0WKJ447698KViKDYZ0nNTI5/LbD1YpxXqFpZ3osS8njkwSh/IqsL7YiSU9IOm5iGNq0iki2jLP7fPM7CDguPD4XNwkZna9mbWYWcuIESOK+xAJiDLIG/DwH9cmGp+90tjtfQ00NTb4qsOpbSro8lsKvbVDFJv/K+0a8n25Jn2SFcoq4BFJbUBX4WYz+26+QWZ2Stw1Sa9Lag5XJ83AGxHd2oETss5HA4+Ec3eEr+9Iuhk4ArgpwWdJnSQG+abGBtZ19lylNDU2VG2l4WWBnV5TQZff3lJKAau4lCpxXmHF9i9E2vNVkiQrlFcI7BeDgV2yjlJYAGS8tmYAd0f0uQ+YJGm30AtsEnCfpEGS9oCu1PqfAMqUbKgwSTyuWs84MDKyvfWMA8sqWxzFZCt2nB5U0OW3t5RihyjWK6zUAMjclRTQZwMqE3l5AUgamla2YUnDgduBMQQK6ywze0tSC/BFM/tC2O98gqzHAN82s59IGgo8CjQQFPl6ALjEzAoWay+Hl1dSj6taWhEU8kpznG5EeXRBTXt5TZg3AcvNTQYIsWxG7YTR1aJHVxRJvbwKKhRJRwM3ADub2RhJBwP/aGZfSkfUylEut+FaUhZJ2HdWW8R/tWgXZafOya0rD0XVf68Wk+6cFBkd3zy0mYXTF1ZBop60rWrj6499ne22vce1WpITkiuUJDaU7xG48C4AMLNnJB1fonz9ir7mceVlgZ3E9FGPrlq3Q2RWJlHKBPqGR1cUibINm1muBa7g9pJTu1w6eTyNDQO7tXlZYCeSPuLRlUs5EjvGeY21rWrjuFuP46B5B3HQvIP46C0fLehRVqh8cF/w6IoiyQrlVUnHACZpMPAV4MXyitW/qfYWWeZefWmbzqkSfcCjK440EzvGeY0teWMJd624iy3bd3hxrt+8nssfu7xLhtx55j49N3I7LkMtraSKJYkNZQ9gLnAKwTb7QuArZvZW+cVLl3LZUIohKg1Lw0AxZ/rB/kB3ao8+akNJmzibzAANiN22yrWDRBngo+b714/+a00Z5CHF1CsEUfLnmdleZrZnmBjygNJFrE+i0rBs2WZ861fPV0kix8lDBZM41jJxNo04ZRI1ptA215CBQ2pSmRRDki2v/wscmqDNSUBcGpa4dsepOnH16OuIuJoq+VYouXaQfIb25qHNzDx0Zp9WJpBHoYTuwscAIyRdknVpV4L4D6cKVNv+4jj1SJzX2NT9p/awoQAM0qAedpA4pdRbF+GMPea1Da+x99C9a0Ih5dvyGgzsTKB0siPk/wZML79o/ZOmxobYa4Ui1T3C3XGqQ5zX2OVHXc5Vx15F005NXX2HDR7G1R+9usfDPc1CXxl7zJoNazCsy0mg2nVTkhjl329mf6mQPGWlVozy/3zb0shrhSLVPcLdcfo2aa0qKh24mWZg40ZJc4ADgS71amb+BOsF0yaOilUohUr4FkpE6dthjlPbpOXKXKsp7pN4ef0C+COwL/At4GXgqTLK1O8ZlSeh5PwlHRw7+yH2ndXGsbMf6radlS8RpW+HOU79UKsp7pMolOFmdgOwxcx+Y2bnA0eVWa5+TVyk+okfGpFXKeSLcI+qyxJXk95xnL5NmvaYNEmiUDLuC2skTZE0kaA2idNL4kr4PvzHtXmVQr7Sv0nqsjiO0z8oR2qZNEhiQ7la0jDg/xDEn+xKUBbYKYGohJIXJ7CtxCWi9ISPjlNdKu3Gm2ZqmbQouEIxs/8xs/Vm9pyZnWhmhwEfqIBsdUfcw7/pfQ2xdpUMnvDRcapHJdx4e1vSuJIkyjYcwSWFuzjFEqUUGgaKdzdtLWhsz7cd5jhOeelNhchiFEStxp3kkmTLKwoV7pJnsLQ7cBswlsBr7Gwzezui370EDgCPmdknstr3BW4FdgeeBj5nZptLkakWiMoCvOG9rT3q0WfsKrnKoq/VZXGc/kKxbrzF1rzPp7BqadurtyuUZHWD45kFPGhm44AHw/Mo5gCfi2i/BrguHP82cEGJ8tQM0yaO4vFZJ/HS7Ck8Pusk1ndG5/hyY7vj1A7FuvEWu6Kp1biTXGIViqR3JP0t4ngHGFnifacC88L384BpUZ3M7EHgnRy5BJwE3FlofH8gX+yJ4zi1QbFuvMUqiFqNO8klVqGY2S5mtmvEsYuZ9XarLMNeZrYmvM8aYM8ixg4H1pnZ1vC8HYjd55F0oaRFkhatXbu21wJXCze2O07tU6wbb7EKolbjTnIpVTHEIukBIOrb+UapU0e0xW7Bmdn1wPUQ5PIq8d4Vx6srOk7foBg33mJr3mfmrbXswrmUTaGY2Slx1yS9LqnZzNZIagbeKGLqvwJNkgaFq5TRwOoSxa1p3NjuOP2L3iiIWow7yaVsCqUAC4AZwOzw9e6kA83MJD1MkEL/1mLH9zc8IaTj9E36goIolt56eZXKbOBUSSuAU8NzJLVI+nGmk6TfAncAJ0tqlzQ5vPQ14BJJKwlsKjdUVPoawRNCOo5TSxSsh9KfqIV6KGni9VEcx6kESeuhVGuF4qSAJ4R0HKeWqJYNxSmCODuJJ4R0HKeWcIWSIuUwkGfsJJm09hk7CQQxKtnXwGNUHMepHr7llRLlMpDnK5zlCSEdx6klfIWSEoUe/L2lkJ3EY1Qcx6kVfIWSEuUykHsuL8dx+gquUFKiXA9+z+XlOE5fwRVKSpTrwe92Esdx+gpuQ0mJuCSOEAQgluL55XYSx3H6Aq5QUiT3wZ/P5dcVhOM4/Q3f8ioj+Ty/HMdx+huuUMqIp0ZxHKeecIVSRtzl13GcStG2qo1Jd05iwrwJTLpzEm2r2iougyuUMuIuv47jVIK2VW20/q6VNRvWYBhrNqyh9XetFVcqrlDKiLv8Ok59UO3Vwdyn53YrJwywadsm5j49t6JyVMXLS9LuwG3AWOBl4Gwzezui373AUcBjZvaJrPafAh8D1odNnzezpeWVune4y6/j9G8yq4PMAz2zOgBKqsjYtqotcYng1za8VlR7uajWCmUW8KCZjQMeDM+jmAN8LubapWZ2SHjUpDJxHKf/U47VQbFbWHsP3buo9nJRLYUyFZgXvp8HTIvqZGYPAu9USqhyMn9JB8fOfoh9Z7Vx7OyHvEyv4/QTyrE6KFZJzTx0JkMGDunWNmTgEGYeOrPXMvSGaimUvcxsDUD4umcv5vi2pGWSrpO0U1wnSRdKWiRp0dq1a3srb0l47XfH6b+UY3VQrJKast8UWo9ppXloM0I0D22m9ZjWkrbcekPZbCiSHgCivtFvpDD9ZcBrwGDgeuBrwJVRHc3s+rAPLS0tlsK9i6Zcqe0dx6k+Mw+d2c2GAqWvDvYeujdrNqyJbI9jyn5TKq5AcimbQjGzU+KuSXpdUrOZrZHUDLxR5NyZb/o9ST8BvlqCqGXHAxwdp/+SeYgnNaAnoRxKqhJUK5fXAmAGMDt8vbuYwVnKSAT2l+fSFzE9vPa74/Rv0l4dlENJVQKZVX4XSNJw4HZgDPAKcJaZvSWpBfiimX0h7Pdb4EPAzsCbwAVmdp+kh4ARgICl4Zh3C923paXFFi1aVJbPlI/cJJEQBDh6TIrjOH0BSYvNrKVQv6qsUMzsTeDkiPZFwBeyzo+LGX9S+aRLn7jU9q5MHMfpT3j6+grhAY6O4/R3PPWK4ziOkwquUBzHcZxUcIXiOI7jpIIrFMdxHCcVXKE4juM4qeAKxXEcx0kFVyiO4zhOKrhCcRzHcVLBFYrjOI6TCq5QHMdxnFRwheI4juOkgisUx3EcJxVcoTiO4zip4ArFcRzHSQVPX9/HmL+kw+uqOI5Tk1RlhSJpd0n3S1oRvu4W0ecQSb+X9LykZZLOybq2r6Q/hONvkzS4sp+gOmQqP3as68SAjnWdXHbXs8xf0lFt0RzHcaq25TULeNDMxgEPhue5bAT+3swOBE4DviepKbx2DXBdOP5t4IIKyFx15ty3vFsZYYDOLduYc9/yKknkOI6zg2oplKnAvPD9PGBabgcz+5OZrQjfrwbeAEZIEnAScGe+8f2R1es6i2p3HMepJNVSKHuZ2RqA8HXPfJ0lHQEMBv4MDAfWmdnW8HI7EGtEkHShpEWSFq1duzYV4avFyKbGotodx3EqSdkUiqQHJD0XcUwtcp5m4GfAP5jZdkAR3SxuvJldb2YtZtYyYsSI4j5EjXHp5PE0Ngzs1tbYMJBLJ4+vkkSO4zg7KJuXl5mdEndN0uuSms1sTagw3ojptyvQBlxuZk+EzX8FmiQNClcpo4HVKYtfk2S8udzLy3GcWqRabsMLgBnA7PD17twOoefWL4GbzOyOTLuZmaSHgenArXHj+yvTJo5yBeI4Tk1SLRvKbOBUSSuAU8NzJLVI+nHY52zgeODzkpaGxyHhta8Bl0haSWBTuaGy4juO4zi5yCzW/NDvaGlpsUWLFlVbDMdxnD6FpMVm1lKon6decRzHcVLBFYrjOI6TCq5QHMdxnFRwheI4juOkgisUx3EcJxVcoTiO4zip4ArFcRzHSQVXKI7jOE4quEJxHMdxUsEViuM4jpMKdZV6RdJa4C8xl/cgyGRci7hsxVOrcoHL1ltctuJJS673m1nB+h91pVDyIWlRklw11cBlK55alQtctt7ishVPpeXyLS/HcRwnFVyhOI7jOKngCmUH11dbgDy4bMVTq3KBy9ZbXLbiqahcbkNxHMdxUsFXKI7jOE4quEJxHMdxUqFuFYqk3SXdL2lF+Lpbnr67SuqQ9INakU3S+yUtlrRU0vOSvlhDsh0i6fehXMsknVMLcoX97pW0TtL/VECm0yQtl7RS0qyI6ztJui28/gdJY8stUxGyHS/paUlbJU2vIbkukfRC+Hf1oKT315BsX5T0bPh/8jFJH64V2bL6TZdkksrjSmxmdXkA1wKzwvezgGvy9J0L3Az8oFZkAwYDO4XvdwZeBkbWiGwfBMaF70cCa4CmassVXjsZOB34nzLLMxD4M7Bf+G/1DPDhnD5fAv4zfH8ucFuF/r6SyDYWmADcBEyvIblOBN4Xvv/fNfad7Zr1/gzg3lqRLey3C/Ao8ATQUg5Z6naFAkwF5oXv5wHTojpJOgzYC1hYIbkggWxmttnM3gtPd6Jyq80ksv3JzFaE71cDbwAFo2zLLVcoz4PAO2WWBeAIYKWZrTKzzcCtoYzZZMt8J3CyJNWCbGb2spktA7ZXQJ5i5HrYzDaGp08Ao2tItr9lnQ4FKuXxlORvDeAqgh9em8olSD0rlL3MbA1A+LpnbgdJA4DvAJfWmmwAkvaRtAx4leAX+epakS1LxiMIfjX9uZbkqgCjCP5dMrSHbZF9zGwrsB4YXiOyVYNi5boA+HVZJdpBItkkfVnSnwke3F+pFdkkTQT2MbOybvUOKufk1UbSA8DeEZe+kXCKLwH3mNmraf9wTEE2zOxVYIKkkcB8SXea2eu1IFs4TzPwM2CGmZX8SzctuSpE1B9M7i/WJH3KQbXuW4jEckn6LNACfKysEmXdMqKth2xm9kPgh5I+A1wOzCi3YBSQLfxhfB3w+XIL0q8VipmdEndN0uuSms1sTfjgeyOi29HAcZK+RGCnGCzpXTOLNXpVULbsuVZLeh44jmDrpOqySdoVaAMuN7MnSpUpLbkqSDuwT9b5aCB3BZnp0y5pEDAMeKtGZKsGieSSdArBj4iPZW371oRsWdwK/L+ySrSDQrLtAnwEeCT8Ybw3sEDSGWa2KE1B6nnLawE7fj3MAO7O7WBm55nZGDMbC3wVuCkNZZKGbJJGS2oM3+8GHAssrxHZBgO/JPi+7qiATInkqjBPAeMk7Rt+H+cSyJhNtszTgYcstJ7WgGzVoKBc4dbNj4AzzKySPxqSyDYu63QKsKIWZDOz9Wa2h5mNDZ9lTxB8f6kqk8zN6vIg2Kt+kOAf/UFg97C9BfhxRP/PUzkvr4KyAacCywg8OpYBF9aQbJ8FtgBLs45Dqi1XeP5bYC3QSfDLbnIZZfo74E8E9qNvhG1XEvxnBhgC3AGsBJ4E9qvEv2FC2Q4Pv58NwJvA8zUi1wPA61l/Vwtq6DubCzwfyvUwcGCtyJbT9xHK5OXlqVccx3GcVKjnLS/HcRwnRVyhOI7jOKngCsVxHMdJBVcojuM4Tiq4QnEcx3FSwRWK06+R9G4RfU+QdEwJ92oKg2Cjrj0iaXJO2z9L+o8i73GPpKYCfSI/s6SfFpM5uJrZkJ2+iSsUx9nBCUCvFQrQRJCuJ4pbCALOsjk3bC+IAgaY2d+Z2boSZCyGC4C3zWx/gtQd11Tovk4fxRWKU3dIOj38xb1E0gOS9gp/fX8RuDisZ3GcpBGS/lvSU+FxbDi+VdKN4apjlaRMEsDZwAfC8XNybnsn8AlJO4VzjCVI7f+YpJ3D2h5Ph/U0pmb6SHoxXMU8Dewj6WVJe4TX5yuoifO8pAtzPuN3wvkelNQj07OkwyT9Jhx/X5iuJpdqZUN2+iqViuT0w49qHMC7EW27QVdQ7xeA74TvW4GvZvW7Gfho+H4M8GJWv98RlA3YgyCSvIGghshzeWRpA6aG72cBc8L3gwhraYTzrSRI+DeWIH38UVlzvAzsEb7PZANoBJ4DhofnBpwXvr+CMMMD8FOCFC8NofwjwvZzgBsj5H0OGJ11/ufMvf3wI+ro18khHSeG0cBt4a/ywcBLMf1OAT6c9aN8V0m7hO/bLEhM+J6kNwhq5hQis+11d/h6ftgu4F8lHU+gQEZlzfcXi0+u+RVJnwzf7wOMI1Bu24HbwvafA3fljBtPkCzw/vCzDSQogpZLrWYldmoUVyhOPfJ/ge+a2QJJJxCsOKIYABxtZp3ZjeFDODvL7TaS/V+aD3xX0qFAo5k9HbafR1CA7DAz2yLpZYI8XxDk0upBKPcpoXwbJT2SNSaXqLT5z5vZ0QXkrVY2ZKeP4jYUpx4ZBnSE77PrVbxDkOo7w0LgosyJpEMKzJs7vhtm9i5BYr4b6W6MHwa8ESqTE4EkddKHERjMN0r6EHBU1rUBBFtbAJ8BHssZuxwYIeloAEkNkg6MuEe1siE7fRRXKE5/532S2rOOSwhWJHdI+i3w16y+vwI+mTHKE1Tca5G0TNILBEb7WMzsTeBxSc9FGOUz3AIcTFAvI8MvwvssIlit/DHB57oXGKSgYudVBCnJM2wADpS0GDiJIOtstpybCRTENZKeIciOG+XddgMwXNJK4BICu4/jxOLZhh3HcZxU8BWK4ziOkwquUBzHcZxUcIXiOI7jpIIrFMdxHCcVXKE4juM4qeAKxXEcx0kFVyiO4zhOKvx/Bz8eEiXESfQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#viz\n",
    "z_em=get_latent(iris.data,iris.data.mean(axis=0),lambda_,psi)\n",
    "\n",
    "for i in np.unique(iris.target):\n",
    "    plt.scatter(z_em[:,0][iris.target==i],\n",
    "                z_em[:,1][iris.target==i],\n",
    "                label=f'Label {i}')\n",
    "plt.legend()\n",
    "plt.title('Factor Analysis by EM')\n",
    "plt.xlabel('Latent Variable 0')\n",
    "plt.ylabel('Latent Variable 1')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXmYFPW19z9nhsGZsI0KRpiBIFGJQSaCo4IYrlvAm4lA3Jd7r4nmmrxRQ2LiG71vohPje8UYNXjz5n2j0WiicUMDMeNC4obLRQUXUHEhuDDjqAQzCAgyDOf9o6qHnp6q7uqluqq7z+d5+unuWk8XQ536/c453yOqimEYhmGkUhW1AYZhGEY8MQdhGIZheGIOwjAMw/DEHIRhGIbhiTkIwzAMwxNzEIZhGIYn5iAMIwURaReRw/PYv1pENonImCz3+4aIPJrrebM4z9Ei8la264zKwxyE0QcReUtEtrg3uMRrVB7H21tEQim2SRxbRK4N4/i5oqo9qjpYVd+JygYRmSgifxGRf7ivZSIyMyp7jNLEHIThxbHuDS7xejcqQ0RkQJrVZwAfAqeKSE2RTIo9IiLAn4H7gD2APYHvAZsitCndv6MRU8xBGIEQkSoRWSAi74lIl4g8KiL7Ja3/lIhcIyLviMgGEVkiIrsAS9z1idHIQe6xLhaRt0XkAxG5SUSGutslRgVfF5F3gMU+9gjwr8BFgAAtSesGuMf4poisdp+gr01av4+IPCIi60Xk7yLyexEZ5nGOBhH5WETqk5Yd4l6DASKyr/s7N7jH+UPK+ce6378iIqtEZKM7ffW9NJe6SkR+5R5zlYgc4R7jVBF5OsW+H4rIAo9jfBoYA1yvqt2q+omqPq6qT/pcy++JyEteI0URaRSRP4rIOhF5U0TOSVo3VUSWun8PnSJybcJRJ12Db4vIauDVTP8uRgxRVXvZq/cFvAUc7bG8CvgaMASoBX4JLEta/2vgIWAkUA0cBtQAezt/Zn2OdTbwOrCXe7xFwG/ddXsDCvwW+BRQ52PnEcAWYBjwf4F7ktYNcI+xyF0/FmekcbS7fl/gKGAgzhP2k8DPk/ZvBw53Py8G/j1p3X8B17if7wJ+6F6bWmBayvnHut/XAYe6n3cDJvv8pm8A24HvuNfuNKALqAfq3M/7JG2/Epjt82/1N/f3zwb2SFl/NPCW+/knwDJguMe6auAF4D/ca7W3+/dxlLv+IOAQ9/eOc/9Nz025Bg8Au7r2p/13sVf8XpEbYK94vdwbwCb3ZtQFLPTZbrj7n32QeyP5BJjgsZ2Xg3gMODvp+wR3/yp2OogxGey8CVjgfv6iu//u7vfEjWhK0vb3AD/wOdYJwLNJ35MdxOnAY0nH/SBxgwf+gOOcGlKOl+og3nVv/kMy/KZvAGsBSVr2HHCq+/l64Cfu5wOAvwM1PscaDfwKWAP0AI8An3XXHQ28DczHGeENTdov2UFMA9akHPfHOCMTr3P+ALgr5RpM97gugf5d7BX9y6aYDC/mqGq9+5oDvZk5PxORNSLyEbDa3XY4zpTGQJyn1iCMwrlBJXjb3X9E0rK1fjuLyCDgeOBWd9ETQCdwasqm7yV9/hgY7O6/p4jcKSId7m+5yf0dXvwR+II4GUnHAOtU9Tl33fdxnvSXichKETnD5xhfBWYB77hTc4f4/TagXd07p8vbONcL4GYchwXwL8AdqtrtdRBVXauq31bVcTgjtW73dybYHcchXaaqH/nY8hlgjDuF1CUiXcD/xIlpICKfE5E2d8rtI+BS+l9Hr39Hz38XI36YgzCC8m/Al4EjcaYH9naXC/A+sA34rMd+XhlM7+LcfBKMcfdf17tT35tkKsfj3FSuE5H3cJzDnq6NQbgCZ8QxUVWH4kydideGqvoxcDfOjflfgd8nretU1W+o6kjgHNeevTyO8bSqzsKZzvozcHsa2xpTvo/BuV6o6hMAIjINxxn+ngCok031K2D/pMV/x3Fat4jIFJ9d1wJvJD0s1KvqEFU91l3/a+AlYG/3Ol5M/+toctEljDkIIyhDcG6q63FiA/87sUJVe3CeTn/hPp1Xi8g0N2D5AaAiMi7pWLcB54vIWBEZ4h7rNlXdEdCWM3CmWybiTLUcAEwHmpMD5xl+y2Zgg4iMxpkaScfvgDNxAuG3JBaKyEki0uB+7cK5GfYk7ygidSJymogMdZ/2N6Zuk8JIETnXDeieguN0H0ha/3ucaa3NqrrU6wAiMlxELhGRceIwAvg60Gd7VX0Ix6kuEpFmj0P9N7BNRL4vIrXuv+tEETnQXT8E2ABsdq/7N9P8LqMEMQdhBOW3OE+y7wIvA0+lrP8esApYjhN4/E+cufSNwOXA0+40RTPOzf0O4HGcOfKNwNwgRrhTPYcDv1DV95JezwB/xXEembgEOBjn5vYnnBFCOpbgxFmeVtX2pOWHAM+KyGacufRz1Lv24QzgbXca5iyckYgfT+HEZD4EWoHjVfUfSet/hzMSSDd6+ATHsTyCE09a6b6fmbqhqj4A/DvwZxE5IGXddpxR48E4sam/44wahrqbfN/9bRvd5XeksckoQST9SN4wDAARWQLcqKo3RWzHIJxR2f6q+maUthjlj40gDCMD7hz9/jhprVFzDvCkOQejGFh1o2GkQURuxYk9nKeqmyO2pR0nG2l2lHYYlYNNMRmGYRie2BSTYRiG4UlJTTENHz5cx44dG7UZhmEYJcXy5cv/rqojMm/Zl5JyEGPHjmXZsmVRm2EYhlFSiMjbmbfqj00xGYZhGJ6YgzAMwzA8MQdhGIZheBJZDEJEanEkDHZx7VigqpdEZY9hGPGnu7ub9vZ2tm7dGrUpsaS2tpbGxkZqagrTYDHKIPUnwJGquskVdXtCRO73EyAzDMNob29nyJAhjB07FhFPAd6KRVVZv3497e3t7LVXP1HhnIhsikkdEj1ya9yXVe0ZhuHL1q1b2X333c05eCAi7L777gUdXUUag3Dlg1/AER/7i6o+nWkfwzAqG3MO/hT62kTqIFS1R1UPwGmScrCI7J+6jYicLSLLRGTZunXr+h/EMAzDCIVYZDGpahfwKE5Lx9R116lqs6o2jxiRdSGgYZQcbWvamLFgBk03NzFjwQza1rRFbZKRxODBwTuktra28vOf/7wgx3/ggQcYP348e++9N/PmzcvqmLkSmYMQkREiUu9+rsNplv5qVPYYRhxoW9NG61OtdG7uRFE6N3fS+lSrOYkKp6enh3POOYf777+fV155hdtuu41XXnkl9PNGOYIYCTwiIiuAZ3FiEH+O0B7DiJz5z81na0/fIOPWnq3Mf25+RBaVNguf72DavIfZ68I2ps17mIXPd4RynnvvvZdDDjmESZMmcfTRR/P+++/3rnvxxRc58sgj2Weffbj++ut7l1955ZUcdNBBNDU1cckl6TP8n3nmGfbee2/GjRvHwIEDOeWUU1i0aFEovyWZyNJcVXUFMCmq8xtGHHlv83tZLTf8Wfh8Bxfds5It3U4L8I6uLVx0z0oA5kxqSLdr1hx22GEsXboUEeE3v/kNP/vZz7jqqqsAWLFiBUuXLmXz5s1MmjSJlpYWXnrpJd544w2eeeYZVJVZs2axZMkSpk+f7nn8jo4ORo8e3fu9sbGRp58OP6enpMT6DKPc2XPQnnRu7vRcbmTHlQ++1uscEmzp7uHKB18ruINob2/n5JNPprOzk23btvWpQ5g9ezZ1dXXU1dVxxBFH8Mwzz/DEE0+wePFiJk1ynpE3bdrEG2+84esgvPr2FCObKxZBasMwHOZOnkttdW2fZbXVtcydPDcii0qXd7u2ZLU8H8477zzOPfdcVq5cya9//es+tQipN3IRQVW56KKLeOGFF3jhhRdYvXo1Z511lu/xGxsbWbt2be/39vZ2Ro0aVfDfkYo5CMOIES3jWmg9tJWRg0YiCCMHjaT10FZaxrVEbVrJMaq+Lqvl+bBhwwYaGpxRyc0339xn3aJFi9i6dSvr16/n0Ucf5aCDDmLmzJnceOONbNrk1Ap3dHTwwQcf+B7/oIMO4o033uDNN99k27Zt3H777cyaNavgvyMVm2IyjJjRMq7FHEIBuGDm+D4xCIC6mmoumDk+r+N+/PHHNDY29n4///zzaW1t5cQTT6ShoYEpU6bw5ptv9q4/+OCDaWlp4Z133uHHP/4xo0aNYtSoUaxatYqpU6cCTmrrLbfcwh577OF5zgEDBvDLX/6SmTNn0tPTw5lnnsmECRPy+h1BKKme1M3NzWoNgwyjclm1ahX77bdf4O0XPt/BlQ++xrtdWxhVX8cFM8cXPP4QN7yukYgsV9XmbI9lIwjDMMqWOZMayt4hhInFIAzDMAxPzEEYhmEYnpiDMAzDMDwxB2EYhmF4Yg7CMAzD8MQchGEYRhZEJfd95plnsscee7D//v3a5oSGOQjDMIwS4Gtf+xoPPPBAUc9pDsIwjPJlxZ1wzf7QWu+8r7gzlNOELfcNMH36dHbbbbdQ7PfDCuUMwyhPVtwJ934Hul1xvg1rne8ATScV9FRhy31HhTkIwzDKk4cu3ekcEnRvcZYX2EGELfcdFeYgDMMoTza0Z7c8D8477zzOP/98Zs2axaOPPkpra2vvunRy39/85jcLbkshsRiEYRjlybDG7JbnQdhy31FhDsIwjPLkqIuhJqX3Q02dszwPEnLfidfVV1/dK/f9xS9+keHDh/fZPiH3PWXKlF657xkzZnDaaacxdepUJk6cyAknnMDGjRvTnvfUU09l6tSpvPbaazQ2NnLDDTfk9TuCYHLfhmGUDNnKfbPiTifmsKHdGTkcdXHB4w9xw+S+DcMwgtB0Utk7hDCxKSbDMAzDE3MQhmEYhifmIAzDMAxPzEEYhmEYnpiDMAzDMDwxB2EYhpEFUch9r127liOOOIL99tuPCRMmMH/+/KyOmSuW5moYhhFzBgwYwFVXXcXkyZPZuHEjBx54IF/60pf4/Oc/H+p5IxtBiMhoEXlERFaJyMsiMjcqWwzDKE/a1rQxY8EMmm5uYsaCGbStaQvlPGHLfY8cOZLJkycDMGTIEPbbbz86OjpC+S3JRDmC2A58X1WfE5EhwHIR+YuqvhKhTYZhlAlta9pofaqVrT1bAejc3EnrU60AtIxrKei5iin3/dZbb/H8889zyCGHFPQ3eBGZg1DVTqDT/bxRRFYBDYA5CMMw8mb+c/N7nUOCrT1bmf/c/II7iGLJfW/atInjjz+eX/ziFwwdOrSgv8GLWMQgRGQsMAl42mPd2cDZAGPGjCmqXYZhlC7vbX4vq+X5UAy57+7ubo4//nhOP/10jjvuuEKZnpbIs5hEZDBwN/BdVf0odb2qXqeqzaraPGLEiOIbaBhGSbLnoD2zWp4PYct9qypnnXUW++23H+eff37B7fcjUgchIjU4zuFWVb0nSluMMqFIPYiN+DN38lxqq2v7LKutrmXu5PzyYaKQ+37yySf5/e9/z8MPP8wBBxzAAQccwH333ZfX7whCZHLf4oy7bgY+VNXvBtnH5L4rmCCyzak9iMHR/z/2WlP0LBOylftuW9PG/Ofm897m99hz0J7MnTy34PGHuFEuct/TgH8FVorIC+6y/1DV8N2iUVoEbT5fxB7ERmnQMq6l7B1CmESZxfQEIBk3NIygN/4i9iA2jEog8iC1YWQk6I2/iD2IjegopS6YxabQ18YchBF/gt74Q+pBbMSH2tpa1q9fb07CA1Vl/fr11NbWZt44ILGogzCMtBx1sXfwOfXGn5huqrAexJVEY2Mj7e3trFu3LmpTYkltbS2NjYUbMZuDMOJPNjd+60Fc1tTU1PSpUjbCxRyEURpEeeMPkmJbIuetxLRPI3fMQRilQxQ36qAptiVw3mKK1xnlQWSFcrlghXIVTFRFcNfs79ycUxk2Gr73Ukmdd8aCGXRu7uy3vH6XeuoG1NmooozJtVDOspiM0iBdLUSYRFVbEcJ5/UTquj7ponNzJ4r2jirC6ptglBbmIIzSIKobdVS1FSGcN6hIXUIS2zDMQRilQVQ36qhqK0I4r5d4nR9hSGIbpYc5CKM0iOpG3XSSE+cYNhoQ570Y4n8hnLdlXAuth7YyctBIBGHkoJEMGzjMc9swJLGN0sOC1EbpEFW6aRmTmtkEjiR266GtFqguI0pRzdUwssOK4ApOwgl41UZYzYRhDsIw4kCEoyMvSWyrmTDAYhCGET2JGo8NawHdWRQXYTe8nzz1kz7TTmDZTZWIOQjDiJqoajx8uGzpZWzp2eK5zrKbKgtzEIYRNTFrdHTX63f5rrPspsrCHIRhRE3MGh3t0B2+6+ZOnltES4yoyclBiMjKQhtiGBVLzBodVYn3bUEQC1BXGL5ZTCJynN8qwMaZhlEoYtbo6MR9T+SO1+7ot/yk8ZZiXGmkS3O9A7gV8KqkK1xPO8Mw/Gs8Ikh//dGUHwFOLGKH7qBKqjhx3xN7lxuVg28ltYgsB85Q1X7awiKyVlVHh21cKlZJbVQUUUmcG2VHGHLf3wU+8ln31WxPZBhlx4o7nb4NrfXOe6HrFmKW/mpUHr5TTKr6eJp1FfkYv/D5Dq588DXe7drCqPo6Lpg5njmTGqI2q/QoB02lYnSai1n6a6ExKY/4Y2muAVn4fAcX3bOSjq4tKNDRtYWL7lnJwuc7ojattIhh1XBOFOPpPmbpr4UkIeVhjYrijTmIgFz54Gts6e7ps2xLdw9XPvhaRBaVKOUybeLVDjTd8lyIOP21bU0bMxbMoOnmJmYsmFHQm/f85+ablEcJUHFifblOE73b5S094Lfc8KFcpk2kGrTHe3mhiDD9NWyxPj/JDpPyiBcZRxAi8mkRuUFE7ne/f15EzirEyUXkRhH5QERC7P6+k2ymiRY+38G0eQ+z14VtTJv3MMPqajyPOaq+znO54UO5TJt4OQe/5WEHs0Mg7Cd8P8kOk/KIF0GmmG4CHgRGud9fx8lwKgQ3AccU6FgZCTpN5OVINm/bTk2V9NmurqaaC2aOD9vs8iJmVcM5M8wnyzt1eT4xlwjjNWE/4Xu1P62trjUpj5gRxEEMV9U7gR0Aqrod8Hl8yg5VXQJ8WIhjBSHoNJGXI+nuUWqqhYb6OgRoqK/j8uMmWhZTtkTVwrPQBHV0+cRcIozXhP2E79X+tFBd7MKMnVQaQWIQm0Vkd9yKahGZAmwI1aqQGFVfR4eHk0idJvLaBuDj7h38p6W25k/cO8MFScMNGh/IJ+ZSxHhNasrp9MbpLFq9qF8r0kI+4Xs1KsoXa3RUWIKMIM4H/gR8VkSeBH4HnBeqVUmIyNkiskxElq1bty6vY10wczx1NX2DiF7TRNXSdyopGctaKnOymdZpOgm+9xK0djnvXk4vn5hLkeI1Ximni1YvYnb9BEb2KKLKyB6ltfGY2N9kLTuqsGR0EKr6HPBPwKHAN4EJqroibMOSzn+dqjaravOIESPyOtacSQ1cftzEjNNEPT7yI2BZS2VPoad18om5FCle43dTXfL+Myx+Zy0r3lrL4nfW0vLk9bEPsFt2VGHJRc11XxFBVe8JyaZQmTOpIeMUUYPPVBRY1lLZkTqd5FvfEHBax2t66thrc0tVLVKaq+9NtTrl+THhKGM8PbjnoD3p3NzpudzInnQxiGPTrFMgbwchIrcBhwPDRaQduERVb8j3uPlywczxXLDgRbp7+o4kaqrEspbKCS+5DARPAeMg0zp+8hvHXutMQeVCEeI1vjfV7R65KDGvV5k7eW6fGARYdlQ+pNNi+nrYJ1fVU8M+Ry4kRhg/ufdl/vFxNwD1dTW0zppgAepywms6CaWfkwg6rZNueiqqp+4AAXfPm6oqc//R1f94Ma9XScRITOOpMGTMYnIzmC4BDsP5X/MEcKmqrg/ZtkgJMhVllDi+T8PqpN8m31TBKXILK2MpDAIKCnreVIcfQkvH9X2PVyL1KmFkR1UqQdJcbweWAMe730/HaSZ0dFhGGUZR8Is5DBvdd0ooqHKr7/Eao1GwzWJE43lT3W1i6avuGnkRxEHspqo/Tfp+mYjMCcugUsRkwEuUoy72bsiTTbFb8g3T73i7jYN7zqZ32ioMaXAv8h3RxL1exQidIHUQj4jIKSJS5b5OAqw00cVkwEuYoFXdQW+0iePV7bZzmVTBm4/RL/BdjIroctG9MiIjXZrrRnZG7M4HbnFXVQGbcOISFU86fScbRZQAQZ6S000debE9aQSxbbP/ccOOTQQdIUWINQ2KN74jCFUdoqpD3fcqVR3gvqpUdWgxjYwzJgNeAWRTsOaZGeVD2E/yEeteZdJEsqZB8SdQPwgR2RXYB+iVX3SF9iqO1HjDsLoaurZ099vOCurKiGwK1gKPCqQ4T/JNJ9E2eNDOp/TXf8PcwYNCf0oPoonkV8F9+dOX2ygiJgRJc/0GMBdoBF4ApgD/DRwZrmnxY+HzHX0K6Dq6tlBdJdRUCd07ds4xmwx4GRI0YJuuGrsXgeYzy6Lxjx/pNJES5/Wr4N6wbQNta9rMScSAIEHqucBBwNuqegQwCchPNa9E+cm9L/erru7ZoQwcUNVP3wno03Aom6B1arMiC3iXEF7TUdUD3cC1O81z3HXwlauLYk5U4nVBNJHSyV+YuF48CDLFtFVVt4oIIrKLqr4qIhX5eJyoqk5l87YeXr5054AqkdmUCF4nMpuAjIHrfPY1YkCEbUK9iEq8Logm0tzJc7nw8Qs99zdxvXgQZATRLiL1wELgLyKyCHg3XLNKm6Cd6wq9rxETgsiAF4moWnsG6RjXMq6F+l3qPfevG2AxvDgQRO77q6rapaqtwI+BG4CKLJSr9+lLnbo8n8wmy4oyConfjXp64/RQu64F7Rh34cEXUi3V/fb/ePvHXLb0soLaZGSPr4MQkaHu+26JF7ASR4tpcJHsixWtsyb060tdUyW0zprQZ5lfBlOQzKZ89jUCsOJOR1Optd55D6O/QTHOERCvG/XsvWezaPWi0NNLW8a1sPiExaw4YwWLT1jsGXRuGddCj3p3ML7r9bsKao+RPeliEH8AvgIsp6/EZeJ9XOjWxYxEDCCTrMYFM8f3iSNA8MymfPY1MhBUUynbYybHG/aZAS/+obDnyJNUnaUZC2ZkzDCKAzt0R9QmVDzp5L6/IiIC/JOqvlNEm2JNEJXXoI6k0PsaGSi0HLeXw1l2I76yGjHRNYpb17UqqfJ0BlXSf4LDKq+LS9osJlVVEfkjcGCR7CkbUm/0iSBzUCdhDiEECi3H7dtPooDnSFBANdi4dV07cd8TueO1OzyXJxNVTUclEySLaamIHBS6JWWGifjFkEKL12Vz089HViMxUtmwFtCd01Y5xjaCZBgVkx9N+REnjz+5d8RQJVWcPP5kfjTlR322i6qmo5IRVZ8nnsQGIq8A+wJvA5txYxCq2hS+eX1pbm7WZcuWFfWcuUp5T5v3sGdf64b6Op68sOKK0ONB6pQQOEVt2egTJT/JSxV4Blg9OtIde63zOZdRwDX7B+tbkQVRTdXkc96mm5tQjxGaIKw4Y0WhTS0rRGS5qjZnu1+QQrl/zsGesiCfojUv55BuuVEE8i1iS3UwXs6hpg6+cBq8sbh/R7pcA+RBp8ZW3An3/xC2fOh8r9sN/vkKz+NH0XUt3ymiuE2NVQIZHYSqvg0gInuQJNZXCeQj5V0tQo/H6KxaxGNro2jk0wTHT6lVqkF3pHc41+yfe4A8iNz4ijth0TnQs23nsi0fwsJvO58D/uYwRxZB9JnS4dk7O8KpsUogiFjfLOAqYBTwAfAZYBUwId1+5YBfcVpH1xamzXs47bSTl3NIt9woAfye5HWHUzWdy75B4hj7zIBlN3gvT/DQpX2dQ4Id3fDQpX0VXX1u/F5P+D9+8sdc/vTlfLTto4wOI5NzyTd7yrN3tmUxhUqQKaaf4ii4/lVVJ4nIEcCp4ZoVD0bV12WcKuro2sIFC14E+k47Nfjs22AFb8WlkL2gs20cVKh931iceXkaR9O2/cNAUzteT/jdO7rZsG1D2v0g2PRRIaaIopgaq2SCZDF1q+p6oEpEqlT1EeCAkO0qKLmqo14wczx1Nf1lAFLp7lF+cu/LGfe1grciU+Dsn6waBxVyX9/Rx9qdFdseNQMJ5u++W6DsnyBP8n5ZQ0EyjOKWPWVkJoiD6BKRwcAS4FYRmQ9sD9eswpFPuumcSQ1cftzEPlLefqQqvXrte/lxE62+oZikK4zLhXw6tOWzr+8oQ3Y6Px+5CqpqeK/aO+6V6hCCPsl7ORKvkUHq8qD6TEZ8CJLmOgjYipO7dzowDLjVHVUUlVzSXAudbjr2Qn+9mrfm2R96rGitx7twTTLHDOKEV3puaipt7+IqJyYCvVlMM17/jecNvH6HUtfTw3sDqtlzYD3Txx3DotWL+o0EUhk5aCSLT+g77fWF333Btxr6xX97sc8yq4YuPgVPcxWRXwJ/UNWnkhbfnItxUVJoddR6nxajfkqvRoTkM+8fJxKjjOQUVr+KbVVo3dBn0dzBg/pl/9Sosgmlq8a5BXR2b2DR6/cwe9/jWNK+hPc2v8ewXYaxadsmtuvOCQO/KSE/3aTU5VYNXVqkm2J6A7hKRN4SkStEpKTiDgkKrY7qp+g6YdQQPnvRfYy9sI3PXnQfP1q4MqfjGwUkn3n/OLI9wEONh/PrN7XTo3yqZwfbq/r+99+q3SxpX9KrwPr4KY9z2WGXBZoSGjlopKc5qcutGrq0SCfWNx+YLyKfAU4BfisitcBtwO2q+nqRbMyLQqujeonpjd29jif/9mHvNj2q3LLU0Te8bM7EPKw38qJY3d0KmSnlh18NRjJpnF+f7J/WeprGeo+iUuMLQbOGgtYo5JvqatNTxSVoodwVwBUiMgm4EbgEyJzekwEROQaY7x7rN6o6L99jphKGOmrqMf1SYW97eq05iKjJpzAuCGFIiHuRtl5CsnNMwxrZc3sPnTX9//vnWpUctEYhn1RXm54qPkGC1DXAMTijiKOAx4DbVHVhXicWqQZeB74EtAPPAqeq6it++0ShxeRFqgRHOixwXeaEoJMU+nlW3EnbXy+gddfBbE2aZqqVGloP+2moN9vUmzw4I40g2UwzFszwdC5eQXOjL2EEqb+EUxDXAjwD3A6craqbc7bbJz9eAAAXwUlEQVSyLwcDq1V1jXu+24HZgK+DiAteEhxeVIvkLPZnlAiFlhD346iLvYUGc4mnNJ1EC8DjlzJ/l+29WUxzp1wU+pN4PtXQcetjUQmkm2L6D5yucj9Q1Q/TbJcrDUDyI1E7cEjqRiJyNnA2wJgxY0IwI3uCZkBNGbdrzmJ/RolQrEypQsdTmk6iJeEoikyu1dAm1ld8fLOYVPUIVb0+JOcATiJ3v9N62HGdqjaravOIESNCMiU7MmVAVYvwL1PG8Nb6Lb5if0aZUMxMqaaTnOmk1i7nPSYd6oqFVWIXnyBaTGHRDoxO+t4IvBuRLVnhlxmVWim9l09RXa41GEYMKVamVAwpdkaRifUVnygdxLPAPiKyF9CBEwQ/LUJ7sqK2pqrXQdTX1dA6a0K/aSNnm/4FRLU1QRROjJIh7EypGFKsjCIvJ2QB6eKR8U4lIlcEWZYtqrodOBd4EEc+/E5VfTn9XtGTyGBK1l76ZLt3FWm2y73OlYvIoGGETTEK3hJOqHNzJ4r2OqG2Nf5yN0ZhCTKC+BLww5Rl/+yxLGtU9T7gvnyPU0yyaSK0wyeDeIdmbmWaTzc7o8IpQuFeWBlFySMGEekn1eHXYMgK6MIhXZrr/wC+DYwTkeSGr0OAJ8M2LK6k03ZKvelXibeTqBIy3vzz6WZnVDBFKtwLI6ModdrKr0Yr9bxWQBce6aaY/gAcC/zJfU+8DlTVfymCbbHEL4Op/lM1/WTF/fTUdhlQlTG7qdAig0aMWXHnzr4O1+zvfPdaFoRCS5y7tK1pY8aCGTTd3MSMBTOY3ji94BlFXtNWXlSl9L4wfafwSJfmukFV31LVU3EyjrpxbnmDRSQeBQkR4NUISICt3T39bvo7gLqaqt4+1In0160egWvoe/MvtMigETK53tC9mhotOsfpJZ1Lo6MQCve8YgGLVi9i9t6zC9rbIej0VOq0kxXQhUeQntTnAq3A+zj3PHAcRVN4ZsWXOZMaWPb2h9y69J3eAYKCZ7YSwNbuHbyZIrfxyKvrPPWbkm/+hRYZNEIkn2kdryd+r97SiVFApuOFULjn94SeUH4tFH7TVqmkKsRaAV14BMm3/C4wXlUnqOpE91WRziHBI6+u85s96ofXE3+QdqTWka6EyGdaJ5sn+yDbhlC4V6wndK9CuFS8prGsgC48gmQxrQU2ZNyqgvBTb03F74k/qMLsnEkN5hBKgXymdfye+P22zUQIhXvFekL3KoSb3ji9t4GRX3aSFdCFRxA11xuA8UAb8EliuapeHa5p/YlSzTU5Q8nvionAqGF1JsxXaeSjtOrVTrR6oNMZbkdS58KauuA9rAtMPgqsRjwouJprEu+4r4Huq+xJTVc94nMjuHt5R0YFV1Vy6nNtlDj5KK36PfF7LYuoWtue0CuXjCOI3g1FBhVQ6jsnijGC8Or14NMevh/1dTUM2mVAnxEEFLZZkRFTitFVzjByJNcRRJAppqnADcBgVR0jIl8Avqmq387N1NwphoOYNu/hwDGGZGqqBAS6e/pez+oqoSepWs5L1M8wDCNMcnUQQbKYfgHMBNYDqOqLwPRsT1QqZFOIVi3Sm2E0uHZAP+cA9HEOYHLfhpEgtfjONJbiRyBZUVVNjcBlbqdWovgVoqU2r6irqebUQ0Yzqt4JSieL92XCqqGNciHXm7wJ8ZUGQRzEWhE5FFARGSgiP8BRXy1L/GoUTp8ypk9NwvEHNnDHs2t7pTWywaqhjXIgn5u8yWOUBkGymL4FzMdpEdoOLMYR8StLgtYoTLp0seeUUiasGtrIixgFw9Pd5DNlOGVTfGdKrdERxEGMV9XTkxeIyDTKWNE1SIFaNlNKu36qhq6Puy2LyciPIim1BiWfCuugxXem1BotQaaY/ivgMsOD+roanr94Bm/Oa+HJC48052Dkjp+kx/0/zE0oME/8KqmDVFgHlcewqahoSdcPYipwKDBCRM5PWjUUqPbeq3Kor6uha0v6UURdTTWtsyYUySKj7PGT7tjyofOCoo4q5k6e61lhHUQDKWjxnSm1Rku6KaaBwGB3myFJyz8CTgjTqFKgddYELrjrRbqT0lircPpNf+wqu1rvaSMwQWILQXWbgiq/5km+FdYt41oybmtKrdHi6yBU9THgMRG5SVXfLqJNJYFXMDshyZHgHx93W5tQIzNBYwtekh5+5NH/IRuC3OTzIZ9RipE/QR5xPxaRK0XkPhF5OPEK3bIS5M8vdmbsFGcY/QgqF950kiPYN2w0IM573W7ex8yj/wOQewOkAtMyroXWQ1sL2pjICE6QLKZbgTuAr+CkvJ4BrAvTqFIgVbMpnTyHFcYZaclGLrzppL6jCi812Dz7P8QtWyrsUYrhT5ARxO6qegPQraqPqeqZwJSQ7Yo9Vz74WkZ11wRWGGekxe9pP2j/h9RRRb6y4CH1tTZKjyAjiESqTqeItADvAnmOX0ufoKOCmiqxwjgjPfnIhUP/UUW+hNTXOuxiNyuoKzxBHMRlIjIM+D5O/cNQnDakFc2o+rpgqq+pIk6GkUoIXeDyosB9rYtR7GYFdeGQcYpJVf+sqhtU9SVVPUJVDwQ+WwTbYo2XZpOXL+juUQtSG5lpOsnpPtfa5bxn4xwKHVAucF/rYhS7WUFdOOSaqH9+5k3KmzmTGrj8uIl9BPz8lJksSG2ERiKgvGEtoDsDyvk4iQLHNYpR7GYFdeEQZIrJC5s4ob9mk1+zIa8gdWpbU9NoMnIiXUA5nymqAsY1ilHsFsY5LKaR+wgiexnTJETkRBF5WUR2iEjWXY7iyhGfGxFoeSJFNiEV3tG1hYvuWcnC5zs89zcMX0IIKBeaoLpLcTqH9atw8HUQIrJRRD7yeG0ERuV53peA44AleR4nVjzyqnd5SOpyrxRZK6gzciKfFNkiUYxit2zPkanRkcU0HNJJbQzxW5cvqroKQKS8Zqr8Yg2py/22y6UXtlHh5JsiWySKUewW9BxBMp4spuFganIFxK8gLnV5uramNs1kZEW6gHJM5DKixGukEGR0kI+UeTkRmoMQkb+KyEser9lZHudsEVkmIsvWrYu3wodf6mtqDOKCmeM9o/wKNs1k9CXITd4rRTaM7CZy70EdBX5xBK9gNvQdHRQjblIKhOYgVPVoVd3f47Uoy+Ncp6rNqto8YoR3EDguzJnUwPEHNvS5+Stw9/KOPiODOZMaLCXWyEw+N/kQ5DJKLXDrN1KoEu/bXvLowEQCHXJNczV8eOTVdf1u/okAdHIaa4NPJbbpNhm95JPC6pvdtNYZieRQsZ1PD+oo8IsX7NAd1FbXZpQQN5HAiGIQIvJVEWkHpgJtIvJgFHaEQdBAtdd0VF1Ntek2GTvJJ4XVN4tJcp52KrXArV+8IDEaqPTRQRAicRCq+kdVbVTVXVT106o6Mwo7wqD+UzWBlntVYl9+3EQrljN2kk8Kq5dcBkK/EqYspp3iFLgNEgtJF0doGdfC4hMWs+KMFSw+YbE5Bx9siqnAqE9wwWt5aiW2YfQhnxRWLwFAv3alAYvq4tLdLagwX74tUaMkLlXc5iAKzIYt3VktNwxf8lV5TZXLuGb/vFRa43LDzSYWUopxhDgp05qDKDB+MuAWfDZyopC9HgpQVBeHG26pxUKyJU7JAFYoV2As+GzEljC6z0VAIWIhca7niJMDtBFEgUnEFEyp1SgIK+4sbCMhrxFJoc8RMvnGQuI0heNFMdRvgyLqF1WNIc3Nzbps2bJIbTCZbqNoJArlUqeECvnUX4xzhEA+QdwZC2Z43oBHDhrJ4hMWR25jqgMDxwHmk4orIstVNWvlbHMQWZCQ6U5WYq2rqbb0VCMcfIPKox1JjVI5R8xourkJ9dAyEIQVZ6woyDnyvckXOospVwdhU0xZkE6m2xyEUXCK0euhBPpJFJpiTOHkG2iOQzIAWJA6K4JWSRtGQShGr4cS6CdRaIohxBenQHM+mIPIgqBy3oZRELyqoQvd66EY54gZxRDi8xuNDB04NLbZU15YDCILLAZhFJ1iZBiVWBZTKeAVgxggAxARunfsLJrNN/gcFAtSF4kfLVzJbU+vpUeVahFOPWQ0l82ZGKlNhmHEj9RA85btW+j6pKvfdoXOnvLCgtRFYOHzHdy9vIMe16n2qHL38g6aP7ObjSAMw+hDaqC56eYmz+3iHJewGEQWpMtiMgzDSEec1HCDYg4iC/yylTq6trDXhW1Mm/ew9ZQ2DMOTUmxjalNMWeAnxAeOyn5H1xYuumclgE05GYbRh7io4WaDjSCywEuILxWbcjIMo1ywEUQWpArx+eV/ZVM4Z9pOhlEZxF0k0AsbQWTJnEkNPHnhkbw5r4WGPAvnEnUVHa6zSUxRWRzDMMqPdPIbccUcRB7k2/vBsqIMo3IoRfkNcxB5MGdSA5cfN5GG+joEaKivy6qq2rSdjNBZcaej2Npa77yvuDNqiyqWUkxztRhEnsyZ1JBzzMDakxqhktrrYcNa5zuYlEYE5NvoKApsBBEh1p7UCJWHLu3bCAic7w9dGo09FU4xRAILjY0gIsTakxqhUoG9HuJOXPo8BMUcRMTkM0VlGGkZ1ujTLa58ez0YhcWmmAyjXKnAXg9GYTEHYRjlStNJcOy1Tn9pxHk/9loLUBuBsSkmwyhnmk4yh2DkTCQOQkSuBI4FtgF/A76uqv07aZQ4JqNhGPEmtalP3MXzik1UU0x/AfZX1SbgdeCiiOwIDZPRMIx4k9BG6tzciaK92khx7xNdTCJxEKq6WFW3u1+XAmWXVmEyGoYRb0pRG6nYxCFIfSZwv99KETlbRJaJyLJ169YV0az8MBkNw4g3paiNVGxCcxAi8lcRecnjNTtpm/8FbAdu9TuOql6nqs2q2jxixIiwzC04fnIZJqNhGPGgFLWRik1oDkJVj1bV/T1eiwBE5AzgK8DpqurXWqFkMRkNw4g3pdgCtNhElcV0DPBD4J9U9eMobAgbk9EwjHhTii1Ai41E8fAuIquBXYD17qKlqvqtTPs1NzfrsmXLQrXNMAyj3BCR5aranO1+kYwgVHXvKM5rGIZhBCcOWUyGYRhGDDEHYRiGYXhiDsIwDMPwxByEYRiG4Yk5CMMwDMMTk/vOgCmyGoZRqZiDSENCkTUhupdQZAXMSRhGhVMJUuE2xZQGU2Q1DMOLSpEKNweRBlNkNQzDi0qRCjcHkQZTZDUMw4tKkQo3B5EGU2Q1DMOLSpEKNweRhjmTGrj8uIk01NchQEN9HZcfN9EC1IZR4VSKVLhlMWVgzqQGcwiGYfShUqTCzUEYhmHkQMu4lrJzCKnYFJNhGIbhiTkIwzAMwxNzEIZhGIYn5iAMwzAMT8xBGIZhGJ6YgzAMwzA8MQdhGIZheCKqGrUNgRGRdcDbIRx6OPD3EI4bJqVmc6nZC2ZzMSg1e6H0bB4ODFLVEdnuWFIOIixEZJmqNkdtRzaUms2lZi+YzcWg1OyF0rM5H3ttiskwDMPwxByEYRiG4Yk5CIfrojYgB0rN5lKzF8zmYlBq9kLp2ZyzvRaDMAzDMDyxEYRhGIbhiTkIwzAMw5OKdBAicqKIvCwiO0TEN/1LRN4SkZUi8oKILCumjR62BLX5GBF5TURWi8iFxbQxxY7dROQvIvKG+76rz3Y97vV9QUT+VGw7XRvSXjMR2UVE7nDXPy0iY4tvZR97Mtn7NRFZl3RdvxGFnUn23CgiH4jISz7rRUSudX/PChGZXGwbPWzKZPPhIrIh6RpfXGwbU+wZLSKPiMgq9z7Rr7VdTtdZVSvuBewHjAceBZrTbPcWMDxqe4PaDFQDfwPGAQOBF4HPR2Tvz4AL3c8XAlf4bLcp4uua8ZoB3wb+n/v5FOCOmNv7NeCXUV7XFHumA5OBl3zWfxm4HxBgCvB0Cdh8OPDnqO1MsmckMNn9PAR43ePvIuvrXJEjCFVdpaqvRW1HNgS0+WBgtaquUdVtwO3A7PCt82Q2cLP7+WZgTkR2ZCLINUv+LQuAo0REimhjMnH6Nw6Eqi4BPkyzyWzgd+qwFKgXkZHFsc6bADbHClXtVNXn3M8bgVVAaq/krK9zRTqILFBgsYgsF5GzozYmAA3A2qTv7fT/IykWn1bVTnD+eIE9fLarFZFlIrJURKJwIkGuWe82qrod2ADsXhTr+hP03/h4dxphgYiMLo5pOROnv9tsmCoiL4rI/SIyIWpjErhToJOAp1NWZX2dy7YntYj8FdjTY9X/UtVFAQ8zTVXfFZE9gL+IyKvuk0UoFMBmr6fa0PKY09mbxWHGuNd4HPCwiKxU1b8VxsJABLlmRb2uGQhiy73Abar6iYh8C2f0c2ToluVOnK5vUJ4DPqOqm0Tky8BCYJ+IbUJEBgN3A99V1Y9SV3vskvY6l62DUNWjC3CMd933D0TkjzjD+9AcRAFsbgeSnxYbgXfzPKYv6ewVkfdFZKSqdrrD2A98jpG4xmtE5FGcJ59iOogg1yyxTbuIDACGEd30Q0Z7VXV90tfrgSuKYFc+FPXvthAk33xV9T4R+ZWIDFfVyET8RKQGxzncqqr3eGyS9XW2KSYfRGSQiAxJfAZmAJ4ZDTHiWWAfEdlLRAbiBFQjyQxyz3uG+/kMoN8ISER2FZFd3M/DgWnAK0Wz0CHINUv+LScAD6sb9YuAjPamzCvPwpmPjjN/Av7NzbKZAmxITE/GFRHZMxGHEpGDce6l69PvFao9AtwArFLVq302y/46Rx19j+IFfBXHm34CvA886C4fBdznfh6HkyHyIvAyzjRPrG3WnZkKr+M8hUdmM84c/UPAG+77bu7yZuA37udDgZXuNV4JnBWRrf2uGXApMMv9XAvcBawGngHGRfy3kMney92/2ReBR4DPRWzvbUAn0O3+DZ8FfAv4lrtegP/j/p6VpMksjJHN5yZd46XAoRHbexjOdNEK4AX39eV8r7NJbRiGYRie2BSTYRiG4Yk5CMMwDMMTcxCGYRiGJ+YgDMMwDE/MQRiGYRiemIMwYo2IbMpi28NF5NA8zlUvIt/2WfeoiMxMWfZdEflVlue4T0TqM2zj+ZtF5CYROSGLc8VKhdYoPcxBGOXE4Ti1FblSj6Pc6sVtOEVpyZziLs+IW5xUpapfVtWuPGzMhrOAf6jq3sA1xL+i2ogZ5iCMkkNEjnWfiJ8Xkb+KyKfdp+NvAd9z9fm/KCIjRORuEXnWfU1z92919f4fFZE1IvId99DzgM+6+1+ZctoFwFeSKr/H4hQpPiEig0XkIRF5Tpz+IbMT24ijz/8rHO2e0eL0GBnurl/oCkG+nCoGKSJXucd7SERGeFyDA0XkMXf/B8VblTNOKrRGKRJ1xaK97JXuhUe/CGBXdvZT/wZwlfu5FfhB0nZ/AA5zP4/BkSFIbPcUsAswHEcioQYYi4/+v7tfGzDb/XwhcKX7eQAw1P08HKfiWtzj7QCmJB3jLdweI+ysLq/DkXHZ3f2uwOnu54txezsAN+FIfdS49o9wl58M3Ohh70tAY9L3vxGT/ib2Ko1X2Yr1GWVNI3CH+9Q8EHjTZ7ujgc8nPTQPTehrAW2q+gnwiYh8AHw6wHkT00yL3Pcz3eUC/KeITMdxCA1Jx3tbHe19L74jIl91P4/GUQNd7x7jDnf5LUCq8Np4YH8chWFwmgh5aeqUokqqESPMQRilyH8BV6vqn0TkcJwRgRdVwFRV3ZK80L2pfpK0qIdg/xcWAleL06qxTt0GLcDpwAjgQFXtFpG3cPSbADZ7Hci1+2jXvo9dJdtar23xlh9/WVWnZrA3Tiq0RgliMQijFBkGdLifz0havhGn3WKCxTiiagCIyAEZjpu6fx9UdRNOy9cb6RucHgZ84DqHI4DPZDhPYp9/uM7hczgtIBNU4UwlAZwGPJGy72vACBGZCo7Ms3g3rImTCq1RgpiDMOLOp0SkPel1Ps6I4S4ReRxI1t+/F/hqIkgNfAdoFqez2is4QWxf1Omj8KSIvOQRpE5wG/AFnFafCW51z7MMZzTxaoDf9QAwQERWAD/FUQRNsBmYICLLcRr9XJpi5zacG/4VIvIijnKnV/bWDcDuIrIaOB8nbmIYgTE1V8MwDMMTG0EYhmEYnpiDMAzDMDwxB2EYhmF4Yg7CMAzD8MQchGEYhuGJOQjDMAzDE3MQhmEYhif/H3v11tXPaLrUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#our implementation looks like a flipped version\n",
    "#kinda like pca\n",
    "z_skl=clf.transform(iris.data)\n",
    "\n",
    "for i in np.unique(iris.target):\n",
    "    plt.scatter(z_skl[:,0][iris.target==i],\n",
    "                z_skl[:,1][iris.target==i],\n",
    "                label=f'Label {i}')\n",
    "plt.legend()\n",
    "plt.title('Factor Analysis by Sklearn')\n",
    "plt.xlabel('Latent Variable 0')\n",
    "plt.ylabel('Latent Variable 1')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimal Number of Factors\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "The optimal number of factors is somewhat similar to PCA. In our case, we already know the number because we merely want to make 2D visualization. For a more general case of latent variable transformation, the optimal could be determined by a couple of different methods. Here, we merely talk about three simple but powerful methods (basically they are related to eigenvalues).\n",
    "\n",
    "* Kaiser Criterion\n",
    "* Scree Plot\n",
    "* Variance Explained Criteria\n",
    "\n",
    "Other methods can be found in Wikipedia\n",
    "\n",
    "https://en.wikipedia.org/wiki/Factor_analysis#Criteria_for_determining_the_number_of_factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get eigenvalues via LAPACK\n",
    "miu=iris.data.mean(axis=0)\n",
    "eigvals,eigvecs=np.linalg.eig((iris.data-miu).T@(iris.data-miu)/iris.data.shape[0])  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Kaiser Criterion\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Don't get overwhelmed by the buzzword. The criterion can be described within 10 words, drop all components with eigenvalues smaller than 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The optimal number of factors suggested by Kaiser Criterion is 1\n"
     ]
    }
   ],
   "source": [
    "print(f'The optimal number of factors suggested by Kaiser Criterion is {eigvals[eigvals>1].shape[0]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Scree Plot\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Scree plot is plotting eigenvalues against the number of features. Then we use elbow method to determine where the steepest drop is.\n",
    "\n",
    "More details of elbow method can be found in the link below\n",
    "\n",
    "https://github.com/je-suis-tm/machine-learning/blob/master/k%20means.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#dual axis plot\n",
    "def dual_axis_plot(xaxis,data1,data2,fst_color='r',\n",
    "                    sec_color='b',\n",
    "                   x_label='',y_label1='',y_label2='',\n",
    "                   legend1='',legend2='',grid=False,title=''):\n",
    "    \n",
    "    fig=plt.figure()\n",
    "    ax=fig.add_subplot(111)\n",
    "    \n",
    "\n",
    "    ax.set_xlabel(x_label)\n",
    "    ax.set_ylabel(y_label1, color=fst_color)\n",
    "    ax.plot(xaxis, data1, color=fst_color,label=legend1)\n",
    "    ax.tick_params(axis='y',labelcolor=fst_color)\n",
    "    ax.yaxis.labelpad=15\n",
    "\n",
    "    plt.legend(loc=3)\n",
    "    ax2 = ax.twinx()\n",
    "\n",
    "    ax2.set_ylabel(y_label2, color=sec_color,rotation=270)\n",
    "    ax2.plot(xaxis, data2, color=sec_color,label=legend2)\n",
    "    ax2.tick_params(axis='y',labelcolor=sec_color)\n",
    "    ax2.yaxis.labelpad=15\n",
    "\n",
    "    fig.tight_layout()\n",
    "    plt.legend(loc=4)\n",
    "    plt.grid(grid)\n",
    "    plt.xticks(xaxis)\n",
    "    plt.title(title)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get distance between a point and a line\n",
    "def get_distance(x,y,a,b):\n",
    "    \n",
    "    numerator=y-x*a-b\n",
    "    denominator=(a**2+1)**0.5\n",
    "    \n",
    "    return np.abs(numerator/denominator)\n",
    "\n",
    "#create line equation from two points\n",
    "def get_line_params(x1,y1,x2,y2):\n",
    "    \n",
    "    a=(y1-y2)/(x1-x2)\n",
    "    b=y1-a*x1\n",
    "    \n",
    "    return a,b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAElCAYAAACxnHbGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4VFX+x/H3Se+BCb0GpLcQwAoqI6N0UJe1DVn1p0tcK0pVdLHTQde1xLWtcbCuSlV0dMAGSAmh9wQIHQYCJEBIcn9/3ImOMZlMQiZ3knxfzzOPM3PbNyHOZ869556jNE1DCCGE8DcBRhcghBBClEQCSgghhF+SgBJCCOGXJKCEEEL4JQkoIYQQfkkCSgghhF+SgBIeKaXuUkr95PZaU0q1MbKmylKZP4tSKlMpZSllWbhSaoFSKlsp9WllHK+6UUo9rZT6wOg6RPUiASWKPlzPKqXOuD3+bXRd8FtAakqp2cXev9H1/nte7mepUupenxRZthFAQyBO07S/VsYOlVJPKKUyXP9WWUqpjytjvxdZ03tKqTxXTU6l1LdKqQ4V2E+pYS9qFwkoUWSopmlRbo8HjS7IzS7gVqVUkNt7fwO2G1RPebUEtmuall/eDYv9zEXv3QkkARZN06KAXsB3lbHvSjDdVVMz4Ajwng+OIWoJCShREYOUUruVUseUUjOUUgEASqkApdSTSqk9SqkjSqn3lVKxrmX/VUqNcT1v6mr93O963cb1jVuVcrxDwAagv2t9E3AVMN99JaXUFUqpX5RSJ5VS6Uqpvq73XwCuBv5dQuvQopTaoZQ6oZR6tagGTz+La3mSa9lxpdSk0n5RSqlngH+iB+wZpdQ9Zfye4l2/m3uUUnuB70vY7aXAEk3TdgFomnZI07Q33Y5pUkq9q5Q64Pq5vnS939fV2pqglDoEvOt6f4hSap3r9/aLUqqb276aKKX+p5Q66mqxPVzaz+pO07RcYC7QpZTfyzCl1CbXMZcqpTq63k8FWgALXL+v8d4cT9RMElCiIm5C/9beAxgO/J/r/btcDzPQGogCisJgGdDX9fxaYLfrvwDXAD9qnsfdeh+91QRwGzAPOF+0UCnVFFgEPA+YgLHA/5RS9TVNmwT8CDxYQutwCPoHfgJwC64Q9PSzKKU6Aa+jt2KaAHHoLYY/0TRtMvAi8LHr2G+X8Xsqci3Q0a0edyuAvymlximleimlAostTwUigM5AA2CO27JGrt9PS2CUUqoH8A6Q7Po5UoD5SqlQ1xePBUA60BToB4xWSpVU0x8opaIAK5BWwrJ2wIfAaKA+sBg9kEI0TUsC9vJ7i356WccSNZimafKo5Q8gEzgDnHR7/N217C7gJ7d1NWCA2+v7ge9cz78D7ndb1h64AAQBl7j2GwC8gf6BmOVa77/AY6XUdhfwExAOHAZi0T+ge6OH0Xuu9SYAqcW2XQLc6Xq+FLi32HIN6OP2+hNgohc/yz+Bj9yWRQJ56KfcSvoZngY+cHvtad/xrrpal/FvZgXsQA5w3K3uxkAhULeEbfq66gxze+914Lli621DD8jLgb3Flj0OvFtKTe8B51z/zofQW7iXFP8dAE8Bn7htFwDsB/q6/T2W+LuUR+16+OIctKiebtQ0ze7luvvcnu9Bb0Xg+u+eYsuCgIaapu1SSp0BuqOfbnsOuEcp1R79w/Bfng6oadpZpdQi4EmgnqZpPyulBrqt0hL4q1JqqNt7wYCjjJ/lkNvzXPTWjMefxbXst9+Bpmk5SqnjZRzHnad9F9mHB5qm2QCbUioYuNH1PA04ATg1TTtRyqZHNU075/a6JXCnUuoht/dCXDUWAE2UUifdlgWit0ZLM1PTtCc91U6xn1/TtEKl1D70VpoQv5FTfKIimrs9bwEccD0/gP6B574sH73lA/ppvhFAiKZp+12v/wbUBdZ5cdz3gTHop7CK24fegqrj9ojUNG2qa3l5h+339LMcxO13oJSKQD89Vhn7LuJVvZqmXdA07VNgPfr1nn2ASSlVp7RNir3eB7xQ7PcWoWnah65lGcWWRWuaNsib2jz4w8/vuu7XHL0VVVKNopaSgBIVMU4pVVcp1Rx4BCjq4vwh8KhSqpXrGkTRtZei3mvLgAeBH1yvlwIPoZ9CLPDiuMuA64FXSlj2ATBUKdVfKRWolApzdQooujZ0GP16j7c8/SyfAUOUUn2UUiHAs5Tv/6Wyfk8eKb3r/WClVLSrw8VA9OtNKzVNOwh8Bbzm+jcKVkpd42F3/wHuU0pdrnSRRfsGfgVOuTpVhLt+r12UUpeW42ctySfAYKVUP1cLcAz69cRfXMvL+28laigJKFGkqNdU0eMLD+vOA9agt3oWAW+73n8HvXXzA5CBfj3C/dTRMiCa3wPqJ/SL+T/gBU33naZpzhKW7UPvsPEEcBT92/84fv8bfxkY4erV5vF0Ylk/i6Zpm4AH0HupHUQ/rZblzc9Q1r69dAr959yLfr1nOvAPTdOKbqhOQr+mtRW9q/fo0nakadpq4O/onTROADvRr/vh+tIwFP20bAZwDHgL/TpghWmatg0Yif5F45jrGEM1TctzrTIFeNLVw2/sxRxLVG9K06Q1LYQQwv9IC0oIIYRfkoASQgjhlySghBBC+CUJKCGEEH5JAkoIIYRfqrYjSQQEBGjh4eFGlyGEENVGbm6upmlatWmYVNuACg8PJycnx+gyhBCi2lBKnTW6hvKoNkkqhBCidpGAEkII4ZckoIQQQvglCSghhBB+SQJKCCGEX5KAEkII4ZckoIQQQvilWhVQmqbx/YHNZJ4+anQpws2xE+dxnjxvdBlCCD9TqwLKeT6HpKWvMX39QqNLES75+YU88uxq7p6wgl17TxtdjhDCj9SqgIoLi+L/2vflk4yV7Mg+ZHQ5Avhu+SEOHDlLQUEhj72wlp17JKSEELpaFVAAD3e6gbCAYGZuWGR0KbVeYaHG3PmZtG4RRcrzlxMaEsCYFyWkhBC6WhdQ9cNjuLdDXz7LWMW27INGl1Or/bjqCHsP5GIdFk/TRhHMebLn7yGVKSElRG1X6wIK4KFO/YkMCmF6ulyLMoqmadjmZ9K8cQTXXt4QgKYN9ZAKCw1gzBQJKSFqu1oZUHFhUfy9g5kv96xh84n9RpdTK61MP86OzNPcPjSewAD12/vFQ2pH5ikDqxRCGKlWBhTAA52uJzI4VHr0GUDTND74MoOGcWFc37vRn5Y3aaCHVHhoIGNelJASorbyq4AypSYHmlKT00ypyT5PDVNoFPd16Mf8vWvZ6Nzn68MJN+lbT7JpRza3DmlJUFDJf4JNGkQw+8keRIQFMebFtWzPkJASorbxq4ACHgG2VNXB7u/Yj5jgcGlFVTHblxnUjQ1hUN8mHtcraklFhAUxdoqElBC1jd8ElCk1uRkwGHirqo5ZJzSSf3Tsx8J960g/vreqDlurbd2VzeqNTv46sAWhIYFlrt+4QbgeUuESUkLUNn4TUMBLwHigsCoP+o+OFmJDIpi2fkFVHrbW+mBeJtGRQQy3NPN6m8YNwnnJFVJjXlzLNgkpIXzCbLW/Y7baj5it9o1lrHep2WovMFvtI3xZj18ElCk1eQhwxJmUssbTekqpUUqp1Uqp1fn5+ZVy7JiQcB7oaOHrrPWkHc+slH2Kku3ee4af1xzl5v7NiQgPKte2jerrIRUZEcTYF9eybbeElBA+8B4wwNMKZqs9EJgGLPF1MX4RUEBvYJgpNTkT+Ai4zpSa/EHxlTRNe1PTtF6apvUKCirfB5wnozpcR92QSKamSyvKl2zzMwgPC+Tm/i0qtH1RSEVF6qf7JKSEqFwOm+UHwFnGag8B/wOO+LoevwgoZ1LK486klGbOpJR44Dbge2dSysiqOn5MSDgPdb6Bb/dvZNXR3VV12Fpl/6Fclq44zLB+zYiJCq7wfhrVD2fOJAkpIYxgttqbAjcBb1TF8fwioPzBve37EhcaJdeifGTugkwCgwK4ZVDFWk/uikIqOjKYMVPWsnVXdiVUKEStEFR0mcT1GFXO7V8CJjhslgJfFFec3wWUMyllqTMpZUhVHzcqOIyHOt/A9wc2s/LIrqo+fI125Pg5vvnxIIP7NsFUJ7RS9tmofjizJ/UgJjKYsVPT2CIhJYQ38osuk7geb5Zz+17AR2arPRMYAbxmttpvrPQqXfwuoIx0T7u+1A+LlmtRleyjhXvQgNuGxFfqfhvV17ugx0QGM27KWgkpIXzMYbO0ctgs8Q6bJR74DLjfYbN86avjSUC5iQwO5eHO/Vl2aAu/HN5udDk1gjP7PIsc+7m+TyMa1gur9P03rBfGnCd7EhsdIiElxEUyW+0fAsuB9marPctstd9jttrvM1vt9xlRj9I0zYjjXrTIyEgtJyen0vebm59Hzy8n0TamEfNvGFPp+69t/vPRTj5amMl7M66keeNInx3nyPFzPPr8GrJP5zF9Yg86tYn12bGEqK6UUrmapvnuf8RKJi2oYiKCQhjdeSA/Hd7Oj4e2GV1OtXY65wJffruPay9v6NNwAmgQ52pJxYQwfupaNu+QlpQQ1Z0EVAnubHc1jcPrMDV9PtW1hekPvvhmH7nnCrhjWHyVHK9BXBhzJukhNU5CSohqTwKqBGGBwTzadQDLj+xk2aGtRpdTLZ09l8//vt7HlYn1aNMyusqO2yAujJcm9aROrISUENWdBFQpktr0oUlEXaamL5BWVAUs+G4/p85cYOTwVlV+7PqukKrrCqlNO05WeQ1CiIsnAVWK0MBgxnQdyK9Hd/H9wc1Gl1Ot5OUV8PHiPSR2rkuntsZ0VqjvuiZVNzaE8VPT2LhdQkqI6kYCygPrJb1pFmliyjq5FlUeX/1wEOfJPENaT+7qm9xCapqElBDVjQSUByGBQYzpOoi1xzOxH/A4+rxwyc8v5MMFmXRqE0tip7pGl/NbSMXVCZWQEqKakYAqwx2XXEXLqHrSivKS/ZdDHD52jpHD41FKGV0OoIfU7Ek9fgupDdskpISoDiSgyhAcEMjYroNY59zL11nrjS7HrxUUasydn8klLaK4IrGe0eX8QVFLql7dUCZMl5ASojqQgPLCra2voFV0fenRV4YfVx1h38FcrMNb+U3ryV29uqHMnqSHlLSkhPB/ElBeCAoIZFzXwWw4sY9F+9YZXY5f0jQN27wMmjeO4JrLGhhdTqmKQqq+SUJKCH8nAeWlEa0uo01MQ6amL6BQKzS6HL+zMv04O/ec4Y6h8QQG+F/ryV29uqHMcQup9VtPGF2SEKIEElBeKmpFbT65nwV704wux69omsYHX2bQMC4MS+9GRpfjlThXSDWIC2XC9HWkb5GQEsLfSECVw83xl9I2phFT0xdQUCitqCLpW06waUc2tw5pSVBQ9fmTcg+piTMkpITwN9Xn08QPBAYEMCFhCNuyD/LlnjVGl+M3PpiXSd3YEAb1bWJ0KeVmquMeUmmsk5ASwm9IQJXTjS170iG2CdPXL5RWFLBlVzZrNjq5ZVALQkMCjS6nQopCqmFcOI9LSAnhNySgyilA6a2oHacO8XnmKqPLMdwHX2YQHRnEsH7NjC7lopjqhDL7yR6/h9Rmp9ElCVHrSUBVwNAWiXSu24zp6xeSX1hgdDmG2b33DL+sPcZf+rcgIjzI6HIuminWFVL1wnl85joJKSEMJgFVAQEqgAndhrDr9BE+y/jV6HIMY5ufQXhYIDf1b250KZXGFBvK7El6SE2csY60TRJSQhhFAqqCBjfvTjdTc2ZsWMSFWtiKyjqUy9IVhxluaUZMVLDR5VSqopBqVF9vSUlICWEMCagKUkoxsdtQMk4f5ePdK4wup8rNnZ9JUFAAfx3UwuhSfMIUq3ecaNxAD6m1ElKiFjBb7e+YrfYjZqu9xOkbzFa71Wy1r3c9fjFb7Qm+rEcC6iL0b9aNxLiWzFxfu1pRh4+d45ufDjLY3ARTbKjR5fhM3dgQZj/RkyYNwnli5jrWbJSQEjXee8AAD8szgGsdNks34DngTV8WIwF1EZRSTEwYyt6c48zd9YvR5VSZjxftAeDWwfHGFlIF6saGMEtCStQSDpvlB6DUP3KHzfKLw2Ypug9jBeDT7rsSUBfJ0qQLPeu1YtaGxeQV5Btdjs85s8+zyLGf6/s0omG9MKPLqRJ1Y0OYPaknTRtKSIlqL0gptdrtMeoi9nUP8FVlFVYSCaiLVNSKyspxYtv1s9Hl+NxnX+0lP7+QO4bGG11KlaoTo4dUs0YRrpA6bnRJQlREvqZpvdweFTpFZ7bazegBNaFyy/sjCahKcF3jTlxW/xJmbfiKcwUXjC7HZ07nXGDet1lce3lDmjeONLqcKlcnJoRZT/RwhVQ6qzdISInax2y1dwPeAoY7bBaf/k8gAVUJlFI8njCUA7knSN3xk9Hl+MznS/aRe64A6/B4o0sxjN6S6kGzxhFMmpXOqvUSUqL2MFvtLYDPgSSHzbLd18dT1XWG2MjISC0nJ8foMn6jaRpDvplJ5uljrL7xOcKDQowuqVKdPZfPbY/8TJd2sbwwprvR5Rgu+3Qej724lqyDuTz/WAKXdoszuiQhyqSUytU0rdTTH2ar/UOgL1APOAxMBoIBHDbLG2ar/S3gL8Ae1yb5Dpull8/qlYCqPD8d2sawb2fzYq9buK9jP6PLqVQfL9rDG3N38Oozl9KpTazR5fiF7NN5jHlxLXsP5vKChJSoBsoKKH8jp/gqUZ9G7bm6YXte2vg1ufl5RpdTafLyCvhk8R56dDZJOLmJjdavSbVoHMGk2en8Kqf7hKhUElCVbGLCUI6cO8W725cZXUqlWbzsAM6TeYysxdeeSlMUUi2bRPDk7HR+TT9mdElC1Bhyis8Hbvr2JTadyCLtpheIDK7eIy3k5xcycswv1KsbyiuTe6GUMrokv5R9Oo+xU9ay50Auzz/ajcsS6hldkhB/YuQpPrPVfgpQgAZEAGddzwEiHTbLnyaUkxaUD0xMGMqx86d5e/tSo0u5aPafD3H42Dmsw1tJOHkQGx3CzMd7EN80kifnrJeWlBDFOGyWGIfNEu2wWWKA9UXPXa/TS9pGAsoHLm9wCdc16cS/Ni3h9IVzRpdTYQWFGnMXZNKmZRRXdJcOAGX5Q0jNTmflOgkpIUpjttrdJ5ErcUoECSgfmdhtGM7zOby11WF0KRX246oj7DuYK62ncoiJCtZDqlkUT81JZ4WElBAl+QH4zGy1J5ut9g+AzSWtJAHlI73qt+KGpl15ZfM3nMo7a3Q55aZpGh98mUHzxhFcfWkDo8upVtxD6p9z0lmeJiElRDFj0cfx6wqsBu4saSXpJOFDaccz6bd4Ck8kDGNst8FGl1Muy9OO8cTMdUwY1YkB1zYxupxq6dSZC4ydspbMrDM8MzqBKxOl44QwltH3QZmt9o5AP/TOEd85bJatntaXFpQPJcbFM7BZAq9usZOdl2t0OV7TNA3bvAwa1gvD0ruR0eVUWzFRwcx6ogetmkcx+aV0lq89anRJQhjGbLX/FfgSaAw8Dkw3W+1WT9v4RUCZUpPDTKnJv5pSk9NNqcmbTKnJzxhdU2WZkDCE7Lxc3tjyndGleG3d5hNs2pHNbUNaEhTkF38i1VZ0pH66r3WLaP750noJKVGbPQH0cdgsk4AjwE3Aw5428JdPn/PAdc6klASgOzDAlJp8hcE1VYpuphYMad6d17bYOXnev09JFvlgXiamOiEMklN7lSI6MpiZExO5pKUeUr9ISInaKcBhsxT98SuHzVJAKb33ftvA9zWVzZmUojmTUs64Xga7HtXz4lgJJiQM5fSFc7y6xW50KWXavDObtZuc3DKoJSEhf7pvTlRQlFtITZaQErVTntlqr+t6Hma22l8FVnrawC8CCsCUmhxoSk1eh970+9aZlOKx8Oqkc91mDG/Zk5St3+M8f6bsDQxkm5dBTFQww/o1NbqUGqcopNq4QurnNRJSolZ5AIhyPZ8LLAfu97SB3/XiM6Um1wG+AB5yJqVsdF/mmp54FEBISEjP8+fPG1BhxWw5eYA+C55ldJf+PJV4k9HllGjX3tPc+/hK7vpLa+68ubXR5dRYZ3IuMG5aGjszTzP54a706SXd+EXVMLoXX3n5XUABmFKTJwM5zqSUmaWtUx26mRd3749vsSRrPWk3vUC9sGijy/mTZ1/ZwMp1x/joX32IjvR4alhcpDO5+YybulZCSlQpPxqLLxwIBHIcNkupH4Z+cYrPlJpc39VywpSaHA5YAI/946uj8d2GcLYgj39v/sboUv5k38Eclq48zPDrm0k4VYGoiCBmTOxBm/honv7XBn5afcTokoTwqWJj8YUAfwVme9rGLwIKvV+8w5SavB5YhX4NaqHBNVW6drGN+Ev8pby1dSlHzp4yupw/+HDBHoKDAhgxsIXRpdQaRSHVTkJK1DIOm0Vz2Cz/A27wtJ5fnuLzRnU8xQew89Rhrpg/mfs69OP5Xn81uhwAfbTyx35mWL+mPHxnB6PLqXXO5OYzYVoa2zJOMfmhrjK0lPAZf7sGZbbaewPLHTZLYUnLJaAM8MDP7/HFntWsvfEFGkUYP0Ptv/67lfnf7cc2uzcN64UZXU6tlJObz3hXSP3zoa5cIyElfMCIgDJb7bvRrz2VyGGztCptWVBpC4TvjO02mE8yVvLypq+ZcumthtbizD7PIscBbujTWMLJQJERQUyfkMj4aWk8+8oGCSlhCLPV/g4wBDjisFm6lLBcAS8Dg4Bc4C6HzbK2jN32cv33AfRBGT5GD6zbgDqeNvSXa1C1Sqvo+tzW+kre2/4DB3JPGFrLp4v3kp9fyB3D4g2tQ/weUh1ax/DsKxtY9utho0sStc97wAAPywcCbV2PUcDrZe3QYbM4HTaLExjssFmmO2yWPQ6bJdNhs0wFzJ62lYAyyNiugyjQCnl549eG1XDqzAXm2bPoe0VDmjWKMKwO8bvIiCCmjS8KqY0SUqJKOWyWHwCnh1WGA++7OjmsAOqYrfbGXu5ema32kWarPcj1SKKMEYMkoAzSMroed1xyFf/d8RNZOZ7+HnzniyX7OHuuAOuwUk8BCwMUtaQ6XuIKqZUSUqLSBCmlVrs9RpVz+6bAPrfXWa73vHE7etfyA8BB4BbgDo/FlrM4UYnGdB3Eh7uXM2fjV8y63OOo85Uu92w+/1uyl6t61KN1i6iyNxBVKiJcD6kJ09N49t8b+Sdw7eUNjS5LVH/5mqb1Knu1UpXU2cGrnnYOm2U3egvMaxJQBmoeFUdSmz6k7vyJ0Z0H0DwqrsqOPf+7LE7n5DPyRmk9+auIcP10X1FIPQX0lZASxsoCmru9bobeIiqT2WqPQ+9g0R8oBOzAIw6bpdQpp+UUn8Ee7TIQhWLWhsVVdszzeQV8sngvPbuY6HiJ8d3cRemKQqpz21ie+/dGHCvkdJ8w1Hzgb2arXZmt9iuAbIfNctDLbV8F1qGfEtzvev2Gpw0koAzWNLIud7btw9xdv7DndKlfJCrV4qUHOJGdx8jh8VVyPHFxIsKDmDquO53bxvL8qxv5fvkho0sSNZTZav8QfZTx9marPctstd9jttrvM1vt97lWWQzsBnYC/6GM0ciL6eiwWWY6bJY89PmgfuGPrbE/kVN8fmB0l4G8v+MnZm5YzCtX/c2nx8rPL+TjhXvo3DaWhI51y95A+AW9JdWdidPX8cJrmwC47spGBlclahqHzXJ7Gcs19PuZKuIPE8yZrXaP4QTSgvILjSPqcFe7a/ho93J2n/LteGz2nw9x+Pg5rMNboVSpN3cLPxQeFsTU8d3p0jaWF6QlJaqfH8xWe4LreRywBBjtaYNyD3VkSk0OAKKcSSmGjnZanYc6Ksnhs9n0+GISw1v25LXed/vkGAWFGnePX05oSABvvnC5BFQ1dfZcPhNnrGPjtpNMeqCLtKSE1/xlLD6z1R7psFnK/AD3KqBMqclzgfuAAmANEAvMdialzLjYQiuqpgUUwFNrPuP1LXZWDHuGNjGV31vLseIwz76ygckPd5XeYNWce0g9cX8X+l0lISXKZvB8UJNLet9hszxT2jbenuLr5Gox3Yh+kawFkFTuCoVHD3e6gbCAYGasr/yZRjRNwzYvg+aNI2S07BogPCyIqeMS6dqhDi++thH7z3K6T/i9026PPOB69CwplbcBFWxKTQ5GD6h5zqSUC3h5c5bwXv3wGO7t0JfPMlaxLdvbnpveWZF2jF17z3DHsHgCA+TUXk0QHhbIlLF6SE15XUJK+DeHzTLb7TEF6Au09rSNtwGVAmQCkcAPptTkloB/zbhXQzzUqT8RQSHMWL+o0vapaRofzMukUf0wLHIqqEYpCqluHeoy5fWNfPtz5X6xEcKHYqmMbubOpJR/Af9ye2uPKTXZ4yi0omLiwqIY1cHMSxuX8FiXgXSq6+0wV6VL23yCzTuzGX13B4KCpONmTRMeFsiLY7szadY6pr6+CTS4vo+343cKUTXMVvt69KGSNPQu5w2BUq8/gZctKFNqckNTavLbptTkr1yvOwF3Xly5ojQPdLqeyOBQplfStagP5mUQVyeEgdfIh1ZNVRRSCR3rMvWNTXz7k7SkhN8ZAgx2/fd6oKHDZnnF0wbefp1+D73PehPX6+2U0X9dVJwpNIrkDtcxf+9aNp3Iuqh9bd6RTdqmE9wyqCUhIYFlbyCqrbBQV0h1qssUCSnhZxw2y170+59uQh/VvHtZ23gbUPWcSSmfoA/whzMpJR+9y7nwkQc6WogODmNa+oKL2s8H8zKIiQpmaL+LP1Uo/F9YaCAvjulOd1dIffOjhJTwD2ar/QngXfRZdOsC75it9kmetvE2oHJMqclxuHrumVKTrwCyL6JWUYY6oZHc39HCwn3rWO/cW6F97NxzmuVpx/jLgOaEh8moVrVFUUgldjIxNWUTS370arBpIXztb8DlDpvlGYfN8jRwGWXcruRtQD2GPortJabU5J+B94GHKl6n8MY/OlqIDYlgWnrFrkXNnZ9JRFggN91Q5pBXooYJCw3khTEJ9OhsYlrKZr7+QUJKGC4TCHF7HQJ4/Pbt9VBHptTkIKA9ei+Mba57oQxTE0eSKMnM9Yt4MX0+3w16nMS4eK+323cwhzvHLee2IS0ZdVtb3xUo/Nr5vAImzUpn7SbIbjCJAAAgAElEQVQn40d1YsA1TcreSNRYBo8k8RZwDVB03WIYsArYBiWPKOFtL76/oU/N2xPoAdzuek/42KgO11E3JJKp5bwWNXdBJsFBAYwY6PFGbVHDhYb83pKa/uZmvl4mLSlhmM3o8z/tdz1eB1bz++gSf+LthYlL3Z6HAf2Atein+oQPxYSE82Cn63lu3ZesPppBr/plz4B76OhZvv3pEMP6NcMUG1oFVQp/VhRST85KZ/p/NqMBA6+VlpSoWg6bZXZ5tyn3aOYAptTkWCDVmZQyrNwbV5LacooP4PSFc/T4YhIJcS34rN8jZa7/8ntbWfj9fmxzetMgLqwKKhTVwfm8Ap6cnc6ajU7G3tuRQX2lZ2dtY/ApvnfQLxH9gcNmudtstT/jsFn+NJhsRYcVyAXkwkYViQ4O46HON/D9gc2sPLLL47rOk+dZtPQAN1zdWMJJ/EFoSCDPP5ZAr65xzHxrC4uX7je6JFG7LES//lT8AbCspA28nW5jAb8PDhsAdAI+cSalTLzIgiusNrWgAHIunKfHl5PoVKcZX1xf+j3Sb8zdwaeL9/D+zKto2iiiCisU1UVeXgFPzlnP6g3HpSVVyxg9H5TZag8GOqDnyTaHzeKxs523LaiZwCzXYwpwjZHhVBtFBofycOf+LDu0hV8Oby9xnVNnLjD/uyz6XtFQwkmUKiQkkOcf7UavrnHM+M8WFklLSlQBs9XeE9gCvIbeQWKr2Wq/1NM23g4WW2LzS1Stu9tdy783f8PU9AXMv2HMn5Z/vmQvZ88VMHJ42R0pRO1WFFJPzVnPzP9sAQ0Gm6UlVduZrfYBwMvog7m+5bBZphZb3gL4L/poEIHARIfNstjL3f8buNVhs6xx7asH8ArQu7QNPAaUKTX5NCXP+6QAzZmUEuNlYaISRASF8EjnATyx+hN+PLSNqxu1/21Z7tl8Pl+yj94969OqeZSBVYrqIiQkkOce7cY/X1rPzLe2oGkw5DoJqdrKbLUHAq+iD+SaBawyW+3zHTbLZrfVngQ+cdgsr5ut9k7oE9jGe3mI8KJwAnDYLGvNVrvH040eA8qZlBLt5YFFFbmz7dW8sukbpqbPp0/DsSild4qZZ8/idE4+I4fHG1ugqFZCQgJ5dnQ3/vnyema9vQWQkKrFLgN2OmyW3QBmq/0jYDj6/UtFNKCoYRILlOfGulyz1R7lsFnOuPYfBZzztEG5BmgzpSY3QL8PCgBnUkrFBokTFRYeFMKjXQcw/tePWHZoK30bd+R8XgGffrWXXl1MdLgk1ugSRTUTEhLIs4/8HlIaGkOva2Z0WaLqNQX2ub3OAi4vts7TwDdmq/0h9AlsLeXYvwU47/b6LPo9taXydiSJYabU5B1ABnp3wEzgq3IUJipRUps+NImoy9T0BWiaxuKlBziRnYdVWk+igvSWVAKXd49j9ttbmf/dxU3zIvxWkFJqtdtjlNuyP92jxJ8v8dwOvOewWZoBg4BUs9XuVY44bJZc4Aqz1W41W+13AiOBEZ628bYX33PAFcB2Z1JKK/TU+9nLbUUlCw0MZkzXgfx6dBffZm3i44V76Nw2loSOdY0uTVRjIcEBPDs6gSu612POOxJSNVS+pmm93B5vui3L4o9TsDfjz6fw7gE+AXDYLMvRz6jV8+bAZqv9A2AacCX6sHk9gV6etvH2FN8FZ1LKcVNqcoApNTnAmZTiMKUmT/NyW+ED1kt6M2fj10ybt5zC460Y/X8dfrseJURFhQQH8Mzobjz98nrmvLOV8LBAru8tMzHXEquAtmarvRX6WHm3oY/B6m4vegPlPbPV3hE9oI56uf8EoJvDZvF6+CJvW1AnTanJUcAPgM2UmvwykO/tQUTlCwkM4tHOAzmbXo+GTYK4PCHO6JJEDRESHMDTD3ele8e6TEvZzIp1x4wuSVQBh82SDzyIPnv6FvTeepvMVvuzZqu9aFi7McDfzVZ7OvAhcFc5AmcjeqvMa96OJBGJfkErALCi996wOZNSjpfnYJWpto0kUZJvfznAi69uRl29g++S75MWlKhUObn5PPrCGvYeyGHGxB50bV/H6JLERTJ4LL7v0U/r/Ypb7z2HzTK0tG28PcU3CvjUmZSShX6TljCYpml8tGAvsfUUP0SsY0nWegY0TzC6LFGDREYEMW18Io88t5onZq7j5ad60bqF3GMnKuzp8m7gbQtqMnAL4AQ+Aj5zJqUcLu/BKlNtb0H9svYok2alMz65I086XycqKIylgydJK0pUukNHz/LQM6vRNI1/Te5FkwYyjFZ1ZfRYfOVVruk2TKnJ3YBbgb8AWc6klPL0ga9UtTmgNE3jgcmrOHEqj9SZV/HZnpXc/8t7vH/tfQxpkWh0eaIGysw6wyPPrSEqMohX/tkLUx2ZZ6w6MvgU325K78quHDbLn8ZoK+90G0eAQ8BxoEG5KxSVYu2mE2zZdYrbh8QTFBTAiFaX0SamIdPWL6RQKzS6PFEDxTeLYsq47jhP5jF+WhpncjwOQi1ESXrxe/dy90fR+3/i1TUoU2ryP9BbTvWBz4C/O5NSNnveynum1OTm6LPzNgIKgTedSSkvV9b+axrbvAzi6oQw4Bq9+29QQCDjug4m+ed3WLA3jeEtS/y3FuKidGoTy7Oju/HEzHU8MSud6RMSCQsNNLosUX2cAQYA2Q6bxasByL1tQbUERjuTUjo7k1ImV2Y4ueQDY5xJKR3Rbwh+wJSa3KmSj1EjbNpxkrTNJ7hlcEtCQn7/cLg5/lLaxjRiWrq0ooTvXNotjifu78LG7Sd59pUN5OfL35rw2pfoN/o+Z7banzZb7XXMVvsXnjbwKqBccz9tMKUmNzGlJrcoelRCwUX7P+hMSlnren4avQ++jFhZAtu8TGKighlabEDPwIAAJiQMYWv2Ab7cs6aUrYW4eOYrGvLIXR1YnnaM6f/ZTGGh99exRa3WxGGzDEcfLf1Gh81ykjLui/J2LL4HgcPAt8Ai12PhxdVa6rHigURgZfFlSqlRRWNI5efXvvuEd+45zfK0Y/xlQHPCw/58dvbGlj3pENuEaekLKSiUb7bCd4ZbmvF/I1rz7U+HeN22nfJ0thK11jaz1d7BYbOcBzBb7WG4DT5eEm/vgxoNtPf1jbmu0Sr+h3468VTx5a5xo94EvRefL2vxR7Z5mUSEBXLTDc1LXB6g9FbU3T+8yeeZq/hr6+IDEQtReUbe2IrsMxf47Ot9xMaEyESZoiz1gTSz1b4C/bLRKvTJEUvlbUDtA7IvrjbPTKnJwejhZHMmpXzuy2NVR3sP5LDs18PcPiSe6MjgUtcb2iKRznWbMWPDIm6K70VQgFzEFr6hlOJ+aztOnb7A25/sIiYqmGH9ZJoOUaqn3Z6fA3Y4bJYTnjbwNqB2A0tNqcmLcJvPw5mUMru8FZbElJqsgLeBLZW1z5rmwwWZhAQHMGKg50t/ASqACd2G8Ldlb/BZxq/cdsmVVVShqI0CAhTjR3XidE4+L727ldioYK69vKHRZQk/5LBZfijvNt4G1F7XI8T1qGy9gST0jhjrXO894UxK8Xau+xrt0NGzfPvzIYZbmlE3tuxf/+Dm3elmas6MDYsY0eoyaUUJnwoKCmDyw10ZN3Utz7+6kciIIHp1lcGLxR+ZrfZT6DfqaujXnoKAHIfNUurM7eUdSSLSmZTiF8M31KaRJF56dyuLHPuZO6c39eM8XlP8zdf70rlj6Wu8cuXfsLbp7eMKhYAzORd45Pk1HDh8ltmTetBRZnf2O/4y1JHZalfAzejTb0wubT1ve/FdaUpN3oze/RtTanKCKTX5tUqpVHh0/MR5Fi87QP+rG3sdTgD9m3UjMa4lM9Yv4kJhgQ8rFEIXFRnM9PGJmGJDmDh9HXv2144vkKL8HDaL5rBZ/gfc4Gk9b2/UfQnojz7EEc6klHTgmouqUHjl06/2UpBfyO1D48u1nVKKiQlD2ZtznLm7fvFNcUIUE1c3lOkTEwkMVIybupbDx86VvZGozcZ4mjLe67H4nEkp+4q9JV/LfSz7dB7z7FmYr2xE00blH0Ha0qQLPeu1YtaGxeQV1L77xoQxmjaMYPqERHLPFTBu6lpOnsozuiThpxw2yy8Om6XUmza9Dah9ptTkqwDNlJocYkpNHovrdJ/wnS++2ce58wVYh8VXaPuiVlRWjhPbrp8rtzghPGjTMpopY7tz+Ng5JkxPI/esfEES5edtQN0HPIA+/FAW0N31WvhI7tl8Pl+yj94969OqecUnibuucScuq38JszZ8xfkCGYFaVJ2u7evw9CPd2LnnDE/OSScvT066iPIpVy8+f1LTe/F9uCCTNz/ayevPXkqHi+wNtezgFm6yv8T0y27j3vbmSqpQCO98+9NBXnx9E1f3qs/kh7sSGFjeWX5EZTF4Pqh3PS132Cx3F3/P2+k2/lXC29nAamdSyjzvyhPeOp9XwKdf7aVXF9NFhxPANY06cGWDNszZ8DXWS3oTHuSLW9mEKNn1fRqTfeYCr6ZuZ/Y7Wxl7b0eZ+dlPma32AejDDwUCbzlslqklrHML+qgQGpDusFnu8HL31wDjXdsVN7OkDby9UTcM6AB86nr9F2ATcI8pNdnsTEoZ7eV+hBcWLz3Aiew8rA9VzthmSikeTxjGsG9n898dP3Jfx36Vsl8hvDViQAtOnb5A6pcZ1IkO4e+3tTG6JFGM2WoPBF5FH208C1hlttrnO2yWzW7rtAUeB3o7bJYTZqu9PBPXnnJ1LS/p2E+V9L63AdUGuM6ZlJIPYEpNfh34Bv0H2VCOAkUZLuQX8tHCTLq0iyWhQ51K22+fRu25umF7Xt70NX9rezUR0ooSVezuEa3JPp3H3AWZxEQHc+vglkaXJP7oMmCnw2bZDWC22j8ChgPu8//9HXi1aAw9h81ypBz799QVObykN709GdwUcD9vGQk0cSalFOA2Np+4eN/+dJAjx88zcnirSj8NMjFhKIfPnuLd7V5NZilEpVJK8fBdHeh7eUPemLuDr5YdMLok8UdN0QcGL5LFn+flawe0M1vtP5ut9hWuU4LeWmG22h8u/qbZah8NrChpA29bUNOBdabU5KXoYyldA7xoSk2OBOzlKFB4UFCo8eGCPbSNj+ayhMofy+zKhm25tlFHXt64hLvaXkNkcGilH0MITwIDFE/c35kzOReY+Z/NREcG0adXec4SiYsUpJRa7fb6Tdc0RqB/thdX/HpRENAW6Is+2eCPZqu9i2vywbI8DjjMVvtfgDTXvhOBxkCJvbe8nVH3beAq9Cl7vwT6OJNS3nImpeQ4k1LGebMPUbZlKw+TdSgX6/B4n11EnpgwlGPnT/P29qU+2b8QZQkOCuDZR7vRrnUMz/57I+u2eJxxQVSufE3Terk93nRblgW4TzbXDCjezM0C5jlslgsOmyUD2IYeWGVy2CwHgM7oDZ79rn3PAjq7lv2Jx27mptTkDs6klK2m1OQeJS0vmqbdCDWtm3lhocbfn1hJQYHGO9OuICDAd72cRnz3MuuO7yHtpheJDvZ+fD8hKlP26TweeXYNR53neOmpnrSNjzG6pBrPUzdzs9UeBGwH+qEHyCrgDofNssltnQHA7Q6b5U6z1V4PvSXU3WGz+GQy27JaUGNc/51VwqPEboGiYpanHWP3vjPcMSzep+EEMLHbMJznc3hrq8OnxxHCk9joEGZMTCQ6Mpjx09LYd7DmfOGsjhw2Sz7wILAEfaSgTxw2yyaz1f6s2Wof5lptCXDcbLVvBhzAOG/DyWy17zZb7Rluj91mqz3DtWxhSdvIjbp+QNM07p+8ipOn8vhg1lVVciPjbd//m1VHd5N20wvEhJTYgUaIKrHvYA4PPbOasNBAXpnci/omadX7isE36ppKet9hszjNVnuMw2Y5VXyZx09CU2ryeLfnfy227MWKFir+aO2mE2zddYrbh8RX2V32ExKGcCIvhze3fl8lxxOiNM0bRzJtQiKnzlxg3NQ0Tp2RIblqqF5AqMNmcQJx6B0tzgGUFE5Q9im+29yeP15sWXm6FwoPbPMyiKsTwoBrGlfZMRPj4hnYLIFXt9jJzsutsuMKUZL2rWJ44bEEDhzO5fEZ6zh7Tsbtq4FmAE6z1V4H+Bp9CqdPPG1QVkCpUp6X9FpUwKYdJ0nbfIJbBrckJKRqp2afkDCE7Lxc3tjyXZUeV4iSJHY28dSDXdm6K5vJL6/nQn6pszCI6qnQYbOcBwahX99K5s/3Wf1BWQGllfK8pNeiAj74MpOYqGCGXtesyo/dzdSCIc2789oWOyfP14zreaJ6u/rSBoy5tyOr1h9nyuubKCiUj5kaJM9stQ8H/gEscr3n8Vt5WTfqJphSk0+ht5bCXc9xvZYrmRdpZ+ZpVqw7xv+NaE14WNW2nopMSBjKwn3reHWLnUndhxtSgxDuBvVtSvbpC7z50U5iooJ55K72MrhszfAP4ElgicNm+clstUcDz3nawGNAOZNSjPnUrCU+mJ9BZHggN93QvOyVfaRz3WYMa9GDlK3f84+O/TCFVnzuKSEqy+1D48k+fYGPF+0hNjqYu0dcYnRJ4iI5bJa1wM1ur0+brfZGnrbxdqgjUcn2Hsjhh1+PcPvQeKIigw2tZXy3ISzYm8arm7/lqcSbDK1FiCLJt7fh1JkLvP9FBrHRwdzcv4XRJYmLYLbaH0QfbLbojmwNaGq22h8DXnLYLC8X30YCyiBz52cSEhzAXwca/z9dp7pNuSm+F29udfCPjhbqhUUbXZIQKKUYc08HTudc4JX3txMdFcz1vauup6uodA8AA4GiS0UasBR9HL4SuxJLQBng0NGzfPvzIW66oRl1Yvxj2otx3QbzReZq/r35G57u8RejyxECgMDAAJ56oAsTpq9jWspmoiODuaJ7PaPLEhVzwGGzZLq/Ybbaj7nuiyqRzL1sgA8X7iFAwa2D/Gc+nPaxjRnR6lLe2rqUo2dLvGdOCEOEhATy/GMJtG4exdMvr2fDNm8Gzhb+xmGz9DNb7d3MVnu8+3uetpGAqmLHT5znq2UH6H9NE+rH+VdHyHHdhnCu8AL/2vyN0aUI8QeREUFMG59IfVMYT8xcx+69Z4wuSZST2Wp/F/gv8K3Zar/PbLXXMVvtr3naRgKqin361V4K8gu5fYj/tJ6KtIlpyC2tLuftbUs5lJttdDlC/EHdWH1w2bDQQMZPW8uBIzICSjVzGdAD6An83TWHVC9PG0hAVaHs03nMs2dhvrIRTRt5mv3YOGO7DuZCYQEvb/ra6FKE+JNG9cOZMTGRC/ka46am4TwpE3pXI7uBBq5x94LMVnsApUz1XkQCqgp9vmQf584XYB0Wb3QppWod04DbWl/Je9t/4GCunOsX/ie+WRRTxnbHefI846elcSZHBpetJs4D6a5TfY2A74DPPG0gAVVFcnLz+XzJPvr0qk+r5v59M+zYroMo0Ap5aeNXRpciRIk6tY3lmdEJ7Nmfw6RZ6ZzPk8Flq4H5wHj0eaQeAh5y2CzPeNpA5oOqIh8uyOTNj3byxnOX0b61/88cOnp5Kh/tXsHqG5+jWWSJ07gIYbjvlx/i+Vc3cmViPZ55pBtBQfKd2xMj54OqCLkPqgqczyvg08V76dXVVC3CCWBM10F8uHs5L238mpmX32F0OUKU6LorG3E6J5+X3t3K9P9sZmJyZ5/PSC0qxmy176aEWTAcNkur0raRgKoCixz7OXEqj5HDS/138DvNo+IY2aY3qTt/4pHO/WkeFWd0SUKUaLilGadO5/HOZ7uJjQrm/pHtZHBZ/+TeYy8UuAlo4GkDCSgfu5BfyEcL99C1fR0SOtY1upxyebTLQGw7f2HWhsW8dGWS0eUIUaqRN7bi5OkLfPb1PmJjQqrVl8HaooQRI14zW+1rgKdL20YCyse++ekgR53nGXtvR6NLKbdmkSbubNuHd7f/wKNdBtIyWoaYEf5JKcUDI9tx6swF3v5kF7HRxsyxVt2ZrfYBwMvo8zS95bBZppay3gjgU+BSh82y2st993R7GYh+P5THDJIrij5UUFDIh/Mzadcqmku7Vc9TZKO7DCRQBTBzw2KjSxHCo4AAxYRRnbiiez3mvLOVZSsPG11StWK22gOBV9EHdO0E3G622juVsF408DCwspyHmOH2eB5IBP7qaQNpQfnQ0pVH2H/4LM+M7lZtz4k3jqjDXe2u4a1tS3m0ywBax3g8ZSyEoYKCApj8cFfGTV3L869uJDIiiF5dq+eXQwNcBux02Cy7AcxW+0fAcGBzsfWeA6YDY8uzc4fNcl15C5IWlI8UFmrY5mXQskkkfXrWN7qcizK6ywBCAgKZJa0oUQ2EhQYyZWx3WjSJ5Kk569myS4bt8lJTYJ/b6yzXe78xW+2JQHOHzbKwvDs3W+2RZqt9itlqX2O22lebrfapZqvdY5d3CSgfWZ52jIysHO4YFl/tu702DI/l7nbX8nHGCnaektMmwv9FRQYzfUIidWNCmDh9HXv2V597Jn0sSCm12u0xym1ZSR9Uv90o6xqaaA4wpoLHfhUIQz+tdysQArziaQMJKB/QNI0P5mXQuH44/a5qaHQ5leLhzv0JCwhmxvpyf3ESwhBxdUOZ8XgigYGKcVPXcvjYOaNL8gf5mqb1cnu86bYsC2ju9roZcMDtdTTQBVhqttozgSuA+War3eOAr256OmyWRx02y26HzbLLYbM8hgwWW/XWbnKyddcpbhvaksDAmvErbhAew70d+vK/zFVsyz5odDlCeKVpwwimT0gk91wB46au5eSpPKNL8mergLZmq72V2WoPAW5DH54IAIfNku2wWeo5bJZ4h80SD6wAhnnbiw/IL+G9Qk8b+MWnpyk1+R1TavIRU2ryRqNrqQwffJlJvbqhDLimidGlVKqHOvUnPDCEGesXGV2KEF5r0zKaKWO7c/jYOSZMTyP3bEmfk8Jhs+QDDwJLgC3AJw6bZZPZan/WbLUPq4RDPGO22usUvTBb7bGA/4/FZ0pNvgY4A7zvTErp4s02/joW38btJ3nomdXcP7Itfx3of3M+Xazn0r7gpY1L+GnoP+lYp2YFsKjZlqcd48nZ6SR0rMPUsd0JCQk0uqQqZ/RYfGarvR76qUENWOGwWY57Wt8vWlDOpJQfgFLnpa9OPpiXQWx0MEPMNfMmwQc6XU9kcCjT5VqUqGauTKzHhFGdSNt0gudf20hBofFfzmsTs9XeF1gN3A68B3xuttqv97SNXwRUTbEj8xQr1x1nxIAWhIfVzG9nptAokjtcx7w9a9h0IsvocoQolxuubswDSe34cdVR5ry9BX84g1SLTAf6OWwWK3p39gGUcYqvWgWUUmpUUffI/Hz/O49sm5dJZHggN15fM1tPRR7oaCE6OIxp6QuMLkWIchsxoAVJN7Zi0dIDvPXxLqPLqU1CHTZL0S9cOWyWs+hdzUtVrQJK07Q3i7pHBgX51yAYe/bn8MOqI9x4Q3OiIoONLsen6oRGcn9HCwv3rWO9c6/R5QhRbnePaM3Q65oyd0EmHy/aY3Q5tYVmttqLpngPNlvt4wGP3xCqVUD5s7kLMgkNCWDEgBZGl1Il/tHRQmxIBNPS5VqUqH6UUjxydweuvbwBb8zdwdfLDpS9kbhYT6LfWwWwHH3Kjbs9beAvvfg+BPoC9YDDwGRnUsrbnrbxp158B4+cZeSYX7j5hmY8kNTe6HKqzMz1i3gxfT7fDXqcxLh4o8sRotzyLhQyaeY61m4+wTOPdKVPr5o91qTRvfjKyy8CqiL8KaDmvLuVxY79zJ3Tm/pxYUaXU2VO5Z0l8YtJXFq/NR9d96DR5QhRIWfP5fPYi2vZtfcM0yck0r2azdtWHtUtoOQU30U6fuI8Xy07QP9rmtSqcAKICQnnwU7X883+Daw+mmF0OUJUSHhYEFPHdadJ/XAmzVzHjsxTRpckXCSgLtIni/dQkF/I7UNr3k253ri3g5m40Cimrp9f9spC+KnY6BBmTEwkOjKY8dPS2HfQP87O1HYSUBch+3Qe87/bz3VXNaJpwwijyzFEdHAYD3W+ge8PbGblEemyK6qv+nFhTJ+YiKbBuKlpHHXK4LJGk4C6CP9bso9z5wu4Y1i80aUY6p52fakfFs1UuS9KVHMtmkQybUIip85cYPy0NE6duWB0SbWaBFQF5eTm88WSfVzdqz6tmkUZXY6hIoNDebhzf5Yd2sLywzuMLkeIi9K+VQwvPJbA/kO5PD5jHWfPFRhdUq0lAVVB8+xZnMnNx3pjK6NL8Qt3t7uWhuEx0ooSNUJiZxNPPdiVrbuymfzyei7ke5wVQviIBFQFnDtfwKdf7eHSbnG0bxVjdDl+ISIohEc6D+DHw9v48dA2o8sR4qJdfWkDHru3I6vWH2fK65tkcFkDSEBVwCLHfk6eusDI4fFGl+JX7mx7NY3D6zA1fb4MwilqhMF9mzLqtjY4Vhzmlf9uk7/rKiYBVU4X8gv5eNEeurWvQ7cONfeGvooIDwrh0a4DWH5kJz8c2mp0OUJUituHxnPr4JbMs2fx3893G11OrSIBVU7f/HiQo87zcu2pFElt+tAkoi5T0hfIt01RYyTf3oaB1zbhv59n8PkSGSC5qkhAlUNBQSFzF2TSrlU0l3Y1GV2OXwoNDGZM14H8enQX3x/cbHQ5QlQKpRRj7ulAn171eeX97Xz780GjS6oVJKDKYenKIxw4fBbr8FYopYwux29ZL+lNs0gTU6UVJWqQwMAAnnqgCwkd6zAtZTMr1x0zuqQaTwaL9VJhoca9j6+gUIN3pl5BQIAElCf/3fEjj674gI+ve5Drm3Y1uhwhKk1Obj6PvrCGvQdymPl4D7q0q2N0SV4ra7BYs9U+AHgZCATectgsU4stfwy4F8gHjgL/57BZfDahlrSgvPTL2qNkZOVgHRYv4eSFOy65ipZR9aQVJWqcyIggpo1PpL4pjMdnrGP33jNGl1QpzFZ7IPAqMBDoBNxutto7FVstDY7aa04AABodSURBVOjlsFm6AZ+hT+PuMxJQXtA0Ddu8TJo0COe6KxsaXU61EBwQyNiug0g7voclWeuNLkeISlU3Vh9cNiw0kPHT1nLgSK7RJVWGy4CdDptlt8NmyQM+Aoa7r+CwWRwOm6Xoh13B7xMQ+oQElBfWbHSydfcpbhvaksBA+ZV569bWV9Aqur706BM1UqP64UyfmEjehULGTU3DefK80SVdrKbAPrfXWa73SnMP8JUvC5JPWy98MC+DenVD6X91E6NLqVaCAgIZ13UwG07sY9G+dUaXI0Sla9UsiqnjEnGePM/4aWmcyfH7wWWDlFKr3R6j3JaVdO2ixG+WZqt9JNALmOGLIotIQJVhw7aTpG85ya2DWxISLL+u8hrR6jLaxDRk2vqFFGoynpmoeTq1jeWZ0Qns2Z/DpFnpnM/z68Fl8zVN6+X2eNNtWRbQ3O11M+BA8R2YrXYLMAkY5rBZfNpslE/cMtjmZRAbHcxgs6eWrihNUStq04ksFuxNM7ocIXzism5xPP6PzmzYfpJnX9lAfvUcXHYV0NZstbcyW+0hwG3AH2YiNVvtiUAKejgd8XVBElAebM84xcr044wY2ILwsECjy6m2bo6/lLYxjZiWLq0oUXNdd2UjHrmrA7+sPcb0/2ymsJoNLuuwWfKBB4ElwBbgE4fNsslstT9rttqHuVabAUQBn5qt9nVmq92nU2nLfVAeTH5pPWs2Huejf11NVESQT49V032euYp7f3yLt66+l5vjLzW6HCF85v0vdvPuZ7sZMaA5949s51c39Zd1H5S/kRZUKfbsz+HH1Ue48YbmEk6V4MaWPekQ24Rp6QspKJRWlKi5km5sxc39m/PZ1/uwzc80upxqTQKqFHMXZBIaEsCIAS2MLqVGCFABTEgYwo5Th/hiz2qjyxHCZ5RSPDCyHZbejXj7k10s+D7L6JKqLQmoEhw8chb7z4cYYm5KnZgQo8upMYa2SKRz3WZMX7+Q/EK/7ukkxEUJCFBMGNWJK7rXY847W1m28rDRJVVLElAl+GhhJoEBcMvglkaXUqMEqAAmdBvCzlOH+SzjV6PLEcKngoICmPxwVzq3jeX5VzeyesNxo0uqdiSgijl24jxfLTtA/2uaUN8UZnQ5Nc7g5t3pZmrOjA2LpBUlaryw0EBeHNudFk0ieWrOerbsyja6pGpFAqqYTxbtoaAQbh8qrSdfUEoxodtQMk4f5ePdK4wuRwifi44MZvqEROrGhDBx+jr27K+6WRiqOwkoN9mn81jwfRb9rmxIkwYRRpdTYw1o1o3uphbM3LCYC9KKErVAXN1QZjyeSGCgYtzUtRw+ds7okqoFCSg3//t6H+fOF/L/7d17dFTVvcDx728mk0lCQiAJ4ZHwUnkrBkVEg4+wEFyl8vCiRiJiWxUUhIC1onb1xmJbtOlVKV7RilehQXFhVapei6EBepUq77eiQgyRR54kM0lmMpPZ9485CZMX5slMkv1ZK2ty5uxz5jfnzJnf2Wef2Xv2tEH+DqVTExGeSJjG9/YC3vpuh7/D0bSLIq53GM89PoZyRxWPrdjDudJKf4cU8HSCMtjL3fxt80luuKYXg+LD/R1Opzep3+VcHTOY9IMfUVnl9nc4mnZRXDYwgt8/eiVnCxw8/txeyiv0Z/9CdIIyfPDpScrK3aRMH+zvULoEEWHZlbeRW1ZExnef+TscTbtoRg/vSdqiK/j2ezu/fn4/lYHduaxf6QQFOJxVbPwkh3Gjoxk2uLu/w+kyJvYdybhel/Kng/+LsyrghynQtDZz3VW9ePzBkew9XMwz/32Iqg7Wb9/F0qn64nO5XOTm5uJwNK8B8p9f2tm4uYSl98Zw2QBrW4bZZYSEhBAfH4/FYmnWcltPH+X2zBe4b8iNJPYeSkxIODEhEURZw4m2hhNs1t1MdVUtPZ47kvIKN7YyN6FWM90jmnfsXEhjx2NH64uvUyWoEydOEBERQXR0dJM7aHS5PcxO/Yx+saG8+Jux7RFqp6eUorCwEJvNxuDBzbtEqpRiztaX+Th3f4Pzu1tCiQmJIDoknBir9zHaGk50SAQx1Y818yIIC9I9f3QWLTmeO6L8IgdF5yqJ6hHcJr+9vNDx2NESVKc6PXU4HAwaNKhZH+bN/zpNQbGTxx4c2Y6RdW4iQnR0NPn5+S1adu3N8ylyllHgsFHotFPgsFHktFPgsFPosFHg9D7mlBWytzCbQqe90dvTQ82WOskrgmjr+VpZTEg40daImppad0top/7y68hacjx3RDE9rVRVKYrOVWI2CVE9WncVpzXHY6DpVAkKaNaHuarKw/q/ZzN0cATXXBHVjlF1fq35EjGJiRgjmTSFUgqby0GBw0aB00ahw16T3AqNZFbgsFPotHGs5AyFDhvlVQ3f0msxmYm2htdLXt4kV7/GFmUNx2zSTbcXS2dPTuB9j71jQvB4FPlFTsxmITKidVcCOst263QJqjmy/p3HqbMV/DZ1dJvtULPZzBVXXFEznZyczLJly7j//vtZunQpI0f6t6YWHh6O3W73awytJSJ0Dw6le3AolxDbpGXK3ZVGrcxWq6ZW6PBJak47B4pyKHDaKaksb/i1EXpaw84nsnqXG+vX2Kzmtmtb0C6u6uPZ5XIRFBTE3LlzSU1NxWQysWvXLtauXcvKlSsbXDY7O5vPP/+c2bNn/+jriAh9eoVS5SnnTL4Dk0mI6KY/N102QXk8ioxNJxgU343Eq3u12XpDQ0PZt29fvedfe+21NnsNrfnCgoIJC4oivlvTasouTxWFDvv5pFZdU/OpsRU57XxbmscX+d9R6LTjaaQ9N9wScj6BGY/exOZNZNU1terHbkHWTnMG3NH5Hs95eXnMnj2bkpISnn76acaOHcvYsY23W2dnZ7N+/fomJSjw9oAe1zuMk6fLOZ1XgbmPEBbaZb+igS6coD7fk092bhlPPjwKk6n9vwxuvvlm0tPTGTt2LGvWrOHZZ5+lX79+DBkyBKvVyqpVq8jPz2f+/Pnk5OQA8MILL5CYmEhaWho5OTkcP36cnJwcUlNTWbRoEY8//jgDBw7k4YcfBiAtLY2IiAjmzZvH9OnTKS4uxuVy8cwzzzB9+vRa8WzdupX09HQ+/PBDABYuXMjYsWO577772L17N0uXLsVutxMTE8Mbb7xB3759WblyJatXryYoKIiRI0fy9ttvt/t28xeLyUyfsEj6hEU2qbxHeThXWV7TblZTQ/OpqRU4bZyuOMfB4lwKHTacnoZ/pBlittS/5Gg8Rhm1M9+aWmRwKCbRlx3bW2xsLK+++irXXHMNaWlpbNu2reYY2rZtG4sXLwa8taHt27ezbNkyjh49SkJCAnPnzmXmzJnMmTOH6pu7Vq1axfXXX8/WrVtJS0sjJiaGQ4cOMeaqq/jdilf44Ww5Z384wq8eW0pZWRlWq5UtW7YQFhbGsmXL2Lp1K06nkwULFjBv3jx/bpp202kT1BM7N3CouOGBwpSCio+GIOFmXrStZ+Xmpq3z8p7x/OGauy5YpqKigoSEhPNxPPEEd911fplTp06xfPly9uzZQ0REBBMnTuTKK68EYPHixSxZsoQJEyaQk5PDlClTOHr0KABfffUVWVlZ2Gw2hg0bxkMPPURycjKpqak1Ceqdd97hk08+ISQkhPfee4/u3btTUFDA+PHjmTZtWpPOyl0uF4888ggffPABvXr1YsOGDTz11FO8/vrrrFixghMnTmC1Wjl37lzTNloXYRITUUYbFZF9frS8Ugq721mrvaxuMquuqR235VHotGN3NXy7tVlMtdvR6tzhGG09n8yq29SCTOa23gTtatW6r/n2e1ubrvOygREsnDOsWctccskleDwe8vLyaj2fnp7OSy+9RGJiIna7nZCQEFasWFHrJLC8vJxPP/2UkJAQvvnmG+6++2527fIO3rl3714OHz5Mv379SExMJOf4XnrHj+Luu+9mw4a3uP668ZSWlhIaGsqaNWuIjIxk586dOJ1OEhMTmTx5crPvoO0IOm2CupCqU+F4CsOwXneStj7xbOwSX7Uvv/ySm266iago76WmO+64g2PHjgGQmZnJkSNHasqWlpZis3kPyqlTp2K1WrFarcTGxnL27FnGjBlDXl4ep06dIj8/n549ezJgwABcLhdPPvkk27dvx2Qy8cMPP3D27Fn69PnxL86vv/6aQ4cOccsttwBQVVVF3759ARg9ejQpKSnMmDGDGTNmtGwDaYD3LDvCEkKEJYRBEU27xOyochltZjafOxxtFDnsNXc6FjjtHCnOpcBhp7iy8V6zewSH1dTCousks+ok1yM4DKspiGCzBaspCIvZjNVkIdgc5J02mbvkpciGfpqTmJjI0qVLSUlJ4fbbbyc+Pr5eGZfLxcKFC9m3bx9ms7nmuAcYN25czTIJCQnk5uYQ0b07sbG96R03CpfbQ/fu3k4ENm/ezIEDB9i4cSMAJSUlfPPNNzpBdSQXquksXr6LUz0ryJg3l2DLxb00cqHfnXk8Hnbs2EFoaGi9eVbr+VtPzWYzbrf38tCsWbPYuHEjZ86cITk5GYCMjAzy8/PZvXs3FouFQYMG1fuxY1BQEB6Pp2a6er5SilGjRrFjR/1OXD/66CO2b9/Opk2bWL58OYcPHyYoqNN+hAJOiNlCXLeexHXr2aTybk8Vxc6yWsmr0LhJpOY2fqeNbHsBuwtPUOiw41aeH1+xD28CCzqfuMxBWExBPs/7PlpqpoMbnO9dj9Xs/X+YO4xzzjJEhHvvGoBJBEEQEUx4k/z5aUGEmun2cvz4ccxmM7GxsTVXNwCWLVvG1KlT+fjjjxk/fjyZmZn1ln3++efp3bs3+/fvx+PxEBJy/jdPDR3fliAT1mAzHo8i90w5A/qGYTabUErx5z//mSlTprT5+0tKybwVeBEwA69lZUxaUWe+FVgLXA0UAndlZUzKbvNADAHz7RK1bl6tDVM055UVP7JIixz8+hwHvjrHgjlDL3pyAu+Z0pIlSyguLiYiIoJ333235q6/yZMns2rVKh577DEA9u3bV+tyYUOSk5N54IEHKCgoYNu2bYD3jCo2NhaLxUJWVhbff/99veUGDhzIkSNHcDqdOBwOtmzZwoQJExg2bBj5+fns2LGD6667DpfLxbFjxxgxYgQnT54kKSmJCRMmsH79eux2Oz169GjjLaS1lSCTmV6h3ekV2rTuu5RSlFSWG5ca7ZS6ynFWuXFWuan0uHFWuYxH73Slx01lrfluKj2uOtNuytxOiqvKGl2Ps8qNov6J29vD52CyFzT7fTecvAQxnj8/LZh8kppvuepkqIC8ilJMIhTmFzD/wQe4f/6DlLoqKHM5cSsPZS4nJ45/x2UjhpE6cjifff45h44cYUD//thsNpRSiAglJSXEx8djMpl48803qaq6cB98w4cP58yZ05w5eZje8aP46ts8hgyOZsqUKbz88stMnDgRi8XCsWPHiIuLo1u31v3+Nikl0wy8BNwC5AI7k1IyN2VlTDriU+wXQHFWxqTLklIyk4FngQu3e7RCQCSoqHXz6m2YqHXzNhXNeeXIhZdsvr++f4LICAs/TYpr61UD9dugbr31VlasOJ9r4+LiePLJJ7n22mvp168fI0eOJDLS2xC/cuVKFixYwOjRo3G73dx4442sXr36gq83atQobDYbcXFxNZfiUlJSuO222xg7diwJCQkMHz683nL9+/fnzjvvZPTo0QwZMoQxY8YAEBwczMaNG1m0aBElJSW43W5SU1MZOnQo99xzDyUlJSilWLJkiU5OnYyI0MPajR7WblzavfdFe12lFFXKUy/Rnfv+DEMi++JRCoVCKW8a8532NPK8QvlM41PO+7xHKZTy1Cmn8BjxKCNlOioquGHceNwuN+YgM1PvuJ07H76f47Z8TpUXU+Zycqz0DM+mP8eu/9uByWzikqFD6D/+cjCZqFBuhl0xkmnJs5h49wwe/dl81r6Vwbgbrie0WxhfnTtFjr2QMreTb0vPIkBpZQUFDhunK0tZ+T+vsPjRR6god2C1hLJ2w/vc+7P7yM7O5qqrrkIpRa9evXj//ffbYleMA77Nyph0HCApJfNtYDrg+z08HUgz/t8IrEpKyZSsjEnt0iVRQHR1FLVu3nVAWtGcV6YY008AFM155Q+NLdNQV0dHjx5lxIgRjb5Oia2Se3+5gzt+MoB7/Nhrud1uJzw8HLfbzcyZM/n5z3/OzJkz/RZPW/mx7a9pzeHvz1PtBNhwoqtOjr6Jzjcx1qyj1nKcT4h1k2Mj6/cohanSgqUijL59QogM+/HeJhrafiJSCRz0eepVpdSrAEkpmbOAW7MyJt1vTM8Brs3KmLSwunBSSuYho0yuMf2dUab5Vd0mCIgaFBAHnPSZzgWurVtIRB4EHgTvmX5zRUYE8/aLiYB/G3bT0tLIzMzE4XAwefJkfcOBpgUgEcGM+PvrohaXy0NQUKsCciulGvvxVkMrrluDaUqZNhMoCapJb9rI9K+CtwbVkhcKDfH/W05PT/d3CJqmdUCW9m03zwX6+0zHA6caKZOblJIZBEQCRe0VkP+/rb2asmE0TdO09rMTGJKUkjkY+AFIBup2g7EJmAvsAGYB/2yv9icInAELdwJDotbNGxy1bl4w3g2zqSUrCoQ2ta5Ib3etPejPVcu0ZLtlZUxyAwuBfwBHgXeyMiYdTkrJ/G1SSuY0o9gaIDopJfNbYCmwrI1CblBA3CQBELVu3k+AF/DeZv560ZxXfneh8m01HpTWeq0ZD0rTGqOP55bpTONBBUyCaq62HFFXa72WjqiraY3Rx3PL6RF1/ayhBKVpmqY1rqMlqEBpg9I0TdO0WnSC0jRN0wKSTlCapmlaQOqwbVAi4gEqWrh4ENDwaHGav+h9Enj0Pgk8rd0noUqpDlMx6bAJqjVEZNcFuvvQ/EDvk8Cj90ng6Wr7pMNkUk3TNK1r0QlK0zRNC0hdNUG96u8AtHr0Pgk8ep8Eni61T7pkG5SmaZoW+LpqDUrTNE0LcF0qQYnI6yKSJyKH/B2L5iUi/UUkS0SOishhEVns75i6OhEJEZEvRWS/sU+e9ndMmpeImEVkr4h86O9YLoYulaCAN4Bb/R2EVosbeFQpNQIYDywQkZF+jqmrcwITlVJXAgnArSIy3s8xaV6L8Q6F0SV0qQSllNpOO47+qDWfUuq0UmqP8b8N78EX59+oujblZTcmLcafbqz2MxGJB6YCr/k7loulSyUoLbCJyCBgDPCFfyPRjEtJ+4A84FOllN4n/vcC8CvA4+9ALhadoLSAICLhwLtAqlKq1N/xdHVKqSqlVAIQD4wTkcv9HVNXJiI/BfKUUrv9HcvFpBOU5nciYsGbnDKUUn/zdzzaeUqpc8BWdNutvyUC00QkG3gbmCgif/VvSO1PJyjNr8Q7lvca4KhS6r/8HY8GItJLRHoY/4cCk4Cv/BtV16aUekIpFa+UGgQkA/9USt3j57DaXZdKUCLyFrADGCYiuSLyC3/HpJEIzMF7RrjP+PuJv4Pq4voCWSJyANiJtw2qS9zWrAUW3ZOEpmmaFpC6VA1K0zRN6zh0gtI0TdMCkk5QmqZpWkDSCUrTNE0LSDpBaZqmaQFJJygt4IiIEpE/+Uz/UkTS2mjdb4jIrLZYVwtee7hxG/1eEbm0zrxsETnoc6v99S1Y/wzd0a7WmegEpQUiJ3C7iMT4OxBfImJu5SpmAB8opcYopb5rYH6SUirB+Pu8hetvVoISkaAWvI6mXRQ6QWmByI13aOsldWfUrQGJiN14vFlEtonIOyJyTERWiEiKMa7RwTo1lkki8i+j3E+N5c0i8kcR2SkiB0Rkns96s0RkPXBQRLqJyEfGWEmHROSuBmJMEJF/G+t5T0R6Gj8+TgXuF5GspmwEEQkXkS0issd4D9N95t1rrH+/iKwzalzTgD8aNbBLG4rDWHariPxeRLYBi0XkDuO97BeR7U2JTdMuBn32pAWql4ADIvJcM5a5EhiBd0iV48BrSqlxxiCIj+BNEACDgJuAS/H2mHAZcC9QopS6RkSswGcistkoPw64XCl1QkT+AzillJoKICKRDcSxFnhEKbVNRH4L/KdSKlVEVgN2pVR6I/FniUgV4FRKXQs4gJlKqVKjNvlvEdmEt5b0FJColCoQkSilVJEx70Ol1EYjtgN14/DZBj2UUjcZ5Q4CU5RSP1R3caRpgUDXoLSAZPRovhZY1IzFdhrjSzmB74DqBHMQb1Kq9o5SyqOU+gZvIhsOTAbuNYaY+AKIBoYY5b9USp3wWdckEXlWRG5QSpX4BmAkrB5KqW3GU28CNzYx/upLfNdWrw74vZFoMvGOk9UbmAhsVEoVACil6o1x1oQ4Nvj8/xnwhog8ALT2MqamtRmdoLRA9gLwC6Cbz3NujM+t0dFssM88p8//Hp9pD7WvFtTt30vhTQaP+LQBDVZKVSe4spqCSh0DrsabqP4gIr9pyRtrohSgF3C1MfTFWSDEiLW1fZT5vqf5wK+B/sA+EYlu5bo1rU3oBKUFLKNm8A7eJFUtG2+CAJiOd7TX5rpDRExGu9QlwNfAP4CHjKE/EJGhItKt7oIi0g8oV0r9FUgHrqoTcwlQLCI3GE/NAbbRMpF4xwByiUgSMNB4fgtwZ3UiEZEo43kbENHcOETkUqXUF0qp3wAFeBOVpvmdboPSAt2fgIU+038BPhCRL/F+UZc1uNSFfY33y7o3MF8p5RCR1/BeBtxj1Mzy8d4VV9cVeG9E8AAu4KEGyswFVotIGN5LiD9rQYwAGcDfRWQXsA9jyAul1GER+R2wzWiz2gvch3ecoL+IyCJgVjPi+KOIDMFbM9sC7G9hvJrWpnRv5pqmaVpA0pf4NE3TtICkE5SmaZoWkHSC0jRN0wKSTlCapmlaQNIJStM0TQtIOkFpmqZpAUknKE3TNC0g6QSlaZqmBaT/B5PLhv+UDGgdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#elbow method\n",
    "a,b=get_line_params(0,eigvals[0],len(eigvals)-1,eigvals[-1])\n",
    "\n",
    "distance=[]\n",
    "for i in range(len(eigvals)):    \n",
    "    distance.append(get_distance(i,eigvals[i],a,b))\n",
    "\n",
    "dual_axis_plot(np.arange(1,len(distance)+1),eigvals,distance,\n",
    "               x_label='Numbers of Factors',\n",
    "               y_label1='Eigenvalues',\n",
    "               y_label2='Perpendicular Distance',\n",
    "               legend1='Eigenvalues',legend2='Distance',\n",
    "               title='Elbow Method for Scree Plot',\n",
    "               fst_color='#11a266',sec_color='#405abb')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The optimal number of factors suggested by Scree Plot is 2\n"
     ]
    }
   ],
   "source": [
    "#apparently different approaches end up with different results\n",
    "print(f'The optimal number of factors suggested by Scree Plot is {distance.index(max(distance))+1}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Variance Explained Criteria\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "This is the same as PCA. The R squared (usually 95% in statistics) exceeding the preset threshold is the optimal number of factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute variance explained\n",
    "def compute_r_squared(data,z,miu,lambda_):\n",
    "\n",
    "    #compute fitted and error\n",
    "    fitted=lambda_@z.T+miu.reshape(-1,1)\n",
    "    error=(data-fitted.T)\n",
    "\n",
    "    #sum of squared errors\n",
    "    #equivalent to np.linalg.norm(error)\n",
    "    numerator=np.diag(error.T@error).sum()**0.5\n",
    "\n",
    "    #sum of squared total\n",
    "    denominator=np.diag((data).T@(data)).sum()**0.5\n",
    "\n",
    "    return 1-numerator/denominator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "#iterate through all possible numbers of factors\n",
    "rsquared=[]\n",
    "threshold=0.95\n",
    "for i in range(1,iris.data.shape[1]+1):\n",
    "    _,_,lambda_,psi=training(iris.data,i)\n",
    "    z=get_latent(iris.data, miu, lambda_, psi)\n",
    "    rsquared.append(compute_r_squared(iris.data,z,miu,lambda_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "#find the optima\n",
    "try:\n",
    "    optima=[i for i in rsquared if i>threshold][0]\n",
    "    ind=rsquared.index(optima)+1\n",
    "except IndexError:\n",
    "    ind=len(rsquared)-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcFNW9///XG4ZFWRVQ2RRUXDAi6khEo06UJJjkuoDrVRNJ/BFcEveF3MRvojfBBeOSaLwmATSaKLhFY9QIwSVqIoPsKItoZEAUVDaRZZjP74+qgWYchhl6mp7l/Xw8+jG1nDr1qW7oT1edU3UUEZiZmW2vJvkOwMzM6jcnEjMzy4oTiZmZZcWJxMzMsuJEYmZmWXEiMTOzrDiRWJ0jaU9JqyU1zUHdP5P0YG3Xu70khaR987Tv/SVNkbRK0o/yEYM1DE4kljVJ50uaIWmNpCWSfiupfQ22f0/SgPL5iHg/IlpHxMbcRLzVOIrSL/a7Kyz/p6Tzd2QsO8g1wIsR0SYi7qq4UtKLktamSb381X97dyapIH1/e2QRs9VBTiSWFUlXAjcDVwPtgCOBvYAXJDXPZ2zb6TPgO/Xty05SwXZsthcwaxtlLkmTevnr9e3YT63YzmO0HcCJxLabpLbAz4EfRsRzEbEhIt4DziD5kjo3LfczSY9KeiS9jPKmpEPSdX8E9gSeTn/xXiOpR/rLtSAt86Kk/5X0WlrmaUkdJD0kaaWkSZlf/JLulLQwXTdZ0jE1OKzlwBjg/23lmLe4NJZtrKlvSlogaZmkWyU1yaj/e5LekvSppOcl7ZWxLiRdLGkeMG8r8Z4kaZak5WlsB6bL/wF8FfhNGud+NXiPkPQbSSUZx3RUxroCST+V9E66vlhSF+DltMisdJ+D0/LDJM2X9LGkJyV1zqgnJF0kaT7wtqQmku6S9JGkFZKmS+pdk9gtByLCL7+26wUMBEqBgkrW3Q/8OZ3+GbABOA1oBlwFvAs0S9e/BwzI2LYHEOX1Ai8C84F9SM56ZgNzgQFAAfAAMDpj+3OBDum6K4ElQMuMWB7cyvEUASXAHsBKYP90+T+B8yvbvhZiDWAisCtJQp0LXJCuOyWt68B0258Ar1XY9oV0250qOZ79SM6wvpa+79ek9TXPiPWCKj7fra4Hzkv3WwBcCywCWqTrhgPTgF4kP1b7ZpQNoEdGPV8HPkrLtATuAf6Rrisv/xywC7AT8C3gjfS9bQL0BvbI9/+Fxv7yGYlloyOwLCJKK1n3Qbq+3OSIeDQiNgC/IvnSOLIG+xodEe9ExArgWeCdiBif7nsccGh5wYh4MCI+jojSiLgNaAHsX90dRcQS4F7ghhrEt12xpm6OiE8i4n3gDuDsdPkPgBER8Va67S+BvplnJen6TyLi80riOBN4JiJeSN/3kSRfxkdVUnZr7krPZpZLerN8YUT8Md1vKXAL0BYo7zRwAfDjiJgXEWURMTUiPtlK/ecAv0/LrAWuA46T1C2jzC8j4tP0GDek+zogjWN2+nlZHjmRWDaWAR23cu26c7q+3MLyiYgoI/nl36UG+/owY/rzSuZbl89IujK9HLRC0nKSX6+ZSa06bga+UX4JroaqHWtqYcb0f9j8vuwF3Fn+RQ58AgjoupVtK+qS1gdset8XVth+W34UEe3T12HlC9NLkG9LWgF8CrRi83vcHXinmvVXjHFlWl+lxxgRfydJ8r8FPpR0r6Q2NTgeywEnEsvG68A6YFDmQkmtgBOBCRmLu2esbwJ0Axani2rtEdRpe8i1JO00u0REe2AFyRdwtUXExyRnBzdWWPUZsHPG/B7bH+0m3TOm92Tz+7IQ+EHGF3n7iNgpIl7LDLWKeheTJCMAJCnd16JsgpX0VeAKYDDQnuSy02o2v8cLSS7tVVRZrBVjbJPWlxnjFttFxB1pUvsSyaWtK7brQKzWOJHYdksv3fwc+LWkgZKapQ3J40jOOP6YUfxwSYPSs5fLSBLQv9J1HwJ711JYbUjabZYCBZKuJ7kUsj1+RXIZ6MCMZVOBY5Xc69KOpD0gW1dL2kVSd+BS4JF0+b3AcEkHAUhqJ+n0GtQ7FviWpBMkNSNpL1oHvFb1ZttU/h4vI2l7+RnJGUm53wP/K2kfJfpK2jWS7twfs+Vn/Wfg+5L6SGoBjABeiYiSynYsqV/6KiBJ6uuBHdpN3L7IicSyEhG3AD8muf6+Evg3yS/SEyJiXUbRv5Bcs/+UpKF2UHrdHpIvj5+kl3CuyjKk50naJeaSXDJZS9WXf7YqvcxyC0lDcfmyF0i+6KcDk4G/ZhkvJO/NZJIk9Qzwh3RfT5BcYntY0kpgJsmZXnXjn0PS8eDXJF/6/wX8V0SszzLevwHjSXqKvUfyuX+Qsf5W4EmSM9KVwH0kbWKQ9Ib7U/pZD4qI50jaop5I69iTpN1ka9qTvD/L031/ANye5fFYlhThga0styT9DNg3Is7NdyxmVvt8RmJmZllxIjEzs6z40paZmWXFZyRmZpaVRvEQtI4dO0aPHj3yHYaZWb0yefLkZRHRaVvlGkUi6dGjB8XFxfkOw8ysXpH0n22X8qUtMzPLkhOJmZllxYnEzMyy4kRiZmZZcSIxM7Os5DSRpE+EnZMOo3ldJev3kjQhHS7zxczBbCRtlDQ1fT2VsbynpH9Lmqdk6Nb6OC64mVmDkbNEIqkpcDfJ00p7A2dXMrbySOCBiOhD8gTQERnrPo+IvunrpIzlNwO3R0QvkifJfj9Xx2BmZtuWyzOSfsD8iFiQPrb6YeDkCmV6s3nwo4mVrN9COjDP8cCj6aL7Sca1zpmioiLGjBkDwIYNGygqKuLBBx8EYM2aNRQVFfHII8nwEStWrKCoqIjHH38cgGXLllFUVMTTTz8NwJIlSygqKuK5554DYOHChRQVFTF+/HgAFixYQFFRES+99BIAc+bMoaioiNdeS4aPmDlzJkVFRUyaNAmAqVOnUlRUxNSpUwGYNGkSRUVFzJw5E4DXXnuNoqIi5syZA8BLL71EUVERCxYsAGD8+PEUFRWxcGHylPXnnnuOoqIilixJRi59+umnKSoqYtmyZKDDxx9/nKKiIlasWAHAI488QlFREWvWrAHgwQcfpKioiA0bkqfDjxkzhqKiok3v5e9+9zsGDBiwaf6ee+7hxBM3PxX9zjvv5KSTNv9mGDlyJIMHD940f9NNN3HWWWdtmr/xxhs599zNDxS+/vrrGTJkyKb54cOHM3To0E3zV111FRdffPGm+csuu4zLLrts0/zFF1/MVVdtfor90KFDGT5883AjQ4YM4frrr980f+6553LjjZvHvTrrrLO46aabNs0PHjyYkSNHbpo/6aSTuPPOOzfNn3jiidxzzz2b5gcMGMDvfve7TfP+t+d/e+Wy/beXa7lMJF3ZchyIEr44xOc0klHWAE4F2kjqkM63lFQs6V+SypNFB2B5xhjhldUJgKSh6fbFS5cuzfZYzMxsK3L20MZ0JLdvRMQF6fx5QL+I+GFGmS7Ab4CewMskSeWgiFghqUtELJa0N/AP4ASSQXJej4h90+27A3+LiIOriqWwsDB8Z7uZWc1ImhwRhdsql8tHpJSw5VjUmWN0AxARi0nH+5bUGhicDt9avo6IWCDpReBQ4DGgvaSC9KzkC3WamdmOlctLW5OAXmkvq+bAWcBTmQUkdZRUHsNwYFS6fJd0/GYkdQSOBmZHcvo0ETgt3ea7JMOUmplZnuQskaRnDJeQjKH9FjA2ImZJukFSeatWETBH0lxgd+AX6fIDgWJJ00gSx00RMTtddy1whaT5JG0mf8jVMZiZ2bY1ioGt3EZiZlZz1W0j8Z3tZmaWFScSMzPLihOJmZllxYnEzMyy4kRiZmZZcSIxM7OsOJGYmVlWnEjMzCwrTiRmZpYVJxIzM8uKE4mZmWXFicTMzLLiRGJmZllxIjEzs6w4kZiZWVacSMzMLCtOJGZmlpWcJhJJAyXNkTRf0nWVrN9L0gRJ0yW9KKlbhfVtJS2S9JuMZS+mdU5NX7vl8hjMzKxqOUskkpoCdwMnAr2BsyX1rlBsJPBARPQBbgBGVFh/I/BSJdWfExF909dHtRy6mZnVQC7PSPoB8yNiQUSsBx4GTq5QpjcwIZ2emLle0uHA7sDfcxijmZllKZeJpCuwMGO+JF2WaRowOJ0+FWgjqYOkJsBtwNVbqXt0elnrp5JUWQFJQyUVSypeunTp9h+FmZlVKZeJpLIv+KgwfxVwnKQpwHHAIqAUuAj4W0Qs5IvOiYiDgWPS13mV7Twi7ouIwogo7NSp0/Yeg5lZ/bRxBZScmvzNsYIc1l0CdM+Y7wYsziwQEYuBQQCSWgODI2KFpP7AMZIuAloDzSWtjojrImJRuu0qSX8iuYT2QA6Pw8ys/ln9FKx+ElY/De3OzemucnlGMgnoJamnpObAWcBTmQUkdUwvYwEMB0YBRMQ5EbFnRPQgOWt5ICKuk1QgqWO6bTPg28DMHB6DmVn9tGLUln9zKGdnJBFRKukS4HmgKTAqImZJugEojoingCJghKQAXgYu3ka1LYDn0yTSFBgP/C5Xx2BmVm+8PwDWTMhY0Dz5s+ZVeDujpWHnE2DP8bW6a0VUbLZoeAoLC6O4uDjfYZiZbbeysmDNho18tq40fW3ks/Xp9PqNtCp9mWN3HkIzrd16JdoZuj0DrYqqtU9JkyOicFvlctlGYmbWKEUE6zeWJV/260pZs34jq9eVsmZ9xSSQJoaMhPDZulLWrNtcfvW6jaxZn9RRtTb03+N6/jDg5+xcsO6Lq2uYRGrCicTMGr1t/drPXJ58uW+ZHFav28iadVuWLy2r3tWeJoJWzQto1aKAnVs0Taeb0rldS1q1SKZbNS9g5xYFtG7RlJ2bF9C6RQE7N2+a/M1Y3qrF12m+/iBYciZExpmJWkLXR3KSRMCJxMzqmYq/9jN/2Wf+gl+d/rKvnV/7m7UoaLLFF3yrFgW0bVlAl3Yt0y/5pumXe/Jl36pFwabkUNl0y2ZN2MrtcNtn/UqSr/YmoBYQ65L5jctrbx8VOJGYWU7VpV/75V/uXdq3TH/BF9Cq/Mt+iy/6LZeXl925eVOaNa3jz7pd8QeINdDiENjtZvjoWlg3Lem9laNuwE4kZpaV9aVl3PvSO7z5/qdfuP6/Pb/2W29xiaeAdjs1o0v5ZZ7mef61Xx80aQe73Qq7XAZqAj2Oh0/ugM9fydkunUjMbLvN/2g1lz0yhZmLVtK7c1va7lRQ41/7yS/9ZFlBXf+1Xx90e3LLeTWFDlcCV+Zsl04kZlZjEcGf3nifG/86m52aNeX/zjucbxy0R77DsjxxIjGzGvl49TqufWwG49/6kGN6dWTk6Yewe9uW+Q7L8siJxMyq7aW5S7lq3DRWrNnAT7/dmyFH9aBJk0bWBmFf4ERiZtu0dsNGbn7ubUa/+h777d6aB77XjwM7t813WFZHOJGYWZXeXrKSS/88lTkfruL8o3pw3YkH0LJZ03yHZXWIE4mZVSoiGPPae4x49m3atmzG6CFH8NX9d8t3WFYHOZGY2Rd8tGotV4+bzktzl3LCAbtx82l96Ni6Rb7DsjrKicTMtvDC7A+59rHpfLaulBtP+RLnfnnPxndTn9WIE4mZAfD5+o387zOzeejf79O7c1vuOrsv++7WJt9hWT3gRGJmzFy0gh89PIUFSz/jB8fuzRVf348WBW5Qt+pxIjFrxMrKgvteWcBtf59Dh1YteOiCL3P0vh3zHZbVMzl9sI2kgZLmSJov6bpK1u8laYKk6ZJelNStwvq2khZJ+k3GssMlzUjrvEu+eGu2XT5Y8Tnn/P7f3PTs2ww4cHeevfQYJxHbLjlLJJKaAncDJwK9gbMl9a5QbCTwQET0AW4ARlRYfyPwUoVlvwWGAr3S18BaDt2swXtm+gcMvOMVppUs55bBfbjnnMPYpVXzfIdl9VQuz0j6AfMjYkFErAceBk6uUKY3UD5a/cTM9ZIOB3YH/p6xrDPQNiJej2Sw+QeAU3J3CGYNy+p1pVw9bhoX/+lNenRsxTM/OoYzjujuXlmWlVwmkq7Awoz5knRZpmnA4HT6VKCNpA6SmgC3AVdXUmfJNuoEQNJQScWSipcuXbqdh2DWcLz5/qd8665XeOzNEn54/L48Oqw/PTu2yndY1gDkMpFU9hOn4rBmVwHHSZoCHAcsAkqBi4C/RcTCCuWrU2eyMOK+iCiMiMJOnTrVLHKzBqR0Yxl3TZjH6fe+TunG4OGh/bny6/vX/ZH+rN7IZa+tEqB7xnw3YHFmgYhYDAwCkNQaGBwRKyT1B46RdBHQGmguaTVwZ1rPVus0s80WfrKGyx+ZSvF/PuXkvl248ZQv0bZls3yHZQ1MLhPJJKCXpJ4kZxpnAf+dWUBSR+CTiCgDhgOjACLinIwy5wOFEXFdOr9K0pHAv4HvAL/O4TGY1VtPTCnhp0/OQsAdZ/bllEMrvQpslrWcJZKIKJV0CfA80BQYFRGzJN0AFEfEU0ARMEJSAC8DF1ej6guBMcBOwLPpy8xSKz7fwPV/mclfpi6mcK9duP3MvnTfded8h2UNmJLOTw1bYWFhFBcX5zsMs5x7491PuPyRqSxZuZbLTujFhUX7eBx0226SJkdE4bbK+c52swZgw8Yy7hw/j3tenE/3XXfm0WH9OXTPXfIdljUSTiRm9dy7yz7jsoenMK1kBWcUduP6/zqI1i38X9t2HP9rM6unIoKxxQv5+dOzada0CfeccxjfPLhzvsOyRsiJxKwe+vSz9Qx/fAbPzVrCUft04LYzDqFzu53yHZY1Uk4kZvXMq/OXccXYqXzy2Xp+/M0DuOAre9OkiR9xYvnjRGJWT6wr3cjI5+fwu1feZZ9OrfjDd4/gS13b5TssMycSs/pg3oer+NHDU3nrg5Wce+Se/M83e7NTcw88ZXWDE4lZHRYRPPiv//C/z7xFqxYF/P47hQzovXu+wzLbghOJWR21bPU6rnl0Ov94+yOO268Tt57eh93atMx3WGZf4ERiVgdNfPsjrn50GivXlvKz/+rNd4/q4TFDrM5yIjGrQ9Zu2MiIv73F/a//hwP2aMNDFxzJ/nu0yXdYZlVyIjGrI2YvXsmlD09h3ker+f5XenL1N/anZTM3qFvd50RilmdlZcGoV9/llufm0G7nZjzwvX4cu58HY7P6w4nELI8+XLmWq8ZN45V5y/ha7925eXAfdm3VPN9hmdWIE4lZnjw3cwnDH5/O2g1l/PLUgzm7X3c3qFu95ERitoOtWV/KjX+dzZ/fWMjBXdtxx1l92adT63yHZbbdnEjMdqDpJcu57OGpvPvxZ1xYtA+XD9iP5gUeeMrqt5z+C5Y0UNIcSfMlXVfJ+r0kTZA0XdKLkrplLJ8saaqkWZKGZWzzYlrn1PS1Wy6Pwaw2bCwL7p44n0H3vMbnGzbypwuO5NqBBziJWIOw1TMSSTOArY7DGxF9qqpYUlPgbuBrQAkwSdJTETE7o9hI4IGIuF/S8cAI4DzgA+CoiFgnqTUwM912cbrdORHhsXOtXli0/HMuf2Qqb7z7Cd86uDO/PPVg2u3cLN9hmdWaqi5tfTv9e3H694/p33OANdWoux8wPyIWAEh6GDgZyEwkvYHL0+mJwJMAEbE+o0wLcnzmZJYrT09bzI+fmEFZWTDy9EMYfFhXN6hbg7PVL+iI+E9E/Ac4OiKuiYgZ6es64BvVqLsrsDBjviRdlmkaMDidPhVoI6kDgKTukqanddyccTYCMDq9rPVTbeV/paShkoolFS9durQa4ZrVnlVrN3DFI1P54Z+nsO9urfnbpcdw2uHdnESsQarOL/1Wkr5SPiPpKKBVNbar7H9MxUtlVwHHSZoCHAcsAkoBImJhevlsX+C7ksofeXpORBwMHJO+zqts5xFxX0QURkRhp06+uct2nMn/+YRv3vUKT05dxKUn9GLcD/qzV4fq/Jcxq5+q02vr+8AoSe1IEsEK4HvV2K4E6J4x3w3IPKsgPcsYBJC2hQyOiBUVy0iaRZI0Ho2IRenyVZL+RHIJ7YFqxGOWU6Uby/j1P+bz63/Mo0v7nRg3rD+H77VrvsMyy7ltJpKImAwcIqktoIpf9FWYBPSS1JPkTOMs4L8zC0jqCHwSEWXAcGBUurwb8HFEfC5pF+Bo4FeSCoD2EbFMUjOSdpzx1YzHLGfe/3gNlz4yhSnvL2fQYV35+UkH0aalG9StcdhmIkkvKf0S6BIRJ0rqDfSPiD9UtV1ElEq6BHgeaAqMiohZkm4AiiPiKaAIGCEpgJfZ3LB/IHBbulzAyIiYIakV8HyaRJqSJJHf1fywzWpHRPD4m4u4/i8zadJE3HX2oZx0SJd8h2W2Qyliqz18kwLSs8Bo4H8i4pD0rGBK2k5RLxQWFkZxsXsLW+1asWYDP35yBs9M/4B+PXfl9jP70rX9TvkOy6zWSJocEYXbKledNpKOETFW0nDYdKaxMesIzeqx19/5mCvGTmXpqnVcM3B/fnDsPjRt4h5Z1jhVJ5F8lnbJDQBJR5I0uJs1OutLy7h9/FzufekdenRoxeMXHUWfbu3zHZZZXlUnkVwBPAXsI+lVoBNwWk6jMquD3lm6mksfnsLMRSs5u193fvrt3uzc3I+rM6vyf4GkJkBLkns89idp+J4TERt2QGxmdUJE8Oc3FnLjX2fTslkT7j33cAZ+aY98h2VWZ1SZSCKiTNJtEdEfmLWDYjKrMz75bD3XPjadF2Z/yDG9OjLy9EPYvW3LfIdlVqdU57z875IGA4/Htrp4mTUgL89dypXjprFizQZ+8q0D+d7RPWniBnWzL6huG0kroFTSWpLLWxERbXMamVmerN2wkVuem8OoV9+l126tuX9IP3p38T93s62pzp3tbXZEIGZ1wZwlq7j04Sm8vWQV3+2/F8O/eSAtmzXNd1hmdVq1upykjynpRdLwDkBEvJyroMx2tIjg/tfe45fPvk3blgWMPv8IvnqAx0wzq47qPCLlAuBSkocuTgWOBF4Hjs9taGY7xker1nL1uOm8NHcpxx+wG7ec1oeOrVvkOyyzeqM6ZySXAkcA/4qIr0o6APh5bsMy2zHGz/6Qax6bzmfrSrnx5IM498i9PGaIWQ1VJ5GsjYi1kpDUIiLelrR/ziMzy6HP12/kF3+bzYP/ep8DO7flrrP60mt3NweabY/qJJISSe1JhsF9QdKnVBhXxKw+mbloBZc+PIV3ln7G0GP35sqv70eLAjeom22v6vTaOjWd/JmkiUA74LmcRmWWA2Vlwe//uYBbn5/Drq2a8+D3v8xXenXMd1hm9V51Gtv3zJh9N/27B/B+TiIyy4EPVnzOlWOn8do7HzPwoD0YMehgdmnVPN9hmTUI1bm09QzJk39F0v23JzAHOCiHcZnVmmdnfMB1j89gw8Yybhnch9MLu7lB3awWVefS1hYDWEk6DPhBziIyqyWfrSvl50/PYmxxCYd0a8cdZx1Kz46t8h2WWYPTpKYbRMSbJN2Bt0nSQElzJM2XdF0l6/eSNEHSdEkvpmO1ly+fLGmqpFmShmVsc7ikGWmdd8k/La0SU97/lG/e9QrjJpdwyVf35dELj3ISMcuR6rSRXJEx2wQ4DFhaje2aAncDXwNKgEmSnoqI2RnFRgIPRMT9ko4HRgDnAR8AR0XEOkmtgZnptouB3wJDgX8BfwMGAs9u+1CtMdhYFtwzcT53TJjHHm1b8vD/dyRf3rtDvsMya9Cq00aS2bm+lKTN5LFqbNcPmB8RCwAkPQycDGQmkt7A5en0RJIuxkTE+owyLUjPnCR1BtpGxOvp/APAKTiRGLDwkzVcMXYqk977lJMO6cKNp3yJdjs1y3dYZg1eddpItvcu9q7Awoz5EuDLFcpMAwYDdwKnAm0kdYiIjyV1J0la+wJXR8RiSYVpPZl1dq1s55KGkpy5sOeee1ZWxBqQJ6cs4qdPziSAO87syymHVvrPwsxyoDqXtp6qan1EnLS1TSsrXmH+KuA3ks4HXgYWkZz1EBELgT6SugBPSnq0mnWWx3UfcB9AYWGhx1FpoFau3cBPn5zJX6YupnCvXbj9zL5033XnfIdl1qhU59LWuyT3jTyYzp8NvAc8v43tSoDuGfPdqHBHfNrmMQggbQsZHBErKpaRNAs4Bng1rWerdVrj8ca7n3D5I1NZsnItV3xtPy4q2oeCpjXuP2JmWapOIjk0Io7NmH9a0ssR8eNtbDcJ6CWpJ8mZxlnAf2cWkNQR+CQiyoDhwKh0eTfg44j4PH2E/dHAryLiA0mrJB0J/Bv4DvDrahyDNSAbNpZx14R53D1xPt133Zlxw/pz2J675Dsss0arOomkk6S9MxrNewKdtrVRRJRKuoTkzKUpMCoiZkm6ASiOiKeAImCEpCC5tHVxuvmBwG3pcgEjI2JGuu5CYAywE0kjuxvaG5HFyz/nwofeZNrC5Zx+eDf+30kH0bpFtYbVMbMc0baGYZc0kKStYUG6qAcwNCL+ntvQak9hYWEUFxfnOwzL0uzFKxky5g3WrNvITYP78K0+nfMdklmDJmlyRBRuq1x1em09J6kXcEC66O2IWJdtgGY18fLcpVz00Ju0aVnAuAv7c8AeHkPdrK7YasukpCMk7QGQJo5DgBuAWyXtuoPiM2PspIUMGTOJbrvsxBMXHe0kYlbHVNXF5f+A9QCSjgVuAh4AVpB2qzXLpYjgVy/M5ZrHpnPUPh0YN6w/e7Rrme+wzKyCqi5tNY2IT9LpM4H7IuIx4DFJU3MfmjVm60vLuO7x6Tz+5iJOP7wbvxx0MM3ctdesTqoykUgqiIhS4ATSu8SrsZ1ZVlau3cCFD07m1fkfc8XX9uOHx+/rx76b1WFVJYQ/Ay9JWgZ8DrwCIGlfkstbZrVu8fLPGTJ6Eu8sXc3I0w/htMO7bXsjM8urrSaSiPiFpAlAZ+DvsbmfcBPghzsiOGtcMrv3jhnSz8PgmtUTVV6iioh/VbJsbu7CscaqvHtv6xbu3mtW37itw/JubPFCfvz4DPbdrTWjhxxB53Y75TskM6sBJxLLm4jg9vHzuGvCPI7p1ZF7zjmMNi09fohZfVPj/pSSmko6JxfBWOOxvrSMq8ZN564J8zj98G6MOv8IJxGzeqqqO9vbShq+o5RiAAASCklEQVQu6TeSvq7ED0meuXXGjgvRGpqVazfwvTGTeOzNEi4fsB+3nNbH94iY1WNVXdr6I/Ap8DpwAXA10Bw4OSJ8Q6Jtlw9WJN1753/k7r1mDUVViWTviDgYQNLvgWXAnhGxaodEZg3O7MUr+d6YSaxeV+ruvWYNSFWJZEP5RERslPSuk4htry269w7rz4Gd3b3XrKGoKpEcImllOi1gp3ReQESEvwmsWty916xhq+rO9qY7MhBreCKCO8bP40537zVr0HLaVUbSQElzJM2XdF0l6/eSNEHSdEkvpmO1I6mvpNclzUrXnZmxzRhJ70qamr765vIYbPuUd++90917zRq8nN2QKKkpcDfwNaAEmCTpqYiYnVFsJPBARNwv6XhgBHAesAb4TkTMk9QFmCzp+YhYnm53dUQ8mqvYLTsr127gogff5J/zl3H5gP340Ql+eq9ZQ5bLO9v7AfMjYgGApIeBk4HMRNIbuDydngg8CVs+zysiFkv6COgELMfqNHfvNWt8cnlpqyuwMGO+JF2WaRowOJ0+FWgjqUNmAUn9SO5feSdj8S/SS163S2pR2c4lDZVULKl46dKl2RyHVdPsxSs59e7XKPn0c8YM6eckYtZI5DKRVHYtIyrMXwUcJ2kKcBywCCjdVIHUmeTGyCERUZYuHg4cABwB7ApcW9nOI+K+iCiMiMJOnTpldSC2ba/MW8oZ//c6AOOG9fc9ImaNSC4vbZUA3TPmuwGLMwtExGJgEICk1sDgiFiRzrcFngF+kvk4+4j4IJ1cJ2k0STKyPHL3XrPGLZdnJJOAXpJ6SmoOnAU8lVlAUkdJ5TEMB0aly5sDT5A0xI+rsE3n9K+AU4CZOTwGq0JEcPsLc7nm0ekcuXcHxg7r7yRi1gjl7IwkIkolXQI8DzQFRkXELEk3AMUR8RRQBIyQFMDLwMXp5mcAxwIdJJ2fLjs/fcbXQ5I6kVw6mwoMy9Ux2NatLy3jx0/M4NHJJZx2eDdGDDrYD140a6S0eQTdhquwsDCKi4vzHUaDsWrtBi5Mu/deNqAXl57Qy917zRogSZMjonBb5TywldVIZvfeW0/rw+mF3be9kZk1aE4kVm1vfbCSIaOTp/eOHnIEx/RybzgzcyKxanpl3lIufNBP7zWzL3IisW0aV7yQ4e7ea2Zb4URiWxUR3DlhHneMn8dX9u3IPeceRls/eNHMKnAisUpt2FjG8MeT7r2DD0u69zYvcPdeM/siJxL7glVrN3DRQ2/yyrxlXHpCLy4b4O69ZrZ1TiS2hczuvbec1ocz3L3XzLbBicQ2yezeO+r8Izh2P3fvNbNtcyIxYMvuvWN/0J/eXdy918yqx4nE3L3XzLLiRNKIuXuvmdUGJ5JGyt17zay2OJE0Qu7ea2a1yYmkkXH3XjOrbU4kjYi795pZLjiRNBL/nLeMCx+czM4tmrp7r5nVqpy2rkoaKGmOpPmSrqtk/V6SJkiaLulFSd3S5X0lvS5pVrruzIxtekr6t6R5kh5Jx3e3KowrXsj5o9+g6y478cRFRzuJmFmtylkikdQUuBs4EegNnC2pd4ViI4EHIqIPcAMwIl2+BvhORBwEDATukNQ+XXczcHtE9AI+Bb6fq2Oo7yKCO8fP4+pHp/PlvXdl7LD+dGnve0TMrHbl8oykHzA/IhZExHrgYeDkCmV6AxPS6Ynl6yNibkTMS6cXAx8BnZR0LToeeDTd5n7glBweQ721YWMZ1zw6ndvHz2XwYd0YfX4/3yNiZjmRy0TSFViYMV+SLss0DRicTp8KtJHUIbOApH5Ac+AdoAOwPCJKq6izfLuhkoolFS9dujSrA6lvVq3dwPfGTGLc5BJ+dEIvRp7ex/eImFnO5PLbpbIbE6LC/FXAcZKmAMcBi4DyJIGkzsAfgSERUVbNOpOFEfdFRGFEFHbq1Hh6Jy1ZsZbT732d19/5mFsG9+GKr+3ne0TMLKdy2WurBMi8SaEbsDizQHrZahCApNbA4IhYkc63BZ4BfhIR/0o3WQa0l1SQnpV8oc7G7O0lSffelZ9vcPdeM9thcnlGMgnolfayag6cBTyVWUBSR0nlMQwHRqXLmwNPkDTEjysvHxFB0pZyWrrou8BfcngM9cY/5y3j9N++TlkEY4f1dxIxsx0mZ4kkPWO4BHgeeAsYGxGzJN0g6aS0WBEwR9JcYHfgF+nyM4BjgfMlTU1ffdN11wJXSJpP0mbyh1wdQ33x6OQSzh/9Bl3aJ917D+rSLt8hmVkjouRHfsNWWFgYxcXF+Q6j1kUEd02Yz+3j53L0vh347bmHu2eWmdUaSZMjonBb5Xxnez21YWMZP358BuMmlzDosK7cNMg9s8wsP5xI6qHMp/f+6IReXO6n95pZHjmR1DNLVqzl/NFvMO+j1dwyuA9nHOGn95pZfjmR1CMVu/ce555ZZlYHOJHUE6/OX8awP6ZP7x3W3z2zzKzOcCKpBx6dXMJ1j01nn06tGT3kCD940czqFCeSOszde82sPnAiqaM2bCzjf56Ywdhid+81s7rNiaQOcvdeM6tPnEjqmCUr1jJkzCTmfrjK3XvNrF5wIqlD3L3XzOojJ5I6wt17zay+ciKpAx6bXMK17t5rZvWUE0keRQS//sd8fvXCXI7aJ+ne224nd+81s/rFiSRP3L3XzBoKJ5I82KJ77/H7crnHVTezesyJZAfL7N578+CDOfOIPfMdkplZVnJ6LUXSQElzJM2XdF0l6/eSNEHSdEkvSuqWse45Scsl/bXCNmMkvVvJELx13ttLVnLqPa/y/sefMer8I5xEzKxByFkikdQUuBs4EegNnC2pd4ViI4EHIqIPcAMwImPdrcB5W6n+6ojom76m1nLoOfHq/GWc/tvX2VgWjB3W3/eImFmDkcszkn7A/IhYEBHrgYeBkyuU6Q1MSKcnZq6PiAnAqhzGt8M8NrmE7456g87tW/LExUf7HhEza1BymUi6Agsz5kvSZZmmAYPT6VOBNpI6VKPuX6SXw26X1KKyApKGSiqWVLx06dKaxl4rkqf3zuPKcdPo13NXxg07iq6+R8TMGphcJpLKuiFFhfmrgOMkTQGOAxYBpduodzhwAHAEsCtwbWWFIuK+iCiMiMJOnXb8ZaQNG8u47rEZ/OqFuQw6tCtjhvTzPSJm1iDlstdWCZD5xMFuwOLMAhGxGBgEIKk1MDgiVlRVaUR8kE6ukzSaJBnVKavWbuDiP03h5blL3b3XzBq8XCaSSUAvST1JzjTOAv47s4CkjsAnEVFGcqYxaluVSuocER8o+WY+BZhZ65Fn4cOVazl/tLv3mlnjkbNEEhGlki4BngeaAqMiYpakG4DiiHgKKAJGSArgZeDi8u0lvUJyCau1pBLg+xHxPPCQpE4kl86mAsNydQw1NWfJKoaMfoMVfnqvmTUiiqjYbNHwFBYWRnFxcU738dr8Zfzgj5PZqXlTRg85wj2zzKzekzQ5Igq3Vc53tteCx99Mnt7bs2MrRg/p555ZZtaoOJFkISL4zT/mc9sLc+m/dwfuPc9P7zWzxseJZDtt2FjGT56YySPFCxl0aFduGuyn95pZ4+REsh1Wryvloofe5OW5S/nh8ftyhbv3mlkj5kRSQx+uXMuQ0ZOY8+Eqbhp0MGf1c/deM2vcnEhqILN77x++W0jR/rvlOyQzs7xzIqmmzO69j/ygP1/q6u69ZmbgRFIt7t5rZrZ1TiRVcPdeM7NtcyLZiojgx0/M5M9vvM+ph3blZnfvNTOrlBPJVkhin06t3L3XzGwbnEiqcMExe+c7BDOzOs/XaszMLCtOJGZmlhUnEjMzy4oTiZmZZcWJxMzMspLTRCJpoKQ5kuZLuq6S9XtJmiBpuqQXJXXLWPecpOWS/lphm56S/i1pnqRHJDXP5TGYmVnVcpZIJDUF7gZOBHoDZ0vqXaHYSOCBiOgD3ACMyFh3K3BeJVXfDNweEb2AT4Hv13bsZmZWfbk8I+kHzI+IBRGxHngYOLlCmd7AhHR6Yub6iJgArMosrOSuwOOBR9NF9wOn1H7oZmZWXbm8IbErsDBjvgT4coUy04DBwJ3AqUAbSR0i4uOt1NkBWB4RpRl1dq2soKShwNB0drWkOTU/BAA6Asu2c1vLDX8mdZM/l7on289kr+oUymUiqeyZIlFh/irgN5LOB14GFgGlFTeqYZ3Jwoj7gPu2HWbVJBVHRGG29Vjt8WdSN/lzqXt21GeSy0RSAnTPmO8GLM4sEBGLgUEAkloDgyNiRRV1LgPaSypIz0q+UKeZme1YuWwjmQT0SntZNQfOAp7KLCCpo6TyGIYDo6qqMCKCpC3ltHTRd4G/1GrUZmZWIzlLJOkZwyXA88BbwNiImCXpBkknpcWKgDmS5gK7A78o317SK8A44ARJJZK+ka66FrhC0nySNpM/5OoYUllfHrNa58+kbvLnUvfskM9EyY98MzOz7eM7283MLCtOJGZmlhUnkq2QNErSR5Jm5jsWS0jqLmmipLckzZJ0ab5jMpDUUtIbkqaln8vP8x2TJSQ1lTSl4qOmapsTydaNAQbmOwjbQilwZUQcCBwJXFzJY3dsx1sHHB8RhwB9gYGSjsxzTJa4lKSzU045kWxFRLwMfJLvOGyziPggIt5Mp1eR/Aep9MkGtuNEYnU62yx9uRdPnqUPwf0W8Ptc78uJxOolST2AQ4F/5zcSg02XUKYCHwEvRIQ/l/y7A7gGKMv1jpxIrN5Jn4LwGHBZRKzMdzwGEbExIvqSPG2in6Qv5TumxkzSt4GPImLyjtifE4nVK5KakSSRhyLi8XzHY1uKiOXAi7h9Md+OBk6S9B7Jk9ePl/RgrnbmRGL1RjqMwB+AtyLiV/mOxxKSOklqn07vBAwA3s5vVI1bRAyPiG4R0YPk8VT/iIhzc7U/J5KtkPRn4HVg//QRLR5AK/+OJhns7HhJU9PXN/MdlNEZmChpOskz9l6IiJx2N7W6xY9IMTOzrPiMxMzMsuJEYmZmWXEiMTOzrDiRmJlZVpxIzMwsK04k1mBICkm3ZcxfJelntVT3GEmnbbtk1vs5PX268cQKy3tI+jyj2/PUdAjrmtZ/maSday9iMycSa1jWAYMkdcx3IJkkNa1B8e8DF0XEVytZ905E9M14rd+OcC4DapRIahi/NUJOJNaQlJKMUX15xRUVzygkrU7/Fkl6SdJYSXMl3STpnHR8jRmS9smoZoCkV9Jy3063byrpVkmTJE2X9IOMeidK+hMwo5J4zk7rnynp5nTZ9cBXgHsl3VqdA5bUT9Jr6ZgTr0naPyOukek+pkv6oaQfAV1Ibh6cuLU4yt8fSTdI+jfQP31fZqd1jaxObNaIRIRffjWIF7AaaAu8B7QDrgJ+lq4bA5yWWTb9WwQsJ7k7uwWwCPh5uu5S4I6M7Z8j+fHVCygBWgJDgZ+kZVoAxUDPtN7PgJ6VxNkFeB/oBBQA/wBOSde9CBRWsk0P4HNgavq6O13eFihIpwcAj6XTF5I8k6x83a7p3/eAjtWII4AzyrcF5rD5Bub2+f6s/apbr4JqZRuzeiIiVkp6APgRyRdvdUyKiA8AJL0D/D1dPgPIvMQ0NiLKgHmSFgAHAF8H+mSc7bQjSTTrgTci4t1K9ncE8GJELE33+RBwLPDkNuJ8J5In7GZqB9wvqRfJl3+zdPkA4N6IKAWIiMrG1qkqjo0kiQhgJbAW+L2kZwA//sS24Etb1hDdQdLW0CpjWSnpv/f04Y+ZDdXrMqbLMubLYIsfWxWfJxSAgB/G5naLnhFRnog+20p8qu6BVMONwMSI+BLwXyRnSeX72Nbzj6qKY21EbARIk1E/ksRyCsmZmdkmTiTW4KS/vseSJJNy7wGHp9Mns/mXe02cLqlJ2m6yN8nlnueBC9PH2yNpP0mtqqqEZDCu4yR1TBuyzwZe2o54IDkjWZROn5+x/O/AMEkFaVy7pstXAW1qEkc6/ku7iPgbSWN9xbMia+ScSKyhug3I7L31O5IvzTeAL7P1s4WqzCH5on0WGBYRa0mGMZ0NvClpJvB/UPUl4/Qy2nBgIjANeDMi/rId8QDcAoyQ9CqQ2bvq9yTtH9MlTQP+O11+H/CspIk1iKMN8Nf06b4vUUlnBmvc/PRfMzPLis9IzMwsK04kZmaWFScSMzPLihOJmZllxYnEzMyy4kRiZmZZcSIxM7Os/P+/X88RBDFXmgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#viz\n",
    "plt.plot(range(1,iris.data.shape[1]+1),rsquared)\n",
    "plt.plot([ind+1],[rsquared[ind]],lw=0,marker='*',\n",
    "         markersize=10,color='gold') \n",
    "plt.hlines(threshold,1,iris.data.shape[1],linestyles=':')\n",
    "plt.xticks(range(1,iris.data.shape[1]+1))\n",
    "plt.title('Optimal Number of Factors')\n",
    "plt.xlabel('Number of Factors')\n",
    "plt.ylabel('R Squared')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The optimal number of factors suggested by Variance Explained Criteria is 4\n"
     ]
    }
   ],
   "source": [
    "#apparently nobody hits the threshold, lol\n",
    "print(f'The optimal number of factors suggested by Variance Explained Criteria is {ind+1}')"
   ]
  },
  {
   "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
}