PowerLawExample-with-diff-Normalizations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find  the relation between\n",
    "\n",
    "The Stable Rank (in log units) and the Power Law exponent (we fit) $\\alpha$\n",
    "\n",
    "$$\\dfrac{\\log_{10}\\Vert W\\Vert^{2}}{\\log_{10}\\lambda_{max}}\\;\\;\\;vs\\;\\;\\alpha$$\n",
    "\n",
    "We do this by \n",
    "\n",
    "- generating a random ($N\\times M$) Pareto matrix $\\mathbf{W}(\\mu)$ for $\\mu\\in[1,5]$\n",
    "\n",
    "- computing the correlation matrix $\\mathbf{X}$ and its eigenvalues $\\lambda_{i}$\n",
    "\n",
    "- fitting the density  to a power law: $\\rho(\\lambda)\\sim\\dfrac{1}{\\lambda^{\\alpha}}$\n",
    "\n",
    "where $Q=N/M1=1$ and we normalize $\\mathbf{X}$ by $1/N$\n",
    "\n",
    "$$\\mathbf{X}=\\dfrac{1}{N}\\mathbf{W}^{T}\\mathbf{W}$$\n",
    "\n",
    "\n",
    "\n",
    "#### We also examine other normalizations\n",
    "\n",
    "- $\\dfrac{1}{N^{2/\\mu}}\\mathbf{X}$\n",
    "\n",
    "- $\\dfrac{1}{N^{1/\\mu}}\\mathbf{W}$\n",
    "\n",
    "- $\\dfrac{1}{\\max|\\mathbf{W}|}\\mathbf{W}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-09T21:53:31.334410Z",
     "start_time": "2018-11-09T21:53:29.413713Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "\n",
    "import powerlaw\n",
    "\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mus = [0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 4, 4.5, 5]\n",
    "Ms = [100, 500, 1000, 5000, 10000]\n",
    "Mcolors=['blue','green', 'teal', 'purple', 'red']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### also create table for paper here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-09T22:54:39.475194Z",
     "start_time": "2018-11-09T22:42:38.351508Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "927fb82c3bd74cf69e8c18e6114d36af",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=16), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating best minimal value for power law fit\n",
      "/Users/charleshmartin/anaconda3/lib/python3.6/site-packages/powerlaw.py:700: RuntimeWarning: invalid value encountered in true_divide\n",
      "  (Theoretical_CDF * (1 - Theoretical_CDF))\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "59775a4ebe6b4c0ab72bd0f6ef060f28",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=16), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "917d50a21b4546e3a333b23bc0858d80",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=16), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating best minimal value for power law fit\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3eef5bb11d2843b587954a9b8a059663",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=16), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating best minimal value for power law fit\n",
      "/Users/charleshmartin/anaconda3/lib/python3.6/site-packages/powerlaw.py:1141: RuntimeWarning: invalid value encountered in log\n",
      "  self.alpha = 1 + (self.n / sum(log(data/self.xmin)))\n",
      "/Users/charleshmartin/anaconda3/lib/python3.6/site-packages/powerlaw.py:1160: RuntimeWarning: invalid value encountered in log\n",
      "  return 1 + len(data)/sum(log(data / (self.xmin)))\n",
      "/Users/charleshmartin/anaconda3/lib/python3.6/site-packages/powerlaw.py:701: RuntimeWarning: overflow encountered in double_scalars\n",
      "  )[1:]\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "/Users/charleshmartin/anaconda3/lib/python3.6/site-packages/powerlaw.py:700: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  (Theoretical_CDF * (1 - Theoretical_CDF))\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "47a79ff90d5048db860a325a7e976642",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=16), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n",
      "Calculating best minimal value for power law fit\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "from IPython.core.display import display, HTML\n",
    "html = '<table><tr><th>mu</th><th>alpha</th><th>logNorm</th><th>formula</th><th>ratio</th></tr>'\n",
    "Mresults = {}\n",
    "l_maxs = {}\n",
    "logNorm2s = {}\n",
    "\n",
    "for M in Ms:\n",
    "    ratios, alphas = [], []\n",
    "    l_maxs[M] = []\n",
    "    logNorm2s[M] = []\n",
    "    \n",
    "    N = M\n",
    "    for mu in tqdm(mus):\n",
    "\n",
    "        # NO W NORMALIZATION \n",
    "        W = np.random.pareto(a=mu,size=(N,M))\n",
    "        normW = np.linalg.norm(W)\n",
    "        logNorm2 = 2.0*np.log10(normW)\n",
    "        logNorm2s[M].append(logNorm2)\n",
    "\n",
    "        # STANDARD NORMALIZATION \n",
    "        X=np.dot(W.T,W)/N\n",
    "        evals = np.linalg.eigvals(X)\n",
    "\n",
    "        l_max, l_min = np.max(evals), np.min(evals)\n",
    "        fit =  powerlaw.Fit(evals)\n",
    "\n",
    "        a = fit.alpha\n",
    "        alphas.append(a)\n",
    "\n",
    "        ratio = logNorm2/np.log10(l_max)\n",
    "        ratios.append(ratio)\n",
    "\n",
    "        Mresults[M]=[ratios,alphas]\n",
    "        \n",
    "        l_maxs[M].append(l_max)\n",
    "\n",
    "            #print(\"mu= {} 1+mu/2 = {}, alpha= {:2f} logNorm2= {:2f}  formula= {:2f} ratio= {:2f}\"\n",
    "            #      .format(mu,(1+mu/2),fit.alpha, logNorm2, formula, ratio ))\n",
    "\n",
    "    #html += \"<tr><td>{:.3f}</td><td>{:.3f}</td><td>{:.3f}</td><td>{:.3f}</td><td>{:.3f}</td></tr>\".format(mu, fit.alpha, logNorm2, formula, ratio)\n",
    "\n",
    "#html += \"</table>\"\n",
    "#display(HTML(html))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/charleshmartin/anaconda3/lib/python3.6/site-packages/numpy/core/numeric.py:553: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return array(a, dtype, copy=False, order=order, subok=True)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x10e3d53c8>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmwAAAJwCAYAAAA0gT3FAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl8VNX9//HXgewKJGETUEgqggqRxRAJW5RVa4n+sIolAWsVbC20YIvlKxUUhKC21mpbq1BogbQNSEWwIiJL4gY1IEqBFtEQBFtQViEsWc7vjztJs5OQydzJ5P18PO5jmHPP3POZZGA+nO0aay0iIiIi4r+auB2AiIiIiFRPCZuIiIiIn1PCJiIiIuLnlLCJiIiI+DklbCIiIiJ+TgmbiIiIiJ9TwiYiIiLi55SwiYiIiPg5JWwi4leMMY8YY7YZY04ZY/5jjPmTMeYyt+MSEXGT0Z0ORMSfGGNeB/4MbAWaAc8BRdbafq4GJiLiIvWwidSBMeZGY4w1xjzmdiwXwxjzR2OM9fz5Mc97iXGzPWvtN621S621u621/wAmAYnGmDY1uH6D/n1UxxjzS2NMkTEmxxgzxu14RMS3lLCJSJWMMbGeBOj1Ks6/6Dn/dhXnf+05X5cEowVQAJysoo2bPW3cXoc2GoIPgJeATsCLxphQl+MRER9SwiYi1TnueWxR/oQxpiUwtprzzYF7gYPA8otp3JOUpAFLrbVnq6iWDJwF1l1MGw2Ftfav1trvA38HLgWudTkkEfEhJWwiUp0TgAWaV3Lu+0A48DmVJGzAfThz0J631ubXtmFjTBNgMRAE/Liaqt8C1ltrT9e2jQbqA89jnKtRiIhPKWETqQfGmGBjzE+NMTuMMWeMMceMMW8YYwZWUf9SY8yvjDFfeOpvM8bcaYz5rme477s+fgsAWGuLgK8pl5AZY4KBHwLbgDcrOd8EZ+7ZaZxhvFoxxhhgIXAdMNxaW9VwaG/gCmBVDa4ZEL8TYL/nsbtL7deIMSbJ83P6TRXn+3nOv1iq7C5jzDvGmK88P/NcY8xKY8wAL8XUwhgz1xjzqef6u4wxozznJhhjCowxV3qjLRFvU8Im4mWeZONl4Gmcv2PPe573AzYaY+4sV78p8DowGWf48FngY5zepW/7LvIqnaBiD9t3gHbAL3HmljX3vO9itwGxwB+ttcdq05jnOvOBgcBQa+2X1VRPxukBXF2Dazb434knUf6J56m/97Bl4SSXdxljgio5XzyvMR3AGPNDIAO4DPgr8GsgE+gF3FjXYDyLW7YCD3liewnoAKQbY64AfoYz9P5pXdsSqQ+V/SUSkboZh5NIvAncaq0tADDGPANkA/ONMW9Ya7/21P8eTnKyHBhtPXvtGGP+gPOF5bbjQAdjjLH/2wdoMnAAWAZ0AQzO8GdxT9gUoAgn0amtF3ASvm8BhaX2YDtqrT1frm4ykG2t/c8Frhkov5OH+V/Pml8nbNZaa4z5C04iNBwnAQbAk8DdhZPQFS9YuQ8nOb7OWptXqq4BouoSiycB/xvOfyIGW2szPeXbgD8Cv8BZzDGiLu2I1Cf1sIl43z2ex6nFiQGAtXY3Ts9RC6D0isYUz+PPSyVEWGvfBt6o51hr4jjOvxWXgrN1Bk6vx6897684SWvhOd8LJ9lZba3dexHtPQC0AjYD/yl1lNmHzRhzuSeOCw6HEgC/E2PMVcDPcX4uH+Ek0XVKZHxgqeex/Crh4UBr4M+lf77AeZwVwSWs42gd4xiF81n5Q3Gy5pHtebwLp3ftYj6vIj6hhE3E+3oCx6y1H1dyblOpOsV6AEestXsqqf++l2O7GCc8j8Xz1KbgzGub73l+spLzAM9cTGPWWlPFsalc1WTPY00StkD4nfweaArcD2z3lPn1PDZr7T9xhpJvN8ZElDpVnBAvLVWWgdMD9k9jzGxjzFBjzCVeCuU7nsfy8+nOeR4LgTleakukXihhE/G+5sChKs79t1SdYs2AquZpHS5fYIwZZYx50xhzxFSx0a1ncv0Bz8TqdcaYb9Q8/AqKt/ZobozpjDNU+QdrbXEiV5ywRXqGL0cD26y1WXVosyaSgdwqkrDyGvTvxLPAYTCQZq3diZMEgZ8Pi3qkA5fgDHPjScJuA7Z73kuxp3B6V0/j9CSuA74yxizyQk/iIOBgNZ+VdGvtJ3VsQ6ReKWET8b6TQNsqzrUtVafY1zjDQ5WpbHf/S4B3gEcqe4ExZhwwC2dyegLOF+BrVUz8ronSe7FNxpnk/+tS50v3sD0IhHCRvWs1ZYy5FGciek1616AB/06MMa1w5ljt5n+9QMWJh1/3sHn8GeczU9yrdjvOzyu9dCXP0OdL1tpeOAsPRuPMF/wuzjyzi2Kc/QJbAv+q5PQtnsdKN4YW8SdK2ES8bzsQZYyp7Ms0qVSdYh8BLY0xXSqpn1i+wFq7xFo7C1hfRfs/wtn7LMNauwNn/taVOPOGLkZxwtYR58tzhbV2X6nzxYlOW5y92Q7iLEa4aMaYnZ6eqtLHllJVbgZCqXnC1pB/J78CooHxpRZdfOR5LNPDZozpaZzbVyUbY7YbY/KMMa8YY0KNMT80xuwzxhw3xvxfudf90hjzb0/v335jzKRS5x7w9BwWz2EMM8ZsNsb8qQaxY609gLMqc7gneUrBWZDyl2pec8hauwz4JvAJcEsd/sNRPBRb5sbZxphoYKbnadFFXlvEZ5SwiXjfYs9jmmd1GgCeL/8JOHPCXi1Vv/iLa3bprTE8e0/dXJuGjXNngJ7AhuIyz9DlB0Df2lyrlOKhz5/i9Iz8stz54oTtBzi9Ur+5mI1yixljwoGuOF/s7UodpVfwJXviqumKzQb5OzHGDAVSgReste+Wev2XOEO85RPQ63Am7d+Hc5eJ23CSnjdw5oeNAB4HnigeZjTGhOG8/1TgGmAG8IwxppvnmgtxfscPen4Wi4EzwPha/BiWAsV79w0DNlprD5Z7r8NL/248InAWu5ynVFJljLnaGHN1Ddv+r+f1CcaY9p7XNwFexOl5g6p7X0X8hrb1EPG+4r26vgV8aIxZg7MtwWicOwOklNsI9g84207cBcQaY9bjJCijcYZqbqXmPQAtcSaml59ndZiL/1Iq7mGLB97x3JC9tJOlzp/G+SKsizic95Bprf1v+ZOeL/VvAm/UIjFscL8TT+L6e5ztU/6vkiofA8OMMZd7erHASdgOAndba894rvMpzgKKn3qef4UzZN0KZyHGWZzh2mJ/9PTAdQF2WmvzjTFPAHNxhirjgH6VbLFSnZdxJvz/HOd7J72SOsuAU8aYd4BcnGTtVpyf+xOeTZyL7fY8Gi7AE/9fcHo13zHGrMIZlk4EpuHc+myqMaaFtVYLD8RvqYdNxMs8Xyz/D+fLoCnObZXuwtmOYbC19q/l6hfg9No8h7Nr/2ScVYrjgI2eal9TMxf8ArsIx0v9ubK5aSdK/bnWG+VWohfwZfkemFL64yRBNR0Obai/kxk4w6YPVnGnh8qGRa8DMoqTNY/2lJ0D1hFnVeTnAMaYK40xC4wx/zLGnDDGnAKuwtlKpdifgHyce8feWtvfsbX2OE6iG4xz39cVlVT7P+BDnETqxzjJ8QGc5PPR2rRXiUk4ye8lOPMsLwPGWmufxJl3GIl/bFItUjVrrQ4dOvz0AJbgzL25tpJznT3nYkqVheIMiY0oV/cdYFYl1/ij88+ABXis/PXq4f1csD2cL9ZC4FSp49elzv8CJ3mICsTfSR1j+y8wqtTzGE88bUqV3QPs8vy5HXAUZ4uWG4GrcfYsKwIuLfWaMThbYBwBmrnxc9eho7Ef6mET8QPGmHaVlA0A7saZdL27wosqYa09hzN5/qZS12kO9MHpTWoIeuHc7aBnqeOxUudHAm/buvfkVauh/U6MMa1xhlhLb13RA/iPtbb0cOx1peokA8etteOttZustf/Cmee211p7ynPd/jjD3N/C2epksrdiFpGa0xw2Ef8w3zMh+h84c8Ku5n/zpH5krS1Z4eZZ3dYRZ6gO4FpjTCT/+5J9DnjBOLfd2QXMBj7FuS2TX/PMT4sDfmmr2HXeWtvVR+E0tN9JD5zeyE/LlX1Urt51/G8BxFHgMmPMEJy5b9/F2WR2red9XQmsBCZba9cZY+YCvzbGPG+dYU4R8RH1sIn4h2U4K+++jXOngH44c36SrLXlb4WUjDPXp3gO1989z+MBrLWLcbYr+BXOrXcuBb5lS92SyY91xVkEUD7JcEND+53EAf8snUjiJGzby9Ur3cP2N5z7pb6KkzwewUlQt3tWkf4dWGit/YOnfjrwFc4N1EXEh0zZv9si0pgYY24HelprHzPOPUJvBJ6tr96TC7VnjEnBGX5rbsuuChQRadSUsImI3zDG/ALob62tsDmtiEhjFnAJW6tWrWxMTIzbYYiIiIhc0NatW7+y1lZ1K7wSAbfoICYmhuzsbLfDEBEREbkgY0xuTepp0YGIiIiIn1PCJiIiIuLnlLCJiIiI+DklbCIiIiJ+TgmbiIiIiJ9TwiYiIiLi5wJuW4+aOHnyJIcPHyY/P9/tUERqLDg4mDZt2tC8eXO3QxERER9rdAnbyZMnOXToEB06dCA8PBxjjNshiVyQtZYzZ85w8OBBACVtIiKNTKMbEj18+DAdOnQgIiJCyZo0GMYYIiIi6NChA4cPH3Y7HBER8bFGl7Dl5+cTHh7udhgiFyU8PFxD+SIijVCjS9gA9axJg6XProhI49QoEzYRERGRhkQJm4iIiIifU8IWIN58801uueUWWrZsSVhYGF27dmXatGkcP37cK9c/cOAAkyZNIjExsWTBxr59++p0zZiYGIwxjBkzptLzN954I8YYBgwYUKd2qlKb93T27FmmTp1Ku3btCA8PJzExkaysrAr1ioqKSEtLIyYmhrCwMHr06MGKFSvqJX4REWk8lLAFgLlz5zJixAjCwsJYsGABa9eu5YEHHmDRokUkJCSUbAVRF3v37mXZsmVERUUxcOBAL0TtaNasGStXruTrr78uU56bm0tWVhbNmjXzWlvl1eY93XfffcyfP59Zs2bx2muv0a5dO0aMGMH27dvL1Hv00Ud57LHHmDhxImvWrKFv377ceeedvP766/X2PkREpBGw1gbUcf3119vq7Nq1q9rzDc2GDRusMcZOnjy5wrnPPvvMRkVF2WHDhtW5ncLCwpI/z58/3wI2JyenTtfs1KmTHTVqlG3evLldtGhRmXOzZ8+2sbGxdtCgQbZ///51aqcqNX1P27dvt4BduHBhSVl+fr7t0qWLHTlyZEnZoUOHbEhIiJ0xY0aZ1w8ePNjGxcV5Le5A+wyLiDRmQLatQX6jHrYG7qmnniI6Opq0tLQK52JjY5k2bRrr1q1j69atdWqnSZP6+aiEh4dzxx13sGTJkjLlS5YsYezYsfW6KrKm72nVqlUEBwczevTokrKgoCDuvvtu1q5dy7lz5wBYu3Yt58+fJzU1tczrU1NT2bFjBzk5Od4LXkREGhUlbHWUng4xMdCkifOYnu67tgsKCsjMzGTYsGGEhYVVWic5ORmA9evXl5QVFhZSUFBwwcNJ/OvfuHHj2LRpEwcOHABg8+bN7Nmzh7Fjx1Zav6bxe8vOnTuJjY0lIiKiTHm3bt04f/48e/fuLakXGhpK586dK9QD2LVrl9diEhGRxsXvEjZjzI3GGFvJ4Z3Z816Ung4TJkBuLljrPE6Y4Luk7ciRI5w5c4aYmJgq6xSfy83NLSkbMmQIwcHBFzwyMzPr+R04kpKSuOKKK1i6dCkAixcvpl+/fhUSn2JXXnlljeKv66KIYkePHiUqKqpCeXR0dMn54sfIyMgKvYLl64mIiNSWP99L9EfAB6Wee6/LxEumT4e8vLJleXlOeUpK/bdfmx6w0sN/L774YoVJ/pXp2rXrRcVVW8YYUlNTWbJkCQ899BAZGRnMnTu3yvqrV68uGYasTvv27b0Sn7W20qHZ8j//mtYTERGpLX9O2HZbaze7HUR19u+vXbm3tWrVivDw8Gp7korPdejQoaSsc+fONUoimjZtWtcQa2zcuHHMmTOHxx9/nNOnT5eZL1betddeW6P4g4K88/GOjo5mfyW/1GPHjpWcL348duxYhcStfD0REZHa8rsh0YakY8falXtbUFAQgwYNYt26dZw9e7bSOqtWrQKcYcdi/jYkCtClSxduuOEG5s2bx8iRI4mMjKyyrq+HRLt160ZOTg555bpTd+3aRUhISMnQbbdu3Th37hyffvpphXrgJJoiIhdUeN7tCMQP+XMPW7oxphVwHFgLTLPW+qjvqmbmzHHmrJX+Ho+IcMp9ZerUqQwdOpRHHnmEZ555psy5nJwcnnzySXr06EFiYmJJub8NiRZ7+OGHWbx4MRMnTqy2nq+HRJOTk5k5cybLly/nnnvuAZwFHxkZGQwfPpzQ0FAAbr75ZkJCQkhPT2fmzJklr1+6dCndu3cnNjbWK/GISAA7+QlsGAp9/wCXDXU7GvEj/piwnQB+CWQCJ4FewCPA+8aYXtbaw+VfYIyZAEwA6Oir7i3+N09t+nRnGLRjRydZ88X8tWJDhgxh1qxZzJgxg3379jFu3DiioqLYtm0b8+bNo6ioiIyMjDKvudhE7OWXXwYo2SJkzZo1tG7dmtatW5fpwYuJiSEmJoZNmzbV6vqjRo1i1KhRF6wXFxdXq+tWpybvqWfPnowePZrJkyeTn59PbGwsL7zwAjk5OaSXWmHSpk0bpkyZQlpaGs2aNaN3795kZGSwYcMGXn31Va/FLCIB6vxxyBoJhafh0ivdjkb8TU02a3P7AHrjLDp44kJ1G9vGucXWrFljhw8fbiMjIy1gARsfH28///xzr7VRfN3yR1JSUpl6rVq1sqNHj77g9Tp16mRTUlKqrZOUlFRvG+daW/P3lJeXZ6dMmWLbtm1rQ0NDbUJCgt24cWOF6xUUFNjZs2fbjh072pCQEBsXF2eXL1/u1ZgD9TMs0qgV5lu7YYS1fw6y9lCm29GID1HDjXONbSAr2Iwxu4DPrbUjqqsXHx9vs7Ozqzy/e/durrnmGm+H53dSU1N55ZVXWL9+PX379vVZu3v27KFr165s2bKFhIQEn7XbmDSWz7BIo7J1Cvz7WUiYD53vdzsa8SFjzFZrbfyF6jWkRQcGp/dDamDhwoX06dOHW2+9ld27d/us3eKNfJWsiYjU0N4FTrLWdbKSNamSP85hq8AYEw90AZa5HUtDERISUus5ZN4wfvx4xo8f7/N2RUQapMNZkP0gtBsBvZ52OxrxY36XsBlj0oEcYBvOCtFewP8BB4HnXQxNRETEe07lwNt3wKXfgP5/hSZ+95UsfsQfPx3/BL4DTAIigP8CfwNmWmu/cjMwERERr8j/GjKTwRbCoNUQUvXekyLghwmbtTYNSHM7DhERkXpRVAjvpcDJ3XDTG9D8KrcjkgbA7xI2ERGRgPbxdDi4GuJ/o81xpcYa0ipRERGRhi1nCex6Ejp/H6560O1opAFRwiYiIuILX22GLfdD25sg/jkwxu2IpAFRwiYiIlLfTn8OWbdDxBUwYDk0CXY7ImlgNIdNRESkPhWchqzboPAMDNkIoS3djkgaIPWwBYg333yTW265hZYtWxIWFkbXrl2ZNm0ax48f98r1Dxw4wKRJk0hMTCQiIgJjDPv27avTNWNiYjDGMGbMmErP33jjjRhjGDBgQJ3aqcqmTZswxlQ4IiMrLq8/duwY999/P61ateKSSy5h6NCh7Nixo0K9s2fPMnXqVNq1a0d4eDiJiYlkZWXVS/wi0gDYInj/Hjj+kbPXWgvdVk4ujhK2ADB37lxGjBhBWFgYCxYsYO3atTzwwAMsWrSIhIQEDh48WOc29u7dy7Jly4iKimLgwIFeiNrRrFkzVq5cyddff12mPDc3l6ysLJo1a+a1tqry3HPP8f7775ccb731Vpnz1lqSk5N54403eP7551mxYgX5+fncdNNNHDhwoEzd++67j/nz5zNr1ixee+012rVrx4gRI9i+fXu9vw8R8UM7ZsHnK6Dn09D+FrejkYasJneIb0jH9ddfb6uza9euas83NBs2bLDGGDt58uQK5z777DMbFRVlhw0bVud2CgsLS/48f/58C9icnJw6XbNTp0521KhRtnnz5nbRokVlzs2ePdvGxsbaQYMG2f79+9epnaps3LjRAnbdunXV1lu5cqUF7IYNG0rKjh8/bqOiouykSZNKyrZv324Bu3DhwpKy/Px826VLFzty5EivxR1on2GRgLUvw9p0rH3/XmuLityORvwUkG1rkN+oh62Be+qpp4iOjiYtreJew7GxsUybNo1169axdevWOrXTpEn9fFTCw8O54447WLJkSZnyJUuWMHbsWIwfrKJatWoV7du356abbiopa9GiBSNHjuTVV18tUy84OJjRo0eXlAUFBXH33Xezdu1azp0759O4RcRFR7Jh8z3QegD0eUErQqXOlLDVUfqOdGKejaHJ402IeTaG9B3pPmu7oKCAzMxMhg0bRlhYWKV1kpOTAVi/fn1JWWFhIQUFBRc8nMS//o0bN45NmzaVDC9u3ryZPXv2MHbs2Err1zT+mkpJSaFp06a0bNmSMWPGsH///jLnd+7cSffu3Su8rlu3buzfv59Tp06V1IuNjSUiIqJCvfPnz7N3794axyQiDVjeF84ig7C2MHAFNA11OyIJAErY6iB9RzoTVk8g90QuFkvuiVwmrJ7gs6TtyJEjnDlzhpiYmCrrFJ/Lzc0tKRsyZAjBwcEXPDIzM+v5HTiSkpK44oorWLp0KQCLFy+mX79+dO7cudL6V155ZY3iv9CiiBYtWvCTn/yEBQsWsGHDBh599FHeeustEhMTOXz4cEm9o0ePEhUVVeH10dHRgLMgoSb1jh49euEfhog0bAVnnO078k/AoFUQ1sbtiCRAaFuPOpi+fjp5+XllyvLy85i+fjopcSn13n5tesBKD2m++OKLFSb5V6Zr164XFVdtGWNITU1lyZIlPPTQQ2RkZDB37twq669evbpGw4vt27ev9nyvXr3o1atXyfOkpCQGDRpEQkICzz33HE888QTg/JwrG5ot//OvaT0RCVDWOhvjHs2GQa9A1HVuRyQBRAlbHew/sb9W5d7WqlUrwsPDq+1JKj7XoUOHkrLOnTvXKIlo2rRpXUOssXHjxjFnzhwef/xxTp8+XWYeWHnXXnttjeIPCqr9x7t379506dKFDz74oKQsOjq60t6x4p614l616OjoCsOppesV97SJSIDaNQ9y/ww95sLlt7kdjQQYDYnWQccWHWtV7m1BQUEMGjSIdevWcfbs2UrrrFq1CnB6j4r525AoQJcuXbjhhhuYN28eI0eOrHQvtGLeGhKtSvmesm7durFz584K9Xbt2kXHjh259NJLS+rl5OSQl5dXoV5ISEiVQ7wiEgA+XwkfPQKdxsC109yORgKQetjqYM6QOUxYPaHMsGhEcARzhszxWQxTp05l6NChPPLIIzzzzDNlzuXk5PDkk0/So0cPEhMTS8r9bUi02MMPP8zixYuZOHFitfW8NSRamezsbPbs2cNdd91VUpacnMyiRYvIzMwsSXxPnjzJ6tWry2z6m5yczMyZM1m+fDn33HMP4CwMycjIYPjw4YSGauKxSEA69hG8nwotE+CGBVoRKvVCCVsdFM9Tm75+OvtP7Kdji47MGTLHJ/PXig0ZMoRZs2YxY8YM9u3bx7hx44iKimLbtm3MmzePoqIiMjIyyrzmYhOxl19+GaBki5A1a9bQunVrWrduXaYHLyYmhpiYGDZt2lSr648aNYpRo0ZdsF5cXFytrluVlJQUYmNj6d27N5GRkXz44YekpaXRoUMHJk2aVFIvOTmZxMREUlNTefrpp4mKiiItLQ1rLQ8//HBJvZ49ezJ69GgmT55Mfn4+sbGxvPDCC+Tk5JCe7rvVwyLiQ2cPQ2YyBEfCoJUQFO52RBKoarJZW0M6GtvGucXWrFljhw8fbiMjIy1gARsfH28///xzr7VRfN3yR1JSUpl6rVq1sqNHj77g9Tp16mRTUlKqrZOUlFRvG+fOnTvXxsXF2ebNm9ugoCB7+eWX2/Hjx9svvviiQt0jR47Ye++910ZFRdnw8HA7ePBgu3379gr18vLy7JQpU2zbtm1taGioTUhIsBs3bvRq3IH6GRZpcArOWvtmf2v/Gm7tkWy3o5EGihpunGtsgK1gi4+Pt9nZ2VWe3717N9dcE/j3cktNTeWVV15h/fr19O3b12ft7tmzh65du7JlyxYSEhJ81m5j0lg+wyJ+zVrYch98tggGLIOOd7odkTRQxpit1tr4C9XTooMAtXDhQvr06cOtt97K7t27fdZu8Ua+StZEJKD961dOstZ9ppI18QnNYQtQISEhtZ5D5g3jx49n/PjxPm9XRMRnDr4O26fCFd+GuBluRyONhHrYREREaurELnjvOxDZAxL/CEZfo+Ib+qSJiIjUxLkjkDkSmkZA0ioIusTtiKQR0ZCoiIjIhRTlw9vfhryDMDQTIi53OyJpZJSwiYiIVMdayJ4EhzdB4hJodYPbEUkjpCFRERGR6nzyO9j7onPLqdhUt6ORRkoJm4iISFX++xZs/TF0GAk9fHfbQZHylLCJiIhU5uQn8Pad0Pwa6JeuFaHiKn36REREyjt/HLJGQpMgSFoNwc3cjkgaOSVsAeLNN9/klltuoWXLloSFhdG1a1emTZvG8ePHvXL9AwcOMGnSJBITE4mIiMAYw759++p0zZiYGIwxjBkzptLzN954I8YYBgwYUKd2qlKb93T27FmmTp1Ku3btCA8PJzExkaysrAr1ioqKSEtLIyYmhrCwMHr06MGKFSsqveb8+fO5+uqrCQ0NpWvXrvz+97/35tsTkYtVVADvjIZTn8HAv8GlMW5HJKKELRDMnTuXESNGEBYWxoIFC1i7di0PPPAAixYtIiEhgYMHD9a5jb1797Js2TKioqIYOHCgF6J2NGvWjJUrV/L111+XKc/NzSUrK4tmzervf7W1eU/33Xcf8+fPZ9asWbz22mu0a9eOESNGsH379jL1Hn30UR577DEmTpzImjVr6Nu3L3feeSevv/56mXrz58/ngQce4I477uA+63lvAAAgAElEQVSNN97gzjvv5MEHH+SFF17w+vsUkVr68Kfw3zehzwvQxnv/3onUSU3uEN+Qjuuvv95WZ9euXdWeb2g2bNhgjTF28uTJFc599tlnNioqyg4bNqzO7RQWFpb8ef78+RawOTk5dbpmp06d7KhRo2zz5s3tokWLypybPXu2jY2NtYMGDbL9+/evUztVqel72r59uwXswoULS8ry8/Ntly5d7MiRI0vKDh06ZENCQuyMGTPKvH7w4ME2Li6uzGtbt25tx40bV6bevffea1u2bGnPnz9fbdyB9hkW8SufvGRtOtZmT3E7EmkkgGxbg/xGPWwN3FNPPUV0dDRpaWkVzsXGxjJt2jTWrVvH1q1b69ROkyb181EJDw/njjvuYMmSJWXKlyxZwtixYzHG1Eu7UPP3tGrVKoKDgxk9enRJWVBQEHfffTdr167l3LlzAKxdu5bz58+Tmlp22X9qaio7duwgJycHgPfff58vv/yyQr2xY8dy5MgR3nnnnbq8LRG5WIcy4YMHod3N0Ospt6MRKUMJWx2l79hBzLPP0uTxx4l59lnSd+zwWdsFBQVkZmYybNgwwsLCKq2TnJwMwPr160vKCgsLKSgouODhJP71b9y4cWzatIkDBw4AsHnzZvbs2cPYsWMrrV/T+L1l586dxMbGEhERUaa8W7dunD9/nr1795bUCw0NpXPnzhXqAezataukHkD37t2rrSciPnQqB965A5p1hv5/dRYbiPgRJWx1kL5jBxNWryb3xAkskHviBBNWr/ZZ0nbkyBHOnDlDTExMlXWKz+Xm5paUDRkyhODg4AsemZmZ9fwOHElJSVxxxRUsXboUgMWLF9OvX78KiU+xK6+8skbx13VRRLGjR48SFRVVoTw6OrrkfPFjZGRkhV7ByuoBFa5Zvp6I+Ej+SeceobYIBq2CkBZuRyRSgf4LUQfT168nLz+/TFlefj7T168nJS6u3tuvTQ9Y6eG/F198scIk/8p07dr1ouKqLWMMqampLFmyhIceeoiMjAzmzp1bZf3Vq1eXDENWp3379l6Jz1pb6dBs+Z9/beoB9TrcKyI1VFQI76bAyX/BTW9C86vcjkikUkrY6mD/iRO1Kve2Vq1aER4eXm1PUvG5Dh06lJR17ty5Rsle06ZN6xpijY0bN445c+bw+OOPc/r06TLzxcq79tpraxR/UJB3Pt7R0dHs37+/QvmxY8dKzhc/Hjt2rELiVlk9cHrS2rVrV1KvuGet+LyI+MBHj8AXr0Gf38Flg92ORqRKGhKtg44tKu82r6rc24KCghg0aBDr1q3j7NmzldZZtWoV4Aw7FvO3IVGALl26cMMNNzBv3jxGjhxJZGRklXV9PSTarVs3cnJyyMvLK1O+a9cuQkJCSoZuu3Xrxrlz5/j0008r1AMn0SyuB/+by1ZVPRGpZ58tht1PwVUPwlU/cDsakWqph60O5gwZwoTVq8sMi0YEBzNnyBCfxTB16lSGDh3KI488wjPPPFPmXE5ODk8++SQ9evQgMTGxpNzfhkSLPfzwwyxevJiJEydWW8/XQ6LJycnMnDmT5cuXc8899wDOgo+MjAyGDx9OaGgoADfffDMhISGkp6czc+bMktcvXbqU7t27ExsbC0BiYiKtWrUiPT2doUOHlqkXHR1N//79vRK3iFTjy/fhH+Oh7WC4/lm3oxG5ICVsdVA8T236+vXsP3GCji1aMGfIEJ/MXys2ZMgQZs2axYwZM9i3bx/jxo0jKiqKbdu2MW/ePIqKisjIyCjzmotNxF5++WWAki1C1qxZQ+vWrWndunWZHryYmBhiYmLYtGlTra4/atQoRo0adcF6cV78+dbkPfXs2ZPRo0czefJk8vPziY2N5YUXXiAnJ4f09PSSa7Vp04YpU6aQlpZGs2bN6N27NxkZGWzYsIFXX321pF5wcDCzZ8/mwQcfpEOHDgwdOpQNGzawcOFCnn/+eUJCQrz2/kSkEqf3w9u3Q8QVMGA5NAl2OyKRC6vJZm0N6WhsG+cWW7NmjR0+fLiNjIy0gAVsfHy8/fzzz73WRvF1yx9JSUll6rVq1cqOHj36gtfr1KmTTUlJqbZOUlJSvW2ca23N31NeXp6dMmWKbdu2rQ0NDbUJCQl248aNFa5XUFBgZ8+ebTt27GhDQkJsXFycXb58eaVt//73v7dXXXWVDQkJsZ07d7a//e1vaxRzoH6GRXwi/5S1r/e0dllza4/r75K4jxpunGusj/ba8pX4+HibnZ1d5fndu3dzzTXX+DAid6SmpvLKK6+wfv16+vbt67N29+zZQ9euXdmyZQsJCQk+a7cxaSyfYRGvs0Xwzp1wYCUk/R3a3+x2RCIYY7Zaa+MvVE+LDgLUwoUL6dOnD7feeiu7d+/2WbvFG/kqWRMRv7PjMfj8b9Drl0rWpMHRHLYAFRISUus5ZN4wfvx4xo8f7/N2RUSqlZsB/5wNV94HXX/sdjQitaYeNhERCWxHsmHzd6H1QIj/HWjTammAlLCJiEjgyvsCsm6DsLYwcAU01SpsaZg0JCoiIoGp4Axk3e7cK3T4exDW2u2IRC6aEjYREQk81sKW++BoNgxaCZG+2x9TpD4oYRMRkcCzcy7k/gV6pMHlyW5HI1JnmsMmIiKB5fNX4OOfQ0wqXPszt6MR8QolbCIiEjiObYf3UqHlDXDDfK0IlYChhE1ERALDmUOQmQyh0c68taZhbkck4jVK2ALEm2++yS233ELLli0JCwuja9euTJs2jePHj3vl+gcOHGDSpEkkJiYSERGBMYZ9+/bV6ZoxMTEYYxgzZkyl52+88UaMMQwYMKBO7VRl06ZNGGMqHJGRkRXqHjt2jPvvv59WrVpxySWXMHToUHbs2FGh3tmzZ5k6dSrt2rUjPDycxMREsrKyKtQrKioiLS2NmJgYwsLC6NGjBytWrKiX9ynSKBSeg7dHwbmvYNCrEH6Z2xGJeJUStgAwd+5cRowYQVhYGAsWLGDt2rU88MADLFq0iISEBA4ePFjnNvbu3cuyZcuIiopi4MCBXoja0axZM1auXMnXX39dpjw3N5esrCyaNWvmtbaq8txzz/H++++XHG+99VaZ89ZakpOTeeONN3j++edZsWIF+fn53HTTTRw4cKBM3fvuu4/58+cza9YsXnvtNdq1a8eIESPYvn17mXqPPvoojz32GBMnTmTNmjX07duXO++8k9dff73e369IwLEW/vEAfPUeJP4Jonu7HZGI99XkDvEN6bj++uttdXbt2lXt+YZmw4YN1hhjJ0+eXOHcZ599ZqOiouywYcPq3E5hYWHJn+fPn28Bm5OTU6drdurUyY4aNco2b97cLlq0qMy52bNn29jYWDto0CDbv3//OrVTlY0bN1rArlu3rtp6K1eutIDdsGFDSdnx48dtVFSUnTRpUknZ9u3bLWAXLlxYUpafn2+7dOliR44cWVJ26NAhGxISYmfMmFGmncGDB9u4uLgLxh1on2GROtv1tLXpWPvxY25HIlJrQLatQX6jHrYG7qmnniI6Opq0tLQK52JjY5k2bRrr1q1j69atdWqnSZP6+aiEh4dzxx13sGTJkjLlS5YsYezYsRg/mDC8atUq2rdvz0033VRS1qJFC0aOHMmrr75apl5wcDCjR48uKQsKCuLuu+9m7dq1nDt3DoC1a9dy/vx5UlNTy7STmprKjh07yMnJqed3JBJADv4dPnwYOt4F3We4HY1IvVHCVkc70nfwbMyzPN7kcZ6NeZYd6RXnNdWXgoICMjMzGTZsGGFhlU+uTU529h9av359SVlhYSEFBQUXPJzEv/6NGzeOTZs2lQwvbt68mT179jB27NhK69c0/ppKSUmhadOmtGzZkjFjxrB///4y53fu3En37t0rvK5bt27s37+fU6dOldSLjY0lIiKiQr3z58+zd+/eknqhoaF07ty5Qj2AXbt21Th2kUbt+E549zsQ1Qv6LtKKUAloStjqYEf6DlZPWM2J3BNg4UTuCVZPWO2zpO3IkSOcOXOGmJiYKusUn8vNzS0pGzJkCMHBwRc8MjMz6/kdOJKSkrjiiitYunQpAIsXL6Zfv34VEppiV155ZY3iv9CiiBYtWvCTn/yEBQsWsGHDBh599FHeeustEhMTOXz4cEm9o0ePEhUVVeH10dHRgLMgoSb1jh49WvIYGRlZofewfD0RqcbZryArGYIugaRXISjiwq8RacB0p4M6WD99Pfl5+WXK8vPyWT99PXEp9X8blNr0gJUe0nzxxRcrTPKvTNeuXS8qrtoyxpCamsqSJUt46KGHyMjIYO7cuVXWX716dcnwYnXat29f7flevXrRq1evkudJSUkMGjSIhIQEnnvuOZ544gnA+TlXNjRb/ufv7XoiUoXC8/DOtyHvIAzNhIjL3Y5IpN4pYauDE/tP1Krc21q1akV4eHi1PUnF5zp06FBS1rlz5xolB02bNq1riDU2btw45syZw+OPP87p06fLzAMr79prr61R/EFBtf949+7dmy5duvDBBx+UlEVHR1fa61Xcs1bcqxYdHV1hOLV0veIetOjoaI4dO1YhcStfT0QqYS1snQSHMyFxKbS6we2IRHxCQ6J10KJji1qVe1tQUBCDBg1i3bp1nD17ttI6q1atApzeo2L+NiQK0KVLF2644QbmzZvHyJEjK90LrZi3hkSrUj6R6tatGzt37qxQb9euXXTs2JFLL720pF5OTg55eXkV6oWEhJQM8Xbr1o1z587x6aefVqgHTkIqIlXY8xvY+xJc+38Qm+J2NCI+ox62OhgyZwirJ6wuMywaHBHMkDlDfBbD1KlTGTp0KI888gjPPPNMmXM5OTk8+eST9OjRg8TExJJyfxsSLfbwww+zePFiJk6cWG09bw2JViY7O5s9e/Zw1113lZQlJyezaNEiMjMzSxLfkydPsnr16jKb/iYnJzNz5kyWL1/OPffcAzgLQzIyMhg+fDihoaEA3HzzzYSEhJCens7MmTNLXr906VK6d+9ObGxsreMWaRT+8yZsmwyX3wY9nnA7GhGfUsJWB8Xz1NZPX8+J/Sdo0bEFQ+YM8cn8tWJDhgxh1qxZzJgxg3379jFu3DiioqLYtm0b8+bNo6ioiIyMjDKvudhE7OWXXwYo2SJkzZo1tG7dmtatW5fpwYuJiSEmJoZNmzbV6vqjRo1i1KhRF6wXF+edn29KSgqxsbH07t2byMhIPvzwQ9LS0ujQoQOTJk0qqZecnExiYiKpqak8/fTTREVFkZaWhrWWhx9+uKRez549GT16NJMnTyY/P5/Y2FheeOEFcnJySE9PL6nXpk0bpkyZQlpaGs2aNaN3795kZGSwYcOGMtuEiEgpJ/8N79wFLbo7Q6FGA0TSyNRks7aGdDS2jXOLrVmzxg4fPtxGRkZawAI2Pj7efv75515ro/i65Y+kpKQy9Vq1amVHjx59wet16tTJpqSkVFsnKSmp3jbOnTt3ro2Li7PNmze3QUFB9vLLL7fjx4+3X3zxRYW6R44csffee6+Nioqy4eHhdvDgwXb79u0V6uXl5dkpU6bYtm3b2tDQUJuQkGA3btxYoV5BQYGdPXu27dixow0JCbFxcXF2+fLlNYo7UD/DIlU6d9TaVV2sfbm1taf2uR2NiFdRw41zjQ2wlWnx8fE2Ozu7yvO7d+/mmmuu8WFE7khNTeWVV15h/fr19O3b12ft7tmzh65du7JlyxYSEhJ81m5j0lg+wyIAFBXApm/C4U0weAO0qZ97C4u4xRiz1Vobf6F66lMOUAsXLqRPnz7ceuut7N6922ftFm/kq2RNRLxi20Pw33XQ5/dK1qRR0xy2ABUSElLrOWTeMH78eMaPH+/zdkUkAH3yIux5Hq5+CK78ntvRiLhKPWwiIuJ/Dm2C7InQ7mbo+ZTb0Yi4TgmbiIj4l1Ofwdt3QLOroP9foYnvNvEW8VdK2ERExH/kn4TMkc6fk1ZBiG82Ihfxd41yDput4l6OIv4u0FZ1i5RRVAjvjoGTe2Dwm9Css9sRifiNRtfDFhwczJkzZ9wOQ+SinDlzhuDgYLfDEKkfH02DL/4O8c9D25vcjkbErzS6hK1NmzYcPHiQvLw89VZIg2GtJS8vj4MHD9KmTRu3wxHxvs/+CLt/AVf9EK76vtvRiPidRjck2rx5cwC++OIL8vPzL1BbxH8EBwfTtm3bks+wSMD48l34xwPQdghc/yu3oxHxS40uYQMnadOXnoiIHzidC2+PgoiOMGAZNNGQv0hlGt2QqIiI+In8U5B5GxSeg6TVEBrtdkQifqtR9rCJiIjLbBG8Pw5O7ICk16HF1W5HJOLXlLCJiIjvfTwDDrwCvZ+F9iPcjkbE72lIVEREfGvfn2HnHLjyfuj6I7ejEWkQlLCJiIjvfPUP2Pw9aD0Q4n8L2sRcpEaUsImIiG/kHYS3b4fwdjBwBTQNcTsikQZDCZuIiNS/gjzIuh3yv3ZWhIa1djsikQZFiw5ERKR+WesMgx7dCoNehcjubkck0uAoYRMRkfq1cw7sz4CeT8LlI92ORqRB0pCoiIjUn/0r4ONHIWYsXDPV7WhEGiwlbCIiUj+OfuhsjtuyL9zwklaEitSBEjYREfG+M/+FrNuc200NegWahrkdkUiDpjlsIiLiXYVnIev/wbkjMOwdCL/M7YhEGjwlbCIi4j3Wwj8egCObYcDLEN3L7YhEAoKGREVExHt2/wJyFkPcLOh4h9vRiAQMJWwiIuIdB1bD9p9Bx9HQ/eduRyMSUJSwiYhI3R3/J7w3BqJ7Q9+FWhEq4mVK2EREpG7OfgmZIyG4mXMng6AItyMSCThadCAiIhev8Dy882048x8YmgURHdyOSCQgKWETEZGLYy1k/xAOZ0G/dGiV4HZEIgFLQ6IiInJx9jwPny6Abo9AzBi3oxEJaErYRESk9r5YC9umwOW3w3Wz3Y5GJOD5fcJmjHnDGGONMU+4HYuIiAAn/gXvjoYW3SFxCRi//yoRafD8+m+ZMeY7QA+34xAREY9zR50VoU1CIGkVBF/qdkQijYLfJmzGmEjgV8BDbsciIiJAUT68cxfk5To3dL+kk9sRiTQafpuwAU8BO621f3E7EBERAbY9BIfWQ8JL0Lq/29GINCp+ua2HMWYAMA4Nh4qI+IdPfg97fgNX/wS+8V23oxFpdPyuh80YEwy8CPzCWvtvt+MREWn0Dm2E7EnQ/pvQ80m3oxFplPwuYQN+BoQDc2r6AmPMBGNMtjEm+8svv6y/yEREGpuv98Lbd0Czq6Dfn6FJU7cjEmmU/CphM8Z0BKYDjwKhxphIz+IDSj2v8K+FtfYla228tTa+devWvgxZRCRwnT8BmcmAgaTVENLC7YhEGi2/StiAbwBhwFLgWKkD4KeeP8e5E5qISCNSVAjvfge+/gQGvgzNrnQ7IpFGzd8WHWwHbqqkfCNOEvcHYK9PIxIRaYy2/wz+swb6/B7aVvbPsoj4kl/1sFlrj1trN5U/PKdzPc9PuRmjiNRcejrExECTJs5jerrbEUmNfLoI/vVL6DIRrnrA7WhEBP/rYRORAJGeDhMmQF6e8zw313kOkJLiXlxyAYffgQ8egMuGQu9fuR2NiHgYa63bMXhVfHy8zc7OdjsMkUYvJsZJ0srr1An27fN1NFIjp/bB2gQIiYQRWyAkyu2IRAKeMWartTb+QvX8akhURALH/v21KxeX5Z+CrGQoOu9ZEapkTcSfKGETkXrRsWPtysVFtgjeT4UTO2HAMmje1e2IRKQcJWwiUi/mzIGIiLJlERFOufiZjx+FA686c9baDXc7GhGphBI2EakXKSnw0kvOnDVjnMeXXtKCA7+z78+wcy5cOR66THI7GhGpglaJiki9SUlRgubXvtoCm78HbZIg/jdOZi0ifkk9bCIijVHeAci6HcLbw4CXoWmI2xGJSDXUwyYi0tgU5EHmbVBwCga/BWGt3I5IRC5ACZuISGNiLWz+Lhz7EJJWQWQ3tyMSkRpQwiYi0pj8czbsXw49n4IO33I7GhGpIc1hExFpLPa/DDtmQuw4uOanbkcjIrWghE1EpDE4ug3eHwetEiHhJa0IFWlglLCJiAS6M/+BrNsgtBUMfAWahrodkYjUkuawiYgEssKzkPX/4NxRGP4uhLd1OyIRuQhK2EREApW1sGU8HNkCA1dAVE+3IxKRi6QhURGRQLX7Kdi3FK6bDVeMcjsaEakDJWwiIoHowCrY/n/Q6W7oNt3taESkjpSwiYgEmuM74L0UiL4eblioFaEiAUAJm4hIIDn7JWSOhODmMOhVCAp3OyIR8QItOhARCRSF5+HtUXD2EAzNgoj2bkckIl6ihE1EJBBYCx/8AL58B/r9BVr2cTsiEfEiDYmKiASCf/8aPlsI3X4OMXe7HY2IeJkSNhGRhu6LN+DDnzhbd1z3uNvRiEg9UMImEmDS0yEmBpo0cR7T092OSOrViX/Bu6OhRRwkLgajf9ZFApHmsIkEkPR0mDAB8vKc57m5znOAlBT34pJ6cu6osyK0aRgkrYKgS9yOSETqif4rJhJApk//X7JWLC/PKZcAU5QP79wJefudG7pf0tHtiESkHqmHTSSA7N9fu3JpwLZOhkMboO+foHU/t6MRkXqmHjaRANKxik6WqsqlgdrzO/jkd3DNw/CNcW5HIyI+oIRNJIDMmQMREWXLIiKccgkQ/90AW38E7b8FPea6HY2I+IgSNpEAkpICL70EnTo5t4/s1Ml5rgUHAeLrvfDOt6H51dA/HZo0dTsiEfERzWETCTApKUrQAtL5486KUNPEWREa3NztiETEh5SwiYj4u6ICePdup4dt8Ftw6TfcjkhEfEwJm4iIv/vwYfjPWkh4CdomuR2NiLhAc9hERPzZp3+Af/8Kuv4YOo93OxoRcYkSNhERf3X4bfjgB3DZcOj1C7ejEREXKWETEfFHp/bB26PgklgYkAFNNINFpDFTwiYi4m/yv3ZWhBYVQNJqCIl0OyIRcZn+yyYi4k+KCuG9FDi5G25cA827uB2RiPgBJWwiIv7k45/DwdVw/fPQbpjb0YiIn9CQqIiIv8hZCrvmQefvQ5cfuh2NiPgRJWwiIv7gq82w5X5ocyPEP+fcW0xExEMJm4iI205/Dlm3Q0QHGPgyNAl2OyIR8TOawyYi4qaC05B1GxTkwZANENrS7YhExA8pYRMRcYstgvfvgeMfwaDV0OJatyMSET+lhE1ExC07ZsHnK6DXL6HDN92ORkT8mOawiYi4IXcZ/PNx+Ma9cPUUt6MRET+nhE1ExNeOboXN34XW/aHPC1oRKiIXpIRNRMSXzvwHMm+D0NYw8G/QNNTtiESkAdAcNhERXyk442zfkX8chr0LYW3cjkhEGgglbCIivmCtszHukX/AwFcgqofbEYlIA6IhURERX9g1D3L/DD3mwhW3ux2NiDQwSthEROrbgVfho0eg0xi4dprb0YhIA+S1hM0YM9hb1xIRCRjHPob3UqBlAtywQCtCReSieLOH7XljTJmdH40xd3nx+iIiDcvZw5CVDMGRMGglBIW7HZGINFDeXHQwDPi7MSYcOAD8CmgKLPNiGyIiDUPhOXj7DidpG/Y2hLdzOyIRacC8lrBZa78wxnwPeBc4Dky21ipZE5HGx1r44Afw5TvQPwOir3c7IhFp4Lw5h+1pYA2QBhzC6V0TEWl8/vUr+GwRdJ8JnTQzRETqzptDos2AOGvtl8aY3wFrjDGXWGsXeLENERH/dvB12D4Vrvg2xM1wOxoRCRBe62Gz1n7fWvul589HcOa0fddb1xcR8XsndsF734HIHpD4RzDaOUlEvMOr/5oYY1p4Fh1grT0BDPfm9UVE/Na5I5CZDE3DYdCrEHSJ2xGJSADx5hy2nwOfAMeNMTuNMX8Evuet64uI+K2ifHjnTsg74GzfcckVbkckIgHGmz1s44COwL+AiUA3YJAXry8i4p+2/hgObYQb5kOrvm5HIyIByJsJ2wlr7VnAWms3At8Eznrx+iIi/mfPb+GTF+Dan0HsWLejEZEA5c2ELc8Y0xQ4Zoxp41mAcJ0Xry8iDVR6OsTEQJMmzmN6utsRecl/33J61zqMdG7qLiJST7y5rccUnK09XgRWGWO2e/HaItJApafDhAmQl+c8z811ngOkpLgXV52d/ATevhOaXwP90rUiVETqlTe39dhmrT1urf0rMBs4CNzmreuLSMM0ffr/krVieXlOeYN1/jhkjYQmQZC0CoKbuR2RiAQ4b/awYYxpC5y01v4d+Ls3ry0iDdP+/bUr93tFBfDu3XDqMxj8Flwa63ZEItIIeHNbj8U4K0QPGWP+aYz5ozFmoreuLyINU8eOtSv3ex9Ohf+shT4vQBsthBcR3/DmpIt+wGVAJDAa2ABc5cXri0gDNGcORESULYuIcMobnL0L4N/PQtfJcOV9bkcjIo2INxO23UChtbbIWrvTWrvYWvtjL15fRBqglBR46SXo1AmMcR5feqkBLjg4nAXZD0K7EdDrabejEZFGxlhrvXMhY3oAPwIetdZ+4ZWLXoT4+HibnZ3tVvMiEohO5cDaPhDaGoZvhpAWbkckIgHCGLPVWht/oXre7GF7DPgWsMMYs9cY81djzE+9eH0REd/LPwmZI8EWwaBVStZExBXeXCV6HXC5tTbfGBMD9AEumDGKiPitokJ4NwVO/gtuWgvNNS1XRNzhzYRtJ2AArLX7gH3Aci9eX0TEtz6eDl+8BvG/hcuGuB2NiDRi3hwSPQ0sM8bov6Ai0vDlLIFdT8JVP4AuD7odjYg0ct5M2Pbg9NhlGWP+a4x5zRjzmBevLyLiG1++D1vuh7aD4fpfux2NiIj3hkSttTOL/2yM6QBc7zlERBqO0/vh7dsh4goYsByaBLsdkYiId29NVcxaexDnXqKr6uP6IiL1ouA0ZN0GhWdhyCYIjXY7IhERoJ4SNhGRBscWwfv3wPGPIek1aHGN2xGJiGFyE4UAACAASURBVJSoc8JmjGkKXI2zrUcP4Dpr7Tfrel0REZ/a8Th8vgJ6PwPtb3E7GhGRMmqVsBljWuFJyko9XgOE4GzpcQ7Y4eUYRUTqV24G/HMWfON7zn1CRUT8TI0SNmPMvcAsoH2p4tPAP4DfAts9x25rbaG3gxQRqTdHsmHzd6H1AOjzO+eGpyIifqamPWzzgGPANODfwE+AvkAmMFdJmog0SHlfOIsMwtrCwBXQNNTtiEREKlXThK018B1r7QbP81XGmB8CacBtxpjvWmv/WS8RiojUh4IzkHW7c6/Q4e9BWBu3IxIRqVJNN879BrCldIG19rc4c9iOAdnGmOmeBQgiIv7NWthyHxzNhn7pEBnndkQiItWqUcJmrd1nrT1dRfkw4EfAVGCzMaa7l2MUEfGuXWmQ+xfoMRcuT3Y7GhGRC/LKramstS8B3YEvgWxvXFNEpF58/gp8NB1iUuDan7kdjYhIjXjz1lQHgG8aY+7x1jVFRLzq2HZ4LxVaJsANC7QiVEQaDG/e/B0Aa+2fvH1NEZE6O3MIMpP5/+3dd3xb9dXH8c9x4kAUiMImDFlsaDEFEjYhgBhdDpRd3Ba6VErLA91DHaRFpZQ+1E9pKbiTtmqBlhU/fWgB0SZhBEhYYs/IrJAQiBJwhhP/nj+undiOnUiypHslfd+vl1+K7r26PpiQHO7vd85hzBZw1C0walO/IxIRyZtGU4lI7VuzEmafAivfhOPvhrET/Y5IRKQgSthEpLY5Bw+eB2/eC0feAFse6HdEIiIFK/mSqIhIoDx9Bbz4B2i+GCKn+x3NWqlMimhblIbpDUTboqQyKb9DEpEAK9kTNjM7EIj0vu10zj1UqnuLiBTl1X/Aw1/zErV9v+t3NGulMiniHXG6ursAyOayxDviALQ2t/oZmogElDnnRnYDs0lACnBAJ94Q+Kbe02c75+aN6BsUaPLkyW7uXHUWEal7S56A2w+DzfeA42fD6JDfEa0VbYuSzWXXO94UbmL+RfMrH5CI+MbM5jnnJm/sulI8YWsHvuicu3NQAMcDvwa0YUREKmvFmzBrGoweB1NvDVSyBtCZ6yzouIhIKfawbTY4WQNwzt0BjCvB/UVE8rdmFdx9GnS96rXvCO3kd0TriYQjBR0XESlFwrbQzD7Zf46omY0ys88Ab5bg/iJSpVIpiEahocF7TZV7X71zMO8CWDgTDvktbH1Imb9hcZKxJKHGgU/9Qo0hkrGkTxGJSNCVImE7B/go8JaZPWVmTwFvAWf2nhOROpRKQTwO2ayXR2Wz3vuyJm3P/gKeb4f3fAt2Ce7m/dbmVtpb2mkKN2EYTeEm2lvaVXAgIsMacdHB2huZbc3AKtGinq6Z2YnAN4D3AFvgzSe9F7jYOffkxj6vogORYIhGvSRtsKYmmD+/DN/w9dvhPx+AHVtgyk1g6lokIsFXyaIDAHoTtAFJmplt65xbWOCttgTmAVfhJWsR4JvAHDNrds4N8VeAiARN5zD754c7PiJLn4G7z4Dwe+GwPylZE5GaU+5JB3NZ99QtL865vwJ/7X/MzB4AngZOA/67ZNGJSNlEIkM/YYuUel/9qre9GaENjXDUDGjcvMTfQETEfyNO2Mxs2gZOl2q68uLe1+4S3U9EyiyZ9PasdXWtOxYKecdLpmc13H0mvPsSHJuGzaIlvLmISHCU4gnbzcBMvIa5gxX9v7q9Vaej8Jrw/hhYAFxX7P1EpLJae/fPJxLeMmgk4iVrraXcV//QV2DBHXDIb2DbKSW8sYhIsJQiYXsO+LRz7qXBJ8zs5RHc935gUu+vnweOLWI/nIj4qLW1xAlaf8+3w7M/h72+BLt9ukzfREQkGEqxM/daYOthzl09gvt+HDgUOBtYCtxhZtGhLjSzuJnNNbO5ixYtGsG3FJGq8MZ/4MEvwMT3wwE/8TsaEZGyK1lbj3IyswnAfOA659x5G7pWbT1Eatw7L8I/D4JNt4UT5sCYsN8RiYgULd+2HlVR++6cW4K3LLq737GIiI+6l8LMFsDB1A4layJSNwrew2Zmd23gdA+Qw+uj9lvn3BvFBjboe24H7A2Ue7CNiARVzxq452yv59oxt8Pm9fX/b6lMikQ6QWeuk0g4QjKW1GQEkTpSTNGBAXsCE4GXgDeA7YBdgNd7338Q+JKZTc1nOsGAm5vdDDwEPIa3d21P4EvAatSDTaR+PfoteO0fcNBVsP2xfkdTUalMinhHnK5ur0dKNpcl3hEHUNImUieKWRK9AlgBTHLO7eacO9w5txtwUO/x6cAeeFMKium4NAc4Ga+Y4R/Al/HahuzvnHu2iPuJSLV78Q/w1OWwx/mwx+f9jqbiEunE2mStT1d3F4l0wqeIRKTSiknYLsGb6/lw/4POuXl4ydolzrlXgMuBowq9uXPuMufcJOfcBOdcyDm3l3Puc865+UXEKiIVlkp5c0QbGrzXEQ97X3QPPPA52C4Gk9pKEGH16cwNPc9ruOMiUnuKWRLdk0EzQ/tZxLrCgBeAccUEJSLVKZUaON0gm/XeQ5H92N7NwuxTIBSBI2/wxk/VoUg4Qja3/pyvSLjUc76Kl8pkSKTTdOZyRMJhkrEYrc3NfoclUjOKecI2H/jMMOfivefB6822eJjrRKQGJRIDR1GB9z5RzMpd9zsw8yRYs9KrCN1ky5LEWI2SsSShxtCAY6HGEMlYKed8FS+VyRDv6CCby+GAbC5HvKODVCbjd2giNaOYJ2w/AP5sZo8BNwILgW2BU4F98RrdAhyHN61AROpE5zArdMMdH5brgfs+AbkMTP0HhPcecWzVrK+wIKhVool0mq7ugaOeu7q7SaTTesomUiIFJ2zOub+a2Zt4+9W+DTTiDWWfC5zgnLuz99IvA2tKFaiIBF8k4i2DDnW8II99D165GQ78Gezw/pLEVu1am1sDk6AN1pnLFXRcRApXVONc59wdzrnDgbHA9sBY59wR/ZI1nHMrnHPdw95ERGpOMgmhgSt3hELe8bzN/ws8kfTmg+51YUnjk/KIhIduYDzccREp3IgmHTjnepxzC51zPaUKSESqV2srtLdDUxOYea/t7QUUHLz5AMz5FGwzBSZf5d1EAi8ZixFqHFgQEmpsJBmL+RSRSO0pZg8bZjYR+AowFdgSr7jgP8AVzrkFJYtORKpOa2uRFaFdr8Lsk2HsRJhyI4waU/LYpDz69qmpSlSkfAoe/m5mewKzgS2Ae4AFeMuihwNvA1Occ8+VOM68afi7SBVa3QV3ToWlT8MJ98IE/UUvIvUh3+HvxTxhuwxvZNQh/ZvZmlkTcHvv+VOKuK+I1CPnvGXQt+bBUbcoWRMRGUIxe9iOAb47ePKAcy4LXNx7XkQkP49fAp3Xw/6Xwk7T/I5GRCSQiknYxgDLhjm3rPe8iMjGdd4Ime9B9GOwz9dHdKtUJkW0LUrD9AaibVFSmZHOxBIRCY5ilkQfAS4ws9v6V4eamQHn954XEdmwtx72muNudQgc8usRVYSmMiniHfG1A9KzuSzxDm8mVlB7l9UrjbASKU4xRQfvB/4Xb1bo9cDreEUHpwN7AB9yzt1e4jjzpqIDkSqwfAH862DAwYkPwtjtR3S7aFt0yFmbTeEm5l80f0T3ltLpG2HVfypCqLGR9pYWJW1St/ItOih4SdQ590/gw3jLnwngl8B3gHeAD/uZrIlUu1QKolFoaPBeU7W4qrdmBcz6CKxcDEfNGHGyBtCZG3r21XDHxR8bGmElIhtWVB+23qTtn2YWwmvv8bZzrmsjHxORDUilIB5fNzw9m/XeQ5F9zYLIObg/DovnwJF/gy0PKMltI+HIkE/YIuFCZ2JJOWmElUjxRjrpoMs596qSNZGRSyTWJWt9urq84zXjqcth/p+geTpETivZbZOxJKHGgTOxQo0hkrFCZmJJuWmElUjx8krYzKzHzNbk+bW63EGL1KLOYVbvhjtedV7pgEe+CZEzYN/vlvTWrc2ttLe00xRuwjCawk20t7Sr4CBgNMJKpHj5Lon+ACisOkFEChKJeMugQx2vekseh3vPhi0OgEN/X5YZoa3NrUrQAk4jrESKV3CVaNCpSlSq1eA9bAChUIHD04NoxSKvInTNCnj/gxDaye+IREQCo6RVomZ22MhDEpENaW31krOmJu8BVFNTDSRra1bB3afB8te9sVNK1kREipLvkujdZrYQmAHcDKSdc90b+YyIFKi1tcoTtP6cg7lfgIWz4PAUbH2I3xGJiFStfKtEd8Tbx9YE3AK8aWbXm9lZZja+bNGJSPV65ufwwm/gvd+G6Nl+RyMiUtXySticcwucc79yzr0f2AY4D68I4WpgoZn9y8zOM7MdyhiriFSL1/4FD38ZdjoJ9vuh39GIiFS9YiYdLHPO/dU5dxZe8vYR4CXgu8DLZna/mX2zxHGKSLXIPQ33nAnhfeGwP4ONqN2jVIlUJkO0rY2G6dOJtrWRymT8Dkmkpoy0cW63c+4259x5zrkdgSOB/wDnlCI4EakyK9+CmS3QMAamzoDGzfyOKLBSmRTRtigN0xuItkVJZap3DlnfjNBsLocDsrkc8Y4OJW0iJVSy//U1s3Odc/c5577hnNunVPcVkSrR0w13nwFdWZhyE4xr8juiwEplUsQ74mRzWRyObC5LvCNetUmbZoSKlN8GEzYz28XM0mb2opldYWab9jv3wKDL/6vfuU+VOE4RCbqHvgxvpOHgdtj2SL+jCbREOkFX98A5ZF3dXSTS1TmHTDNCRcpvY0/YrgJuAk4HtgTSZrZ577nGQdf2b13+xdKEJyJV4bmr4dlfwN5fgV3P9TuawOvMDT1vbLjjQacZoSLlt7GEbTvn3C+dc/Occ+cC/4uXtIVZf1RV//elnzsjIsH0xr9h7gUw8QOw/2V+R1MVIuGh540NdzzoNCNUgiCTytAWbWN6w3Taom1kUrW1h3JjCdsm/d845y4FbgDSwOaDrt2rtzfbV4GQmYVKF6aI+CmVgmgUGhq811TfVqtlz8Ps02DzPeCIv0LDKB+jrB7JWJJQ48A/IkONIZKxpE8RjUxrczPtLS00hcMY0BQO097SohmhUjGZVIaOeAe5bA4c5LI5OuIdNZW0bWzSwbNmdrxz7o6+A865n5pZD/DTQdceDUzq/XoXeMvMOoFHgcecc2rGJFKFBs84zWa9943kOCM8zTs4tQPGaPkrX31D6hPpBJ25TiLhCMlYsqqH17c2NytBE9+kE2m6uwYWvnR3dZNOpGlurY3flxtL2M4a6qBz7gozu37QsQeBB/vem9kYYD/WJXEiUoUSiYED6QFWLF/DNs9+FPZ+Do69HTbfzZ/gqkAqkxoyMev7EpGRy3UOXeAy3PFqtMGEzTm3cgPnXt3IZ1cBc3u/RKRKdQ6xD/6yj36DY/a6DSZfDdsdU/mgqkRf+46+itC+9h2AkjWREgpHwt5y6BDHa4VakIvIBkUG7YM/96jf89UP/Td/uO+LsMfn/AmqStRa+w6RoIolYzSGBha+NIYaiSVrp/BlY0uiAJjZS6xfFZo359yuxX5WRPyVTK7bw3b4nvdwzac/x11PHseYQ3/md2iBV6n2HalMhkQ6TWcuRyQcJhmLaT+Z1JW+fWrpRJpcZ45wJEwsGauZ/WuQZ8IG/KGcQYhIcLX2rtxddXmWmz/3EV5ZEmXxPjdwdmu+f3zUr0g4QjaXHfJ4qfSNheqbNNA3FgpQ0iZ1pbm1ecQJWiaVCWzSl++fuC+VNQoRCbTWM96hdYtp8O4qtj1hBruGt/A7pKqQjCUH7GGD0rfv2NBYKCVsIvnraw3SV23a1xoECETSVoknbA744wg+LyJ+cj1w38ch9zgcfRuE9/Y7oqpRifYdGgslUhpBbw2SV8LmnFNxgki9euy78MotMOl/YOIJfkdTdcrdviMSDpMdIjnTWCiRwgS9NUheiZiZ7VzuQEQkgOb/BZ74Eez2WdjzAr+jKYtUJkW0LUrD9AaibVFSmdTGPxQgxYyFSmUyRNvaaJg+nWhbG6lM7XSDFynWcC1AgtIaJN8nZy+a2Q1mdkRZoxGR4HjzAZjzKdj2KJj8C7DaGxHc1yctm8vicGv7pFVT0lboWKi+IoVsLodjXZGCkjapd0FvDWLObbxbh5l9H4gD2+ONmmoDruttjhsokydPdnPnqlevyIh0vQL/PAhGjYUTH4BNt/Y7orKItkWHrOJsCjcx/6L5lQ+oAqJtbUMuoTaFw8y/6CIfIhIJDj+qRM1snnNu8sauy3cP23QzSwKnARfgFSH8xMyuAX7lnFswkmBFJEBWd8Gsk2H1O3DsnRtM1oYbu1QtKtUnLUhUpCDVqFKJVClag5RL3sUEzrnVzrnrnHNHAJOB24CvAfPN7E9mdlC5ghSRCnEO5nwS3noIjvgrTHjvsJfWwnLicP3QStknLWiGK0ZQkYIEVV+7jVw2B25du41Mqr6W8Yuq/nTOPeSc+ySwMzAdOAqYY2b3mdmQA+NFpAo8fgl03gD7XwY7fniDlwZ57FK+m+qTsSShxtCAY6XukxY0xRQpiPhpQ+026klRCZuZjTOzicA2QBo4D/gbcBDw59KFJyIV03kjZL4Hu3wC9vnqxi8P6HJiIZvqW5tbaW9ppynchGE0hZtob2mvqmXdQhVapCDit6C326iUfGeJPgKEe782Z12i179sbDWQ6/0SkWry1sNw3ydg68Pg4GvyqgitxNilYhTa+b/cfdKCqLW5WQmaVI1wJOwthw5xvJ7k+4RtP6AJr0K0FTgOOBjYC5gIjHPOjXHObaVB7yJVZvkCmDUNNtkKptwEozbN62NBXU7UpnqR2hL0dhuVkm/Ctj9eZeghwJXA0cCrzrnnnHNvOOeWlyc8ESmrNStg1kdg5VswdQaM3T7vjwZ1OVGb6kVqS3NrMy3tLYSbwmAQbgrT0t4S2GrOcskrYXPOPeac+xRekcH/AJ/Eqw79o5kdWM4ARaRwqRREo9DQ4L2mhircdA7uj8PiOXD4n2CL/Qv7HgFt6aFN9SK1p7m1mYvmX8T3e77PRfMvqrtkDQosOnDOLXbO/QjYBfhE7+tcM7vbzE41M80cFfFZKgXxOGSzXk6WzXrv10vanroc5v8J9vsh7HxKYd8jwC09tKlepL5kUhnaom1Mb5hOW7StZtt95DXpYIM3MDsAuBA4E3gD+IVz7qcliK0omnQg9SiVgkQCOju9p2pr1qx/TVMTzJ/f++aVDph1EjSdCYf/peCxU/U4IWAkUpkMiXSazlyOSDhMMhZTAilSAn092vq3/WgMNVbVkmlJJx2Y2fZ4FaLjh/l6A/gn0AJcBviWsInUm74nal29LdGGStbAS+YAWJKBe8+GLSfBIb8rakZoUFt6BFFfm5G+ytW+NiOAkjaREdpQj7ZqSdjylVfCBrwG9D2KG/ynuwPewWvn8SRq6yFSUYnEumRtQyIRYMUimDkNGjeHo26B0WOL+p5BbekRRIW2GREJIj9mbOajnnq05ZuwXQi8zbo+a0v6/TrnRrquKiJFW/vkrDkDsTSEc5ALQzoGGe8P1FAILk2ugtmnwooFcNwsCO1Y9PdMxpLEO+IDJh0EoaVHEKnNiFS7wcuOfaOhAN+Ttnrq0ZZvleiVzrk/O+c6nHOzeqtGs865JUrWRPwVieAlay0dMCHnPQOfkPPeN2doaoL2dsdHdz8fFs2GQ35P6rVnibZFaZjeQLQtWnCxQFBbegSR2oxItQvyaKh66tGW7xM2EQmoZBI+Pi+NGzPwD1TGdNP0qTTzL2qGp38OD/0W3vsdUsvWDHg61lfhCRSUcNXjhIBiJGOxAXvYQG1GpLoEedmx7wlfEJdrS23EVaJBoypRqUc2ffrQx4Gezx4KMz8IO54EU/7O1pfHWLx8f7w6ohzeOOCMKjzLSFWiEkT57ktri7YNvezYFOai+RdVItSaVuoq0ZdYV3RQMI2rEimvpnCY7BB7oo7ZcgXccyaEm+GwP5J6/AkWLz8cGNN7xQS84m7ozD1esXjrjWZ3StAUsi8tlowN2TqjFpcdgyzfJdE/lDMIERmZoZbddtykm5sn/gUaNoGpt0LjZiTSadYla33GADEi4aWVDFlEfFRIO4x6WnYMsnwTtpfKGoWIjEjf05u+Zbddw5sxe5e/MX75Aoj9G8Y1ARuqTAyrwrOGBLUFgwRHofvSmlub9XvIZ5V4wuaAP47g8yKShwHLbg9+EZ67Hw79A2xz+NprIsMsnW41tlEFBDUiyC0YJDjqqR1Grci3rUfDCL5GlfsfQkT6ee5X8NwvYZ+vwa7nDDg13GD0//nAtEpGKGUU5BYMEhz11A6jVqith0gtWXAXzL0AdvgQvO/S9U4PXjodXLGopbTqF+QWDBIc2pdWfZSwidSKZc/D3afB+L3giL9Aw9APt4erWNRSWm3QUpfkS/vSqkteS6L5MLMDzezk3q8DS3VfkaBIpSAahYYG7zVV2HCA8lqVg5ktYA0wtQMaxxd8i3ItpWVSGdqibUxvmE5btI1MKjOi+8mGaalLpDaN+AmbmU0CUnjFBZ14vTqbzAzgbOfcvJF+DxG/pVIQj68bsp7Neu8BWv3eq9+zBu45y3vCduwdsFlxbQ/LsZSmp3aVp6UukdpUiiXRduCLzrk7+x80s+OBXwN62iZVL5FYl6z16eryjvuesD3ydXj9n3DwNbDd0UXfphxLaYX0epLS0VKXSO0pxZLoZoOTNQDn3B3AuBLcX8R3nZ2FHa+YF34HT18Be14Au8dJZTJE29pomD6daFsbqUz+y4/lWErTBngRkdIoRcK20Mw+aWZrdzib2Sgz+wzwZgnuL+K7SKSw4xWxcDY8eB5sfzwceAWpTIZ4RwfZXA4HZHM54h0deSdtza3NtLS3EG4Kg3lzAlvaW0b0pGa4p3PaAC8iUphSLImeA1wNtJnZa73HdgAe6D0nUvWSyYF72ABCIe+4L96ZD7NPgXG7wJHXQ8NoEun0gNFUAF3d3STS6bznWJZ6KU0zCEVESmPECZtz7kXgBDPbGuh73tDpnNPTNakZffvUEglvGTQS8ZI1X/avdS+DWdOgZ7VXETpmC2D4sVPDj6MqP22AFxEpjZL1YetN0AYkaWa2rXNuYam+h4ifWlsDUGDgeuDej0HuSTj6Nhi/59pTw42dioT9XX7UBngRkZErWR+2Ycwt8/1F6suj34FXZ8CBP4OJxw84NdzYqWRMy48iItWuFH3YNjSEcNOR3l9Eer2Ugicvhd3jsOcX1zu9sbFTIiJSvcw5N7IbmK0BZuI1zB3sUOfc2BF9gwJNnjzZzZ2rB3tSY968H+6cClsfCsfcDqPG+B2RiIiUgJnNc85N3th1pdjD9hzwaefcS0ME8XIJ7i9S3959GWadBKEdYcqNStZEROpQKfawXQtsPcy5q0twf5H6tfpdL1lb3QVHzYBNtvI7IhER8UEp2npcuoFzfnWpEql+rgfuOxfefsRr3zHhvX5HJCIiPilZWw8RKbHHfwgv/x0OuBx2/JDf0YiIiI8KTtjM7K4NnO4BcsA84LfOuTeKDUykGqQyKRLpBJ25TiLhCMlYktbmEjRr6/wbZC6GXc6Bvb8y8vuJiEhVK+YJmwF7AhOBl4A3gO2AXYDXe99/EPiSmU11zj1ZolhFAiWVSRHviNPV7c2ryuayxDviACNL2t6aB/edA1sfDgdfAzZUAbaIiNSTYooOrgBWAJOcc7s55w53zu0GHNR7fDqwB7AI0B42qVmJdGJtstanq7uLRDpR/E2Xvw4zT4JNtoEpN8GoTUYYpYiI1IJiErZLgIudcw/3P+icm4eXrF3inHsFuBw4auQhigRTZ66zoOMbtXo5zDoZVr0NU2fA2O1GEJ2IiNSSYhK2PRk0M7SfRcDuvb9+ARhXTFAi1SASjhR0fIOcgwc+C4sfgMP/DFu8b4TRiYhILSkmYZsPfGaYc/He8+D1ZltcxP1FqkIyliTUGBpwLNQYIhkrYifAk5fB/BTsdwns/JESRSgiIrWimKKDHwB/NrPHgBuBhcC2wKnAvsDZvdcdB9xfiiBFgqivsGDEVaKv3AqPfhuazoL3frsMkYqISLUrapaomR2Pt19tEtAIdANzge875+7svWZTYI1zrrt04W6cZolKVXn7MbjjcBi/Dxw3C0ZXdPSuiIj4rKyzRJ1zdwB3mFkD3tLnm865nkHXrCjm3iJ1Y8VCmDUNGsNw1K1K1kREZFgjmnTQm6QtLFEsIvVjzUqYfSqseAOOmw2hHfyOSEREAqyohM3MJgJfAaYCW+IVF/wHuMI5t6Bk0YnUIufgwc/DorvhiOtgq40+CRcRkTpXcJWome0JPAL8F/AO8ADwLnAh8IiZ7VHSCEVqzTNt8OLvYd/vQtOZfkcjIiJVoJgnbJcBS4FDnHPz+w6aWRNwe+/5U0oSnUitee02ePirsPMp0Hyx39GIiEiVKKYP2zHAd/snawDOuSxwce95ERks9xTccxaEm+GwP4IV85+fiIjUo2L+xhgDLBvm3LLe8yLS38rFMLMFRm3qjZ0arSEgIiKSv2IStkeAC3pbeqxlZgac33tepOqlUhCNQkOD95pKFXmjnm64+wzoehmm3ALjihhdJSIida3YSQf/CzxlZtcDrwPbA6cDewAfKl14Iv5IpSAeh64u7302670HaC1wkAHzLoQ37oJDr4VtDitpnCIiUh8KfsLmnPsn8GG85c8E8EvgO3gVox92zt0+koDM7DQzu9HMsma23MyeMbNLzWzzkdxXpBCJxLpkrU9Xl3e8IM9eBc/9Cvb5Ouz6iZLFJyIi9aWo0VRrP2wWArYA3nbOdW3s+jzvOQfoBG4FXgEOwCtmeBo4fPBEhcE0mkpKoaHBa5c2mBn0bPB3YD8L0vDvE2HiB+CoW6BhVEljFBGR6lfW0VR9epO0kiRq/bQ45xb1ez/TzN4CrgWOBu4q8fcTWU8k4i2DDnU8L0ufg7tPh/F7P4UI9AAAHo9JREFUwxEpJWsiIjIieS2JmlmPma3J82v1SAIalKz1ebD3dceR3FskX8kkhEIDj4VC3vGNWrUEZrV4bTumzoDG8WWJUURE6ke+T9h+ABS/djpyU3tfn/IxBqkjfYUFiQR0dnpP1pLJPAoOelZ7vdaWvQDH3gmb7Vr2WEVEpPaNaA9bJZjZjsDDwKPOueOHuSYOxAEikcik7FBrWSKVMO9L3uipg9th98/6HU3eMqkM6USaXGeOcCRMLBmjubXZ77BERGpeRfawlZuZbYZXfLAa+ORw1znn2oF28IoOKhOdyCDP/8ZL1va6cIPJWtCSo0wqQ0e8g+6ubgBy2Rwd8Q4AJW0iIgER2Nk4ZrYpMAPYFTjROfeKzyGJDG/hLJh7Pmx/Ahzw02Ev60uOctkcuHXJUSaVyevbpDIpom1RGqY3EG2LksoU2813nXQivTZZ69Pd1U06kR7xvUVEpDQCmbCZWSNwI3Aw8EHnXH5/m4n44Z2XYPapMG4XOPJ6aBj+wfVIkqNUJkW8I042l8XhyOayxDviI07acp25go6LiEjlBS5h6x15lQJiwEnOuTk+hyQyvO5lMHOaV2wwtQPGTNjg5SNJjhLpBF3dA7vodHV3kUgX2s13oHAkXNBxERGpvMAlbHiTE04Hfgq8a2aH9vvayefYRNbpWQP3tsLSp2DK32D8nhv9yEiSo85cZ0HH8xVLxmgMNQ441hhqJJaMjei+IiJSOkFM2D7Q+5oA7hv09Rm/ghJZz2MJeLUDJv0PbH9cXh8ZSXIUCQ/dtXe44/lqbm2mpb2FcFMYDMJNYVraW1RwICISIIGrEnXORf2OQWSjXvoTPHkZ7H4e7HF+3h/rS4KKqRJNxpLEO+IDlkVDjSGSsXy6+W48LiVoIiLBFbiETSTw3pwD938Gtj0aJv/cGzBagGKTo9Zmr2tvIp2gM9dJJBwhGUuuPS4iIrUr8I1zC6Xh71JW774M/zoIRo+DEx+ATbbyOyIREaliNdE4VyRQVr8Ls6bB6i6I3aVkTUREKiaIRQciFZVXM1rXA/edA0segyOug/B7Kh+oiIjULT1hk7rW14y2byN/XzNaYODesMx0ePlGOOC/YccP+hGqiIjUMT1hk7qWVzPa7PXw+A9g10/C3l+qcIQiIiJK2KTObbQZ7eK5MOdc2OYIOOhXBVeEioiIlIISNqlrG2xG2/Ua3f/8EEvfHMvlpx9G226/yntIu4iISCkpYZO6lowlCTWGBhwLNYb48THfp+vmE3Erl5C67Ey6cpuRy+boiHcoaRMRkYpTwiZ1rbW5lfaWdprCTRhGU7iJ9g9fw1nv3MFY9wQ3XXUKC1/efu313V3dpBNpHyMWEZF6pCpRqXutza0DK0Kf+BFk/8pdf4vxzLy917s+15mrYHQiIiJK2EQGevlmeDQBTWeTeXQysHS9S8KRcOXjEhGRuqYlUZE+bz8K930ctjoYDvkNseRxNIYaB1zSGGokloz5FKCIiNQrPWETAVixEGZOg8YJcNQtMHrs2gHt6USaXGeOcCRMLBkranD7cFKZDIl0ms5cjkg4TDIWo7W5dPcXEZHaoIRNZM1KmH0KrFwEx8+GsRPXnmpubS5pgtZfKpMh3tFBV3c3ANlcjnhHB4CSNhERGUBLolLfnIMHz4NF98Bh18KWkyr2rRPp9NpkrU9XdzeJtKpQRURkICVsUt+evgJe/APs+32InF7Rb92ZG7radLjjIiJSv5SwSf169f/g4a/BzqdB8/cq/u0j4aGrTYc7LiIi9UsJm9Sn3JNwz1mwxf5w2B/AKv+fQjIWI9Q4sAo11NhIMqYqVBERGUhFB1J/Vi6GmS0wOgRH3Qqjx/kSRl9hgapERURkY5SwSX1Zswpmnwpdr8Jx/4FxO/saTmtzsxI0ERHZKCVsUj+cg3kXwMKZcNifYOtD/Y5IREQkL9rDJvXj2V/C8+3wnm/CLh/zOxoREZG8KWGT+vD6HfDQRbBjC7wv6Xc0IiIiBVHCJrVv6bNw9xkwfh84POVLRaiIiMhI6G8uqW2r3vYqQhtGw9QOaNzc74hEREQKpqIDqV09q+Hus+Ddl+DYNGwW9TsiERGRougJm9Suh74CC26Hg35F6o4pRKPQ0ADRKKRSfgcnIiKSPz1hk9r0fDs8+3PY60uk5nyaeBy6urxT2SzE496vW1v9C1FERCRfesImteeNmfDgF2Di++GAn5BIrEvW+nR1QSLhT3giIiKFUsImteWdF+HuU2Hz3eGI66BhNJ2dQ1863HEREZGgUcImtaN7KcycBq4HjpoBY8IARCJDXz7ccRERkaBRwia1oWcN3NMKS5+GI/8O4/dYeyqZhFBo4OWhkHdcRESkGihhk9rw6Lfgtf+FyVfC9scOONXaCu3t0NQEZt5re7sKDkREpHqoSlSq34vXwlOXwx7nwx6fH/KS1lYlaCIiUr30hE2q26J74YE4bHcsTGrzOxoREZGyUMIm1evdTpj9EQhF4Mi/QUOj3xGJiIiUhRI2qU7d73gVoWtWejNCN9nS74hERETKRnvYpPq4HphzDuQyMPUfEN7b74hERETKSgmbVJ/MxfDyTXDgz2CH9/sdjYiISNlpSVSqy/zr4PEfwm6fhr0u9DsaERGRilDCJtVj8YNw/ydhmykw+SqvqZqIiEgdUMIm1aHrVZh1Emy6PUy5EUaN8TsiERGRilHCJsG3ejnMOhm6l8HUGbDpNn5HJCIiUlEqOpBgcw7u/xS8NQ+OugUmNPsdkYiISMXpCZsE2xM/gux1sP+lsNO0jV6eSkE0Cg0N3msqNcQ1mQzRtjYapk8n2tZGKpMpedgiIiKlpCdsElwv3wyPfQeiH4N9vr7Ry1MpiMehq8t7n81672HdHNFUJkO8o4Ou7m7vmlyOeEeHd02znt6JiEgw6QmbBNPbj8C9H4OtDoFDfp1XRWgisS5Z69PV5R1fe006vTZZW3tNdzeJdLoUUYuIiJSFnrBJ8Cx/wxs7tcmW3r61UZsOeVkqkyGRTtOZyxEJh8mOjwHrPyXr7Oz361xuyHsNd1xERCQIlLBJsKxZ6Q10X7kYjr8bxm4/5GVDLW3aSR04gEwzzWSIkSZMjncbwmRSMZpbm73EbojkLBIOl/EfSkREZGS0JCrB4Rw8EIc374PDroUtDxj20qGWNt3obuy4NM1kaKGDCeQwYLM1OTriHWRSGZKxGKHGxgGfCzU2kozFyvFPJCIiUhJK2CQ4nv5veOmP0DwdIqdt8NLhljBdOMeJo9KMYWAy193VTTqRprW5mfaWFprCYQxoCodpb2lRwYGIiASalkQlGF79Bzz8dYicAft+d6OXD7e02RQOs1nP0MlcrtM73trcrARNRESqip6wif+WPAH3fBS2OAAO/X1eFaEbWtoMR4bejzbccRERkaBTwiYVM2RT2xVvwswWGD0Opt4Ko0N53WtDS5uxZIzG0MBkrjHUSCypfWoiIlKdzDnndwwlNXnyZDd37ly/w5BBBje1BQhvvopnfnUC242aA8fNhK0P8a7NpEikE3TmOomEIyRjSVqbWwv6fplUhnQiTa4zRzgSJpb0qkRFRESCxMzmOecmb+w67WGTili/qa3jJ2d8ke1sJhyaGpCsxTvidHV7F2dzWeId3riCQpK25tZmJWgiIlIztCQqFdG/eS3AF0/4BfFjf82Pbv02RM9eezyRTqxN1vp0dXeRSCcQERGpV3rCJhURiXizPQGOb76dto9fxC1zT+LXD/yQb/e7rjPXSfNjzcTSMcK5MLlwjnQszeP7Pe5L3CIiIkGgJ2xSEckkhEKw58RnuOGCM3j85X2J/+HPXHKJ91swlUkRbYuy72P70tLRwoTcBAxjQm4CLR0tTH1uqs//BCIiIv7REzapiNZWGMPbHPjmNFatHsP518/gZ1duRmvrwH1rF6UvYkz3mAGfHdM9huPuOs6nyEVERPynhE0qo2c1p088A0a9BMfexT2fb1p7qv++tXBu6F5pqxesrkiYIiIiQaQlUamMh74MC+6Eg66GbY8ccKozt64iIRceekqBmt6KiEg9U8Im5ffcNfDslbD3l2G3T613OhKOrP11OpZmVeOqAefV9FZEROqdEjYprzf+A3O/CBPfD/v/ZMhLkrEkoUZvwkFmvwwdLR3kJuTAINwUpqW9RT3VRESkrmkPm5TPOy/C7FNh8z3giOugYdSQl/U1xO2bbrB0ylImXTyp4OkGIiIitUqjqaQ8upfC7YfB8gVw4v2w+e55fzSVyZBIp+nM5YiEwyRjMVqb9YRNRERqj0ZTiX961sA9H4Wlz8KxtxecrMU7Oujq7gYgm8sR7+gAUNImIiJ1S3vYpPQe/Sa89n8w+UrY7piCPppIp9cma326urtJpNOljFBERKSqKGGT0nrxD/DUT2GPL8Ae5xX88c7c0G09hjsuIiJSD5SwSeksugce+BxsF4NJPyvqFpHw0P3WhjsuIiJSD5SwSWm8m4XZp0AoAkfeAA2NAGRSGdqibUxvmE5btI1MKrPB2yRjMUKNjQOOhRobScbUh01EROqXig5k5LrfgZknwZqVEOuATbYEvGStI95Bd5e3Jy2XzdER9woIhuur1ldYoCpRERGRdZSwyci4HrjvE5DLwNR/QHjvtafSifTaZK1Pd1c36UR6g41wW5ublaCJiIj0o4RNRuax78ErN8OBP4Md3j/gVK5z6EKB4Y6LiIjI0LSHTYo3/y/wRBJ2+zTsdeF6p4cb2K5B7iIiIoVRwiZF+WfqAVbO/BSznp7CHmdfReovtt41sWSMxtDAAgINchcRESmclkSlYDelXuWw3Mm8tmoip7bdyJvLxhCPe+da+43/7Nunlk6kyXXmCEfCxJIxDXIXEREpkJ6wSWFWd7Hbqyez2SbLmHbFDN5ctg0AXV2QSAy8NJXJ0LIozZc/meP3V4TZpkPJmoiISDH0hE3y5xzM+RTNO8zj5J/dwuMvD0y+OjvX/VozQUVEREpHT9gkf49fAp3X85PbL6XjoWnrnd5yaopoWxSbbnzspt9pJqiIiEiJKGGT/HTeCJnvQfTj7Hzi1wmFBp5unJRi2TFxsrls75GhK0E1E1RERKRwSthk49562GuOu9WhcEg7ra1Gezs0NYGZ9zr+IwlWua5+Hxo6MdNMUBERkcIpYZMNW74AZp0Em2zJjdv8iOiVV9MwfTqJRW0kOzL09MD8+fDW6s5BH0wDqwYc0UxQERGR4phzzu8YSmry5Mlu7ty5fodRG9asgDuPgSWP8X+7/J7T73x+wL60UGMjP214Hyuvfo4l2SXkwjnSsTSZ/foGvDcDMSDMKHuXaz/yCRUciIiI9GNm85xzkzd2napEZWjOwf1xWDwHjvw759/08npFBLvN6+a1jrmM7gbDmJCbQEtHC0Bv0uZ92eoQ157RrmRNRESkSFoSlaE9dTnM/xM0/wAipw5ZLBBLw+iBORxjuscQS8dgzShwwJIm3K3ttDa3rvd5ERERyY+esMn6XumAR74JkTNh3+8AXrFAdlDSFh6m4DOcmwA/XL32fVNT2SIVERGpC3rCJgMteRzuPRu2PBAO/Z1XBgokYzFCjQPngi6dMPQtcv1aeoRCkEyWLVoREZG6oIRN1lmxCGa2QOPmcNStMHpds7XW5mbaW1poCocxoCkcZo9vTF5vuLuNaeSRrWJr2320tw+cLyoiIiKF05KoeNasgrtPg+Wvw3GzILTjepe0NjevVziQ2Tmy3nD372leqIiISEkpYROvInTu+bBwFhyegq0Pzvujza3NGuguIiJSZloSFXjm5/DCb+G9CYie7Xc0IiIiMogStnr32r/g4S/DTifDfj/wOxoREREZghK2epZ7Gu45E8L7wmF/AtNvBxERkSDS39D1auVbXkVowxiYOgMaN/M7IhERERmGig7qUU833H0GdHVC7C4Yp862IiIiQRa4J2xmtpOZXWlm95lZl5k5M4v6HVdNmfcleCMNB18D2xzhdzQiIiKyEYFL2IDdgTOAt4HZPsdSe577FTz3S9jnq7DruX5HIyIiInkIYsI2yzm3nXPug8Df/A6mpiy4C+ZeADt8EN73Y7+jERERkTwFLmFzzvX4HUNNWva8N8lg/F5wxF+hYZTfEYmIiEieApewSRmsysHMaYDBUTOgcbzfEYmIiEgBaiJhM7O4mc01s7mLFi3yO5xg6VkD93wUlj0HU/4Om+/md0QiIiJSoJpI2Jxz7c65yc65ydtss43f4QTLI1+H12+Dg34J2x3jdzQiIiJShJpI2GQYL/wOnr4C9rwAdo/7HY2IiIgUSQlbrVp4Nzx4Hmx/HBx4hd/RiIiIyAgoYatF78yH2afAuCgceQM0aKCFiIhINQvk3+RmdlrvLyf1vn7AzBYBi5xzM30Kqzp0L4NZ07zxU1M7YMwWfkckIiIiIxTIhI31G+Ze1fs6Ezi6sqFUEdcD934Mck/C0bd5PddERESk6gUyYXPOmd8xVKVHvwOvzoBJP4eJx/sdjYiIiJSI9rDVipdS8OSlXjXonl/0OxoREREpISVsteDN++H+T8O2U2HSlWB6QCkiIlJLlLBVu65XYNbJMHYHOPLvMGqM3xGJiIhIiQVyD5vkafW7MPMk7/XYO2HTrf2OSERERMpACVu1cj1w37nw9sNe+44J7/U7IhERESkTJWzV6vEfwst/hwMuhx0/5Hc0IiIiUkbaw1aNOv8GmYthl3Ng76/4HY2IiIiUmRK2avPWQ3DfObD1YXDwNaoIFRERqQNK2KrJ8tdh5jTYZGuYcjOM2sTviERERKQCtIetWqxe7rXvWPU2nHAPjN3O74hERESkQpSwVQPn4IHPwuIHYMpNsMX+fkckIiIiFaQl0Wrw5GUwPwX7XQI7f8TvaERERKTClLAF3Su3wqPfhqaz4L3f9jsaERER8YEStiB7+zG4txW2nASH/E4VoSIiInVKCVtQrVgIs6ZBYxiOuhVGj/U7IhEREfGJig6CaM1KmH0qrHgDjpsNoR38jkhERER8pIQtaJyDBz8Pi+6GI66DrSb7HZGIiIj4TEuiQfNMG7z4e9j3u9B0pt/RiIiISAAoYQuS126Dh78KO58CzRf7HY2IiIgEhBK2oMg9BfecBeFmOOyPYPpXIyIiIh5lBUGwcjHMbIFRY2HqDBg9zu+IREREJEBUdOC3nm64+3Toehli/4FxEb8jEhERkYBRwua3eRfCG/+GQ6+FbQ7zOxoREREJIC2J+unZq+C5X8E+X4ddP+F3NCIiIhJQStj8siAN8/4LdvgwvO9HfkcjIiIiAaaEzQ9Ln/P2rY3fG45IQcMovyMSERGRAFPCVmmrlsCsFq9tx9QZ0Dje74hEREQk4FR0UEk9q71ea8tegGPvhM129TsiERERqQJK2Crp4a/B6/+Cg9thu6l+RyMiIiJVQkuilfL8b7w5oXtdCLt/1u9oREREpIooYauEhbNg7vmw/QlwwE/9jkZERESqjBK2cnvnJZh9KozbBY68Hhq0Ci0iIiKFUcJWTt3LYOY0r9hgageMmeB3RCIiIlKF9LinXHrWwL2tsPQpOPo2GL+n3xGJiIhIlVLCVi6PfQde7YBJV8LE4/2ORkRERKqYlkTL4aU/wZM/ht0/B3t+we9oREREpMopYSu1N+fA/Z+BbY+GyVeCmd8RiYiISJVTwlZK774Ms06G0E4w5e/Q0Oh3RCIiIlIDtIetVFa/C7NOgtVdELsLNtnK74hERESkRihhKwXXA/edA28/AlP/F8Lv8TsiERERqSFK2Eoh8wN4+UZvisGOH/Q7GhEREakx2sM2Utkb4PHpsOu5sPeX/Y5GREREapAStpF4ax7MORe2PhwOuloVoSIiIlIWStiKtfx1mHkSbLINHHUzjNrE74hERESkRmkPWzFWL/eSte4lcPw9sOm2fkckIiIiNUwJW6Gc8xrjvvUgTLkZtnif3xGJiIhIjdOSaKGe/DFk/wLvS8LOJ/sdjYiIiNQBJWyFcA5yT0HTR+E93/I7GhEREakTWhIthBkcdi30dKsiVERERCpGT9gKZQajxvgdhYiIiNQRJWwiIiIiAaeETURERCTglLCJiIiIBJwSNhEREZGAU8ImIiIiEnBK2EREREQCTgmbiIiISMApYRMREREJOCVsIiIiIgGnhE1EREQk4JSwiYiIiAScEjYRERGRgFPCJiIiIhJwSthEREREAk4Jm4iIiEjAKWETERERCTglbCIiIiIBp4RNREREJOCUsImIiIgEnBI2ERERkYBTwiYiIiIScErYRERERALOnHN+x1BSZrYIyPodR0BtDbzpdxB1Qj/rytHPunL0s64c/awrx++fdZNzbpuNXVRzCZsMz8zmOucm+x1HPdDPunL0s64c/awrRz/ryqmWn7WWREVEREQCTgmbiIiISMApYasv7X4HUEf0s64c/awrRz/rytHPunKq4metPWwiIiIiAacnbCIiIiIBp4RNREREJOCUsNUwM9vJzK40s/vMrMvMnJlF/Y6rFpnZaWZ2o5llzWy5mT1jZpea2eZ+x1ZrzOxEM7vLzBaY2Uoze8XMbjCz9/gdW60zs3/2/jlyid+x1BozO7r3Zzv4a4nfsdUqM/ugmc0ys3fMbKmZzTWzY/2Oazij/Q5Aymp34AxgHjAbOMHfcGraV4FO4NvAK8ABwMXAMWZ2uHOux8fYas2WeL+nrwIWARHgm8AcM2t2zqlxdhmY2UeB9/kdRx34L+DBfu9X+xVILTOzzwG/6P36Id4DrP2BkJ9xbYgStto2yzm3HYCZfQYlbOXU4pxb1O/9TDN7C7gWOBq4y5eoapBz7q/AX/sfM7MHgKeB04D/9iOuWmZmE4CfAV8C/uJzOLXuKefcHL+DqGW9K01twNecc239Tv3Ll4DypCXRGqanOpUzKFnr0/d/yTtWMpY6tbj3tdvXKGrXT4AnepNlkWr3KaAHuNrvQAqhhE2kfKb2vj7laxQ1ysxGmdkYM9sDuAZYAFznc1g1x8yOBD4BnO93LHUiZWZrzGyxmf3FzCJ+B1SDjsR7In+Wmb1gZqvN7Hkz+4LfgW2IlkRFysDMdgR+ANzpnJvrdzw16n5gUu+vnweOdc4t9DGemmNmjXjJ8E+dc8/4HU+Ny+Et588EluLtg/02cJ+ZHaDf2yW1Q+/X5Xg/4xeA04FfmNlo59z/+BnccJSwiZSYmW0G3Iq3WfiTPodTyz4OjAd2xSv6uMPMjnTOzfc1qtryDWAskPQ7kFrnnHsYeLjfoZlmNgt4AK8Q4Tu+BFabGoDNgXOdczf1Hrurd2/bt8zs5y6AUwW0JCpSQma2KTADL4k40Tn3is8h1Szn3FPOuft791XFgM3wqkWlBHqX4hLAd4FNzGxCb/EB/d6P8i/C2uecewh4FjjI71hqTN+e1zsGHb8d2A6YWNlw8qOETaREepePbgQOBj7onMv4HFLdcM4twVsW3d3vWGrIrsCmwJ+Bt/t9gfdE822g2Z/Q6ooBgXvaU+WeGOa49b4GsmBPCZtICZhZA5DCe9JzksryK8vMtgP2xtuLIqXxCHDMEF/gJXHH4CXJUiZmNhnYE2+/ppTOzb2vJw46fiLwinNuQYXjyYv2sNU4Mzut95d9m7M/YGaLgEXOuZk+hVWLfom3aTUJvGtmh/Y794qWRkvHzG4GHgIew9ucvSdef7DVqAdbyfQ+tfzP4ONmBpB1zq13TopnZingJbzf20vwig6+BbwKXOljaLXo/4B/A9eY2dbAi3g9HE8gwPuOLYD76qSEzGy4f8EznXNHVzKWWmZm84GmYU5Pd85dXLloapuZfQNvgsduwBjgZbzE4lIVHJRf758pSeecNsGXkJl9C/go3p8jIbw2NbcB33fOve5nbLXIzMYDl+Ilalvgtfn4sXMusI2hlbCJiIiIBJz2sImIiIgEnBI2ERERkYBTwiYiIiIScErYRERERAJOCZuIiIhIwClhExEREQk4JWwiIiIiAaeETURERCTglLCJiAzBzBrN7Dtm9qKZLTez+8xsLzObbGZdZraD3zGKSP3QpAMRkUHMrBFvLNB+ePMc3wCuwhvCvSnwtHPua/5FKCL1RsPfRUTWdz5wLHCkc+5eADM7CLgQMOBc/0ITkXqkJVERkfWdB9zel6z1WgKEgSucc4v9CUtE6pUSNhGRfsxse2Bv4P8GnRoDvAVcUfGgRKTuKWETERlo997Xl/oOmNko4BPAC865Zb5EJSJ1TQmbiMhAPb2vW/Y79nngPcCoyocjIqIqURGRAcxsHPAqsBj4CrAT8BO8JdIW4GTg3865Fb4FKSJ1RwmbiMggZnY88DNgT2AR8HXgX3hJ2yRgvHPuXf8iFJF6o4RNREREJOC0h01EREQk4JSwiYiIiAScEjYRERGRgFPCJiIiIhJwSthEREREAk4Jm4iIiEjAKWETERERCTglbCIiIiIBp4RNREREJOD+H7dSMoL+F0ouAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [10,10]\n",
    "matplotlib.rcParams.update({'font.size': 16})\n",
    "\n",
    "plt.title(r\"$\\log_{10}\\Vert W\\Vert_{F}^{2} / \\log_{10}\\;\\lambda_{max}$ vs. $\\alpha\\;$\")\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "plt.ylabel(r\"$\\log_{10}\\Vert W\\Vert_{F}^{2} / \\log_{10}\\;\\lambda_{max}$\")\n",
    "\n",
    "for im, M in enumerate(Ms):\n",
    "    [ratios,alphas] = Mresults[M]\n",
    "    plt.scatter(alphas, ratios, label=\"Q=1, M= {}\".format(M), color=Mcolors[im])\n",
    "\n",
    "plt.plot(mus, mus,  color='orange')\n",
    "plt.legend()\n",
    "#plt.savefig(\"img/Alpha-LogNorm-Relations.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [10,10]\n",
    "matplotlib.rcParams.update({'font.size': 16})\n",
    "\n",
    "plt.title(r\"$\\log_{10}\\Vert W\\Vert_{F}^{2} / \\log_{10}\\;\\lambda_{max}$ vs. $\\alpha\\;$\")\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "plt.ylabel(r\"$\\log_{10}\\Vert W\\Vert_{F}^{2} / \\log_{10}\\;\\lambda_{max}$\")\n",
    "\n",
    "M=1000\n",
    "im=4\n",
    "[ratios,alphas] = Mresults[M]\n",
    "plt.scatter(alphas, ratios, label=\"Q=1, M= {}\".format(M), color=Mcolors[im])\n",
    "\n",
    "plt.plot(mus, mus,  color='orange')\n",
    "plt.legend()\n",
    "#plt.savefig(\"img/Alpha-LogNorm-M1000.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [10,10]\n",
    "matplotlib.rcParams.update({'font.size': 16})\n",
    "\n",
    "plt.title(r\"$\\log\\Vert W\\Vert_{F}^{2} $ vs. $\\alpha\\;$\")\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "plt.ylabel(r\"$\\log\\Vert W\\Vert_{F}^{2}$\")\n",
    "\n",
    "for im, M in enumerate(Ms):\n",
    "    [ratios,alphas] = Mresults[M]\n",
    "    plt.scatter(alphas, logNorm2s[M], label=\"Q=1, M= {}\".format(M), color=Mcolors[im])\n",
    "    \n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heavy Tailed  normalization on X\n",
    "\n",
    "$$\\dfrac{1}{N^{2/\\mu}}X$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-09T22:54:39.475194Z",
     "start_time": "2018-11-09T22:42:38.351508Z"
    }
   },
   "outputs": [],
   "source": [
    "Mresults = {}\n",
    "l_maxs = {}\n",
    "logNorm2s = {}\n",
    "\n",
    "for M in Ms:\n",
    "    ratios, alphas = [], []\n",
    "    l_maxs[M] = []\n",
    "    logNorm2s[M] = []\n",
    "    \n",
    "    N = M\n",
    "    for mu in tqdm(mus):\n",
    "\n",
    "        W = np.random.pareto(a=mu,size=(N,M))\n",
    "        normW = np.linalg.norm(W)\n",
    "        logNorm2 = 2.0*np.log10(normW)\n",
    "        logNorm2s[M].append(logNorm2)\n",
    "\n",
    "        # CHANGED NORMALIZATION HERE\n",
    "        X=np.dot(W.T,W)/N**(2/mu)\n",
    "        evals = np.linalg.eigvals(X)\n",
    "\n",
    "        l_max, l_min = np.max(evals), np.min(evals)\n",
    "        fit =  powerlaw.Fit(evals)\n",
    "\n",
    "        a = fit.alpha\n",
    "        alphas.append(a)\n",
    "\n",
    "        ratio = logNorm2/np.log10(l_max)\n",
    "        ratios.append(ratio)\n",
    "\n",
    "        Mresults[M]=[ratios,alphas]\n",
    "        \n",
    "        l_maxs[M].append(l_max)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [10,10]\n",
    "matplotlib.rcParams.update({'font.size': 16})\n",
    "\n",
    "plt.title(r\"$\\log\\Vert W\\Vert_{F}^{2} / \\log\\;\\tilde{\\lambda}_{max}$ vs. $\\mu;$\"+\"\\n\"+r\"$(1/N^{2/\\mu})\\mathbf{X}$ normalization\")\n",
    "plt.xlabel(r\"$\\mu$\")\n",
    "plt.ylabel(r\"$\\log\\Vert W\\Vert_{F}^{2} / \\log\\;\\tilde{\\lambda}_{max}$\")\n",
    "\n",
    "for im, M in enumerate(Ms):\n",
    "    [ratios,alphas] = Mresults[M]\n",
    "    plt.scatter(mus, ratios, label=\"Q=1, M= {}\".format(M), color=Mcolors[im])\n",
    "\n",
    "plt.legend()\n",
    "#plt.xlim(1,2)\n",
    "plt.ylim(0,5)\n",
    "\n",
    "plt.savefig(\"img/LogNorm-Lmax-Scaled.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [10,10]\n",
    "matplotlib.rcParams.update({'font.size': 16})\n",
    "\n",
    "plt.title(r\"$\\log\\Vert W\\Vert_{F}^{2} $ vs. $\\alpha\\;$\"+\"\\n\"+r\"$1/N^{2/\\mu}$ X normalization\")\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "plt.ylabel(r\"$\\log\\Vert W\\Vert_{F}^{2}$\")\n",
    "\n",
    "for im, M in enumerate(Ms):\n",
    "    [ratios,alphas] = Mresults[M]\n",
    "    plt.scatter(alphas, logNorm2s[M], label=\"Q=1, M= {}\".format(M), color=Mcolors[im])\n",
    "    \n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Max Value  normalization on W\n",
    "\n",
    "$$\\dfrac{1}{N^{1/\\mu}}\\mathbf{W},\\;\\;\\;\\dfrac{1}{N}\\mathbf{X}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mresults = {}\n",
    "l_maxs = {}\n",
    "logNorm2s = {}\n",
    "\n",
    "for M in Ms:\n",
    "    ratios, alphas = [], []\n",
    "    l_maxs[M] = []\n",
    "    logNorm2s[M] = []\n",
    "    \n",
    "    N = M\n",
    "    for mu in tqdm(mus):\n",
    "\n",
    "        # CHANGED NORMALIZATION HERE\n",
    "        W = np.random.pareto(a=mu,size=(N,M))\n",
    "        W = W/N**(1/mu)\n",
    "        \n",
    "        normW = np.linalg.norm(W)\n",
    "        logNorm2 = 2.0*np.log10(normW)\n",
    "        logNorm2s[M].append(logNorm2)\n",
    "\n",
    "        # STANDARD NORMALIZATION HERE\n",
    "        X=np.dot(W.T,W)/N\n",
    "        evals = np.linalg.eigvals(X)\n",
    "\n",
    "        l_max, l_min = np.max(evals), np.min(evals)\n",
    "        fit =  powerlaw.Fit(evals)\n",
    "\n",
    "        a = fit.alpha\n",
    "        alphas.append(a)\n",
    "\n",
    "        ratio = logNorm2/np.log10(l_max)\n",
    "        ratios.append(ratio)\n",
    "\n",
    "        Mresults[M]=[ratios,alphas]\n",
    "        \n",
    "        l_maxs[M].append(l_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [10,10]\n",
    "matplotlib.rcParams.update({'font.size': 16})\n",
    "\n",
    "plt.title(r\"$\\log\\Vert W\\Vert_{F}^{2} / \\log\\;\\lambda_{max}$ vs. $\\alpha\\;$\"+\"\\n\"+r\"$(1/N^{1/\\mu})\\mathbf{W}$ normalization\")\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "plt.ylabel(r\"$\\log\\Vert W\\Vert_{F}^{2} / \\log\\;\\lambda_{max}$\")\n",
    "\n",
    "for im, M in enumerate(Ms):\n",
    "    [ratios,alphas] = Mresults[M]\n",
    "    plt.scatter(alphas, ratios, label=\"Q=1, M= {}\".format(M), color=Mcolors[im])\n",
    "\n",
    "plt.ylim([-1,100])\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.log10(np.array(l_maxs[500]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [10,10]\n",
    "matplotlib.rcParams.update({'font.size': 16})\n",
    "\n",
    "plt.title(r\"$\\log\\Vert W\\Vert_{F}^{2} $ vs. $\\alpha\\;$\"+\"\\n\"+r\"$\\dfrac{1}{N^{1/\\mu}}$ W normalization\")\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "plt.ylabel(r\"$\\log\\Vert W\\Vert_{F}^{2}$\")\n",
    "\n",
    "for im, M in enumerate(Ms):\n",
    "    [ratios,alphas] = Mresults[M]\n",
    "    plt.scatter(alphas, logNorm2s[M], label=\"Q=1, M= {}\".format(M), color=Mcolors[im])\n",
    "    \n",
    "plt.plot(mus, mus,  color='orange')\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Max Value  normalization on W\n",
    "\n",
    "$$\\dfrac{\\mathbf{W}}{\\max|\\mathbf{W}|},\\;\\;\\;\\dfrac{1}{N}\\mathbf{X}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mresults = {}\n",
    "l_maxs = {}\n",
    "logNorm2s = {}\n",
    "             \n",
    "for M in Ms:\n",
    "    ratios, alphas = [], []\n",
    "    l_maxs[M] = []\n",
    "    logNorm2s[M] = []\n",
    "    \n",
    "    N = M\n",
    "    for mu in tqdm(mus):\n",
    "\n",
    "        # CHANGED NORMALIZATION HERE\n",
    "        W = np.random.pareto(a=mu,size=(N,M))\n",
    "        W = W/np.max(W)\n",
    "        \n",
    "        normW = np.linalg.norm(W)\n",
    "        logNorm2 = 2.0*np.log10(normW)\n",
    "        logNorm2s[M].append(logNorm2)\n",
    "\n",
    "        X=np.dot(W.T,W)/N\n",
    "        evals = np.linalg.eigvals(X)\n",
    "\n",
    "        l_max, l_min = np.max(evals), np.min(evals)\n",
    "        fit =  powerlaw.Fit(evals)\n",
    "\n",
    "        a = fit.alpha\n",
    "        alphas.append(a)\n",
    "\n",
    "        ratio = logNorm2/np.log10(l_max)\n",
    "        ratios.append(ratio)\n",
    "\n",
    "        Mresults[M]=[ratios,alphas]\n",
    "        \n",
    "        l_maxs[M].append(l_max)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [10,10]\n",
    "matplotlib.rcParams.update({'font.size': 16})\n",
    "\n",
    "plt.title(r\"$\\log\\Vert W\\Vert_{F}^{2} $ vs. $\\alpha\\;$\"+\"\\n\"+r\"W/max|W| normalization\")\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "plt.ylabel(r\"$\\log\\Vert W\\Vert_{F}^{2}$\")\n",
    "\n",
    "for im, M in enumerate(Ms):\n",
    "    [ratios,alphas] = Mresults[M]\n",
    "    aa = np.array(alphas)[0:14]\n",
    "    ln2 = np.array(logNorm2s[M])[0:14]\n",
    "    plt.scatter(aa, ln2, label=\"Q=1, M= {}\".format(M), color=Mcolors[im])\n",
    "    \n",
    "#plt.plot(mus, mus,  color='orange')\n",
    "plt.legend()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [10,10]\n",
    "matplotlib.rcParams.update({'font.size': 16})\n",
    "\n",
    "plt.title(r\"$\\log_{10}\\;\\lambda_{max} $ vs. $\\alpha\\;$\"+\"\\n\"+r\"W/max |W| normalization\")\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "plt.ylabel(r\"$\\log_{10}\\;\\lambda_{max}$\")\n",
    "\n",
    "\n",
    "for im, M in enumerate(Ms):\n",
    "    [ratios,alphas] = Mresults[M]\n",
    "    aa = np.array(alphas)[0:15]\n",
    "    llm = np.log10(l_maxs[M])[0:15]\n",
    "    plt.scatter(aa, llm, label=\"Q=1, M= {}\".format(M),  color=Mcolors[im])\n",
    "    \n",
    "\n",
    "#plt.plot(mus, mus,  color='orange')\n",
    "plt.legend()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}