examples/gics_fit_cov.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6d72e5d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/parshakova.tanya/opt/anaconda3/envs/routing/lib/python3.9/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import copy, time\n",
    "import random\n",
    "import pickle\n",
    "\n",
    "import mlrfit as mf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4759ee57",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1001)\n",
    "random.seed(1001)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "573d6bef",
   "metadata": {},
   "source": [
    "#  Matrix definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b435696c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "level=0, 1\n",
      "level=1, 11\n",
      "level=2, 24\n",
      "level=3, 69\n",
      "level=4, 157\n",
      "level=5, 5000\n"
     ]
    }
   ],
   "source": [
    "rank = 30\n",
    "mtype = 'assetcov_gics'\n",
    "\n",
    "with open(\"../data_finance/crsp_gics_5k_hpart.pkl\", \"rb\") as f:\n",
    "       gics_hpart = pickle.load(f)\n",
    "\n",
    "with open('../data_finance/crsp_gics_5k_C.npy', 'rb') as f:\n",
    "       C = np.load(f)\n",
    "\n",
    "m = n = gics_hpart['rows']['pi'].size\n",
    "for level in range(len(gics_hpart['rows']['lk'])):\n",
    "    print(f\"{level=}, {gics_hpart['rows']['lk'][level].size-1}\")\n",
    "\n",
    "Z = mf.demean_clip(C)\n",
    "A = Z @ Z.T / (Z.shape[1]-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d07fff32",
   "metadata": {},
   "outputs": [],
   "source": [
    "U, Vt, sigmas = mf.frob_low_rank(A)\n",
    "frob_losses = np.zeros(m)\n",
    "normalization = np.linalg.norm(A, ord='fro')\n",
    "# assert np.allclose(normalization, ((sigmas**2).sum())**0.5)\n",
    "sigmas2 = (sigmas**2)[::-1]\n",
    "frob_losses = np.concatenate([(np.cumsum(sigmas2)**0.5)[::-1][1:] / normalization, np.array([0])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7efe9192",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAAFfCAYAAABZUFYHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAABJ0AAASdAHeZh94AAAwrUlEQVR4nO3dfWxbZ4Lf+x+pF8uWbFFUrGQm49kV6TTjSWZ3TduabbvbNBm6vi3ajhtL8f61KNpY8rhF2wVqKipwMddAAVnMon9c4DqWnBZoey/uWmaC6aItOiGdINvO7kCWmMyrZybSkWMnM45s0ccyZVlvZP+QSIsWbYnUochz9P0AAsnz8pxHfiDp5+d5znNc6XQ6LQAAABtzl7sCAAAAG0WgAQAAtkegAQAAtkegAQAAtkegAQAAtkegAQAAtkegAQAAtkegAQAAtkegAQAAtldd7gpUmng8rosXL+rQoUMyDEOhUKjcVQIAAGtw8eiDXAcOHNDIyIgkKRaLKR6PE2oAAKhwDDmtEI/Hcz77fD719/eXqTYAAGC9HDnkZJqmBgcH1d/fn+1tWSkcDsvj8UiSxsbG1NPTI4/Ho+HhYfl8vuxxXq9XiUSi6Dp8+OGH2rNnj7Zt21ZUGQAAbEWzs7O6ceOGXnrppezf67U4LtDEYjEZhiHTNGWa5qr94XBYktTZ2SlpqVemo6ND0WhUkladk6+M9fjwww919OjRos4FAADS9773PX37299e17GOCzTBYFCSFIlE8u7v7e3N6bUJBAIaHh6WYRjy+Xw5PTKJRCKnx6YQe/bskbTUGHv37i2qjEclk0kNDQ2pra1NDQ0NlpSJ8qAtnYF2dAbasfKMjo7q6NGj2b+l6+G4QPMk8XhcpmnK6/XmbPd6vYpEIqsm/xqGofb29qKulRlm2rt3r1544YXiKvyIqakp3bx5U/v27dOuXbssKRPlQVs6A+3oDLRj5SpkysaWCjSZ3pdHx+M8Ho8mJyclSRcuXFB3d7cOHTqkK1euqK+vb81yJyYmdOvWrZxto6OjkpaS/9TUlAW1l6anp3NeYV+0pTPQjs5AO1aeZDJZ8DlbKtA8aT5MZl8gEFAgEJCkdffOnDt3TmfOnMm7b2hoSDdv3iyonmsZGhqytDyUD23pDLSjM9COleP69esFn7OlAk2mZ8Y0zZxemmIn/macOnVKHR0dOdsy439tbW3at2/fhsrPmJ6ezo7z1tfXW1ImyoO2dAba0Rlox8pz9erVgs/ZUoEmM3cmkUisGnby+/1Fl9vS0qKWlpa8+xoaGiwfk62vr2ec1yFoS2egHZ2BdqwcxUzO3lIL6wUCAXk8nlU9MoZhZO+OAgAA9uPoQJNvUbyenh7FYrHs53g8njNvBgAA2I/jhpwyD5eMxWIyTVNdXV3y+/3ZW7JDoZDC4bAGBgYkLa0UfPny5XJWGQAAbJDjAs16elt42CQAAM7i6CEnAACwNTiuh6bS/LP/L66dX7pjSVmLi4u6d69Kb3/6sbbVVMvtdqna7VLVitel9+5V+x7/2f3IuS65XS5VV+X77H742e1SVdXyq2v52KrMPvcjnx/uzx7rdqsq5/NSnQAAKAaBpsR+9cU91abuWliiSzem71lYXuVwubQq4DwudFXl+bwySD0ujOUev7y9avX2fNfOCWZ5y3Pnr1fV6rJm7s8pOS9NzczLvW0hJygS7ACgcASaEvumz6umZ/OvUVOohYUF3b59Wx5vs1xutxYW00ql01pIpbWYSq/+nEoplZIWUiktZrc9fE0tv1aKdFpaWK7/bLkrsymqpeEfrtrqcml1SFsZkqpytz/ao5b3vDwhbHUPnzvPddY6z73q+Gq3W9tq3NpW7da26qql15oV76vdcrkIbQCsRaApsX979BuWPpzygw8+0Msvv2jZ4k/pdFqptFaFn5WBZ2UQyheU1jp2MefzcuhazPc5tXRuOq3FxRXl5HxOaTGtpWOXA9z84qNlph4b3lbtX3wkAFZAvkunpfnFpe9LSpW7OiVRW50/8NStCkIPQ9CO2iptr63Wjtoq1a94v/RVvby/SvW11dq+vL2mimmCwFZBoNniXC6XqlxSlbtqeUvVE493ulQmQK3q0Uple8Hybs8b6lJ5jl/anpye0dVf/lL+vc+ppnZbTi9bNtjluc6qslYEvdXnLAfAlHL3Lz5SxzxBtNTmFlKaW0jpnhZKep2aKlc27DRsq9bOumrtrKvRru01y++rtauuJud1Z93DfTvrarRzWzXDgIANEGiAFdxul9xyqabEuW5qakof3P2FXj70bMUttb6y1+5xPVk5wSlnqDOl2eWwMpv5ml98+H5hUbPzK94vpJY/Lz7x+Afzi7o/t6iZ+cWCvpf5xbTuzszr7sx80f8eLpfk2V6jpvpaNe1Y+vLWL33e7k7piwmX9KtJfWX3gpp21Oqpndu0c1s1w2rAJiPQAMixuteucqRSaT1YWNT07KJm5hZ1f37h4fu5Bc3ML+27P7egmblFTc8tamZuQdNzi0o+WNC92Xnde7Cw/DWvqZkFzS0+eVgvnZbu3J/XnfvzkqbzHFGl/3/s5zlbttdUqWXXNrXs3KaWXXVq2blNTz/y+mXPdtVv41cwYBV+mgDYhtudGUKy7lfXg/nFbMDJhJ2pB/PZz1MzS2EmcX9Od6bnlJie0537c7ozPf/YMDQzv6hPJ+/r08n7T7y2t75WX2navvy1I/t+T9MO/VZzvWqrmQMErBeBBsCWVldTpbqaKu3eua2g89LptL6YNPXfL/+F9v3eQc25apWYntWte7P6YmpWE/dm9cXUg+XPD3R/bvVwWWI5IP34s9VLO1S5Xfqqd4f8u+vl390gf0uD/Lsb9NzTDdpVV1P09ws4FYEGAIrgcrm0o7ZKzXXSC1/aueZcqOTsgr6YeqCJqaWA87k5o8/uzOizO/f12Z0ZfX5nJqfHZzGV1vjtaY3fnlbs6kROWb/dvEMvPtuoF59t1DeebdSLX25U4w5CDrY2Ag0AbIKGbdVq2L3Uy5JPKpXWreSsPrtzX9du35dxO6mxiWmN3krq2u3pnDWjrk3e17XJ+/qvP/5NdttzLQ36fV+zft/XrLZWb8E9ToDdEWgAoAK43S49vatOT++q04Hf8ubsm19M6Ubivj6ZSOrqb6b008/v6ief39UXUw+XoPxkIqlPJpL6zz/8VJL0tWd26u+88Iz+ztef1gtf3sVdV3A8Ag0AVLiaKrd8uxvk292gIy88k90+ce+Bfvb5lD66YWpofFLx66bmFpaGrX5x855+cfOe/u/Ln+hZz3b9/d/9kv7o0FfV+lR9ub4NoKQINABgUy0769TytTq9/LWlx6s8mF/Uj26Y+suxSV3+xRf66edTkqTPzRn1f2io/0NDv+/z6uRLfr3013bTawNHIdAAgEPU1VTpm75mfdPXrD85/Nf02Z37iv38C/3XH/9Gw5/ekST90Ejoh0ZCv7vHo3/zd7+mb/qay1xrwBoscgAADvWVph36x3+zVZHv/A1F/+Rv6Z/+Qavqa5cWTPzRDVPHB36o7siPlZwt7SMogM1AoAGALeC5p3fq//z7X9cP3nhF//zlvaqrWfr1f3H4hv7R//MDjd/OtwoyYB8EGgDYQjw7avWvjzyv6J+8pL+5d2m46ZOJpDrO/6XGbiXLXDugeAQaANiC9nh36D/9k2/qO3/bL0m6nZzTH//7Id29X/yDPIFyItAAwBZV5Xap+//4mv5V8DlJS3dDfffPf1rmWgHFIdAAwBb3L7/1nF5ZvvX7ex//Wj+6YZa3QkARCDQAsMW5XC6d+YcvqKZqaV2a/r8YK3ONgMIRaAAA2uPdoX/wu1+WJEV//gVzaWA7BBoAgCTp27/3rCRpfjGtH45Plrk2QGEINAAASdKh327KDjv90CDQwF4INAAASdKO2mo9/8xOSdLoBGvSwF4INACArN9uXnoaNysHw24INACArK96d0iSfm3OKJVKl7k2wPoRaAAAWc0N2yRJqbR0d4Y7nWAfBBoAQFbTjprse5NAAxupLncFKk0sFpMkGYahaDSqCxcuyOPxlLdSALBJdtRWZd/PzC2WsSZAYeiheURHR4cOHjyozs5Oeb1e9fb2lrtKALBpttU8DDQPFgg0sA8CzSPGx8ezPTLDw8M6fPhweSsEAJto+8pAQw8NbMRWQ06maWpwcFD9/f0aGRlZtT8cDmfDyNjYmHp6egoeLvJ4PIrH47p48aKOHz+uYDBoQc0BwB7q6KGBTdkm0MRiMRmGIdM0ZZrmqv3hcFiS1NnZKUmKx+Pq6OhQNBot+FqBQEA+n08dHR0KBAKEGgBbRrXblX2/sMht27AP2wSaTKiIRCJ59/f29ub02gQCAQ0PD8swDPl8PkUikceGG4/Ho76+PknKHu/xeNTV1aWuri6NjfHkWQAAKpltAs2TxONxmaYpr9ebs93r9SoSiSgUCqm9vV3t7e1PLCcWi6m7uzsbjBKJBHc4AQBgA44INIlEQpJWhQ+Px6PJyfU/YC0YDKqrqys7vBWNRnXp0qU1z5uYmNCtW7dyto2OjkqSksmkpqam1l2HJ5mens55hX3Rls7gxHZc+b3MzMxY9vurkjmxHe0umSz8WWKOCDT55tSsZ18+mTk4j75/knPnzunMmTN59w0NDenmzZsF1WEtQ0NDlpaH8qEtncFJ7fjZtJT50/Djn/xE+nzrzKNxUjva3fXr1ws+xxGBJtMzY5pmTi9NoWGmWKdOnVJHR0fOttHRUR09elRtbW3at2+fJdeZnp7W0NCQ2traVF9fb0mZKA/a0hmc2I6/+CIp/fgjSdI3XnxRLz//VJlrVHpObEe7u3r1asHnOCLQZObO5Jvz4vf7S379lpYWtbS05N3X0NCgXbt2WXq9+vp6y8tEedCWzuCkdqxPPuyR2b5jh2O+r/VwUjvaXUNDQ8HnOGJhvUAgII/Hs6pHxjAMbrkGAGALsGWgyUwCXqmnpyf7HCZp6c6nQCCgQCCwmVUDAMdIb53pM3AA2ww5ZVbvjcViMk1TXV1d8vv9CoVCkqRQKKRwOKyBgQFJSysFX758uZxVBgDbccm19kFABbJNoFlPb0sm3AAAgK3FlkNOAAAAKxFoAACPwSQa2AeBBgCQ5WIKDWyKQAMAAGyPQAMAAGyPQAMAyIt1aGAnBBoAQBZzaGBXBBoAAGB7BBoAAGB7BBoAQF5MoYGdEGgAAFk8ywl2RaABAAC2R6ABAAC2R6ABAAC2R6ABAOTFwnqwEwINACCLhfVgVwQaAABgewQaAABgewQaAEBeaZbWg40QaAAAWUyhgV0RaAAAgO1VW1lYT0+PmpubFQgE9Morr1hZNAAAwGNZGmii0aiGh4etLBIAUCasQwM7sXTI6fjx45KkqakpK4sFAGwS1qGBXVkaaDwejyQpnU7r5MmTOnLkiN5//30rLwEAALCKpYHGtRztGxsbdf78eQWDQebSAACAkrM00IRCIR0/flxvv/22rl27pqampuy+a9euWXkpAECJMYUGdmLppODXXntNBw4c0HvvvadQKKS7d+8qGo3q8OHDGhkZ0VtvvWXl5QAAlmMSDezJ0kDT19enxsZGnThxQpI0Pj6uWCym9957T++88w6BBgAAlISlQ06NjY05n1tbW3XixAkNDg7q/PnzVl4KAAAga9NWCs702gAA7CHNQjSwER59AADIYh0a2FVRc2g2snDerl27ij53s4TD4eyaOgcPHlQgEChvhQAAwBMVHGjGx8fV19cnl8tVVHdkU1OTent7Cz5vs3R3d+vQoUNqb29XJBLRxYsXCTQAAFS4ggNNa2urYyf4mqapSCSivr4+xeNxBYNBtbe3l7taAABgDbaaQ2OapgYGBnTgwIG8+8PhsAYGBjQwMKDu7m6ZpllQ+ZkHa0YiEfl8PnV3dysej2+02gBgG0yhgV0VNeTU3d2dfcxBIdLptJqbm4tajyYWi8kwDJmmmTeohMNhSVJnZ6ckKR6Pq6OjQ9FotKDrJBKJbK9MV1eXTpw4oZGRkYLrCwAANk9RQ06Dg4OlqMsTBYNBSUu9J/n09vbmBI9AIKDh4WEZhiGfz6dIJPLYcOPxeNTX1yefzyev15uz3TAMC78LAABQCpauFFwu8XhcpmnmhBFJ8nq9ikQiCoVCam9vX3M+jM/ny97dJEmGYWSDFAAAqFwlDzRvv/22Xn/99ZJeI5FISFJOGMl8npycLKisS5cuqbu7W36/X2NjY7pw4cKa50xMTOjWrVs520ZHRyVJyWRyQ7e5rzQ9PZ3zCvuiLZ3Bie2YTN7Pvr9/f8ay31+VzIntaHfJZLLgcywNNFNTU0okEtm5LoZhqL+/v+SB5kmTfwudGOzz+dTX11fQOefOndOZM2fy7hsaGtLNmzcLKm8tQ0NDlpaH8qEtncFJ7TgxI2X+NPz85z/Xjls/K2t9NpOT2tHurl+/XvA5RQWaN998U4ZhrAovGSvXpylm8nChMj0zpmnm9NIUGmaKderUKXV0dORsGx0d1dGjR9XW1qZ9+/ZZcp3p6WkNDQ2pra1N9fX1lpSJ8qAtncGJ7fhpYkb6eOmOz69//et6+cWWMteo9JzYjnZ39erVgs8pKtD4fD719vbK7/ertbVVPp9Pzc3N2QnDfX198nq9amxs1MmTJ4u5REEyc2cSicSqYSe/31/y67e0tKilJf8PfUNDg+WrI9fX19tixWWsjbZ0Bie1Y8NcVfb99u3bHfN9rYeT2tHuGhoaCj6nqEBz7NgxBYNBDQwMyO/369VXX83uGx8fV2tra/ZzV1dXMZcoSCAQkMfjWdUjw6ReACheWjycEvZR9MJ6jY2NOn36tPbv36+enh5du3ZN0uohpv3792+ogvlkJgGv1NPTo1gslv0cj8cVCAR4bAEAFICF9WBXG54U3Nraqt7eXr3zzju6dOlS3rBhhXg8rosXLyoWi8k0TXV1dcnv9ysUCkmSQqFQdqVgSRobG9Ply5dLUhcAAFBZLLvL6dixY5KkN954Q++++27OMJQV1tPbkgk3AABga7H8WU5nz57NDkO9//77OnLkiNWXAABsgjRTaGAjJXk4ZWYYamRkJGdeCwCgsm3CShtASZR0peDTp0/rypUrpbwEAABAaXpoVip01V0AAIBClTzQrFyTBgBgH8yhgZ0UNeS0kYeVsQojAFQuFyvRwKYKDjTj4+Pq6+uTy+XKeWbTejU1Nam3t7fg8wAAAB6n4EDT2tqq8+fPl6IuAAAARSn5HBoAgD0xhQZ2UtSQU3d396pnNq1HOp1Wc3Oz3nrrrYLPBQCUHuvQwK6KGnIaHBwsRV0AAACKwpATAACwvYIDzd27d/Xxxx+XoCoAAADFKaqH5s/+7M905MgR/emf/umG1qQBAFSuYpbmAMql4Dk0jY2NOnv2rCTpo48+UigU0p07d3T8+HG9+uqrllcQAABgLRuaQ7N//36dP39eFy9eVDqd1muvvabvfOc7DEkBAIBNZdmk4GPHjmlwcFBnz55VNBrNDkldu3bNqksAAADkZfldTo2NjTp9+rS+//3v61vf+pbOnj2r48eP6+2337b6UgCAEmIGDeykpLdtrxySampq0muvvaaenp5SXhIAsAEsrAe7Kupp28U4duyYjh07tlmXAwAAWwgL6wEAANuztIemp6dHzc3NCgQCeuWVV6wsGgCw2ZhEAxuxNNBEo1ENDw9bWSQAYBMV8+BhoBJYOuR0/PhxSWL1YAAAsKksDTQej0fS0nLZJ0+e1JEjR/T+++9beQkAAIBVLA00ma7KxsZGnT9/XsFgkLk0AGBTaSbRwEYsDTShUCi7iN61a9fU1NSU3ceKwQBQ+ZhBA7uydFLwa6+9pgMHDui9995TKBTS3bt3FY1GdfjwYY2MjOitt96y8nIAAACSLA40fX19amxs1IkTJyRJ4+PjisVieu+99/TOO+8QaAAAQElYOuTU2NiY87m1tVUnTpzQ4OCgzp8/b+WlAAAllmYKDWxk01YKzvTaAAAqF8vQwK549AEAALC9oubQbGThvF27dhV97maIx+OKxWLy+XwyDEOhUKjcVQIAAGsoONCMj4+rr69PLpdL6SIGWJuamtTb21vweZvBNE11dHRobGxMkjQwMKCBgQF1dnaWuWYAsPmYQgM7KTjQtLa2OnaCb6ZnJuPgwYM6ceIEgQYAgApn6W3bpWaapgYHB9Xf36+RkZFV+8PhcPbxC2NjY+rp6cl+Xg+fz6dEIpH9nEgkZBjGRqsNALbhYmk92FRRQ07d3d1FPZE1nU6rubm5qPVoYrGYDMOQaZoyTXPV/nA4LEnZ3pR4PK6Ojg5Fo9F1XyMQCOjgwYOKxWI6ePCgotGovF5vwXUFAACbq6ghp8HBwVLU5YmCwaAkKRKJ5N3f29ub02sTCAQ0PDwswzDk8/kUiUQeG248Ho/6+vokSf39/YrH4zIMQ4cOHaKHBgAAG7DVkNPjxONxmaa5qjfF6/UqEokoFAqpvb1d7e3ta5bl9/uzk4J7e3uzQQcAthoW1oOdOCLQZOa9PDpfxuPxaHJysqCyurq6FIlEZBiGurq6ciYJP87ExIRu3bqVs210dFSSlEwmN3Sb+0rT09M5r7Av2tIZnNiOyeRs9v2DBzOW/f6qZE5sR7tLJpMFn+OIQJNvTs169uVTzLoz586d05kzZ/LuGxoa0s2bNwsu80mGhoYsLQ/lQ1s6g5Pa8e6clPnT8Mtf/lIfmL8oa302k5Pa0e6uX79e8DklCzRvvPGGLl++rCtXrpTqElmZnhnTNHN6aQoNM8U6deqUOjo6craNjo7q6NGjamtr0759+yy5zvT0tIaGhtTW1qb6+npLykR50JbO4MR2nLg3K40s/WF//vnn9fL+L5W5RqXnxHa0u6tXrxZ8TskCjd/vV09PT6mKz5GZO5NIJFYNO/n9/pJfv6WlRS0tLXn3NTQ0WL46cn19fcWvuIz1oS2dwUntOKMH2ffb6uoc832th5Pa0e4aGhoKPqdkz3Ly+Xy6dOlSqYrPEQgE5PF4VvXIGIaRvTsKALA2VqGBXZUs0Fy6dElnz57Vc889p56eHn388ceWlb1y8buMnp4exWKx7Od4PK5AIKBAIGDZdQEAQGUqWaA5cOCAotGoPvnkEx08eFD9/f0bKi8ej6u7u1u9vb0yTVNdXV3ZxfSkh5N5M89funjxoi5fvryhawIAAHso2RyaEydO6KOPPtLU1JSOHTumY8eObai89fS28GRsALAO69DATkrWQyNJ+/fvL8uqwgCAIjGJBjZlaQ/N1NRU9oGOpmnKMAz19/fr9ddft/IyAAAAOYoKNG+++aYMw1gVXjLSK/opi3mIJQAAQCGKCjQ+n0+9vb3y+/1qbW2Vz+dTc3Nz9sGVfX198nq9amxs1MmTJ62uMwBgEzCFBnZSVKA5duyYgsGgBgYG5Pf79eqrr2b3jY+Pq7W1Nfu5q6tr47UEAGwKF5NoYFNFTwpubGzU6dOntX//fvX09OjatWuSVg8x7d+/f0MVBAAAWMuGJwW3traqt7dX77zzji5dupR30TsAAIBSsuy27WPHjun06dNKp9N69913rSoWAFAuLEQDG7F8HZqzZ89mh6Hef/99HTlyxOpLAABKhBtTYVclWVgvMww1MjKS83wlAACAUijZow8k6fTp07py5UopLwEAAFDaRx9IUl9fX6kvAQAAtriSB5qVa9IAAOyDKcGwk5IHGgCAfTAnGHZFoAEAALZHoAEAALZHoAEA5MW6erATywPNxx9/bHWRAIBN8ujz+AC7sDzQ9Pb2Wl0kAADAE1keaNL0UQIAgE1meaChuxIAnIH/oMJOmBQMAMjiv6SwKwINAACwPQINAACwPQINACAvZtDATgg0AIAs7uuAXXHbNgAAsD3LA83x48etLhIAAOCJLA80x44ds7pIAEAZ0OEOO2EODQAgy8VKNLApAg0AALC96nJXoFzi8bi6u7sVjUZztsdiMRmGobGxMZmmqf7+/jLVEAAArNeW7KGJxWJKJBIyDCNnu2maikaj6uzsVF9fnySpq6urHFUEAAAFKHmgefvtt0t9iYIFg0EFg8FV2xOJhGKxWPbz4cOHNTw8vJlVA4CKwZxg2ImlQ05TU1PZng/TNGUYhvr7+/X666+vuwzTNDU4OKj+/n6NjIys2h8Oh+XxeCRJY2Nj6unpyX7eKJ/Pl3PNK1eu6ODBg5aUDQC2wJxg2FRRgebNN9+UYRirwkvGysX1XAUsO5mZv2KapkzTXLU/HA5Lkjo7OyUtzYPp6OhYNQ/GCoZhKB6Pl6RsAABgraICjc/nU29vr/x+v1pbW+Xz+dTc3KzW1lYNDg6qr69PXq9XjY2NOnny5LrLzQwDRSKRvPt7e3tzelACgYCGh4dlGIZ8Pp8ikchjA4jH48nOi1lLpmeJMAMAgD0UFWiOHTumYDCogYEB+f1+vfrqq9l94+Pjam1tzX62alJtPB6XaZryer05271eryKRiEKhkNrb29Xe3r6h6xiGoVgslg0/4XBYoVBoQ2UCgB3xKBvYSdGTghsbG3X69Gnt379fPT09unbtmqTVQ0z79+/fUAUzEomEJK2aL+PxeDQ5OVlQWbFYTN3d3TIMQ+FwWPF4XNLS/J0DBw6oq6tLLpdLLpdLY2NjltQfAOyAh1PCrjY8Kbi1tVW9vb165513dOnSpWzwsFq+OTXr2ZdP5i6nR4egPB6P7ty5U3DdJiYmdOvWrZxto6OjkqRkMqmpqamCy8xneno65xX2RVs6gxPb8d6Dhez72dlZy35/VTIntqPdJZPJgs+x7C6nzDOc3njjDb377rs5w1BWyPTMmKaZ00tTaJgphXPnzunMmTN59w0NDenmzZuWXm9oaMjS8lA+tKUzOKkdZxakzJ+GTz75RB9M/6qs9dlMTmpHu7t+/XrB51i+UvDZs2c1Pj6unp4eHT58WH19ffr+97+/4XIzc2cSicSqYSe/37/h8jfi1KlT6ujoyNk2Ojqqo0ePqq2tTfv27bPkOtPT0xoaGlJbW5vq6+stKRPlQVs6gxPb8d6DBenKX0mS9u7dq5e/+ZUy16j0nNiOdnf16tWCzynJow8yw1BvvvlmzkJ1GxEIBOTxeFb1yBiGkXeRvM3U0tKilpaWvPsaGhq0a9cuS69XX19veZkoD9rSGZzUjq7a+ez7uro6x3xf6+GkdrS7hoaGgs8p6UrBp0+fzg5FFSrfXJyenp6cgBSPxxUIBBQIBIquIwAAsL+SP5xyvWu/SEsB5eLFi4rFYjJNU11dXfL7/dnbpkOhkMLhsAYGBiQtrRR8+fLlktQbAADYR8kDzco1adaynt4W1oQBgM3BMjSwk6ICzUZu42N8EgAqVyGPqwEqScGBZnx8XH19fXK5XEWtItnU1KTe3t6CzwMAAHicggNNa2urzp8/X4q6AAAAFKWkdzkBAOwrLSbRwD6KGnLq7u4uapw1nU6rublZb731VsHnAgBKjxk0sKuihpwGBwdLURcAAICiMOQEAABsj0ADAMiLdWhgJwQaAEAWy9DArgg0AADA9gg0AADA9gg0AADA9gg0AIC8mBMMOyHQAACyXCytB5si0AAAANsj0AAAANsj0AAA8mJhPdgJgQYAkMXCerArAg0AALA9Ag0AALA9Ag0AIK80K9HARgg0AADA9gg0AADA9gg0AADA9gg0AIC8WIcGdkKgAQBksQ4N7IpAAwAAbK+63BUAAFSm2NUv9LVndsrlWnoKt8sluVwuuXM+L713L+9zuSS3S9KKbW6X5Ha5lr7cUpXLJZfLpSr3in1ul6qW97tdy+9XfnYvlV2VfU9XEnIRaAAAWS49DAofXTf1T//jcBlr82Q5gWhl6FkOR0uhKXd/dZVLNW63qqtcqq5yq8btkiud0t27bl2a+KnqamtVs2Jfdc57t6oz29zuh8dVubWtevmrpkrbqt2qW37NeV9TpboVx9RUMUhiJQINACCrttqtP9j7lP7X6O1yV2VNi6m0FiXJkgUA3frV3TsWlLN+VW5XNvDsqK1Sw7Zq1S9/NWyrUn1t5v2KbduqtbOuRp4dNWraUaOmHbVq3F6jasIRgQYAkOs//ZM2fTKR1NxCSmmllUpL6fTSusHpdFrptHK2pdJpKbNNmf0Pj0+llj6n0pnXtBZTS8ctptLZbanlz+nl/bnHL79PpbWYKWv53MX0I2Utn7u0/WFZi6m05hdTWlhMayGV0vzy64PZeU3euav6nTuVklsLiykt5Dk2e35q6f1GLabSuj+3qPtzi0pMb6ysXXXVaqqvlWdHrZp21Oiphm16Zlednm6s0zO76vSlxjo9vatOzfW1crudOVxHoHmEYRiSJJ/Pp3g8Lp/PJ4/HU95KAcAmcrtdev6ZneWuxqaZmprSBx98oJdf3q9du3at65xUKq2F1HLYWUhrdnFRs/MpzS4s6sHy69LnlB7ML2p2YfW+B8uvM/NLoSY5u6Dp5a/k7ELOtrXy09SDBU09WNCnk/efeFxNlUt7mnao9an6pa/d9fI91aCvPbNTTfW16/0nq0hbNtDE43F1d3crGo3mbO/r69PAwIA8Ho86OzvV19dXphoCACqV2+1SrdulWrmlWkmqKdm10um0HsynsuHm3oMF3bk/t/Q1Pac79+dl3l96zWy/fW9OE/cerApC84tpGbenZdxe3SW0x7tdv/OsR7/zlUb9vq9Z33i20Va9OVsy0MRiMUkPe2NWOnDggNKsJgUAqBAul0vba6u0vbZKu3duW/d5C4sp3U7O6Td3Z/TF1APdvPtAv7n7QNcmpzV+e1rXbt/X3GIqe/yNxIxuJGb0337yG0lS044a/eFzu/X3vvElvfK1FtVWV/Y8nS0ZaILBYLmrAABASVVXufVMY52eaazLu38xldavzRmN3krqZ5/f1Y8+u6uffHZXN6ceSJLu3J/Xn//o1/rzH/1a3vpadRz8ijr/0KfmhvWHqs1UcYHGNE0NDg6qv79fIyMjq/aHw+HsnJaxsTH19PRYOsfFNE1FIhFJ0pUrVywvHwCASlDldmmPd4f2eHfo5edbsts/N2f0P391S3/xyS39xa9uKzm7oMT0nPo/NPSf/+pT/bOX96rrb/kq7s6qigo0sVhMhmHINE2ZprlqfzgcliR1dnZKWpoH09HRsWoezEaEQqHs+0QioRMnTujSpUuWlQ8AQCV71rNdf9T2Vf1R21c1M7eo7//spv7fH36q4U/v6P7cot78/i/1/i8mdOGPD8pbQROJKyrQZIaCMj0kj+rt7c3ptQkEAhoeHpZhGPL5fIpEIo8NNx6PZ80JvoZhqKurK1uGz+fLO88GAICtYHttlY7uf1bf/r0v6wejk/run/9UY7emNfLpHb3W/1eKnPzr8uyojFBTUYHmSeLxuEzTlNfrzdnu9XoViUQUCoXU3t6u9vb2oq/h9Xp1+PDh7OdoNKrjx48XXR4AAE7gcrn0B889pf/2L/5Qb7zzY33v419rdCKpP7n4sf7DPz5UEY+isE2gSSQSkrRqPovH49Hk5GRBZcViMUWjURmGoXA4rGAwqEAgII/Ho0AgoIGBgeyxK4egHmdiYkK3bt3K2TY6OipJSiaTmpqaKqh+jzM9PZ3zCvuiLZ2BdnQG2rEw/9ff9Sk5M6vYLyf1wS9v6b8Mj+uV55+y9BrJZLLgc2wTaPLNqVnPvnyCwaCCwWDeIahi7oA6d+6czpw5k3ff0NCQbt68WXCZTzI0NGRpeSgf2tIZaEdnoB3X71u7pL+sqtL9RZf+3f/4uVy/XrS0/OvXrxd8jm0CTaZnxjTNnF6aQsNMKZw6dUodHR0520ZHR3X06FG1tbVp3759llxnenpaQ0NDamtrU319vSVlojxoS2egHZ2BdizOeO01Dfzghq5Pu9T6O9/UbzfvsKzsq1evFnyObQJNZu5MIpFYNezk9/vLUKOHWlpa1NLSkndfQ0PDupfSXq/6+nrLy0R50JbOQDs6A+1YmPa2Vg384IYk6cdfzOp3Wp+xrOyGhoaCz6msm8ifIDPH5dEeGcMwWCgPAIBNtnd3gxq3Lz3yIX59c59Unk/FBprMJOCVenp6so8tkJbufAoEAgoEAptZNQAAtryVDzFd66GYm6Gihpzi8bguXryoWCwm0zTV1dUlv9+fvdMoFAopHA5n70IaGxvT5cuXy1llAAC2rK94tmtI0ud3ZspdlcoKNOvpbVnPbdQAAKD0nl5+TtSt5KzS6XRZ16Op2CEnAABQ2XbWLfWLLKbSml1IrXF0aRFoAABAUXZuezjQM/Vgvow1IdAAAIAi1a8INPdnrV1cr1AEGgAAUJTqqocxYiGVLmNNCDQAAKBI1e6Hk4AXCTQAAMCOqlYEmoUUk4IBAIANreyhWVikhwYAANgQPTQAAMD2VgaadHk7aAg0AABg48qcZwg0AACgOC6V71EHjyLQAACADWPICQAA2NLKZ1Gmy5xoCDQAAKAolTPgRKABAAAWYFIwAACwp5whp/JVQyLQAACAInGXEwAAcJR0mQedCDQAAKAoK+9yKvckGgINAAAoSuUMOBFoAACABbjLCQAA2JLLVTl9NAQaAABgewQaAABgewQaAABgewQaAABgewQaAABgewQaAABgewQaAABgewQaAABgewQaAABgewQaAABge9XlroBTzc7OSpJGR0ctKzOZTOr69eu6evWqGhoaLCsXm4+2dAba0Rlox+IZn93V3K1Pl95/0ijPbJMl5Wb+dmb+lq4HgaZEbty4IUk6evRoeSsCAMAm+OP/YH2ZN27cUCAQWNexrnQ6Xe4HZDqSaZr68MMPtWfPHm3btk2S9OKLL+qnP/3pqmPXu310dFRHjx7V9773Pe3du7d0lV+nx9V7s8sr5Lz1HLvWMU/an29fvm2V1JZObce1jitmXyX/TFZKOxZ67kZ/Jp3WjlLltGU5f7fOzs7qxo0beumll+TxeNZVBwLNJnK5XMr3z73e7T/72c+yDf7CCy+UtK7r8bh6b3Z5hZy3nmPXOuZJ+/Pty7etktrSqe241nHF7Kvkn8lKacdCz93oz6TT2lGqnLaspN+t68Gk4E303e9+15LtlcLq+hVbXiHnrefYtY550v58+2hH689b77GFttVa+yq5LSulHQs9d6M/k05rR6ly2rKSfreuBz00NlJp/4tA8WhLZ6AdnYF2dAZ6aAAAgO0RaGxk9+7d+u53v6vdu3eXuyrYINrSGWhHZ6AdnYEhJwAAYHv00AAAANtjYT0HiMfjunjxog4dOiTDMBQKhcpdJRQpHo+ru7tb0Wi03FVBkWKxmAzD0NjYmEzTVH9/f7mrhCLFYjFJkmEYikajunDhwrrXRMHmY8jJAQ4cOKCRkRFJSz+A8XicUGNDmV+eXV1dGhsbK3NtUAzTNNXb26u+vj5JS20piVBjU01NTRofH5fH41FXV5c8Hk+2bVF5GHKyuXg8nvPZ5/Pxy9OmgsGggsFguauBDUgkEtlgKkmHDx/W8PBwGWuEjciEGUkaHh7W4cOHy1shPBFDThXCNE0NDg6qv78/29uyUjgczv5gjY2NqaenRx6PR8PDw/L5fNnjvF6vEonEZlUbjyi2HVFZim1Hn8+Xc/yVK1d08ODBzao28tjIz6TH48kO6R8/fpz/cFQ4Ak0FyIy5m6Yp0zRX7Q+Hw5Kkzs5OSUu9Mh0dHdl5Fo+ek68MlN5G2xGVwap2NAxD8Xic9i0jK9oyEAjI5/Opo6NDgUCAUFPJ0qgYly5dSvt8vlXbPR5PemxsLO+2aDSaDgQC2e1jY2N5y8DmKaYdV6L9KsNG2nFsbCwdCoVKXkesT7FtuXLf48pA5WAOTYWLx+MyTVNerzdnu9frVSQSWfW/BcMw1N7evplVxDqs1Y6wh/W0o2EYisVi2cmjmV4AVJa12jIWi6mjoyO7PZFIMDxc4RhyqnCZ+TCP/iB5PB5NTk5Kki5cuKDu7m4dOnRIV65cYRZ+BVpPO8ZiMUWjURmGoXA4rGAwqEAgsNlVxROs1Y6maerAgQMyTTN7h1NmOAOVZa22DAaD6urqyg5bRaNRXbp0qQw1xXoRaCrck+bDZPYFAoHsHz56ZyrTetoxc5cTgbRyrdWOHo9Hd+7c2bwKoWjr+ZlcGUYJppWPIacKl/nfAxN/7Y12dAba0TloS+ch0FS4zPhuvlux/X7/ZlcHRaIdnYF2dA7a0nkINBUuEAjI4/Gs+l+DYRjcPmgjtKMz0I7OQVs6D4GmwuT730JPT0/O6qPxeDxn3gwqD+3oDLSjc9CWzseznCpAZiXKzHOYOjs75ff7c57HxAqzlY92dAba0Tloy62FQAMAAGyPIScAAGB7BBoAAGB7BBoAAGB7BBoAAGB7BBoAAGB7BBoAAGB7BBoAAGB7BBoAAGB7BBoAAGB7BBoAAGB7BBoAAGB7BBoAAGB7BBoAAGB71eWuAABspkgkomg0qu7ubhmGoXg8rsnJSfX19ZW7agA2gB4aAFuGaZryeDw6fPiwurq6ZBiG2tvbNTAwINM0y109ABtAoAGwZcRiMQWDQV25ckWS1NnZKa/XqwsXLsjj8ZS3cgA2xJVOp9PlrgQAbKbDhw+ro6NDnZ2d5a4KAIsQaABsOS6XS2NjY/L5fOWuCgCLMOQEYEuJx+PyeDyEGcBhCDQAtpTMPBoAzkKgAbCljI2N6fjx4+WuBgCLMYcGAADYHj00AADA9gg0AADA9gg0AADA9gg0AADA9gg0AADA9gg0AADA9gg0AADA9gg0AADA9gg0AADA9gg0AADA9gg0AADA9gg0AADA9gg0AADA9gg0AADA9v43ouyyIyp25jMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 600x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mf.plot_frob_error_rank(frob_losses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a6bc7fe2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5000, 300),\n",
       " (5000, 5000),\n",
       " -0.0018394775345592505,\n",
       " 0.014836883129622934,\n",
       " 0.0002737358427138843,\n",
       " 0.00025728376347159727)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C.shape, A.shape, A.min(), A.max(), A.mean(), A.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8e62bfa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "eps_ff = 1e-2\n",
    "eps = 1e-3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "80d79e96",
   "metadata": {},
   "outputs": [],
   "source": [
    "method = 'bcd'\n",
    "svds_v0 = True\n",
    "max_iters_ff = 2\n",
    "max_iters = 3*10**3\n",
    "PSD = True\n",
    "symm = True\n",
    "\n",
    "filename = \"%s_r%d_%d\"%(mtype, rank, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49200753",
   "metadata": {},
   "source": [
    "# Rank allocation + spectral partitioning + greedy refinement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "fc600efe",
   "metadata": {},
   "outputs": [],
   "source": [
    "info = {} "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "0d3f0f76",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.16178478761767232]\n"
     ]
    }
   ],
   "source": [
    "# LR\n",
    "hpart = mf.random_hpartition(m, n, num_levels=1, symm=symm, perm=False)\n",
    "B1, C1 = mf.single_level_factor_fit(A, np.array([rank]), hpart, level=0, symm=symm, PSD=PSD)[:2]\n",
    "info['LR'] = {'loss':[mf.rel_diff(B1 @C1.T, den=A)]}\n",
    "print(info['LR']['loss'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2be59ee8",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/parshakova.tanya/Documents/mlr_fitting/mlrfit/utils.py:404: NumbaPerformanceWarning: np.dot() is faster on contiguous arrays, called on (array(float64, 2d, A), array(float64, 2d, A))\n",
      "  hat_A_except_level[r1:r2, c1:c2] += np.dot(B_level[r1:r2], C_level[c1:c2].T)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fa_losses[-1]=0.1537922868665804 6 time=55.460\n"
     ]
    }
   ],
   "source": [
    "# LR+D fitting\n",
    "hpart = mf.random_hpartition(m, n, level_list=[0, int(np.log2(m))+1], symm=symm, perm=False)\n",
    "assert hpart['rows']['lk'][1].size == min(m,n)+1\n",
    "ranks = np.array([rank-1, 1])\n",
    "hat_A = mf.MLRMatrix(hpart=hpart, ranks=ranks)\n",
    "start_time = time.time()\n",
    "fa_losses = hat_A.factor_fit(A, ranks, hat_A.hpart, method=method, eps_ff=eps_ff/10000, \\\n",
    "                                  freq=100, svds_v0=svds_v0,\\\n",
    "                                  printing=False, max_iters_ff=max_iters, symm=symm, warm_start=False)\n",
    "fa_time = time.time()-start_time\n",
    "info['LR+D'] = {'loss':fa_losses, 'time':fa_time}\n",
    "print(f\"{fa_losses[-1]=} {len(fa_losses)} time={fa_time:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1e847bb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# GICS hierarchy\n",
    "num_levels = len(gics_hpart['rows']['lk'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7ff51aa4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5 5 5 5 5 5]\n",
      "itr=0, t=3, losses[0]=0.18304971186124508, losses[-1]=0.18115165653212797, [5 5 5 5 5 5]\n",
      "itr=1, t=6, 0.17760319373154995, [6 5 5 5 5 4]\n",
      "itr=2, t=9, 0.17526203538344337, [7 5 5 5 5 3]\n",
      "itr=3, t=12, 0.17320062948232995, [8 5 5 5 5 2]\n",
      "itr=4, t=15, 0.17123781655475356, [9 5 5 5 5 1]\n",
      "itr=5, t=18, 0.16976666828550335, [10  5  5  5  4  1]\n",
      "itr=6, t=21, 0.1683443592560673, [11  5  5  5  3  1]\n",
      "itr=7, t=24, 0.16699033223213045, [12  5  5  4  3  1]\n",
      "itr=8, t=27, 0.1658319270936777, [13  5  5  4  2  1]\n",
      "itr=9, t=30, 0.16473492603709036, [14  5  5  3  2  1]\n",
      "itr=10, t=33, 0.16369949161116465, [15  5  5  3  1  1]\n",
      "itr=11, t=36, 0.16271825209289337, [16  5  5  2  1  1]\n",
      "itr=12, t=39, 0.16183771490312615, [17  5  4  2  1  1]\n",
      "itr=13, t=42, 0.1609730679688495, [18  5  3  2  1  1]\n",
      "itr=14, t=45, 0.16005469382701504, [19  5  3  1  1  1]\n",
      "itr=15, t=48, 0.1592639628654469, [20  5  2  1  1  1]\n",
      "itr=16, t=51, 0.15860377442339196, [21  4  2  1  1  1]\n",
      "itr=17, t=54, 0.15789430844667762, [22  3  2  1  1  1]\n",
      "itr=18, t=57, 0.15714543741668827, [23  3  1  1  1  1]\n",
      "itr=19, t=60, 0.156527399872881, [24  2  1  1  1  1]\n",
      "itr=20, t=63, 0.1557835304779305, [25  2  1  1  0  1]\n",
      "itr=21, t=66, 0.15523730467736668, [26  1  1  1  0  1]\n",
      "itr=22, t=69, 0.15463389148281434, [27  1  1  0  0  1]\n",
      "itr=23, t=72, 0.1541334597291659, [28  1  0  0  0  1]\n",
      "itr=24, t=75, 0.1537924062023245, [29  0  0  0  0  1]\n",
      "quit: new rank allocation is worse\n",
      "itr=25, t=78, 0.15379228615265847, [29  0  0  0  0  1]\n",
      "td+ra time=1188.407\n"
     ]
    }
   ],
   "source": [
    "# uniform rank initialization + GICS hierarchy + rank allocation\n",
    "ranks = mf.uniform_ranks(rank, num_levels) \n",
    "print(ranks)\n",
    "hat_A = mf.MLRMatrix(hpart=gics_hpart, ranks=ranks)\n",
    "start_time = time.time()\n",
    "ra_losses, epochs, ranks_history_uni = hat_A.rank_alloc(A, ranks, hat_A.hpart, freq=1,\\\n",
    "                        eps=eps, eps_ff=eps_ff, max_iters=max_iters, PSD=PSD, symm=symm,\\\n",
    "                        max_iters_ff=max_iters_ff, warm_start=False)\n",
    "ra_time = time.time()-start_time\n",
    "info['MLR uniform'] = {'loss':ra_losses, 'ra_time':ra_time, 'epochs':epochs}\n",
    "print(f\"td+ra time={ra_time:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b8ab3ca7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0  0  0  0  0 30]\n",
      "itr=0, t=3, losses[0]=0.9964034275153691, losses[-1]=0.9964034275153691, [ 0  0  0  0  0 30]\n",
      "itr=1, t=6, 0.25789623933139616, [ 1  0  0  0  0 29]\n",
      "itr=2, t=9, 0.22336450751609954, [ 2  0  0  0  0 28]\n",
      "itr=3, t=12, 0.20445675337736138, [ 3  0  0  0  0 27]\n",
      "itr=4, t=15, 0.1959904016189144, [ 4  0  0  0  0 26]\n",
      "itr=5, t=18, 0.1904434327633549, [ 5  0  0  0  0 25]\n",
      "itr=6, t=21, 0.18578546720178915, [ 6  0  0  0  0 24]\n",
      "itr=7, t=24, 0.18289954291436586, [ 7  0  0  0  0 23]\n",
      "itr=8, t=27, 0.18053266891609757, [ 8  0  0  0  0 22]\n",
      "itr=9, t=30, 0.17831097103325927, [ 9  0  0  0  0 21]\n",
      "itr=10, t=33, 0.1763781102590341, [10  0  0  0  0 20]\n",
      "itr=11, t=36, 0.17466271344630813, [11  0  0  0  0 19]\n",
      "itr=12, t=39, 0.17297096622362212, [12  0  0  0  0 18]\n",
      "itr=13, t=42, 0.17158598680520312, [13  0  0  0  0 17]\n",
      "itr=14, t=45, 0.17027252180523023, [14  0  0  0  0 16]\n",
      "itr=15, t=48, 0.1689720247460934, [15  0  0  0  0 15]\n",
      "itr=16, t=51, 0.1677186255495341, [16  0  0  0  0 14]\n",
      "itr=17, t=54, 0.16647850615838616, [17  0  0  0  0 13]\n",
      "itr=18, t=57, 0.1653154334415647, [18  0  0  0  0 12]\n",
      "itr=19, t=60, 0.1641668753724694, [19  0  0  0  0 11]\n",
      "itr=20, t=63, 0.16302524230347018, [20  0  0  0  0 10]\n",
      "itr=21, t=66, 0.16194041223166816, [21  0  0  0  0  9]\n",
      "itr=22, t=69, 0.16086034070097954, [22  0  0  0  0  8]\n",
      "itr=23, t=72, 0.15979156242581022, [23  0  0  0  0  7]\n",
      "itr=24, t=75, 0.15872549835483032, [24  0  0  0  0  6]\n",
      "itr=25, t=78, 0.15771415700360494, [25  0  0  0  0  5]\n",
      "itr=26, t=81, 0.1567084990408996, [26  0  0  0  0  4]\n",
      "itr=27, t=84, 0.15571110814277317, [27  0  0  0  0  3]\n",
      "itr=28, t=87, 0.15475142386771168, [28  0  0  0  0  2]\n",
      "itr=29, t=90, 0.15379228799036154, [29  0  0  0  0  1]\n",
      "quit: new rank allocation is worse\n",
      "itr=30, t=93, 0.15379228605054865, [29  0  0  0  0  1]\n",
      "td+ra time=1235.205\n"
     ]
    }
   ],
   "source": [
    "# bottom rank initialization + GICS hierarchy + rank allocation\n",
    "ranks = np.array((num_levels-1)*[0]+[rank]) \n",
    "print(ranks)\n",
    "hat_A = mf.MLRMatrix(hpart=gics_hpart, ranks=ranks)\n",
    "start_time = time.time()\n",
    "ra_losses, epochs, ranks_history_bottom = hat_A.rank_alloc(A, ranks, hat_A.hpart, freq=1,\\\n",
    "                        eps=eps, eps_ff=eps_ff, max_iters=max_iters, PSD=PSD, symm=symm,\\\n",
    "                        max_iters_ff=max_iters_ff, warm_start=False)\n",
    "ra_time = time.time()-start_time\n",
    "info['MLR bottom'] = {'loss':ra_losses, 'ra_time':ra_time, 'epochs':epochs}\n",
    "print(f\"td+ra time={ra_time:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "0a2f723a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[30  0  0  0  0  0]\n",
      "itr=0, t=3, losses[0]=0.16178478761767176, losses[-1]=0.16178478761767173, [30  0  0  0  0  0]\n",
      "itr=1, t=6, 0.15379235716161116, [29  0  0  0  0  1]\n",
      "quit: new rank allocation is worse\n",
      "itr=2, t=9, 0.15379228616060095, [29  0  0  0  0  1]\n",
      "td+ra time=226.239\n"
     ]
    }
   ],
   "source": [
    "# top rank initialization + GICS hierarchy + rank allocation\n",
    "ranks = np.array([rank]+(num_levels-1)*[0])\n",
    "print(ranks)\n",
    "hat_A = mf.MLRMatrix(hpart=gics_hpart, ranks=ranks)\n",
    "start_time = time.time()\n",
    "ra_losses, epochs, ranks_history_top = hat_A.rank_alloc(A, ranks, hat_A.hpart, freq=1,\\\n",
    "                        eps=eps, eps_ff=eps_ff, max_iters=max_iters, PSD=PSD, symm=symm,\\\n",
    "                        max_iters_ff=max_iters_ff, warm_start=False)\n",
    "ra_time = time.time()-start_time\n",
    "info['MLR top'] = {'loss':ra_losses, 'ra_time':ra_time, 'epochs':epochs}\n",
    "print(f\"td+ra time={ra_time:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "8fd5fe6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# mf.plot_loss_all_info(info, figsize=(7, 5), dpi=120, logscale=False, hline=True, ylim=[0.14, 0.3], filename='')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f13cc3fe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAHZCAYAAAB0C6brAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAABJ0AAASdAHeZh94AABuMklEQVR4nO3de3zT9b0/8Nc3Sdu0aZNQ7tIgN61tCoKCk1R3s6OpbnPWI8HdpBtlO9sZOTuW/bb9RhzFs4utv63s7AJhKzu7ETzrzm42rfU2R6qAwqShsim6fkFABZtvb/SSfn9/lISm1zTNva+nDx6ab/L9fj/1i/jqp+/P+yPIsiyDiIiIiIj8FLEeABERERFRvGFIJiIiIiIagSGZiIiIiGgEhmQiIiIiohEYkomIiIiIRmBIJiIiIiIagSGZiIiIiGgEVawHQERERMlHlmV0dXVBkiT09vaC2zJQLAmCgLS0NGi1Wmg0GgiCMPk53EyEiIiIwkmWZbz11lu4dOkSACAlJQUKBX94TbEzODiI/v5+AEB2djbmzZs3aVDmTDIRERGFVVdXFy5duoSMjAwsXLgQqampsR4SEfr6+nDu3DlcunQJGo0GmZmZE36e39YRERFRWEmSBAAMyBRXUlNTsXDhQgBXf49OhCGZiIiIwqq3txcpKSkMyBR3UlNTkZKSgt7e3kk/y5BMREREYSXLMmuQKW4JghDUQlL+DiYiIiKiGSOYzhYAQzIRERER0SgMyURERETTIIoicnNzkZubO+GCsPE+53K5kJubC6fTOeF93G63//zhv0pLS+FwOIIeb1VVFYqKioL+/EzFkExEREQUBlqtdsKweuDAAWi12mnfp6amBqdOncKpU6dw5MgRbN26FXa7HaWlpdO+Nl3FkExEREQUBiUlJROG5IMHD6KkpGTa9xketLVaLcxmM+rq6uB2u2G326d9/alyOp3jzqBP9F68Y0gmIiIiCgOz2QxRFCGK4qj3XC4XJEmCyWSKyL21Wi2MRiNOnDgRkeuPR5IkWK1WeDyeKb2XCBiSiYiIiMJAp9PBZDLhwIEDo947cOAALBZLWMotxiJJEtxuN1auXBmR689EDMlEREREYWKxWHDw4MFRxxsaGmCxWMJyD1EUIUkSJEmCKIpwOBwoLS2F0WhEeXl50NdxuVwoLS31L/5zu92jPuN0OgM+M7ycZNu2bVi3bh0AoKioyP+Zyd4L5toAYLfbUVZWBpfLFXANURThdrtRVlYWcCzcVGG/IhEREdE4Lv3+t+h8oTnWwxgl8z3rkX33vdO+jtlshtVqhcvl8pdWOBwOGAwGGI1GuFyuad/DZrPBZrP5X2u1WmzdunVKAVkURTidTuzatQsejwfV1dUoLS3FkSNH/LPddrsde/fuRUVFBUwmE9xuN3bs2IG2tjZs374du3fvhiiKKCoqQl1dHYxGo//6E70XzLUBoL29HS0tLbDb7aipqYHH44HVakVZWRkMBgMqKir8x6xWK+rq6qbzr3UUhmQiIiKiMLJYLDhw4IA/JDudThQXF4ft+rW1tf5rV1VVYd++fVOepdZqtaisrPS/3r9/P9atWweHw4Hy8nJIkoTq6mrU1NTAbDYDAAwGAwwGA0pLS7Fp0yYYDAbodDr/9UYa771gr+37bGVlpf+1b3FkTU2NP3hv3LgR+/btm9LXH4yEDcl2ux3t7e3Q6/UQRREWi2XUdykTkSQJe/bsQUdHh3+Kvry8fMKCet89fQoLCyNWgE9ERJSMsu++NywztvHMYrGgtLTU39XB5XIFBNJw2r59Ow4ePIjq6upp3WPkwj/fjLcvxPr4spbL5Qq5fGSq1/YF5OH/PDzzLV68OKRxTCYhQ/K2bdtQWFgY8GOFoqKigO8qJuL7Dmb4byan04mysjJUVFSM+nGFb3WmxWLxv+d0OmG1WnHkyJEwfVVERESUDIxGIwwGA+rr6wNeR8quXbv8OWUqE4bBkCRpzFnicLR1C+bakVroGIyEW7jncrnQ3Nw86ruX8vJyWK3WoK6xZ8+eUd9tmc1mFBcXo7q6elTx9+bNm2EwGAK+4xFF0f9jBCIiIqLhiouL4XQ64XA4wrZgbzxmsxlGoxE7duwI+Roju2P4wrYv6Pv4ZoFH/iR9otA88r2pXjtWEi4kO51OFBQUjDpeUFDgX+05mYaGhjF3pSksLASAgNWdLpcLbrcbFRUVAZ8tLy9HU1PTVIdPREREM8CmTZv8GSLYkOzxePx9ln2/xuo4MZZdu3bB7XZPurW1jyRJsNlsEEXR3+VCq9X6f2JuMBiwZcsW2Gw2OBwOuN1uOByOUTPWvple32eGb2Yy3nvBXjvWEq7cor6+fszdanw/xqivr5/0N6PBYBizsbVvZnj4e3a7HQaDIabT/URERJRYfEFQr9cHfc7wjhXDDe84MR6j0Yji4mLs2LFjVK3vSIsXL4bJZILRaITVaoUoiigoKBjVHWL79u1YuXIl9u7dC5vNBqPRiIqKilE5y2KxwOFwwOVyjVqgON57wV47lgRZluVYD2IqcnNzsWXLFn97kGDfC4ZvhejwViXr1q1DQUEBKisr4XQ6odVq0dbWxkV7RERE4zh9+jQAYNmyZTEeCdFowf7+TLiZ5Ml0dHSEfG5DQwOMRmPANL8kScjKyoLL5QpY0FdaWgqLxRLydzyDg4MYGBgAACgUCgiCEPK4iYiI4olv/i3B5uEogfh+bwmCELEMlVAhOZh64+Et2qbCbrfD4/GgtrZ21P0aGhrw8MMPB3x+69atsFqtKCkpCakUY2BgIOr7qxMREUVLVlYWuru7Yz0MSnLp6ekRC8kJt3APmHi2eCq1Pz6iKGLv3r3Yv3//mC1ajEbjqCA83spMIiIiIkp8CTWT7AuqY80W+2Z9s7Kypnxdq9U6Zo9l3/3GavXmOxbsqtORFIqr35+sXLkSSqUypOsEq7W1FQCQl5cX8jWk557BO7/cj4XWB5GevzJMI6NQheOZUnzhM00+M/WZvvHGGwCAjIyM2A4kQi5fvgwAUKvVMR7JzCTLMnp6egAgouWqCRWSgfGbSvs6Ukx11xWbzYatW7eOuwhvvE4YYx2biuEPValURjwk++43nfsoBAHCQD8ECBEfL00uHM+U4gufafKZqc/U93Un+3qbZP/6ZrqEK7coKSkZs9zCN5M8lY4TdrsdRqNxVKuU4T0GTSbThLXQ8dLLj4iIiIjCJ+FCstls9u/IMpzL5YLBYAh620ffNUZ2p3C5XAGz1RaLZcxNSnyfG6tnc7Lyf7/M1cpERESU5BKu3MJkMsFkMo3a5tHhcIzaahoAioqKYDAYArpWiKIIm82G4uJiVFVVBXy+ubkZNTU1/tdGoxEWiwXf+MY3sHv3bgBDs9Z2ux27du2aUZuM9La9AQDoecUN9fLrgjpHoVZDmGE/ZiQiIqLEl3AhGQBqa2tht9vhcDgADC2eq6ysHLPUQqvVjlrMV1ZWBlEUsW/fvjGvP3I2urKyEg6HA9u2bYNer0d7e/uYC/2SWdexo/A82QgAaH/8j2h//I9BnZcyfwEWf+d7rNsiIiKihJKQIRlAwMYeExm5xSIANDU1Tfl+09k4JBl0twZ28UhdZIB6xcSzyT2nWtF//hzkvl4IaVwBTERERIkjYUMyRVdGnhGexqs9oWf/iwWaNWsnPOftX/wMnvPnMNhzGQqGZCIiIkogDMkUFM2atVhorUD3KyeRcUP+pAEZAAR1OgBg8HIPAH1kB0hEREQURgnX3YJiR7NmLebe/+mgAjIA/+zx4JWm60RERMlIFEXk5uYiNzd3wrax433O5XIhNzc3oAXtWNxut//84b9KS0v967Tigc1mw7p167Bu3boxO5IlCoZkihhF+lBIli/3xHgkREREkafVaicMqwcOHAhLV6yamhqcOnUKp06dwpEjR7B161bY7XaUlpZO+9pTsW3bNuTm5o46Jooi9u/fj5qamqBb88YjhmSKGIW/3IIzyURElPxKSkomDMkHDx4My/4Kw4O2VquF2WxGXV0d3G437Hb7tK8frMLCwlFNDRoaGlBRUQGj0QiTycSQTDQWhZrlFkRENHOYzWaIoghRFEe953K5IEnSlHYGngqtVguj0YgTJ05E5PpjsVgsAXtU+L7uZNlDgiGZIiZw4R4REVFy0+l0MJlMOHDgwKj3Dhw4AIvFErEAKUkS3G43Vq5cGZHrz0TsbkERw5lkIiIa6beXfo/mzhdiPYxR1me+B/dm3z3t61gsFuzYsQPbt28PON7Q0IC6ujp4PJ5p30MURf/CP4/HA5fLBbvdDqPRGNQ+ElVVVTh48CCOHDniP2a327F3717/MbvdDpfLhfLyclRXV8PtdsNoNAbUGVdVVaGhoQFNTU2oqqryb9JWVFQEYGivCt/Ga06nE3v37vVfZ+T+E7771dTU4Bvf+Aaam5tx5MiRgHHYbDaIooj8/Hx85zvfQUdHB374wx/C5XKNGls4cCaZIsZXk8yFe0RENFOYzWZIkhTQ1cHhcMBgMIRtp97h3SOKiopQXV0Ni8Uy5gZqoWpvb/eH7127dqG2thaiKMJms435+e3bt/vvX1tbi1OnTvm/Xrvdjh07dsBisaCpqQlbt25FdXU1qqqqAu4niiJKS0uh1+tRU1PjP97S0gK73Y6amhrU1tbizJkz+MIXvoDdu3fjwQcf9I/NarWG7esHOJNMEXR1JpkhmYiIhtybfXdYZmzjmcViwYEDB/z1x06nE8XFxWG7fm1trf/avhncSO0KXFlZ6Z+dLSkpQX19/bif9ZWS6HQ6/zFJklBdXY2amhqYzWYAgMFggMFgQGlpKTZt2uS/viiKqKioGDUbLklSwDjMZjMOHjyI7373uzAajRAEARs3bvTPZIcLZ5IpYhTp7G5BREQzj8ViQUNDAyRJ8s8qb9q0KSL32r59O7RaLaqrqyNy/eHlC6GUMvhm1H0B2cc3yzyyj/J45SLD77148WIAQF5e3qhj4cSQTBGj4MI9IiKagYxGIwwGA+rr61FfX+9/HSm7du2Cw+GA2+0O63XDuchwvE1Whh8f799RrLplMCRTxAgqFaBSYbCHM8lERDSzFBcXw+l0wuFwRKwUwsdsNsNoNGLHjh0hX2OstnXh4JsxHlmm4ZtBjlRLvHBgTTJFlCJNzZlkIiKacTZt2uTv8hBsSPZ4PKPCqiRJQS3427VrF0pLS+F0OkeVNoyk1+v9LeOMRiNcLhccDkdEZmwNBgO2bNniX/BXUFCAlpYW/2LDcC1mjASGZIoohVoNmTXJREQ0w/jCoV6vD/qc8TpHHDlyZNIAazQaUVxcjB07dkwaki0WC+rr61FaWgqDwYDi4mJUVlZOuFvgdGzfvh0rV67E3r17YbPZYDQaUVFREfEZ9ukSZFmWYz2Imcjr9eL48eMAgNWrV0OpVEb0fr46pWh/x9b2je2Q+/tx7Xe/H9X7zgSxeqYUOXymyWemPtPTp08DAJYtWxbjkURGT8/QT0jTryxQp+iSZRnd3d0AgIyMDAiCMKXzg/39yZpkiiiFWo3B3t5YD4OIiIhoShiSKaIEdTprkomIiCjhMCRTRPlqkuXBwVgPhYiIiChoDMkUUf6tqVlyQURERAmEIZkiihuKEBERUSJiSKaIUqjVALg1NRERESUWhmSKKM4kExERUSJiSKaIEjiTTERERAmIIZkiylduIXMmmYiIiBIIQzJFlL/cooczyURERJQ4GJIpoq4u3ONMMhERESUOhmSKKP9Mci9nkomIKDmJoojc3Fzk5uZCkqQpf87lciE3NxdOp3PC+7jdbv/5w3+VlpbC4XCE7euhIQzJFFFCOhfuERHRzKDVaicMqwcOHIBWq532fWpqanDq1CmcOnUKR44cwdatW2G321FaWjrta9NVDMkUUWwBR0REM0VJScmEIfngwYMoKSmZ9n2GB22tVguz2Yy6ujq43W7Y7fagruF0Oiec9SaGZIowf3cLLtwjIqIkZzabIYoiRFEc9Z7L5YIkSTCZTBG5t1arhdFoxIkTJyb9rCRJsFqt8Hg8ERlLsmBIpojiTDIREc0UOp0OJpMJBw4cGPXegQMHYLFYwlJuMRZJkuB2u7Fy5cqIXH8mYkimiBKUSggpKaxJJiKiGcFiseDgwYOjjjc0NMBisYTlHqIoQpIkSJIEURThcDhQWloKo9GI8vLyCc/dtm0b1q1bBwAoKiryL/wbzul0orS0dNxFgQ6HA9u2bYPL5UJZWRnWrVuH0tJSuN3usHx98UIV6wFQ8lOo0zmTTEREAIDf//YSXmjujPUwRnnP+kzcfW/2tK9jNpthtVrhcrn8pRUOhwMGgwFGoxEul2va97DZbLDZbP7XWq0WW7dunTQgA8Du3bshiiKKiopQV1cHo9EY8L7dbsfevXtRUVEBk8kEt9uNHTt2oK2tDdu3bwcAtLW1oaGhAQBQXl4OSZJQXV2N0tJSHDlyJGKz5dHGmWSKOEGt5kwyERHNGBaLJaDkwul0ori4OGzXr62t9Xe32LJlCyRJmtIstU6nA4BRYdYXdnft2gWLxQKDwQCz2Yz9+/dj3759AbXWWq0Wu3fvhslkgtlsRm1tLQAkVSs6ziRTxCnU6dyWmoiIAAB335sdlhnbeGaxWFBaWurvHuFyuVBZWRmRe23fvh0HDx5EdXX1tO/hm+U2m80Bx32zzS6Xa9ww7pspD2bhYKLgTDJFnIIzyURENIMYjUYYDAbU19ejvr7e/zpSdu3aBYfDEbaa4PFaw820lnEMyRRxQyGZM8lERDRzFBcXw+l0wuFwhG3B3njMZjOMRiN27NgxpfNGhl7fjHF9fX3Acd8M8/D2dSPPFUUx6bprMCRTxCnU6ZD7+iAPDsZ6KERERFGxadMmuFwuuN3uoEOyx+Px91n2/Qp2dnjXrl1wu92Tbm0NXK1F9s0++zYgMRgM2LJlC2w2m/89h8MBq9UKi8UyapGfzWaDKIr+LhdarTaoxYOJgjXJFHGC+urW1MqMjBiPhoiIKPJ8gVOv1wd9zvCOFcMF0zHCaDSiuLgYO3bsGFVTPBaLxQKHwwGXyxWwqHD79u1YuXIl9u7dC5vNBqPRiIqKilFB37d5idVqhdvthslkQl1dXRBfZeIQZFmWYz2Imcjr9eL48eMAgNWrV0OpVEb0fr7vREd+FxgNb//q5/A8UY8l/++HUGXPjvr9k1UsnylFBp9p8pmpz/T06dMAgGXLlsV4JJHR0zNUQpienh7jkcRGVVUVDh48iCNHjsTk/rIso7u7GwCQkZEBQRCmdH6wvz9ZbkERpxg2k0xERESUCBiSKeIU6dyamoiIiBILQzJF3NWZZIZkIiIiSgwMyRRxCvWVmeQellsQEREluu3bt8esHjmaGJIp4gQ1yy2IiIgosTAkU8T5yi1kLtwjIiKiGAu2sRtDMkUca5KJiGYWQRAwyA2kKE7JshxU2ziGZIo4f00yZ5KJiGaEtLQ09Pf3o6+vL9ZDIQrQ19eH/v5+pKWlTfpZhmSKOPZJJiKaWXy7w507d45BmeJGX18fzp07BwCT7mAIcFtqigIu3CMimlk0Gg2ys7Nx6dIlvPbaa0hJSYEgCFPeGS1eeb1eAIj4brk0NlmWA57BZL+vZFmGLMvo7+8HAGRnZ0Oj0Ux6H4Zkijgu3CMimlkEQcC8efOg0WggSRJ6e3uDXiyVCDo7OwEAOp0uxiOZuTo6OgAAer1+0s8KggCFQoGMjAxotVpoNJqgvmFL2JBst9vR3t4OvV4PURRhsVhgNBqDPl+SJOzZswcdHR0QRREAUF5eDpPJNOm5oiiiqqoKu3fvDnn8M4mgUEBITeNMMhHRDCIIAjIzM5GZmRnroYSd2+0GACxbtizGI5mZvF4vjh8/DgBYsmRJxGb0EzIkb9u2DYWFhSgvL/cfKyoqQk1NTVBBWZIkVFdXo7Ky0n/M6XSirKwMFRUVAdcdi81mg8fjCf0LmIEUajVrkomIiChhJNzCPZfLhebmZlgsloDj5eXlsFqtQV1jz549AQEZAMxmM4qLi1FdXe2fWR6L0+mc8H0amyI9nTPJRERElDASLiQ7nU4UFBSMOl5QUABRFCFJ0qTXaGhoQGlp6ajjhYWFAK7+GGUkSZIgiiLy8/OnOGriTDIRERElkoQLyfX19TAYDKOO+47V19dPeo2xzgeuFuCPV0qxZ8+eSUsxaGyCmjPJRERElDgSriZZkiRkZWWNOu7rd9fW1jbpNWpra8c8fuLECQAYc6ba5XL5Z5pp6hRqNeQeziQTERFRYki4kDwZX0uQUDQ0NMBoNI65+M/tdkdsFrm1tTXivSN7e3sBjF9KEmmq3l4oB/rhfvllgH0lwyLWz5TCj880+fCZJic+19iKVjvBhArJwdQbt7e3h3Rtu90Oj8cz5iyz3W5nmcV0pV7Z/rGvD0hPj+1YiIiIiCaRUCHZZ6LZ4mCaSo8kiiL27t2L/fv3j6pXdrvdU+q/HIq8vLyI79rj+2430l/LeN5++UV4Xj6G65cuQcrsOTEZQ7KJ9TOl8OMzTT58psmJzzW2hvdJjqSEWrjnqzsea7bYN8s8Vr3yZKxW67g9lh9//PGgNhihiSm4NTURERElkISbSfYF5ZF8HSkWL148pevZbDZs3bp1zCDscDhw8uRJ2Gy2gOMnT56Ex+PxHx/Zc5lG829NzcV7RERElAASLiSXlJSMuZmHbyZ5KrO+drsdRqMRZrM54LjT6YTZbIbFYhm1aQkAlJWVQavVMhxPAWeSiYiIKJEkVLkFMLQznsvlGnXc5XLBYDCM2wN5rM8DGBWCXS7XuLPVFDrfTDI3FCEiIqJEkHAzySaTCSaTCQ6HIyDgOhyOMWd2i4qKYDAYArpWiKIIm82G4uJiVFVVBXy+ubkZNTU1E46B21JPHWeSiYiIKJEkXEgGhjYDsdvtcDgcAIZWmVZWVo5ZaqHVakct5isrK4Moiti3b9+Y1x9vNrqqqgrNzc3+kFxaWor169dj+/bt0/lyZgTBP5PMkExERETxLyFDMoCg+xbX1dWNOtbU1BTSPRmGQ6dI980ks9yCiIiI4l/C1SRTYvJ3t+BMMhERESUAhmSKCkXalXILtoAjIiKiBMCQTFHBhXtERESUSBiSKSqEtDQArEkmIiKixMCQTFEhKBQQ1GoM9jIkExERUfxjSKaoUajTuXCPiIiIEgJDMkWNQq1muQURERElBIZkihqFOh2DPZxJJiIiovjHkExRw5lkIiIiShQMyRQ1gjqdLeCIiIgoITAkU9Qo1GrA64Xc3x/roRARERFNiCGZooYbihAREVGiYEimqFGouaEIERERJQaGZIoaziQTERFRomBIpqgR/CGZM8lEREQU3xiSKWoUajUAhmQiIiKKfwzJFDW+cgtuTU1ERETxjiGZooYzyURERJQoGJIpahTpXLhHREREiYEhmaLGP5Pcw5lkIiIiim8MyRQ1AlvAERERUYJgSKao8c0ky6xJJiIiojjHkExRw81EiIiIKFEwJFPUCKmpgCCwuwURERHFPYZkihpBEKBQqzmTTERERHGPIZmiSlCnMyQTERFR3GNIpqhSqNWQL/fGehhEREREE2JIpqhScCaZiIiIEgBDMkXVUE0yF+4RERFRfGNIpqjyzSTLshzroRARERGNiyGZokqhVgODg5D7+2M9FCIiIqJxMSRTVAnp3FCEiIiI4h9DMkXV1a2pGZKJiIgofjEkU1T5t6bu4eI9IiIiil8MyRRVvplkllsQERFRPGNIpqhSpPlCMmeSiYiIKH4xJFNUceEeERERJQKGZIqqqwv3OJNMRERE8YshmaLKv3CPIZmIiIjimCoaN2lsbERWVhYMBgNycnKicUuKU1dDMsstiIiIKH5FZSa5qqoKACBJUjRuR3HsancLziQTERFR/IrKTHJxcTHWr18fjVtRnONMMhERESWCqIRkvV7v/+d9+/ZBEASYTCbk5eVF4/YURwTOJBMREVECiEq5hU6n8//zli1b8PLLL6OlpQVnzpyJxu0pjggpKYBSyW2piYiIKK5FJSS3tLQEvF65ciXuu+8+LuKbgQRBgEKt5kwyERERxbWohGSHw4G8vDzce++9ePTRR9HS0oLOzk7/+83NzdEYBsULhQJ9Z0V0HTsa65EQERERjSkqNckVFRUoLi5Gc3MzDh06hObmZqxbtw75+fkwmUwQRZEL+2aIrmNHMdjRAQA4V1MNITUVglI58UkKBeZ8/AFoC98bhRESERERRSkkb9myBQBgMBiwceNGAIAoiv7Q/Pzzz0djGBQHulvdAa8VWVlIXXDNhOf0uE+gp+VlhmQiIiKKmqiE5LEYDAZ/aN63b1+shkFRlpFnhKex3v963ifLoFmzdsJzXvt8GQYkT6SHRkREROQXs5A8nMViifUQKEo0a9ZiobUC3a+cRMYN+ZMGZABQ6XTwehiSiYiIKHriIiRnZWXFeggURZo1a4MKxz5KnQ79585FcEREREREgeIiJIfCbrejvb0der0eoijCYrHAaDQGfb4kSdizZw86OjogiiIAoLy8HCaTaVqfpfBTanW4/I+/Q/Z6J1/kR0RERBQG0w7Jra2taG9vD8NQhnbmC2YXvm3btqGwsBDl5eX+Y0VFRaipqQkqKEuShOrqalRWVvqPOZ1OlJWVoaKiIuC6U/ksRYZKpwdkGd4OCSr9rFgPh4iIiGaAaYfkEydOQJKkcIwFOp1u0pDscrnQ3NyM3bt3BxwvLy+H1WpFU1PTpPfZs2dPQOgFALPZjOLiYlRXV8NsNsNgMEz5sxQZSu3Qjo1ej4chmYiIiKJi2iHZ19ItWpxOJwoKCkYdLygogCiKkCQJWq12wms0NDSgubkZdXV1AccLCwvR0NAAt9vtD75T+SxFhlKnBwAMeNqRFtuhEBER0QwRlR33wqm+vn7MUOo7Vl9fP+q98T47kk43NGPpGdZJYSqfpchQXQnJXk97TMdBREREM8e0Z5IbGhrCWm6xYcOGCT8jSdKY3TB8s8dtbW2T3qe2tnbM4ydOnACAgJnqqXyWIkN55RsSL3slExERUZRMOyRPVtowFeFoBddxZcvjUDQ0NMBoNAa1+G8qn51Ma2srBEGY9nUm0tvbCwBwu92TfDIOSR6kAbhw+jTeTMTxR0hCP1MaE59p8uEzTU58rrEly3JU7jPtkLx+/fpwjCMowcxYh9ppw263w+PxjDtzHOpnKQwyNAAAobszxgMhIiKimSIh+yRPNFus1+unfD1RFLF3717s379/0kV4U/lssPLy8qCMcP9f33e74Zj5joXTGRqkA1iUoOOPhER/pjQan2ny4TNNTnyuseX1enH8+PGI3yehFu75SjvGmi32zTKHUrJhtVqD7rE8lc9S+Ci5NTURERFFUcxC8qOPPhrSeePVQPu6TCxevHhK17PZbNi6dWtQu+dN5bMUXkqtDgPsbkFERERRErWQ3NnZidbWVjQ3N+Oxxx5DQ0NDSNcpKSkZs9zCN5M8lQBrt9thNBphNpsDjjudzml9lsJPpdNjsKsT8sBArIdCREREM0BYapJbW1shiiI8Hg/a2togSRLOnDmD9vZ2nDlzZtSCO1mWQ+7mYDabUVZWNuq4y+WCwWAIuk7Y5XIBACwWy6jjI2erp/JZiozhbeBU2bNjPBoiIiJKdtMOyaIo4oEHHvDP7vracmi1WuTk5MBoNCInJwctLS3YtGkTDAYDcnJyYLVaQ7qfyWSCyWSCw+EICK0Oh2PU9tEAUFRUBIPBENCJQhRF2Gw2FBcXo6qqKuDzzc3NqKmpCemzFDnDd91jSCYiIqJIm3ZINhgMOHz4MKqqqtDR0YGvfOUryMzMHPU5m80WsIX1nXfeGfI9a2trYbfb4XA4AAytMq2srByz1EKr1Y5azFdWVgZRFLFv375xv6ZQPkuRc3XXPS7eIyIiosgLWwu47du3QxRFVFVVobCwcNTOeSPLK7Zs2TKt+5WXlwf1ubq6ulHHmpqagr7PVD5LkaPUXim34OI9IiIiioKwLtwzGAzYuXMnsrKy8NBDD6G1tdX/XrR2R6HkdLXcgjPJREREFHkR6W6xfv167Ny5E4cOHcJDDz2Ezs7OiG+7TMlNNWzhHhEREVGkRXTHvS1btqCjowOPPPIIzp49G8lbUZJTZg11EWG5BREREUVDxPskZ2VlobKyEg8++GBACUZjY2Okb01JRFCpoMjKwgBnkomIiCgKoraZSH5+Pnbu3IkTJ07AarVix44d0bo1JQmVVseZZCIiIoqKqG9LvXHjRmzdunXUBiNEk1Hq9AzJREREFBVRD8kAYDQasX79+ljcmhKYUqvDYE8PBvv6Yj0UIiIiSnIxCckAsHPnzljdmhLU1Q1F2mM6DiIiIkp+MQvJ3KmOpkrJNnBEREQUJTELyURTxQ1FiIiIKFqm3Se5tbUV7e3tYRgKoNfrkZeXF5ZrUfLxbyjCcgsiIiKKsGmH5BMnToStU4VOp2NIpnEptXoALLcgIiKiyJt2SN64cWM4xkE0KV9N8gBnkomIiCjCWJNMCUOZpQUEAV7WJBMREVGETXsmuaGhIazlFhs2bAjLtSj5CAoFlFla1iQTERFRxE07JGu12nCMAwCQlZUVtmtRclLq9KxJJiIiooibdkjmznkUTUqdDpf/cT7WwyAiIqIkx5pkSigqnR5yby8GL1+O9VCIiIgoiTEkU0JRan277rXHdiBERESU1KYdkqurq9Hc3ByOsRBNirvuERERUTRMOyRbLBYcOnQIpaWleOihh9Da2hqOcRGN6equewzJREREFDnTXrhnMBhQUVEBADh58iR+85vfwO12w2QywWKxICcnZ9qDJPJRcmtqIiIiioJph+Th8vPzUVlZCQBobm5GVVUVOjs7YTabUVJSgszMzHDejmagq+UW7TEdBxERESW3sIbk4davX+9vD9fQ0IBt27ZBEARYLBZuGEIhU/kX7rHcgoiIiCInYiF5uOLiYhQXF6OjowP19fX4zGc+A61WC4vFwj7LNCUKTSagVLLcgoiIiCIqKiHZJysrCxs3bsTGjRshiiIaGhpQVVWFlStXYtOmTcjLy4vmcCgB+bamZncLIiIiiqSY9Uk2GAzYsmUL6urqYLFY8Oc//xn//u//HqvhUAJR6fScSSYiIqKIiupM8njy8/ORn58f62FQglDq9Og7ewayLEMQhFgPh4iIiJIQd9yjhKPU6SAP9GOwpyfWQyEiIqIkFZWZ5MbGRmRlZcFgMLBvMk2balivZGVGRoxHQ0RERMkoKjPJVVVVAABJkqJxO0pySraBIyIiogiLWgs4tnqjcOGGIkRERBRpUQnJer3e/8/79u2DIAgwmUxs+UYhUV0JyV62gSMiIqIIiUq5he5KDSkAbNmyBS+//DJaWlpw5syZaNyekoy/3IIzyURERBQhUQnJLS0tAa9XrlyJ++67j4v4KCQstyAiIqJIi0pIdjgcyMvLw7333otHH30ULS0t6Ozs9L/f3NwcjWFQklBkZAAqFRfuERERUcREpSa5oqICxcXFaG5uxqFDh9Dc3Ix169YhPz8fJpMJoihyYR8FTRCEK7vuMSQTERFRZEQlJG/ZsgXA0FbUGzduBACIougPzc8//3w0hkFJRKnTwdv+bqyHQUREREkqZttSGwwGf2jet29frIZBCUqp1aP3n29AHhyEoODGkURERBRecZEuLBZLrIdACUal0wFeLwa7u2I9FCIiIkpCcRGSs7KyYj0ESjBX28CxLpmIiIjCLy5CMtFU+drAXaw7iK5jR2M7GCIiIko6065Jbm1tRXt7exiGMrQzH3fho2AMXLoIAOh68TC6XjwM3YY7oV6ydMJzlDo9MowrozE8IiIiSnDTDsknTpyAJEnhGAt0Oh1DMgXFF5J9PI2PI5jCi8XfehSp1yyKzKCIiIgoaUw7JPtauhFFU+Ytt6Lz+UP+1/oPfwzqZSvG/Xzv66/h3T/+Dpdf/TtDMhEREU0qZi3giKYj86Z1WGitQPcrJ5FxQz40a9ZO+Pm0a5cMheQ3TkP73g9EaZRERESUqKYdkhsaGsJabrFhw4awXIuSn2bN2knDsY8qezaUWh16X38twqMiIiKiZDDtkKzVasMxDgBsBUeRIwgC0pYsQ7f7Zcj9/RBSUmI9JCIiIopj0w7J69evD8c4iCIubekydL98DL1n2qBeujzWwyEiIqI4xj7JNGP4gjFLLoiIiGgyMQvJjz76aKxuTTNU2tJlAIDLr5+O8UiIiIgo3oVUbtHa2jrlfsadnZ0QRRHt7e04c+YMGhoa8OCDD4ZyewCA3W5He3s79Ho9RFGExWKB0WgM+nxJkrBnzx50dHRAFEUAQHl5OUwmU0TuR7Gn0umhyp6NXoZkIiIimsSUQ3JjYyM8Hk9ASG5tbYUoivB4PGhra4MkSThz5ow/EI/sfiHLMgRBCHnQ27ZtQ2FhIcrLy/3HioqKUFNTE1RwlSQJ1dXVqKys9B9zOp0oKytDRUVFwHXDcT+KH2lLl6PrpSMY7L0MRZo61sMhIiKiODWlkNzY2IgTJ04EzACLoogHHngAHR0dAIYCMDDU9SInJwdGoxE5OTloaWnBpk2bYDAYkJOTA6vVGtKAXS4XmpubsXv37oDj5eXlsFqtaGpqmvQae/bsCQjIAGA2m1FcXIzq6mqYzWYYDIaw3Y/iR9rSZeh68TB6//kG0q+/IdbDISIiojg1pZD817/+FSUlJQHHDAYDDh8+jKqqKnR0dOArX/kKMjMzR51rs9kCdue78847Qxqw0+lEQUHBqOMFBQUQRRGSJE3alq6hoQHNzc2oq6sLOF5YWIiGhga43W5/SA7H/Sh+DF+8x5BMRERE45nSwr3KykocOnQIra2to97bvn07ysvLUVVVhcbGxlHvjyyv2LJlyxSHOqS+vt4fYIfzHauvr5/0GmOdDwxtZgIAHo8nrPej+JG2hIv3iIiIaHJT7m5RUVGBtra2Md8zGAzYuXMnsrKy8NBDDwWEaV8ZxnRJkjTmpiO+2dzxxjZcbW3tqFlkADhx4gQABMwch+N+FD+UGg1S5i9gGzgiIiKaUEjdLYqLiyd8f/369Vi/fj327duHAwcOYPv27dNaqDcVvtroUDQ0NMBoNE5pMd507ufT2toa8X8/vb29AAC32x3R+yQCVfZsKFvdcL/4IqBO3MV7fKbJh880+fCZJic+19gK18TrZKa9495EtmzZgo6ODjzyyCM4e/bstK83skvGWNrb20O6tt1uh8fjQW1tbVTuR7EjL7gGaHVDuHAO8rVLYz0cIiIiikMRDckAkJWVhcrKSpw8eRIPPfQQNm3ahLy8PDQ2NmLDhg0hXXOi2Vu9Xj/l64miiL1792L//v1j1h+H+34j5eXlQalUTvs6E/F9t8uWdUBPihJnn34CCwe9mJXA/z74TJMPn2ny4TNNTnyuseX1enH8+PGI3ydqO+7l5+dj586dOHHiBKxWK3bs2DHla/jqgMeavfXN+o5VPzwZq9U6Zs/jSN2PYivt2iWAIOAy65KJiIhoHBGfSR5p48aNMBqNY3bACMZ47dZ8HSkWL148pevZbDZs3bp13J32wn0/ij1Fmhqpi3LQ+wY7XBAREdHYojaTPJzRaMT69etDOrekpGTM8gffzO54YXcsdrsdRqMRZrM54LjT6YzI/Sh+pC1djoGL72BA8kz+YSIiIppxYhKSAWDnzp0hnWc2m+FyuUYdd7lcMBgM4/ZAHuvzAGCxWEYdHz57HK77UXxRX+mX3HuaJRdEREQ0WtTLLXxCDZcmkwkmkwkOhyMg4DocjlFbTQNAUVERDAZDQNcKURRhs9lQXFyMqqqqgM83NzejpqYm5PtRYkjz7bz3xmloVt8U49EQERFRvIlZSJ6O2tpa2O12OBwOAEOrTCsrK8csfdBqtaMW15WVlUEURezbt2/M648M8FO5HyWGNMNiQKnk4j0iIiIaU0KGZAAoLy8P6nNj7azX1NQUsftRYhBSUpBmuBa9r5+GLMtR2+yGiIiIEkPMapKJYi1t6TJ4JQ8GLl2M9VCIiIgozjAk04yl9tUlv85WcERERBSIIZlmrLSlQx0uWJdMREREI0UtJD/66KPRuhVRUFKvyQFUKnQ0/xVdx47GejhEREQUR6K2cM/lcuHBBx+M1u2IJtX98jFgYADeSxdxrqYaisxMKFJTJz1P+/4PIfuj90RhhERERBQrUQvJsixH61ZEQeludQe8FpQqKGfNnvCc/vNvor3hT5h110chKJWRHB4RERHFUNRCMltsUbzJyDPC01jvfz1v8xZo1qyd8JxLv/8tLv3uMfScakVGfkGkh0hEREQxwoV7NGNp1qzFQmsFdMV3YqG1YtKADACam9YBALpeYg0zERFRMmNIphlNs2Yt5t7/6aACMgCk5higmjsPXS8dYQkRERFREmNIJpoCQRCQefM6DFy6iN5/vh7r4RAREVGEMCQTTZG/5OLFIzEeCREREUUKQzLRFKlXXA9llhZdLzEkExERJSuGZKIpEhQKaNasRd/ZM+g7fy7WwyEiIqIIYEgmCoHmpqGFfuxyQURElJyiFpLZCYCSSXp+AQS1miUXRERESSpqIbmmpiZatyKKOEVqKjQrV+Pya//AQHt7rIdDREREYRa1kGwwGKJ1K6Ko0Ny0FpBldB1nyQUREVGyYU0yUYgybrwJUCpZl0xERJSEGJKJQqTMyEB6nhHdJ1sw2NMd6+EQERFRGDEkE01D5k3rgIEBdL18PNZDISIiojBiSCaaBs2amwFw9z0iIqJkw5BMNA2qWdlIW34dul4+Drm/P9bDISIiojBRxerGjz76KB588MFY3Z4obDJvWouLr/0DZx/9NjIKViH9+hsmPUeRlobUxUsgCEIURkhERERTFbWQ3NnZCVEU0d7ejjNnzqChoYEhmZKCkJ4OALj8yklcfuVk0OcttFZAs2ZtpIZFRERE0xCWkNza2gpRFOHxeNDW1gZJknDmzBl/IJYkKeDzsixzBo2SRv/5cwGv1bl5yLx53bifl72DuHjwV/A8+zRDMhERUZyadkgWRREPPPAAOjo6AFzdflqr1SInJwdGoxE5OTloaWnBpk2bYDAYkJOTA6vVOt1bE8WFjDwjPI31/tezzHdNGn57WlvQ/fIxDHjaodLpIzxCIiIimqpph2SDwYDDhw+jqqoKHR0d+MpXvoLMzMxRn7PZbNi4caP/9Z133jndWxPFBc2atVhorUD3KyeRcUN+ULPDWbe9D90vH0fn84egL74rCqMkIiKiqQhbTfL27dshiiKqqqpQWFiIDRs2BLw/srxiy5Yt4bo1Ucxp1qydUumEZvXNUGRoID33LHQb7mT5ERERUZwJaws4g8GAnTt3IisrCw899BBaW1v97/nKMIgIUKSmIvNWE/rOtKGv7Y1YD4eIiIhGiEif5PXr12Pnzp04dOgQHnroIXR2dnKmjGgE7W3vAwBIf302xiMhIiKikSLaAm7Lli3o6OjAI488grNnz0byVkQJJ23pcqQsvAYdzx/CHMsnIahi1raciIiIRoj4jntZWVmorKzEgw8+GFCC0djYGOlbE8U1QRCgve19GOzoQNffjsV6OERERDRM1Lalzs/Px86dO3HixAlYrVbs2LEjWrcmiltZptsBQUAHSy6IiIjiStRCss/GjRuxdevWURuMEM1EqlnZyChYha6Xj8HL/yaIiIjiRtRDMgAYjUasX78+FrcmijtZhe8DvF50PP/XWA+FiIiIrohJSAaAnTt3xurWRHFFc9NaKNIzIP31L7EeChEREV0Rs5BsMBhidWuiuKJITUXme0zoa3sDvf98I9bDISIiIsQwJBPRVb6eyRd++hN0HTsa49EQERHRtBuztra2or29PQxDAfR6PfLy8sJyLaJEMiC1AwD62t7AuZpqZNx4E1LmzpvwHOWli5CvXQYYjVEYIRER0cwy7ZB84sSJsHWq0Ol0DMk0I/W8cjLgdfffXpr0HBUA+dhR9N50E9IM10ZoZERERDPTtEPyxo0bwzEOohktI88IT2O9//W8LZ9HRsGNE57z92YXUhy/wNv//TMs+vo3ufU7ERFRGHEfXKI4oFmzFgutFeh+5SQybsiHZs3aSc+RFy/B4Oq1uHz8KDoO/cVf10xERETTN+2Q3NDQENZyiw0bNoTlWkSJRrNmbVDheLiB29+P1Nf+jouOX0GzZi2UGk2ERkdERDSzTDska7XacIwDAJCVlRW2axHNCOp0zN74cbz105/gUp0Dcz/1mViPiIiIKClMOyRz5zyi2MoqfC+kvzwNz1NPIOv2D0C9ZGmsh0RERJTw2CeZKMEJCoV/Bvnt//4p5MHBGI+IiIgo8TEkzxCnFK+iUfkUjnYdi/VQKALSFl8LXZEZvadfRcdzz8R4NERERImPIXkGONp1DAdTfocXVC+i+lwNg3KSyr7nPih1erz9m1/grZ/v4859RERE08AWcDOAu7s14PXJ7lewVrMmRqOhSFFmZCDzPYXwNP4Z0tNNkJ5ugkKrgyItbcLzFCmpmL3pk9CsWh2dgRIRESUAziTPAMaMwF0M56dMvN0xJbLAemRBlqHM0Ez4q//SO7jwkx+g/523YzRmIiKi+MOZ5BlgrWYNNvbfg2OKl/Gq8jT+990/YlWGEQtTF8R6aBRmo3bu+8zWSXsvdx1/Cee+/wgu/OQHWPRVGwQV/1ggIiJK2P8b2u12tLe3Q6/XQxRFWCwWGI3GKV9HFEVYrVbU1NTAYDBMej8A6OjogMlkgtlsDnX4UZc7uAK5gysgXdOF3Rd+jF1vPoKHFn2Vs8pJJpSd+zSrb4Ku+E54Gh7Hpf/9H8z+l01RGCkREVF8S8iQvG3bNhQWFqK8vNx/rKioCDU1NUEH5W3btkGv16O9vR1ut3vCXQNtNhsqKioCNk6x2WzweDywWCyhfyExsD7rFnjhxQ8v7MWus9/FQ4u+hrkpc2I9LAqjUHbum3Pfx3H51Ct498+/R3qeERnGlREaHRERUWJIuJpkl8uF5ubmUeG0vLwcVqs16Ovs3r0blZWV2LRp4lkzl8sFg8EwamfByspK2O324AceR27LWo/Pz/ssLg5cwq6z34Wr43m83N2C1y6/jgv9b6HT24VBmb12ZxJBpcKCL2yDkKbGhT3/hQFPe6yHREREFFMJN5PsdDpRUFAw6nhBQQFEUYQkSWHdKtvtdvvLLJLJ+7S3wQsv9r5Vi90XfjLq/TQhFdYFX8BNmtXRHxzFRMq8BZi3uRwXfrIbF+w/wjX/8VUIioT7PpqIiCgsEi4k19fXo6SkZNRxXz1xfX19WEsgtFotqqursXLlyoAaZFEUJ6xhTgQf1L4Py9OW4c2+c+ga7EKntwtdg0O/mjsO47/f+Q1uzFgJpaCM9VApSrJuNaHn5AlIf3kaZ7+9E7Pu/MiUSzeIiIiSQcKFZEmSkJWVNeq4b/a4ra0trPezWCyw2+2wWq0oLi7Gww8/DI/HA5vNhsrKyrDeKxauTTPg2rTRYX+WchZ+++7v8VyHC+/X3h6DkVGspBesgvSXp3H5H6dwruYUFOkZEFJSJj0vdVEOFvzbf0Cp0URhlERERJGVcCF5Mh0dHWG/Zl1dHTZv3oyGhgY0NzejoKAAtbW1Yb9PPLlLXwyn5wnUXfoDbstaD5WQdL9VaByXX/17wGshLQ0p8ybugiL396On1Y0Le3+IhdYKlmkQEVHCS6jkM1EHCp9I1A9rtVqUlJQgJycHDQ0NcLlcqKqqwvbt28Ny/dbWVgiCEJZrjae3txfAUI11sNYpb8Izqr/iN68cxE2DN0ZqaBSiUJ5pMBSaLAyfN778gQ+he0XuxCfJMlSP/x7df3sJp366B17Te8M6ppkiUs+UYofPNDnxucaWLMtRuU9CTvdMNFus1+vDfr+ysjKYTCbs3r0bTU1NMBqN2LdvH8rKysJ+r3hyi/dmpMvpeE7VjAEMxHo4FCWDK3LRf89GDNz8HvTfsxGDkwVkABAEDGy4C4Nz50N56FkoTr8a+YESERFFUELNJPvqjseaLfbNMo9VrzwdNpsNZrPZ33/ZYDCgrq4ONpsNDocDTqdz2puK5OXlQamM7OI433e7U91w5Z53P4JfXzyICznvYIPujkgMjUIU6jMNSojX7L/m6xB3fh1K5x9geOhbSJk3P8wDS24RfaYUE3ymyYnPNba8Xi+OHz8e8fsk3EzyeO3dPB4PAGDx4sVhvZ/D4RizW0ZlZSWKi4vhcrnCer94s0F3B3RKLX536Y/oG+yL9XAozqXMm4/5n/s3DHZ349x/fQ+DV34kSURElGgSLiSXlJSMWW7hm0k2mUxhvd9EPZcLCwvDPnMdb9SKNNw96y68621Hk/RMrIdDCUCzag2y774XfW1v4Nz3HsHbv/45uo4djfWwiIiIpiShyi0AwGw2j1kL7NsZL9y9i3U6Hdxu95g/UnG73dMutUgERdoP4I/v1uP37/4Jd2jfhzRFWqyHRHFu1kdL0Xn8JfS84kbPK254GuuRdfsHkJYzyX+fCgGZN98CVfbs6AyUiIhoHAkXkk0mE0wm06gyCIfDMWbf4qKiIhgMhnFbtk3WMaOmpgY7duxATU1NQAB3Op3IysoK+8x1PEpVpOJjsz6M2nd+iW+/WY0Pz7oTazVrYj0simOCQgH10mXoe+O0/1jHc08jmAaN7fV/Qs6OXVDNyo7cAImIiCaRcCEZAGpra2G32+FwOAAMzehWVlaOGVi1Wu2YJRE2mw2iKKKlpQUAsHnzZhQUFCA/Pz+gtZvRaERNTQ3sdnvA+SaTKWwt4BKBXqUDALxy+R945VwNNujugDH9BmQoMpChSIdaoYZSUEEFJRSCAkoooRQUSFekc8e+GUqzajWkp5v8r7Pvux/p198w4Tm9b5zGO7/6Od589NtY9LWHoNRkRnqYREREY0rIkAwA5eXlQX2urq5uzONT2S3PYDAkxe560/FKzz8CXjd6nkSj58lJz5unmotv5nwd2apZkRoaxSnNmrVYaK1A9ysnkXFDflDbW6dflwu5fwAXD/4K52qqcU3F16FITY3CaImIiAIlbEim6DJm5KHe0+h/XTrrbixKXYjuwW50D3ajZ/AyvPIgBuH1/73b24O/djaj5vyPsWPRV7hr3wykWbM2qHA8nL7kwxjwvAtPw+O4sOcHWPDFL3MHPyIiijqmFgrKWs0aVCy04mT3K8jPuCHomuSstzNR73kCBy7+Dz45Z1OER0nJQBAEzLF8El5POzqfd+Ht//4Z5j7w2YjvSklERDQcQzIFba1mzZQX7H1ijgWv9r6OP7U7cb36OtySeXOERkfJRFAoMH/LF+Dt6ID0TBMGL/dAqdUiI8845ZlpIiKiUDAkU0SpBBX+fcEX8NW2h/Djt/bBkLoIC1MXxHpYlAAElQoL/+3LaLN9FZ3PHwIAeBrrIaSmQlBO/EeXUqvF/M9/Ceqly6MxVCIiSkIs9KOIm63KxpcWfA6XBy/je+d/yJ37KGiK9Axk5BcEHsvKQtqy5RP+Grh0EW8+8p+4fPrVGI2ciIgSHWeSKSpWZRTgX7I/hscu/Q6PnPs+DKk5MGbksd8yTUqz+iZIzz7lfz3vk2WTllx0n2zBue8/gjer/hMLH/wa0ldcH+lhEhFRkmFIpqi5Z9ZHcLTzJbT0nERLz0nUexpxl64YuenXIUVIQaqQCqWghAABCkEBxdA/QSUokZO6CAqBP/iYiUJpJZeRX4Br/uOrePN738Wb1d/CNf/x1Ul7NBMREQ3HkExRoxAUWK5ehtf7/uk/9mdPA/7saZj03DUZN6Ji4TZuTDJDhdJKLv2G/KtB+dFvDwXl3LwIjZCIiJINQzJF1WrNKjRJT/tf362/C4vTDOiX+9En913psTwIWR7EIGQMyoP4x+VX8WL3cdS+/Qt8du4DbAVGQUvPzcM1D34Nbz76HZyt+k9kFKxCRsEqpOcZgzo/df5CCCr+MUlENBPxT3+KqlD6LQ/IA/jum99Dk/QM5qbMxd2z7orCSClZpF+Xi+yP3oOLB3+N7uMvofv4S0Gfm7LgGiz6P9+AalZ2BEdIRETxiCGZom6q/ZZVggpfXvBFfPPst/Cbi49hrmo2TFm3RnCElGwG2t8NeJ22fAUyCm6c8BxvhwTpqSdw5lvfxKL/swMpc+ZGcohERBRnGJIpIWQoM/CVhV/GjjO78KML+zBLNQt56bmxHhYliIw8IzyN9f7X2R/+WFA1zmmLr8XbP/8pzn5751BQnjc/ksMkIqI4wnYBlDDmpMzG/7nmy1AJKlSf243G9ifx87d/jaNdx2I9NIpzvg4ZuuI7sdBaEfQiQN37izDvs5/HwKWLOPPtneg7/2aER0pERPGCM8mUUJakXYt/X/AFfPfc9/Czd34BAKj3NGKD7g5cp16OFEEFlZDibx+ngABBGPqnbNUsXJO6MMZfAcVKKB0yAEB72/sgqFS4sPeHOPvtSsz6yD3of+s8t8gmIkpyDMmUcFZrVsGYnoeWnpP+Y42eJ9HoeXLScz8/77N4v/b2SA6PklDWrYUQlCqc/3EN3vllLYChLbIz19+G1GsWTXiuIi0NWYXvg1KjicZQiYgoTBiSKSGZ9R8KCMn3zrobS9OuRb88gH65HzLkK63k5KFWchjEH979M/a+VYssZSZu5k5/NEWZ694DzaE16Dr+ov9YZ/NfgzrX82QjFmx7EGmLDJEaHhERhRlDMiWkUFrJFaTn46Ez/4nvn/8RvnHNV5Cbfl0URkrJRPu+DwSE5Dn3fxrpN+RPeE7PKyfxzsFf4UzlNzB/y78icx07sxARJQKGZEpYU20ld03qAnz1mi+j8ux38ci57+Gbi74OQ1pOBEdIySaULbLTrl2CtCVLcf6H38f5H34f+rvuxux7LRAUXDdNRBTP+Kc0zSjL1cvw4MIv4fJgL771ZjXe7n8n1kOiBKNZsxZz7//0lBbtpefmwbDz20hbtgLtf/49xG9+DW/9fB+6jh2N4EiJiGg6OJNMM86qjAJ8cX45fnBhDx468zDWaFZjjebGKc1KE02ValY2cr72EN78/iPocZ9AX9s/IT3dBNXsOVCkZwAAUi5fBgC0qdWjzk81LMbcTzwAZWZWVMdNRDRTMSTTjGTKuhUt3a14quNZPCk9gyelZ7AoZSH0Kj1UghJKKKEUlACEoWZyggBAgAICVmeswnu1hbH+EigBCSkpSF2Ugx73Cf+xwZ4eCCkpQy+8AwAAeaA/4DzZ60Vn819x+VQr5n/BivQV10dtzEREMxVDMs1YaYq0gNdvD1yExyvBK3sxgAF45UHIV/4aztX5Av5++VVsnvsJqAT+J0RTM3L3v/nl/+ov3XC73QCA643GgHNkWUbHc8/g7V/8DGe/vROz/+V+6M13QRCEaA2biGjG4f/hacYyZuSh3tPof71twb+OW3Ihy0NhuXOwE7vP/wRN0tMQ+87iywu+CL1KF60hUxIIZfGfIAjQvvcDSFu6HOd/+H1cdPwSl//eCs0t69H7+mvc2ISIKAIYkmnGmkobOUEYKrvQKrX42jUP4lfvHMTjngZ8/cxOVCzYhmXqJdEbOCW8UHf/SzMshuGb38JbP9+Hzua/ouvYUDs6T2M9tB/4ENIWXzvh+YJSCc2am6HM0oY0biKimYQhmWa0qbaRAwCloMSn596PJWmLYX+7Fg+d/U9s0H4Qg5BhzMjjAkCKKIVajflbvwivpx09J1v8x6WnnwjqfKV+Fhb86zak5+ZFaohEREmBIZkoRO/VFmJR6kJ8+81H8WdPAwCg3tOIbGU2tMpMKAUVUgQVFIJiaPEfhv4uCAIyFOkonfVR9mmmkAiCAP2HzAEhOfue+6C+buIFff3nz+Edx69w9juVyL7nPsz68MfYr5mIaBwMyUTTsFy9DLdkrsVT0rP+Y31yL3pkFfoHB+CVBzAoD2LwykbZvm2y++Q+vNT1N3x+3mexPuuWGH4FlKhCqW1G/kqk32DE+R/V4FLdQfScasX8rV+ESqeP+HiJiBINQzLRNN2kWR0Qkj8/f8ukJRev9Pwd3z//Q9Rc+BFe7T2Nj8++70rLOaLghVLbnHrNIuTYHsY7v/45pGeehGj7KnQf/BAGJAnpuXnQrL5p0msIKhVnoIko6TEkE03TVBYA+tyQfj2+Zfgmvn/+h/hzuxOv974B6/wvQKfigiqKPEVqKuZtLkf6DUZc+OmPcel3jwEApCcbgjpfqdVh7gNbkHnzukgOk4gophiSicIglAWA2apZsC36Kn7xzm/Q4HkSD7Z9Hcb0PKzTrMHazJugElRQQsleuBQxWbea0H3iODoO/cV/LGXhNVAvWTbuOTKA7peP4/wPHoX2A0WYs+lTUKSljft5IqJExZBMFEMqQYWyuZ+CAkrUexrxQtcRvNB1BHjr6md8u/8Jvr+EoZ3/vKmDWCjPx3/0fwlzU+bE7oughJa59paAkDxn48cnLeEYePcSLuz9IaSnm3D51CuY/69fQpph4vZzRESJhiGZKA4tSrkGS9IWY0AegBdeDMjeob3/ZNm/C6Cny4PXFf/EV0UbPjfvs7gl8+ZYD5sSUCgLAFWzsnHN9v+L9sf/iIu/O4gzO7+BzMLboEhTc2MTIkoaDMlEcWDk7n/3z7lv0vINt9uNfyhew5/TGvH/zv8Axbo78InZFqQqUiM9XEoyoSwAFBQKzPrw3UjPN+JcTRU6nn0awNDGJqo5c6BIz5jsCshce8tQGzolF60SUfxhSCaKA6Es/gOA6waX4zuLK/Ff5/egwfMkTvX8Ax/Uvh/n+s9zYxOKCvWyFdDc/B5IT139Jm+wuxuCcuL/vQz29ODS7x5D18vHseBzX0TKvAWRHioR0ZQwJBPFiVAW/wHAbFU2diz6P/jtpd/jt+/+Hj97578BDG1s8h7NWixOMyBVSIFKSIHyysYmAhRQCEMbnKgEFdZkrIJGqQn3l0QzhGblqoCQPL/8C5POTMv9/bj4WwfanX9Cm+2rmPupzyDLdDsXqhJR3GBIJkoCCkGB+2bfgzN9Z/FC11H/8Re6jga8Ho9eqcPmuZ/EezRrGVJoykKpaxZSUjBn0yeRsfJGXLD/CG/Zf4Tuvx1Dxpqb0fv6a6xtJqKYY0gmSiK3awsDQvFn5nwKuenXo1/uR7/cH7D736A89E8X+t/GY5d+h++f/yHWatbgM3M/jWzVrBh+FZSIQqlrBoAM40osfvgRvFW7F52Hm9F5uBnAUG2z+voboMqePeH5ygwN9Hd9FCmz2eGFiMKLIZkoiYRa23xr5jr89zu/hqvzBbi7W1GYdStShBQYM/JZ10wRp8zMwoJ/+w+8WfWf6DnZ4j9++e+vBHW+5HoOczZ9Etr3fZA/CSGisGFIJkoyodQ261U6bFvwr7itaz1+fGEfmqRnAAD1nidwY3oBFqTOhxIqqITAns0KYWhrYgUEzFbNRmHWrdxem0IiCAL0HzIHhOQF//ZlaFZP3Nqwu9WNt2v34O39dnQefQHzyrZyVpmIwoIhmYj8btKsxvrMW9AoPeU/9reeFvytp2WCs676U7sTn5n7KdyQfn2khkhJLJTaZs3KG6F+uBoXD/wC0l+eRtv/3T60AFClZF0zEU0LQzIRBVilWRkQkr84bytWZuRjQPZiAF545YErG5sAsq/CWR7E4a4X8cd36/HNs9/Ce7MK8fHZG6FX6WL4lVAiCqW2WZmRgXmf+Rw0a9+DC3v/y99pw9NYj1TDYiiztJNeQ73iesz68MegSGWfcSIawpBMRAFCrWteql6C92YVYv87v8RfOg7hSNdLuFVzC9IUqShgbTNFgWbVamTesh7SU0/4j/WdexOKixcnPE/2DqDnZAs6Dzdj3mc/j/TrciM9VCJKAAzJRDRKqD2bF6YuwFcXPogjXS/B/tZ+PN3xLADA6XkCGkUG0hXpUAlDtc0KXOnZLCggAEO9myFgjeZG3KU3Q61IC/NXRTOBZuWNASF54Rf/ffKezV4v2hv+jEt1j+Hst74J3YfMmH2vBYo0daSHS0RxjCGZiMJKEATcknkzWrpPolF60n88VUiFTqmDVx5AvzyAfvjKNoZ+DWIQlwcv47FLv0OT52lsmv0vuD3L5F8cSBSMkHo2K5WYdedHoVl9M9762U/gaaxH9/GXkHX7++HtkFjbTDRDMSQTUUSs0hQEhOTPzntg0tnpAXkAT0nP4rGLv8OP39oHp+cJfHL2JvTIl+HubuVW2xSUUHs2p16zCIu+vhOeJ5x457Ff49JvHQCGapsVWVoo0ib+6YaQkgJ98V1sRUeUJBiSiSgiQqltVgkqbNDdgcLMW/G/7/4J9e1PYNeb3/W/X+9pxB1Z78Ny9TKoBBWUvrINQQgo2dAoNbhBfT2DCk2ZoFBAX3wnesV/ouOvz149LghQZmVNeO7AO+8MtaI73DzUim7uvEgPl4giiCGZiCIm1NpmjVKDT8yx4A7dB/Cds9U4P/CW/70nO57Fkx3PTnD2kGVpS3H/7PuwMiN/yvcnyrx5XUBInldWPunstLe7GxcP/grSM0+i7f9ux+z7NgELDQC/WSNKSAzJRBS3FqTMwyfn3o/qczX+YxuzS7E8bSkGMFTbDMgYlGUM/+v05dfxhOdp/Oebj2BluhEfn30flqqXxOzroMQTSm2zMiMD8zaXI/OW9Xirdi/e+dXPkZI9G/KCa9DR1QHN6puCuregUkFQsBafKNYYkokoroVStnF7lgkl+g147NLv8NeOZnztzDdxg/p6ZKtmIT/9BqzOWAWVoPSXbCgxtEugQlD4dxMc6rzBGcCZLNTa5oz8Aix++BGc/9FudP/tJeDSRVw4eSLo81Wz52DOJzYj8yYuFiSKJYZkIop7oZRtzEuZiy/O34oP60uw58I+vHL57wAAV+cLQV9jSdq1uFt/F96TuZZdNmhKFGlqpMyfH3AsZeEiqJcum/A8WZbRffwlnN9dDc2atZjzyc3cZpsoRhiSiSipXZtmQG56Lk73/dN/bEnqtbg+fQW8shcD8gC88AIyAko2BmQvXu4+gZoLP8LCSwtw96y7cFvWeqgE/rFJwcnIM8LTWO9/PWfj/UHNTA+0v4t3fv3f6DzcjO6TJ5B9z0akzJ2LnlOtbEdHFEX8036GePWUAm3/FNDX04U1azWxHg5RVBkz8lDvafS//pfZHwtqZlrydqC+vRENnib85K2f4n8u/S9WpRegT+7DderlWJVRAKWvbAPKYeUagCAMddpQCkqkCCkR/OooXmnWrEX/PRshtP0ThttuDzrcqvSzsOALVnTd9j68/Yuf4eKBX/jf8zTWQ3fHBqQtXjLxRRQKZKxaDZVOH/oXQDTDCbIsy7EeRCjsdjva29uh1+shiiIsFguMRuOUryOKIqxWK2pqamAwGIK6p09hYSFMJtOU7wkAXq8Xx48fBwCsXr0aSqUypOsE49jRLtRUn/O/1uoUSEub+EfHKSkCPvP5+VhxHXecildutxsAQvp9PxMd7To25a22fboHe/CE50n877t/Rs9gz5TOFSDgPZnrcLf+zkkXD/KZJp/pPtPB3l6c/c5O9L5+esrnKtIzMHvjx4f6NnMhYFjxv9XYilaGSsiZ5G3btqGwsBDl5eX+Y0VFRaipqQn6N+y2bdug1+vR3t4Ot9sNSZLG/awkSbBarbBYLP57Op1OWK1WHDlyZHpfTBS0ursDXsuygIyM8X9DDQ7KENv68IKrgyGZkkao7egAIEORjrtnfRgX+98N2CBledpS3JCeC6/shVf2YmjfQBm+jhsA8PbAO3i+8zCe7zyMlelGfHTWnShIz+eiQAqKIi0N2R8txbmaav+xWXffC/Xy6yY8z9v+Li7+7jG8/fN9kP76LOZt3oI0w7WRHi5RUkm4kOxyudDc3Izdu3cHHC8vL4fVakVTU1NQ1/Gd73K50NDQMOFnN2/ejIKCApjNZv8xURSh0+mmOPrYyDNmoLHe43/9ma3zJiy5kGUZ/1b+Ol77x+VoDI8oYYzcRfCe7I8GFbz/2duGP7z7OFydL+BEjxvzVfMwL2UurlMvR156LlSCCipBiXPCBSiggLb3LBSCv3ADCgjQq3RQK/hN60wUSjs6AMhcdysu1h2Ep8kJ8aGvQb/hTqQtW47Lr/6dtc1EQUi4kOx0OlFQUDDqeEFBAURRhCRJ0Gq1Ybufy+WC2+3G/v37A46Xl5cHzGTHszVrNbhnYz/a/ingttsNk9YkC4KA5depcfJEN/r7ZaSkcMaLCAitHR0AXJu2GF9a8HlY+u/Fz97+BY53v4wLA2/hRI8beHfYB1Ov/F0cfQ21oEax7g7cNcsMrXLind8o+YTSjk6Rno65n3gAWYXvxds/t6Pd+Sf/e57GeqQsyoEyc+LfSyqdDtkfuw+p1ywKadxEiSzhQnJ9fT1KSkpGHffVE9fX18NisYTtfna7HQaDIazBOxZW5A5iRS5gNAa3aG/5CjVePtaNtjd6sZwlF0R+0ynbmJcyFwtTFuA4XvYfK0jPx1rNGgzIAzh74U144cXcuXP9XTYGIWNQ9uKlrr/h9+1/htPThA26D+LD+hLoVIn95xJFh3rJUuTseBhnH3kYl1856T8+8NZbGJyg1BAALp9qReeLRzCr5COY9ZF7oEhLi/RwieJGwoVkSZKQlTX6O19fiG1rawvr/VpaWvyz1E6nE1qtFm1tbdNatJcIfMH4tX9cZkgmCqORnTbM+g/5Q7f7zSuLgWaPXltx/+z78ELnUdS9+wf8sb0eDZ4nsSrDCBVUWKFePvTPgsr/SwklBGGoVGNoYxQFlFAgTcGQMxMJCgVmFd+Jc8NC8oIvbJt0drrn1b/j7Z//FO/+6X/R8fwhzP1kWdA7BxIluoQLyZPp6OgI6/V8odzlcgWUV5SWlsJisYRl1rq1tTXii3h6e3sBXF2RO+nn+wEgDS+9eAGLrj0buYFRyKb6TCk+pCMVGxX34J9CG66VFyP9jVS4MfQMJ3umWmjwADbhFcXf8YTyGRztOgYAeL7rCHAxuPsvHszBLd6bkDt4HRRgx4NIi6v/TlPTobjSkk5efC3eSE0HghnXxk9C+dJhyH99Fue+/wi8CxcBulkYzMmBPFkrOgCAAFk/C4hgF6doi6vnOgNFqzFbQoXkiTpQ+Axv0Rau+zU0NODhhx8OeG/r1q2wWq0oKSlJ+FKMsajVwOw5gzh3lv8TJQq33MEVyMWKkM4VICBvMBdtwlkcVrzoP75o8BosGVyMQXgxAC9kYXDY1ihDG6X0oAf/UJxGW8oZaOUsrPWuwRrvKoiKs/7QnjsY2rgoMQyuyAVW5E7tJIUC3rW3wpubD9Uffgvlm2eAc2ehfKUl6EvI+lnoLzJDXsrfX5Q4Eiok+0w0W6zX68N+P6PROCoI+1rNhaMGOi8vL6J9koHQejrmF1zAc890YNGiXOj1CflbJamxT2fymcoz7enqw+FzV0Py/YvuC6pWun2gHU9Kz+IJz1N4SvgLnlUdGtpxEMALeBG3Z5qwOM0ApaCEEgoIguJKj42rf2UpM3FjxkqkKlInvhkl3X+nb5/+OzxvnvG/Tlu2AhkFqyY8Z7CnB9JfnkLq//wGmrXvwdxPPADVrOxIDzWiku25JprhfZIjKaGSjy+ojjVb7Jv1Hateebr3G6vVm+9YMv+oZfkKNZ57pgOnX72Mm9Zmxno4RDRMqJ029Co97s2+G3fPugvPdx7BL975DTzeqz+le67TBXROfp10QY33ZK7D7Vkm5KXnQiHwp04zwcittrM/8rGgum7ozXfhnV//N7qOvoDulr8h+2P3cattinsJFZIBjFva4PEM9QFevHhxWO9nMBj81x7rfsnMt2Dv9D96GZKJ4tB0Om2oBBVuy1oPtUKN6nM1/uOfmr0J16df598gBcDVDVIgA7KMs/1v4rkOF57peA7PdDyH2apsLE9bBkDG0rQlQ6EZCigFBRRQQBjW89k/Ky0ImJ8yj1t2J5hQezanzJ6DhV/6D3T97Rje/mXtqK22M1atQcrceRNeQ5GuhvYDH0LK7DnT+hqIgpVwIbmkpASiOLqJqG8mOdwdJ0wmE1wu17jvJ/OPWhYZUpGWJuC1V7mpCFGyCmVG+kasxJ36YpzpO4vnJBeekp7F4a6jAIDDXS9OcvZVGoUGt2Xdivdr34uladwNLlGE0rPZf+6Na5CeZ8TZ71Si9/Sr/uPdLx8L6vz2Ridm3fVR6Es+AkUqy30oshIuJJvNZpSVlY067nK5YDAY/P2Sw8ViscDhcIzapMTlckGr1Y7ZszlZKBQCli5X4/SrlzE4KEOh4KYiRMko1BnpnNRFuH/OfeiT+1DvecJ/PE+di5s1qzEI+cps9NXFg7jy9355AC92HUOD50k0eJ7EktTFeL/2dmgUGpzufQPGjLyQZ8kpvilSU5H9kY8FbLU9b8u/IqPgxgnP633jNN458Etc+t1jkP7yNOZs+iQ0a9/DLd4pYhIuJJtMJphMJjgcjoAFcw6HA5WVlaM+X1RUBIPBgNra2jGvN1nHDKPRCIvFgm984xv+rawlSYLdbseuXbuSsrPFcMuvS8MrJ3twVuyD4Vr2VyWi0YwZ+QEh+a5Z5qAC7qbse/Fa7+t4RvoLDnW+gP3v/Mr/Xr2nEVpFFjKUGVBCCaWghAICcKX3s28pYbpCjfWZt6Aw61Zu251AQinbUK2+CRkFq9De1IB3f/8/OP/D7yM1x4CUhdcg/fo8pOcF95PdlNlzoEhPn+6XQDNAwoVkAKitrYXdbofD4QAwtHiusrJyzFILrVY75mI+m80GURTR0jLUwmbz5s0oKChAfn4+tm/fHvDZyspKOBwObNu2DXq9Hu3t7aipqUnqUguf4ZuKMCQT0VhCXUQoCAJWqJdhhXoZPjXnfvy/8z/A37qHtRUTgFQhFV7Ziz65b6id3aCM4X+JfZ1o6TmJX75zALdlmVCkez/eHrgId3crZ6PjXChlG4JKhVnmu5Blug0X9vwXetwn0HdGRNeRF4K+hkKjQfbH7oPuA0UQVAkZgyhKEvZ3x/CNPSZSV1c35vGxZp0nEq6NQxLNshVXQvKrl/H+otFdPoiIgOktIgSANEUaPqS7IyAkb533mUmv2T3Yg0MdzWjyPI0npKfwhPSU/716TyNKdBtwnXoZlIISKkF1ZQOVwAWEAoYWEc5LmRvy+Cm6VFodUhfloMd9wn8sbfl10KxaPeF58sAAOg49i3d+tR/S002Y8/FPT9rCjmauhA3JFB16vQpz5qrw2qu9sR4KESW5UGakMxTp+JDugyjSfgCv9r6Gn1z4Gc72v+l/v97TiPogmxEtT1uKWzNvwa2Z6zA3hR0U4t2odnQfvjuomelZH74b7z7+R7Q//ge8Wf0tZKy+GRmrVqP//JtsR0cBGJJpUstWqHG4uRNdXV5oNMmzrSgRxZ9QZ6QFQcB16hW4f859AS3t7pn1ESxJWwyv7MWA7IUX3qEFhDL8uxEOwou/X34VL3Yew68uOvCriw4sT1uGnNRr0DXYhSVp1+IG9fVX6qIVUAgKBG6xAkAQkCqkYGHKAvaMjpJQ29Ep0tSYfc990N7+flw8+Gt0Hm5G9/Ghriyexnqo5s2HUqOZ8Bopff3wrloDOS8PgoLPO1kxJNOkll83FJJff60XBasyYj0cIqJxhVofvUF3B/rm9uHlHjee7ziMw51H8VrvaQDA0a7g2pMBwIKU+bhD+368T3sbtMrwbW5FY5tOO7qUOXOx4AtWnJMHA2qavZIH8HonPFeQJKTU/wFnWlsw5+OfRvr1N4Q0BopvDMk0qeGL9xiSiSjehTobnapI9Z9b+1YGGqQn/e8Z0/OwTnMTvBiEV/ZebWonX21w5/F24FDH8/jVRQccF3+LWzLXIk2Rgi6hGz1dfVxEGKe0ptsDQvKCz/3bpMHb/eKLUDY/Bxw7grPf+iYyb7kV6cZV6DsrsmQjiTAk06SuXZIGlQo4zU1FiGiGWKkpCAjJJfoNQYXcT8zZiMOdL6LJ8zRcnc8DVzYUPHzuRSxUzYdWpYNKUEIlKKGA8uo+hMLVxYQ6pRarNTeiID0fagW7CkVaSGUbajW8H/gQlv2LBRcdv0Tn4efRefh5AEMlG1m3fwBpOZPv26CaOw+a1TexZCNOMSTTpFJSBCxekobXXr0MWZbZuJ2Ikl6oZRspQgoKs25FYdat+NEFO/7Sccj/3iVvO7rlnqHa6Cv10YGz0UN/H8QgmqRnkCKokJ+ehzUZNyJFUOFM35tsaxchoZZtpC5YiIXW7Tj3g0fR9eIR//GO555GR5DXSFu6DLM3fgIZQfZ5puhhSKagLF+hxulXe/HWhX7MX8CtQIko+U23rd0tmWsDQvKXFnw+qOud77uAl7r/hmNdf0NL90n8rftqm7N6TyPShXSkKVKHLSIcmoX0tbMDAAUUWK5eirWaNViVUYB0BTfPiCTtbe8LCMnZ990/eZ2yLKPzxcPwNDXgze/uQsaNN2HOxo8jdVFOhEdLwWJIpqAsv06NJ5wevPaPXoZkIqIgrNWswcb+e/BPoQ23G24LOnAvSJ2PO1M34E79BvQM9uAH53+Cl7r/5n8/XaHG3JQ58MpeDGIQg5ABWfZ36wBk9Mn9eK7Dhec6XFBBhYKMfNysWQ0llGjrE2FMz8e6zJsi8nXPRKF22ki//gbo7yjGxf85gM7DzWh7+RjSjSuhyNBAvXQZ0nPzJr+IQoG0nMXcGCUC+G+UgjJ88Z7pdq7YJiIKRu7gCuRiBYya0H6Unq5Ixwd17w8IyZ+Z9+mgAvc7/RdxtOsYXuw6hhPdbhzvftn/ntPTBCUUSFWkIVVIgVJQQQEhYHZaISigUWhgTL8BqzIKsEK9DCqBsWE8oZZspMybjwVfsOKy+S5c2Pdj9LQMPaeuw81BX0M1dx5m/8smZN6yniWRYcTf7RSUOXNV0OqUXLxHRBRlodZHz0mZDbO+CGZ9Ebq93dh94ScBQTlblY0FKfPRL/ejXx4YqoaWB+H/S/airU/Eqct/R927f0C6oIYxIw96pQ5d3m4sVy+DMSMPVyI1lILCX+7h+zuuBLZspR6pCv4UciLqZSuQUbAKnjfPXj12/Q3QrL55wvO8nR2Qnn4CF368G+0Nf8YcyyeDm4GmSTEkU1AEQcCcOSq8froXh5/vwC23hj6bfOxoF1rd3cgzZmDN2okbthMR0fTrozOUGSjSfSAgJD8w9xOTXtMre/Ha5dfxck8LXu5uwYtdx/1FHc1dh4GLwd1fAQUWpV6DJWmLsTTtWixNW4JLA+/i1cunp7UY8WjXMbi7W5NmQePIXQRnlXw4uF0ESz6MS3+og+epJ3D22zuhWbMW6QWr0H/h3LRa0nUdO4ruVveMbWsnyLIsx3oQM5HX68Xx48cBAKtXr4ZSGdmd7NxuNwDAaAztR37HjnahpvpcwLGRP9EJ5ic88lDpXMA5k52nUAAPbJmH29+vDXK0M8N0nynFHz7T5BNvz/Ro17Epz0gP99O3fo4npKf9r69PW4EbNSsxKA9e2c0Q/j/kfWF6EIO40P82Xu99A+8MjJ+q/bPPI/55PMOrsH3nBHPeB7XvxZZ5myf93EQi+Vy7jh2dcm2zT9+F87j42G/QdfSFwDeC+Z/tSKH8D1sQkH3Pfcj+8Memdq8pilaG4kwyBaXV3R3wet58FQzXDuvfOc63WiMPn2nrxVsXBvyv585TIWfxxH1AFQKwKIc/piMimq7pzkjfqFkVEJI/mn3XlK7X4e3EG71t+O2l3+GVy//wH5+vmo9r04a6Oow/cxf4TlvvGVwYeMv/ep5qLhanTd4Z4ob064MebyxMZxfB1PkLsPDfvozzP6pB57CaZtXceUjLWTyla/WeacPAWxemdg0BUC9ZOqX7xDOGZApKnjEDjfUe/+v7PzU3pFKJkTPSoV6HiIiiL9T6aJ8sZSZWZuSjV+7FK+dq/Mc/NXfTlK91tOsYqgOucX9SlFyEQ9b6woCQPPf+T005eHcdO4pzNdXTukaiY0imoKxZq4G1YiFeOdmNG/JDryUO13WIiCg2pjsb7bvGdMJ2uK6RrEJtSRfuayQ6hmQK2pq1mrCE2nBdh4iIEle4wjbD8dimU7YRzmskMm4WTkREREQ0AkMyEREREdEIDMlERERERCMwJBMRERERjcCQTEREREQ0AkMyEREREdEIDMlERERERCOwT3KMyMP2Q/d6vVG7XzTuRdHBZ5p8+EyTD59pcuJzja3h/96H56lwE+RIXp3G1dfXhxMnTsR6GEREREQJa+XKlUhNTY3ItVluQUREREQ0AmeSY2RwcBADAwMAAIVCAUEQYjwiIiIiovgnyzIGBwcBACqVCgpFZOZ8GZKJiIiIiEZguQURERER0QgMyUREREREIzAkExERERGNwJBMRERERDQCQzIRERER0QgMyUREREREIzAkExERERGNwJBMRERERDQCQzIRERER0QgMyUREREREIzAkExERERGNwJBMRERERDQCQzIRERER0QiqWA+AIs9ut6O9vR16vR6iKMJiscBoNMZ6WDQBSZKwZ88edHR0QBRFAEB5eTlMJtOYn+czTlyiKKKqqgq7d+8e830+28The1Y+hYWFY/43y2eaGIY/z46ODphMJpjN5gk/y2eaZGRKal/60pfkAwcOBBy744475JaWlhiNiCbj8XjkHTt2BByrr6+Xr7/+ennv3r2jPs9nnNg2b94s33PPPWO+x2ebGDwej7x582a5vr7ef6y+vl5eu3btqM/ymSaGHTt2yB6PZ9Sxkc9OlvlMkxnLLZKYy+VCc3MzLBZLwPHy8nJYrdYYjYoms2fPHlRWVgYcM5vNKC4uRnV1tX9mGeAzTnROpzPgeQ7HZ5s4Nm/eDIPBEDDLKIoidDpdwOf4TBODy+WCwWCAVqsNOF5ZWQm73T7qs3ymyYshOYk5nU4UFBSMOl5QUABRFCFJUgxGRZNpaGhAaWnpqOOFhYUAALfb7T/GZ5y4JEmCKIrIz88f830+28TgcrngdrtRUVERcLy8vBxNTU0Bx/hME4Pb7Q4om5kIn2lyY0hOYvX19TAYDKOO+47V19dHe0gUhLGeGQD/rJTH4/Ef4zNOXHv27EF5efm47/PZJga73T7mrONY+EwTg1arxb59++B0OgOOi6I46vnxmSY3LtxLYpIkISsra9Rx3x/mbW1t0R4SBaG2tnbM4ydOnACAgFkLPuPE5HK5/D8ZGA+fbWJoaWnxzxo6nU5otVq0tbWNuWiPzzQxWCwW2O12WK1WFBcX4+GHH4bH44HNZhtVCsdnmtwYkmewjo6OWA+BpqChoQFGo3FKK6b5jOOT2+2ecBY5GHy28cEXklwuV8AzLS0thcViGVWrOhE+0/hRV1eHzZs3o6GhAc3NzSgoKBh3AmMifKaJjeUWSSqYOqhga64o9ux2OzweD2pqavzH+IwTk91unzQg89kmBt9zamhoQElJScB7W7duhc1m83+GzzSxaLValJSUoLi4GJIkweVyoaqqKuAzfKbJjyE5yU30Xaxer4/eQChkoihi79692L9//5i1b3zGicPtdoftJwF8tvHDaDSOqkn2PeeRNal8pomhrKwMJpMJu3fvRlNTE4xGI/bt24eysrJRn+UzTV4MyUnK9wf2WN/F+r77HauOiuKP1WpFTU3NqHDFZ5x4Hn/88XE3hBmOzzYx+J7TyFZvw4/5utHwmSYOm80Gs9ns/zPXYDCgrq4OFosFLpfLv6CPzzT5sSY5iY232trXHWHx4sXRHA6FwGazYevWreMGKz7jxOFwOHDy5EnYbLaA4ydPnvQvCgLgXxjEZ5sYDAZDQMcZn7GO8ZkmBofDgVOnTo06XllZifb2drhcLn9PbD7T5MaQnMRKSkrG3KjA9x1uMDNaFDt2ux1Go3HUNqhOp9N/jM84cYy3iKusrAxarXbUqnk+28RgMpngcrnGfX/4T4D4TBPDRO38CgsLAzpW8JkmN5ZbJDGz2TzmH96+3YTG68dLsed7biNDlcvlCvgDnM84efHZJgaLxTLmphG+/1aHL+jjM00MOp0uYNOm4dxud0D7Rj7T5MaQnMRMJhNMJhMcDkfAcYfDMWrWiuKHKIqw2Wxob29HVVVVwK/q6uqAP3T5jBPfeLty8dkmBqPRCIvFgm984xv+Y5IkwW63Y9euXQHf1PKZJoaamhrs2LFj1Ayx0+lEVlZWwOwwn2lyE2RZlmM9CIosu93u/4Pa7XbDbDbzR0BxrKioaMwf3/mMVSvHZ5x4qqqq0Nzc7J+xMhqNWL9+PbZv3x7wOT7bxOBwOHDo0CHo9Xq0t7fjc5/73LidTPhM458oirDb7QHHTCbTqPI3Hz7T5MSQTEREREQ0AsstiIiIiIhGYEgmIiIiIhqBIZmIiIiIaASGZCIiIiKiERiSiYiIiIhGYEgmIiIiIhqBIZmIiIiIaASGZCIiIiKiERiSiYiIiIhGYEgmIopDDocDZWVlkCQp1kMJSqKNl4hoMgzJRERxyO12w+VywePxxHooQUm08RIRTUaQZVmO9SCIiCg4DocDJpMJBoNhRt6fiChaOJNMRJRAbDYbXC7XjL0/EVG0MCQTEREREY3AkExEFIfsdjtyc3P9r7dt2+Z/bbPZkJubi9zcXIii6P+My+VCaWkpcnNzUVRUBIfDEXBNh8OBoqIiAEBVVRWKiorgdDoBAKIowmazoaioCLm5uSgtLfW/F8z9R47XRxRFbNu2DevWrUNRURFsNtuoz/jGJUlSwGdHjp+IKJoYkomIEsDDDz+MpqYmAEBFRQWamprQ1NTkrw12Op0oKytDSUkJamtrUV5ejurq6oCgK0kSPB4PysrKcPDgQRQXF8NkMvnPb29vR0VFBerq6rB+/XpYrVZ/acVk9x+L2+32h/KamhpUVFSgpaXFH4iHj0sURZSWluLOO+9ETU0NtFotbDZbwDcBRETRpIr1AIiIaHJarRZardb/zyPD6Y4dO1BRUYHy8nL/MYPBAJvNBrPZ7D/mC8pHjhwJOH/4eQBgNBrR0NCAAwcOwGQyTXr/sVitVhQXF2P37t3+Y2azGUVFRaiurkZlZWXA5ysqKvxj3bVrF0pLS+FyuWCxWCa9FxFRuHEmmYgowbndbkiShOrqan8ZRG5uLsrKysacia2oqAjqugaDAWfOnAl5TKIoYtOmTaPeKy4uRn19/ajjvllt372JiGKJM8lERAnO15u4qakJOp1u0s8PD6PDud1uOBwOuFwuf7g2Go0hjcl3/lhhd/HixZAkCZIk+WenAQT8MxFRrHEmmYgowRUUFAAYCqa+sojhv4Jhs9mwefNmGAwG1NTU4MiRI+OG6WD4wvFYM9ltbW0AGIqJKL4xJBMRJZiRWz9rtVqYTCZUV1eP+mwwC98kSYLD4fDXNBuNRmi12nF3zwtm62nfNex2+6hzfYsGiYjiGcstiIgSiMFgQH19PYxGIw4dOoTCwkKYTCZUVlaitLQUZWVl/kV4TqcT9fX1oxbpjeSbcbbb7dDpdDAYDHA4HHC73aPKLca7/1hqampQVlaGbdu2YdOmTf66aWCoWwYRUTzjTDIRUQIpLy/39zQWRdFf1mAwGFBXV4esrCxYrVZYrVa0t7dj//79QV23pqYGwFBHih07dsBgMMBisSAnJyeo+4/FZDKhrq4OHR0dsFqtqK6uhslkwpNPPslSCyKKe4Isy3KsB0FEREREFE84k0xERERENAJDMhERERHRCAzJREREREQjMCQTEREREY3AkExERERENAJDMhERERHRCAzJREREREQjMCQTEREREY3AkExERERENAJDMhERERHRCAzJREREREQjMCQTEREREY3AkExERERENAJDMhERERHRCAzJREREREQj/H8+w8mNd0APVgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 720x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mf.plot_loss_all_info(info, labels=['MLR bottom', 'MLR uniform', 'MLR top'], figsize=(6, 4), dpi=120, \\\n",
    "                        logscale=False, hline=False, ylim=[0.15, 0.26], filename=\"plots/loss_%s.pdf\"%filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "75ec43f6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LR             0.16178478761767232\n",
      "LR+D           0.1537922868665804\n",
      "MLR uniform    0.15379228615265847\n",
      "MLR bottom     0.15379228605054865\n",
      "MLR top        0.15379228616060095\n"
     ]
    }
   ],
   "source": [
    "for i, algo in enumerate(info.keys()):\n",
    "    print(f\"{algo:14s} {info[algo]['loss'][-1]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "996f34fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('outputs/info_%s.pickle'%filename, 'wb') as handle:\n",
    "    pickle.dump(info, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "64b4a692",
   "metadata": {},
   "outputs": [],
   "source": [
    "rh_top = np.stack(ranks_history_top, axis=0)\n",
    "rh_bot = np.stack(ranks_history_bottom, axis=0)\n",
    "rh_uni = np.stack(ranks_history_uni, axis=0)\n",
    "\n",
    "with open('outputs/ranks_%s.pickle'%filename, 'wb') as handle:\n",
    "    pickle.dump([rh_bot, rh_uni, rh_top], handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "47fa0fd2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlUAAAMPCAYAAAD8ZCCXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABXmklEQVR4nO3db2wbd57n+Q+dONOdlorq+dMzPePSLDDYc0TKjcvsqtcq4bC4iW3VGLsHtDA20092YnccP5O9sHRPOqZs2XnSkoGRnjnq2JgAhwmTHu7i7hDTiXOHfaBiI8Z240akNB7c4K5d6ruZoG/aLGszvUnHvAcaVYsmZZPFkqpIvl9AA2aRVfWNOrA+qfrVh4lqtVoVAAAA2rIv6gEAAAC6AaEKAAAgBIQqAACAEBCqAAAAQkCoAgAACAGhCgAAIASEKgAAgBAQqgAAAEJAqAIAAAjBs1Ge3PM83bp1S5JULpeVTqeVyWRqPrO0tCTTNOW6rkzTlG3bUYwKAADwRJGGqtdff13T09MyTVOSdPDgQZmmKcuyJEnZbFa2bfuvJycnZZqm0ul0ZDMDAAA0Euntv/X1dTmO4782DEPlctl/ncvl/EAlSWNjY8rlcns6IwAAQDMivVKVz+f9P3ueJ8/z/BDlOI4Mw6jbZ3sIa8WPf/xjVatV7d+/P9iwAAC06PPPP1cikdCLL74Y9SjYA5GGqi2e5+n111/X7Oysf2vP87y6zyWTSVUqlUDnqFarqlarbc0JAEAr+L3TWyIPVVuL1QcGBpRMJv3tlUql5rW0eXuwUdhqxtYVqt8u3NZn7nrT+z3/L/5Qv/6n/05XPvwr/eTnG03t868Gf0tnDv/zlvZhv+jPxX7Rn6vb9+uEGTtlv06YUZL+9J/9mga+zB2SXhF5qDIMw3/ib2RkRJJk23bDq1Ke5zW8JdiKz9x1ffa3f9v05/cfOCBJ+snPN/Q3P3vY1D6DA19peR/2i/5c7Bf9ubp9v06YsVP264QZJemXg881/Vl0vlj1VA0PD+vNN9+U1PiqVKOrVwAAAHEQWagql8saGRmR67r+tv7+fv/19qf+tmxfyA4AABAnkYUqwzBkmmbNlafV1VWdPHnSfz0+Pl7ztJ/jOHXloAAAAHEQ2Zoq0zT12muvKZfL+f1U4+Pjmp6e9j+zuLioubk5ua4rz/OUyWQiK/78/a/2Nf3ZrxvP7+IkAAAgjiJdqN7MV85sD1mRePYZffGoqotHvxHtHAAAINYif/ov9n75hZ7Zl9A//Pnb+vzvP2lqly8NvaDkv/03uzwYAACIE0JVkz79zz9qqYqBUAUAQG+JVaUCAABApyJUAQAAhIBQBQAAEAJCFQAAQAhYqL6LWum2kui3AgCgkxGqdgPdVgAA9BxC1W4I0G0l0W8FAEAnI1Ttola7rST6rQAA6FQsVAcAAAgBoQoAACAEhCoAAIAQsKYqhqhiAACg8xCq4oQqBgAAOhahKk6oYgAAoGMRqmKIKgYAADoPC9UBAABCQKgCAAAIAaEKAAAgBKyp6iJUMQAAEB1CVTegigEAgMgRqroBVQwAAESOUNVFqGIAACA6kYYqz/N0/fp1SdLq6qoymYxs2/bfd11XnucpnU7L8zyVSiVZlhXVuAAAADuKNFTNz89rdnZW0mbAGhkZUT6fVzqdliQVCgXNz8/7n5+amiJUAQCAWIosVLmu61+JMgxDhmHIsixdv35di4uL/ufu3LnjX60CAACIq0ivVJVKJVUqFRmGIUkyTVOlUqnmM6ZpRjFaT6GKAQCA9kUWqkzT1N27d2u2OY5Td3uvUChIklZWVnTo0KGaNVdoE1UMAACEJjZP/5XLZVUqFU1NTfnb0um0hoeHZRiGbNvWkSNHZJomtwLDQhUDAAChiU2ounjxovL5vH8rUFLdVatUKqVcLucvbkc4qGIAAKB9sfjuv7m5OV25cqVu/VQul6t53WjNFQAAQBxEHqpyuZyOHz/u39JzHEfS5tOB2WxWruv6n3348KEOHDgQyZwAAABPEmmo2gpQhmHIdV2Vy2WVy2VJm1elpqamaq5eOY6js2fPRjIrAADAk0S2psrzPJ06dapu+/aF6rZt+7cAy+WyZmdnWaQeI1QxAADwK5GFKsMwdO/evSd+xjRNZTKZPZoITaOKAQCAOrF5+g8dhCoGAADqEKoQGFUMAAD8SuRP/wEAAHQDQhUAAEAICFUAAAAhIFQBAACEgIXq2HOt9FvRbQUA6BSEKuwd+q0AAF2MUIW9E6Dfim4rAECnIFRhz7Xab0WoAgB0AhaqAwAAhIBQBQAAEAJCFQAAQAhYU4WO0EoNg0QVAwBg7xGqEG/UMAAAOgShCvEWoIZBoooBALD3CFXoCK3WMEhUMQAA9lbgUFUsFjUzMyPXdRu+v7a2FngoAACAThM4VGWzWaVSKV2+fFnJZDLMmQAAADpO4FBVqVT0xhtvqK+vtaeyAAAAulHgUGXbtkqlkg4fPhzmPECoqGIAAOyVwKFqcHBQp06dUiqV0vDwsAzDqHn/woULbQ8HBEYVAwBgjwUOVSsrKxodHZWkusXqiUSivamAdlHFAADYY4FD1cLCQtsn9zxP169flyStrq4qk8nItu2azywtLck0TbmuK9M0694HnoQqBgDAXmm7p6pYLMpxHEnSN77xDR09erTpfefn5zU7OytpM2CNjIwon88rnU5L2nzC0LZtWZYlSZqcnJRpmv77AAAAcdHWFyqfPn1a586dU7lcVrlc1ne/+139yZ/8SVP7uq4r13XleZ4kyTAMWZblX7mSpFwu5wcqSRobG1Mul2tnZAAAgF0R+ErVzMyM+vv79fHHH9dsn5yc1KVLl3Tp0qWnHqNUKqlSqfiL3E3TVKlUkiQ5jlO3+H1rOwAAQNwEDlWFQkEfffRR3farV6/q6NGjTw1Vpmnq7t27Ndscx/GvTG1dwdoumUyqUqkEHRloGlUMAIBWBQ5VhmGoUqnUlX9WKhVVq9WWj1cul1WpVDQ1NeUf5/GmdsMwGoYtIDRUMQAAAgocqo4dO6ZsNquFhQU/WG1sbOj8+fPKZDItH+/ixYvK5/P+Lb9GV6U8z2t4SxAIDVUMAICAAoeq6elpTU5OamRkRKZpStpcfG5ZVsvFn3Nzc7py5Yp/HKnxValGV6+A3UAVAwCgVW1VKiwuLvpP/nmeJ8uylEqlWjpGLpfT8ePH/ZqErXVV25/627J1DgAAgLhpu6cqnU4H7o3aepLPMAy/XqFcLvvBaXx8vGbxuuM4/porAACAOGk6VM3MzMi2bf+raa5du/bEzz/tFqDneTp16lTd9u2haXFxUXNzc37gymQyFH8CAIBYajpUraysaGxszH9dLpd3/Gwz3/1nGIbu3bv31M9NT083NyAQA1QxAEDvajpU5fP5mtc3btwIfRigY1HFAAA9L/Caqo2NjbqOKklaX1+XJB04cCD4VECnoYoBAHpe4FA1MjKitbW1uu2u6+r73/++3nrrrbYGAzoRVQwA0LsCf6HyTq3pw8PD/vf3AQAA9IqWr1QdPXpUiURCiURCx44dq3vfdd2Wu6oAAAA6XcuhanZ2VtVqVadPn25Ym2CaJqEKAAD0nJZD1VZP1fj4uMbHx0MfCAAAoBMFXqi+sLAQ5hxAT6PfCgA6X1tfU7O2tuZ/1cx2hmHoxIkT7Rwa6A30WwFA1wgcqm7fvq1z584plUppbW1NQ0NDkqTV1VWNjY0RqoBm0G8FAF0jcKi6du2abt68qdHRUU1MTOjtt99WX1+fCoWCisVimDMCXY9+KwDofIF7qu7fv+8vWjdN0++msm1bhUIhnOkAAAA6ROBQtXXbT5Isy9LS0pIk6d1335XneeFMBwAA0CEC3/47e/asSqWShoaGlMlklMvlNDQ0pGq1qldffTXMGQEAAGIvcKh6vKMqn8/LdV0NDAyov7+/7cEAPB1VDAAQH4FDVbFY1MOHD2u+qsY0TV27dk3Hjx/3nwYEsAuoYgCA2Akcqubn5zU1NVW3fXh4WPPz83rrrbfaGgzAE1DFAACxEzhUlctlHTp0qG67ZVk6f/58OzMBaBJVDAAQH4Gf/jNNU+vr63XbHzx4oAMHDrQ1FAAAQKcJHKpOnjypyclJ/fSnP/W3ra+v6/z587JtO5ThAAAAOkXg239nzpyR67p66aWXlEwmJUme5+nkyZO6cOFCaAMCAAB0gra+UHl2dlZnzpzR6uqqpM1CUNM0QxkMwO5ppYqBGgYAaE5boUravOW39RU1+/btI1QBcUYVAwDsmrZC1enTp1UqlTQ8PCxJyuVyun79un7wgx+EMhyAkAWoYqCGAQCaEzhUzczMqL+/Xx9//HHN9snJSV26dEmXLl166jFc19U777yjwcFBZTKZuvc8z1M6nZbneSqVSrIsK+i4ALZptYqBUAUATxf46b9CoaA33nijbvvVq1d169atp+7vOI7K5bJWV1cbfgFzoVDQxMSEDh48qJGREZXL5aCjAgAA7LrAV6oMw1ClUlFfX+2C10qlomq1+tT9t646vf/++zt+5s6dO/7VKgAAgDgLfKXq2LFjymaz2tjY8LdtbGzo/PnzdbfygjJNk0AFAAA6QuArVdPT05qcnNTIyIj/xJ/rurIsK7SeqkKhIElaWVnRoUOHKBUFItRKDYNEFQOA3tPW03+Li4sql8sql8vyPE+WZSmVSoUyWDqd1vDwsAzDkG3bOnLkCFeugChQwwAATWm7pyqdTu9K0Hn8Sb9UKqVcLqfZ2dnQzwXgCQLUMEhUMQDoPU2HqmvXrrV04HZvAeZyuZq1WaZpqlgstnVMAMG1WsMgUcUAoLc0HapaqTRIJBKBhtniuq6y2awsy/LXaz18+FAHDhxo67gAAAC7pelQdePGjd2co4Zpmpqamqr5yhvHcbSwsLBnMwAAALSi7TVVQTmOo+XlZRWLRa2ururBgwd6+eWX/SBl27ZyuZykzatks7OzLFIHAACxFdmaKsuyZFmWpqenG75vmmZofVcAokMVA4BeEcs1VQC6AFUMAHpMLNdUAegCVDEA6DGRrakC0BuoYgDQK9oKVWtra3Icp267YRg6ceJEO4cGAADoKIFD1e3bt3Xu3DmlUimtra1paGhIkrS6uqqxsTFCFQAA6CmBQ9W1a9d08+ZNjY6OamJiQm+//bb6+vpUKBRoPgcAAD0ncKi6f/++RkdHJW3WH5RKJR0+fFi2bWtmZkaXL18ObUgAvYcqBgCdJnCo2n7bz7IsLS0t6fDhw3r33XfleV6YMwLoJVQxAOhQgUPV2bNnVSqVNDQ0pEwmo1wup6GhIVWrVZ05cybMGQH0EqoYAHSowKFqfHy85nU+n5fruhoYGFB/f3/bgwHobVQxAOg0ofZUbf8CZAAAgF6yL+oBAAAAugGhCgAAIASEKgAAgBDw3X8Augr9VgCisiuhqlgs+sWgALAn6LcCELHAoWpmZkbf/va39cILL/jbNjY29L3vfU/vvfee1tbWQhkQAJpCvxWAiLXVqP6tb31LN2/e1OHDh1UsFnXu3DkdOnRIH3zwQZgzAkDT6LcCEJXAoSqTyUiSXnnlFY2NjclxHF2+fFknT54MbTgAAIBO0daaqkwmI9M0dfr0aU1PTxOoAABAz2o6VL311ls7vmdZlubn53X//n0NDg5Kkr7zne+0Px0AAECHaDpULS8vP/H90dFRua4r13WVSCQIVQA6ClUMANrVdKi6cePGbs4BANGgigFASCj/BNDbqGIAEJLAoapYLGpmZkau69Zsr1arSiQSTfVUua6rd955R4ODg/7ThNstLS3JNE25rivTNGXbdtBxAeCJqGIA0K7AoSqbzSqVSuny5ctKJpMt7+84jjzP0+rqqgYGBhoe37ZtWZYlSZqcnJRpmkqn00FHBgAA2DWBQ1WlUtEbb7yhvr7WFndu2QpL77//fsP3c7mcZmdn/ddjY2N12wAAAOJiX9AdbdtWqVQKcxaf4zgyDKPhdgAAgDgKfKVqcHBQp06dUiqV0vDwcF0IunDhQuChPM+r25ZMJlWpVAIfEwB2QytVDNQwAN0tcKhaWVnR6OioJNUtVk8kEm0NValU6tZpGYbRMGwBQCSoYgDwmMChamFhIcw5ajS6KuV5XsNbggAQiQBVDNQwAN0tlj1Vja5KNbp6BQBRa7WKgVAFdK+2QtXa2lrDxeOGYejEiROBj7v1ZOB2nuc13A4AABAHgUPV7du3de7cOaVSKa2trWloaEiStLq6qrGxsbZClSSNj4/LcRw/SDmOo6mpqbaOCQAAsFsCh6pr167p5s2bGh0d1cTEhN5++2319fWpUCioWCw+dX/HcbS8vKxisajV1VU9ePBAL7/8skzTlCQtLi5qbm5OruvK8zxlMhmKPwEAQGwFDlX379/3n/4zTVOlUkmHDx+WbduamZnR5cuXn7i/ZVmyLEvT09M7fuZJ7wFAp2qlhkGiigHoFIFD1fbbfpZlaWlpSYcPH9a7775L9QEANEINA9DVAoeqs2fPqlQqaWhoSJlMRrlcTkNDQ6pWqzpz5kyYMwJAdwhQwyBRxQB0isChanx8vOZ1Pp+X67oaGBhQf39/24MBQLdqtYZBoooB6ASBv/uv0WJ00zQJVAAAoCcFDlXZbFYffvhhmLMAAAB0rMCh6tVXX9Xc3Jw2NjbCnAcAAKAjBV5TlUgk1NfXp5deekmjo6N+v9SWCxcutD0cAOBXqGIA4i1wqFpeXtbAwIAGBgbkeZ7K5bL/XiKRCGU4AICoYgA6ROBQtbCwEOYcAICdUMUAdIS2vlAZALB3qGIA4i3wQnUAAAD8CqEKAAAgBIQqAACAEBCqAAAAQsBCdQDocvRbAXuDUAUA3Yp+K2BPEaoAoFvRbwXsKUIVAHQ5+q2AvcFCdQAAgBAQqgAAAEJAqAIAAAgBa6oAAA1RxQC0hlAFAKhFFQMQCKEKAFCLKgYgkFiHKtd15Xme0um0PM9TqVSSZVlRjwUAPYEqBqA1sV6oXigUNDExoYMHD2pkZETlcjnqkQAAABqK9ZUqSbpz545/tQoAACCuYh+qTNOMegQAAICnin2oKhQKkqSVlRUdOnRItm1HPBEA4EmoYkCvinWoSqfTGh4elmEYsm1bR44ckWma3AoEgDiiigE9Ltah6vEn/VKplHK5nGZnZyOaCACwI6oY0ONiHapyuZwymYz/2jRNFYvFCCcCADwNVQzoVbGtVHBdV9lsVq7r+tsePnyoAwcORDgVAABAY7ENVaZpampqqubpP8dxdPbs2QinAgAAaCzWt/9s21Yul5Mklctlzc7OskgdAADEUqxDlWmaNWuqAADdq5UqBmoYEEexDlUAgB5AFQO6BKEKABCtAFUM1DAgjghVAIBYaLWKgVCFuInt038AAACdhFAFAAAQAkIVAABACFhTBQDoWK3UMEhUMWB3EaoAAJ2HGgbEEKEKANB5AtQwSFQxYHcRqgAAHavVGgaJKgbsHhaqAwAAhIBQBQAAEAJCFQAAQAgIVQAAACFgoToAoOfQb4XdQKgCAPQO+q2wiwhVAIDeQb8VdhGhCgDQc+i3wm5goToAAEAICFUAAAAhIFQBAACEgDVVAAA0qdUqhmf3JXZpEsQRoQoAgKcJWMWwsrKySwMhjghVAAA8TcAqhkf/3Zj29bV2dQudK/ahamlpSaZpynVdmaYp27ajHgkA0KNarWKofvNfSoSqnhHrUJXNZmXbtizLkiRNTk7KNE2l0+mIJwMAAKgV66f/crmcH6gkaWxsTLlcLsKJAAAAGottqHIcR4ZhNNwOAAAQN4lqtVqNeohGCoWCLl68qLt37z5xW7N+9KMfqVqtan8iIbXyj7xvnxLPPqvqZ59Jjx41t88zzyixf39r+7Bf9Odiv+jP1e37dcKMnbJfJ8wo6fN9+5RIJPSHf/iHzZ8LHSu2a6oqlYqSyWTNNsMw5HleoOMlEptdIYn9+4Pt/9xze7IP+0V/LvaL/lzdvl8nzNgp+8V9xsTnn/u/f9D9YhuqksmkKpVKzTbP8xreEmzGiy++GMZYAAAADcV2TVWjq1KNrl4BAADEQWxD1fan/rZ4ntdwOwAAQNRiG6okaXx8vOZpP8dxlMlkIpwIAACgsdg+/bdlbm5Og4OD8jyPRnUAABBbsQ9VAAAAnSDWt/8AAAA6BaEKAAAgBIQqAACAEBCqAAAAQkCoAgAACAGhCgAAIASEKgAAgBAQqgAAAEJAqAIAAAgBoQoAACAEhCoAAIAQEKoAAABCQKgCAAAIAaEKAAAgBIQqAACAEBCqAAAAQkCoAgAACAGhCgAAIASEKgAAgBAQqgAAAELwbJQn9zxPt27dkiSVy2Wl02llMpmazywtLck0TbmuK9M0Zdt2FKMCAAA8UaSh6vXXX9f09LRM05QkHTx4UKZpyrIsSVI2m5Vt2/7ryclJmaapdDod2cwAAACNRHr7b319XY7j+K8Nw1C5XPZf53I5P1BJ0tjYmHK53J7OCAAA0IxIr1Tl83n/z57nyfM8P0Q5jiPDMOr22R7CWvHjH/9Y1WpV+/fvDzYsAAAt+vzzz5VIJPTiiy9GPQr2QKShaovneXr99dc1Ozvr39rzPK/uc8lkUpVKJdA5qtWqqtVqW3MCANAKfu/0lshD1dZi9YGBASWTSX97pVKpeS1t3h5sFLaasXWF6j/+XVU/+flG0/v9q8Hf0pnD/1xXPvyrpvcLsg/7RX8u9ov+XN2+XyfM2M5+v//VPl08+o2mP98LVlZWoh4BeyjyUGUYhv/E38jIiCTJtu2GV6U8z2t4S7AVP/n5hv7mZw+b/vzgwFda3i/IPuwX/bnYL/pzdft+nTBjO/sBvS5WPVXDw8N68803JTW+KtXo6hUAAEAcRBaqyuWyRkZG5Lquv62/v99/vf2pvy3bF7IDAADESWShyjAMmaZZc+VpdXVVJ0+e9F+Pj4/XPO3nOE5dOSgAAEAcRLamyjRNvfbaa8rlcn4/1fj4uKanp/3PLC4uam5uTq7ryvM8ZTIZij8BAEAsRbpQvZmvnNkesgAAAOIqVgvVAQAAOhWhCgAAIASEKgAAgBAQqgAAAEJAqAIAAAgBoQoAACAEhCoAAIAQEKoAAABCQKgCAAAIAaEKAAAgBIQqAACAEBCqAAAAQkCoAgAACAGhCgAAIASEKgAAgBAQqgAAAEJAqAIAAAgBoQoAACAEhCoAAIAQEKoAAABCQKgCAAAIAaEKAAAgBIQqAACAEBCqAAAAQvBslCf3PE/Xr1+XJK2uriqTyci2bf9913XleZ7S6bQ8z1OpVJJlWVGNCwAAsKNIQ9X8/LxmZ2clbQaskZER5fN5pdNpSVKhUND8/Lz/+ampKUIVAACIpchCleu6/pUowzBkGIYsy9L169e1uLjof+7OnTv+1SoAAIC4ivRKValUUqVSkWEYkiTTNFUqlWo+Y5pmqOf8/a/2tfT5rxvPh3p+AADQnSILVaZp6u7duzXbHMepu71XKBQkSSsrKzp06FDNmqsgLh79Rlv7AwAANBLplartyuWyKpWKpqam/G3pdFrDw8MyDEO2bevIkSMyTbOtW4H/6f/9P/Tgv240/fkDX/kt/Yvf+m8Cnw8AAPSG2ISqixcvKp/P+7cCJdVdtUqlUsrlcv7i9iD+1vt/9Pf/+POW9iFUAQCAp4lFT9Xc3JyuXLlSt34ql8vVvG605goAACAOIg9VuVxOx48f92/pOY4jafPpwGw2K9d1/c8+fPhQBw4ciGROAACAJ4k0VG0FKMMw5LquyuWyyuWypM2rUlNTUzVXrxzH0dmzZyOZFQAA4EkiW1PleZ5OnTpVt337QnXbtv1bgOVyWbOzs/RVAQCAWIosVBmGoXv37j3xM6ZpKpPJ7NFET9ZKvxXdVgAA9J7YPP0XV/sSCX3xqEq/FQAAeCJC1VM8qlb1zL5ES/1WdFsBANB7CFVNarXfilAFAEBvibxSAQAAoBsQqgAAAEJAqAIAAAgBa6p2USs1DBJVDAAAdDJC1S6ghgEAgN5DqNoFQWoYJKoYAADoZISqXdRqDYNEFQMAAJ2KheoAAAAhIFQBAACEgFAFAAAQAtZUxRBVDAAAdB5CVYxQxQAAQOciVMUIVQwAAHQuQlUMUcUAAEDnCRyqisWiZmZm5Lpuw/fX1tYCDwUAANBpAoeqbDarVCqly5cvK5lMhjkTAABAxwkcqiqVit544w319bX2pBoAANg9X3zxhT7//POox+gK+/fv1zPPPNP05wOHKtu2VSqVdPjw4aCHQMioYgCA3lWtVvV3f/d3evDgQdSjdJWBgQH9zu/8jhKJxFM/GzhUDQ4O6tSpU0qlUhoeHpZhGDXvX7hwIeih0SKqGAAAW4Hqa1/7mp5//vmmQgB2Vq1W9emnn+qTTz6RJH39619/6j6BQ9XKyopGR0clqW6xOv9H7i2qGACgt33xxRd+oPqN3/iNqMfpGl/+8pclSZ988om+9rWvPfVWYOBQtbCwEHRX7BKqGACgN22toXr+eZZ1hG3rZ/r555/vXqjaUiwW5TiOJOkb3/iGjh492vS+nufp+vXrkqTV1VVlMhnZtl3zmaWlJZmmKdd1ZZpm3fsAAGATd4rC18rPtK1Qdfr0aZVKJQ0PD0uScrmcrl+/rh/84AdN7T8/P6/Z2VlJmwFrZGRE+Xxe6XRa0mZtg23bsixLkjQ5OSnTNP33AQAA4mJf0B1nZmbU39+vjz/+WDdu3NCNGzf08ccf63d/93d16dKlp+7vuq5c15XneZIkwzBkWZZ/5UraDGlbgUqSxsbGlMvlgo4MAEDP+eJRtafOG6XAV6oKhYI++uijuu1Xr17V0aNHmwpWpVJJlUrFf3LQNE2VSiVJkuM4dU8Ubm1vx29+qbWi0oFf6/4eLqoYAKB7PbMvoSsf/pV+8vPmH2Rq1+9/tS/wE+mu6yqbzapSqSifz4c8We153nnnHQ0ODiqTyYRyzMChyjAMVSqVuvLPSqWiavXp6dQ0Td29e7dmm+M4/pWprStY2yWTSVUqlaAjS5L+h9+3nv6hHkEVAwD0hp/8fEN/87OHUY/RlK310+VyedfO4TiOPM/T6uqqBgYGQjtu4FB17NgxZbNZLSws+MFqY2ND58+fD5T4yuWyKpWKpqamJG2Gs8e//sYwjIZhqxWf/ONdff5F8/9iffnZ39avf6k713BRxQAAiKPl5WUdP358146/dQHn/fffD/W4gUPV9PS0JicnNTIyItM0JW1eSrMsK1Dx58WLF5XP5/1bfo2uSnme1/CWYCv+y+fr+sUX/1+Le3VnqNpCFQMAIE6KxaKuXr0a9Rgta+vpv8XFRZXLZZXLZXmeJ8uylEqlWj7O3Nycrly54oczqfFVqUZXrwAAQPdwXVfJZLLtiyhRaLunKp1Ot1VxkMvldPz4cf8YW+uqtj/1t2UruAEAgO60fX31TrLZ7FOPY5qmzpw5E9ZYTWk6VM3MzMi2bf+raa5du/bEzzdzC3DrST7DMPx6hXK57P8wx8fHa364juP4a64AAED3aWY91VbHZdw0HapWVlY0Njbmv37Sqvxm2kc9z9OpU6fqtm8PTYuLi5qbm/MDVyaTofgTAIAuViwWNT09HfUYgTQdqh7virhx40ZbJzYMQ/fu3Xvq5+Lyg22l36oXuq0k+q0AAOFyXVeSar6erpGOv/33uI2NjbqOKklaX1+XJB04cCD4VLGyT4+qj+i32oZ+KwDoLK3+R3BU53Nd1//qu3K5vGOo6vjbf48bGRnR2tpa3XbXdfX9739fb731VluDxccj7Uvsa6nfqpu7rST6rQCgk0T1H8FfPNr8XdEKy7JUKBSUy+VCazlvxHEcLS8vq1gsanV1VQ8ePNDLL7+8Y4hrVuBQtVNr+vDwsP9VM92k9X6r7g1VW+i3AoD4azXYRH3evbgKtdUyEPYSo5ZD1dGjR5VIJJRIJHTs2LG6913XDdRVBQAA0MlaDlWzs7OqVqs6ffp0w9oE0zQJVQAAoOe0HKq2eqrGx8c1Pj4e+kAAAACdKPCaqoWFhTDn6Eqt1DBIVDE0Qg0DAKBTtPU1NWtra34r+naGYejEiRPtHLrDUcPQCFUMAIBuFjhU3b59W+fOnVMqldLa2pqGhoYkSaurqxobG+vxUNV6DYNEFUMj1DAAADpF4FB17do13bx5U6Ojo5qYmNDbb7+tvr4+FQoFFYvFMGfsWK3XMEhUMdQjVAEAOsG+oDvev3/fX7RumqbfTWXbtgqFQjjTAQAAdIjAoWrrtp+0WaK1tLQkSXr33XfleV440wEAgLY82qGsu1vPG6XAt//Onj2rUqmkoaEhZTIZ5XI5DQ0NqVqt6tVXXw1zRgAAENC+REL/808c/ewXlT07529+KRn4YS3XdZXNZlWpVJTP50OebJPnebp+/bqkzbXgmUxGtm23fdzAoerxjqp8Pi/XdTUwMKD+/v62B+tlVDHUa/XLOaliAIBf+dkvKi1/rVhUTNOUbdsql8u7do75+Xn/63A8z9PIyIjy+bzS6fbWNQcOVcViUQ8fPqz5qhrTNHXt2jUdP37cfxoQraCK4XHUMABA71leXtbx48d35diu68p1XXmeJ8MwZBiGLMvS9evXtbi42NaxA4eq+fl5TU1N1W0fHh7W/Py83nrrrbYG601UMTwuSA2DRBUDAHSyYrGoq1ev7trxS6WSKpWKDMOQVPvAXTsCh6pyuaxDhw7VbbcsS+fPn29npp5HFUO9VmsYJKoYAKATua6rZDLpB56wmaapu3fv1mxzHEeW1f5dosChyjRNra+v64UXXqjZ/uDBAx04cKDtwQAAQO9pJuBks9mnHsc0TZ05c+apnyuXy6pUKg3vvrUqcKg6efKkJicndfPmTf3e7/2eJGl9fV3nz58PZQU9AADoPc2sp9paZB6GixcvKp/Ph3JlLHCoOnPmjFzX1UsvvaRkcvNpNc/zdPLkSV24cKHtwQAAQO8pFouanp7ek3PNzc3pypUrMk0zlOO19YXKs7OzOnPmjFZXVyVtFoKGNRhaRxVDPaoYAKBzuK4rafPWneu6O2aKMG7/5XI5HT9+3K9RCGNdVVuhStq85be1Yn7fvn2EqkhQxfA4qhgA4Fda/Y/uqM7nuq6Gh4clba512ilTtHv7z3EcSZJhGH69QrlcjjZUnT59WqVSyf8B5HI5Xb9+XT/4wQ/aGgqtoorhcVQxAMCmR9VqJP/R/aha1b5EoqV9LMtSoVBQLpdTJpPZlbk8z9OpU6fqtke6UH1mZkb9/f36+OOPa7ZPTk7q0qVLunTpUruzoUVUMdSjigFAr2s12ER93jAXoTdiGIbu3bu3K8cO/IXKhUJBb7zxRt32q1ev6tatW00dw3Vdzc3NKZfLNXxvq6Le8zz/Uh0AAEAcBQ5VhmGoUqn/csZKpaJqE99M7TiOyuWyVldX5Xle3fuFQkETExM6ePCgRkZGdvU7gAAAANoV+PbfsWPHlM1mtbCwoL6+zSesNjY2dP78+abug24tBnv//fd3/MydO3fkeV7bX3AIAACw2wKHqunpaU1OTmpkZMRfne+6rizLCq2najeeJHzumYGWPr//mf7QZ4gbqhjqUcUAAGhVW0//LS4uqlwuq1wuy/M8WZalVCoV1mwqFAqSpJWVFR06dCiUpvbf+8p/3/YxugdVDI+jigEAEFTbPVXpdHpXbs+l02kNDw/LMAzZtq0jR47INM22z1X9+/9N+qyFp8GeN5X4jW+2dc74oorhcVQxAACCajpUXbt2raUDt3sL8PECrlQqpVwu1/6jlhv/p/SLv2ttn64NVZuoYqhHFQMAoFVNh6pWnr5LhNCJ8Xjxl2maKhaLbR8XAABgNzQdqm7cuLGbc9RwXVfZbFaWZfmL1R8+fKgDBw7s2QwAAACtaHtN1W4wTVNTU1M1T/85jqOFhYUIpwIAoPNUq4+USASupdzz825dWKlUKsrn87sw2Wap+FZReblcVjqdDuVrcSJbU+U4jpaXl1UsFrW6uqoHDx7o5Zdf9oOUbdt+03q5XNbs7GykfVWtVDH0Qg2DRBVDI1QxAIibRGKffvpf/nd99sWDPTvnc88MBH7a3jRN2ba9q6Xfr7/+uqanp/3McfDgQZmmuXdfqBz2mirLsmRZlqanpxu+b5rmrn2ZYksS+1StPqKKoQZVDI+jigFAnH32xYMADyRFZ3l5WcePH9+146+vr8txHD9nGIahcrm8d6FqL9dUxco/Xb5sqYqhq2sYJKoY6lHFAADhKRaLunr16q4df/ttRc/z/K7NdsVyTVUstVrF0NWhahNVDPWoYgCA9riuq2QyKcMwdv1cnufp9ddfD22JUVuham1tTY7j1G03DEMnTpxo59AAAKAHOY7z1KtG2Wz2qccxTVNnzpzZ8f2txeoDAwNKJltbI7yTwKHq9u3bOnfunFKplNbW1jQ0NCRJWl1d1djYGKEKAAC0rJn1VG0XgWvzAtDWmqqRkRFJavvr8AI/Y3nt2jXdvHlT+XxeQ0NDevvtt5XP5/Vnf/Zn9EkBAIBAisXinj/tPzw8rDfffLPt4wS+UnX//n2Njo5K2rzEViqVdPjwYdm2rZmZGV2+fLnt4QAAQO9wXVfSZq5wXbemr3K7dm7/lctlvfLKK8rn8/7x+/v7VSqV2ph8U+BQtf22n2VZWlpa0uHDh/Xuu+/K87y2B+sGrXRbSfRbNdIL3VYS/VYAdlerv4+iOp/ruhoeHpa0GX52ClXt3P4zDEOmadaso1pdXdXJkycDH3NL4FB19uxZlUolDQ0NKZPJKJfLaWhoSNVq9YkLw3oC3VY7oN/qcfRbAdhtUf0+CtKoblmWCoVC3ff/hsk0Tb322mvK5XJ+P9X4+PiOvZmtCByqxsfHa17n83m5rquBgQH19/fGFZcdBem2kui3aqCbu60k+q0A7L4ovqKmnfOGsQj9adpdkL6TUHuqdrpM17Na7baS6LdqqHtD1Rb6rQCg80UTXwEAALoMoQoAACAEhCoAAIAQ8N1/MUQVQ71Wahgkqhh2QhUDAOyeXQlVxWLRLwZFC6hiaIAahkaoYgCA+AkcqmZmZvTtb39bL7zwgr9tY2ND3/ve9/Tee+9pbW0tlAF7ClUMDbRewyBRxbATqhgAYPe01aj+rW99Szdv3tThw4dVLBZ17tw5HTp0SB988EGYM/YeqhjqtF7DIFHF0BihCgB2R+BQtdV0+sorr2hsbEyO4+jy5cuh1LwDAAB0mrbWVGUyGZmmqdOnT2t6eppABQBAzAT5uphOPm+Umg5Vb7311o7vWZal+fl53b9/X4ODg5Kk73znO+1PBwAA2pJI7FN1PS/915/t3Ul/7TeVODARaFfXdZXNZlWpVJTP50MerPH53nnnnb397r/l5eUnvj86OirXdeW6rhKJBKEqAlQx1KOKobFWqhioYQC6wH/9WetrdSNimqZs21a5XN6T82Wz2dC+Zq/pUHXjxo1QTohdQBVDA1QxNEIVA4BOsLy8rOPHj+/6eRzHCfV7iyn/7AZUMTRAFUMjQaoYqGEAsNeKxaKuXr266+dxXVemacp13VCOFzhUFYtFzczM1A1SrVaVSCToqYoCVQx1qGJorNUqBkIVgL3iuq6SyaQMw9jV8+RyOWUyGS0tLYV2zMChKpvNKpVK6fLly0omW1u3smVrcdjg4KBf0bDd0tKSnyC37rECAIDu5TiOLOvJSzey2exTj2Oaps6cOdPwPc/zQr3ttyVwqKpUKnrjjTfU1xdsYa/jOPI8T6urqxoYGKh7P5vNyrZt/wc7OTkp0zSVTnf/VQQAAHpVM+upZmdn2zrHrVu3Gl7MaVfgAgnbtlUqlQKf2LIs2bat/v7GT6DlcrmapDo2NqZcLhf4fAAAIP6KxeKuXkBp5kpYUIGvVA0ODurUqVNKpVIaHh6uu/d54cKFwEM5jtPwXqrjOIGP6fu132zt8899tf1zxhxVDPWoYqjXSg2DRBUDgNZtrdPevvSnkXZv/xUKBf/Pt27dkrS55CiTybS1litwqFpZWdHo6Kgk1S1WTyQSgQeSNu91Pi6ZTKpSqbR1XEmBy8i6ElUMDVDF8DhqGIAu0OoFhYjO57quhoeHJUnlcnnHUNXO7T/LsmquVD148ECu6+4YwFoROFQtLCy0ffKdVCqVusXvhmE0DFut+uXy/6Rq5ZOmP5/43YN69r/d/a6MSFDF0ABVDI8LUsMgUcUAxEW1+iiSCwpBvqbGsiwVCgX/ybzdtrS0pNu3b0uS5ubm2m5Vj2VPVaOrUp7nhfJ4ZfX//rGqn/xfre3UraFqC1UMdahiqNdqDYNEFQMQB1F9/17Q87a7CL0VZ86cCeUK1Za2QtXa2lrDdU6GYejEiROBj9voqlSjq1cAAABxEThU3b59W+fOnVMqldLa2pqGhoYkSaurqxobG2srVDVale953q6t1gcAAGhX4GuC165d082bN5XP5zU0NKS3335b+Xxef/Znf6YDBw60Pdj4+HjNVTDHcfbk/ioAAEAQga9U3b9/33/6zzRNlUolHT58WLZta2ZmRpcvX37i/o7jaHl5WcViUaurq3rw4IFefvllf6X/4uKi5ubm5LquPM9TJpOJtvizlScZeqCGQaKKoRGqGOpRxQCgVwQOVdtv+1mWpaWlJR0+fFjvvvtuU0/pbT3S+KSV9u2uwg/Fvmcie3IitqhiaIAqhsdRxQCg1wQOVWfPnlWpVNLQ0JAymYxyuZyGhoZUrVZDXUkfuUdfKJHY11IVQ1fXMEhUMTREFcPjqGIA0GsCh6rx8fGa1/l8Xq7ramBgYMevnulkLVcxdHOo2kIVQx2qGOpRxQCgVwReqF4sFuu2mabZlYEKAADgaQKHqmw2qw8//DDMWQAAQMiqjx711HmjFPj236uvvqq5uTmNjo6qr6/7n2ACAKATJfbt0y9v/Zmq/7C+d+f89QN69o/PB9rXdV1ls1lVKhXl8/lwB9t2Ds/zlE6n5XmeSqVSKF2YgUNVIpFQX1+fXnrpJY2OjtZ96eGFCxfaHg4AALSv+g/rrX9FW0RM05Rt2yqXy7t2jkKhoPn5ef/11NRUtKFqeXlZAwMDGhgYkOd5Nf/wiUSi7cG6Qqvf0k2/VZ1e6LaS6LdqhH4roHctLy/r+PHdfeDrzp07/tWqsAQOVQsLC6EN0XXotmqMfqsG6Ld6HP1WAIrFoq5evbqr53j8DlsY2vpCZewgQLeVRL9VQ13dbSXRb1WPfiugt7muq2QyKcMwdvU8hUJBkrSysqJDhw7Jtu22j0mo2kUtd1tJ9Fs10tWhahP9VvXotwJ6k+M4T13flM1mn3oc0zR3LCNPp9MaHh6WYRiybVtHjhyRaZpt3wokVAEAgNhoZj3V7OxsW+d4PLSlUinlcrm2jxu4pwoAACBsxWIx1MXjjeRyuZrXpmmqVCq1fVyuVAEAgFhwXVfSZshxXXfHxeTt3P7b6sGyLMs//sOHD3XgwIE2Jt9EqIojqhjqtFLDIFHFsBOqGOpRxYBekPj19gPDXpzPdV0NDw9Lksrl8o6hqp3bdKZpampqqubYjuOE0mpAqIoTqhjqUcOwA6oYHkcVA9BY9dGjwO3m7Z43sa+1VUaWZalQKCiXyymTyezSZJJt2/4twHK5rNnZ2VBuORKq4oQqhnpBahgkqhh2QBVDPaoY0O1aDTZRn7fdxeLNME1zV0IboSqGqGJooNUaBokqhh11Z6jaQhUDgKjw9B8AAEAICFUAAAAhIFQBAACEgDVV3YQqhjpUMTTWShVDL9QwSFQxAGgfoaobUMVQjyqGHVDF8DiqGACEhVDVDahiqEcVww5ar2Lo5hoGiSoGAOGJdahyXVee5ymdTsvzPJVKpad+c3Uvo4qhAaoYGmq9iqF7Q9UWqhgAtCvWC9ULhYImJiZ08OBBjYyMqFwuRz0SAAAdpfroUU+dN0qxvlIlSXfu3PGvVgEAgNYk9u3TJ3Pz+sxd37NzPmce0NempwLtu/WFx5VKRfl8PuTJas/zzjvvaHBwUJ7nybbtHb9rsFmxD1Xt/gMCANDrPnPX9dnf/m3UYzTFNE3Ztr2rd6dc19W5c+f80DY5OSlJOnPmTFvHjX2oKhQKkqSVlRUdOnRItm23dbxWvzk7kfxaW+frCFQx1KGKoV4rNQzSr6oY9nq/vdZKFQM1DEBzlpeXdfz47q35nZubq/nuv7Nnz4ZyESfWoSqdTmt4eFiGYci2bR05ckSmabZ1KzCKb+qOLaoY6lHF0EDwGoa93m8vUcUA7J5isairV6/u2vFv376ts2fP+q/DWmIU61D1+JN+qVRKuVyurW+w/oc/f1uf/33ztQNfGnpByX/7bwKfL9aoYqhHFUMDrdcwSL+qYtjr/fZKkCoGahiAp3NdV8lkUoZh7Mrxt24ruq6rUqkkz/NkGEbNlaugYh2qcrlczT+kaZoqFottHfPT//yjlu8rd22o+idUMTRAFUOd1msYJCm95/vttVarGAhVwJM5jvPU+qRsNvvU45im2XCNlOu6/p+3MsapU6eUTCbbXmIU21C1tfrfsiz/PufDhw914EBra6IAAEDnaGY9VTt3rLaugG2/5WdZlt588822Q1Vse6pM09TU1FTNwjHHcWrugQIAgO5SLBZ3tUZpK1c8vjB9+xWsoGJ7pUqSbNtWLpeTtHkPdHZ2lr4qAAC61FawMU1Truvu+EReO7f/to75+PG7/uk/0zRDWTgWhlaqGHqihkEKXsXQyn4dVt8QtIphr/frZp1Q4dBKDYNEFQPa95y5t0tngp7PdV0NDw9L2ryYslPQaef2nyS9+uqrchzHzxiO4+i1115r65hSzENVLDz7jKqPHlHFsF0bVQxdW+HQRhXDXu/Xvfa++qFV1DAgCtVHjwK3m7d73sS+1lYZWZalQqFQ96Ba2KanpzU3N6dcLqf79+/Ltu2211NJhKqn++UXSuzb11IVQ1fXMEhtVzG0sl/H1De0WcWw1/t1p2iqH1qaMEANg0QVA9rTarCJ+rztXoVq1vT0dOjHJFQ1qdUqhq4OVf8kaBVDy/t1QqjaErSKYa/362KdUOHQag2DRBUD0Ali+/QfAABAJyFUAQAAhIBQBQAAEAJCFQAAXaJarUY9Qtdp5WfKQvVd1Eq3ldRD/VZB7EUn1vb9ulwr/VZRdWl1yn6t9Fu124lFvxV2sn//fknSp59+qi9/+csRT9NdPv30U0m/+hk/CaFqN9BtFR46scIVsE+rUzq49nq/IP1WQTux6LfCkzzzzDMaGBjQJ59s1tU8//zzSiQSEU/V2arVqj799FN98sknGhgY0DPPPPPUfQhVuyFAt5XUA/1WQexhJ9b2/bpWkD6tiLq0OmG/fb/xzZb6rdrtxKLfCk/yO7/zO5LkByuEY2BgwP/ZPg2hahe12m0l9Ua/VRB71on1T/t1vVb7raLo0uqQ/VrvtwreiUW/FZ4kkUjo61//ur72ta/p888/j3qcrrB///6mrlBtIVQBANBFnnnmmZaCAMLD038AAAAhIFQBAACEgNt/MUQVQwzsdYVDnPfrhBkj3G8v6ymoYgDijVAVJ1QxRC+CCodO2K8TZoxqv72qp6CKAYg/QlWcUMUQvYgqHOK8XyfMGOV+e1VP0Wp9g0QVA7DXCFUxRBVD9Pa6wiH2+3XCjBHtt5f1FFQxAPHGQnUAAIAQEKoAAABCQKgCAAAIAWuqukjQKoY479cJM3b7fp0wY5T77WU9RdyrGFr9PNBtEtVqtRr1EHthZWVFkvTVpbdaWgT+lX/9r/Xb/+OU1ifPN71fkH3a2u+lP9LXzp9TYl/rFx6rjx7Ffr9OmLHb9+uEGSPZ75++oHq392lnv0fVR9oXYL8vHlX1zL7Enu3XrbZ+9xw6dCjiSbAXuFLVDdqsYojzfp0wY7fv1wkzdsp+nTBjO/vt/+2v6df/9N81/Xmg28Q+VC0tLck0TbmuK9M0Zdt21CPFVtAqhrjv1wkzdvt+nTBjp+zXCTMG3e+5P/gDQhV6WqxDVTablW3bsqzNFuHJyUmZpql0Oh3xZAAAALVi/fRfLpfzA5UkjY2NKZfLRTgRAABAY7ENVY7jyDCMhtsBAADiJrZP/xUKBV28eFF379594rZm/ehHP1K1WtUz//iP0i9/2fR+iV/7Ne3r69MXlUrT+wXZh/2iPxf7RX+uTtzv0caGql980dw++/dr3/PPt7RPzX6ffio9etT8fs8+q8SXvqTqL36haqv7Pfecqp991tL5tG+fEs891/zne8Bnn32mRCKhP/zDP4x6FOyB2K6pqlQqSiZru1UMw5DneYGOl0hsPuL7TLK1vpYtQfbby3N1+36dMGO379cJM0ax376+1ruZguwjSfuefz7QfokvfUlBSg4ISO1LJBL+7x90v9iGqmQyqUqlUrPN87yGtwSb8eKLL4YxFgAAQEOxXVPV6KpUo6tXAAAAcRDbULX9qb8tnuc13A4AABC12IYqSRofH6952s9xHGUymQgnAgAAaCy2T/9tmZub0+DgoDzPo1EdAADEVuxDFQAAQCeI9e0/AACATkGoAgAACAGhCgAAIASEKgAAgBAQqgAAAEJAqAIAAAgBoQoAACAEhCoAAIAQEKoAAABCQKgCAAAIAaEKAAAgBIQqAACAEBCqAAAAQkCoAgAACAGhCgAAIASEKgAAgBAQqgAAAEJAqAIAAAgBoQoAACAEhCoAAIAQPBvlyT3P061btyRJ5XJZ6XRamUym5jNLS0syTVOu68o0Tdm2HcWoAAAATxRpqHr99dc1PT0t0zQlSQcPHpRpmrIsS5KUzWZl27b/enJyUqZpKp1ORzYzAABAI5He/ltfX5fjOP5rwzBULpf917lczg9UkjQ2NqZcLrenMwIAADQj0itV+Xze/7PnefI8zw9RjuPIMIy6fbaHMAAAgLiINFRt8TxPr7/+umZnZ/1be57n1X0umUyqUqkEOsePf/xjVatV7d+/v61ZAQBo1ueff65EIqEXX3wx6lGwByIPVVuL1QcGBpRMJv3tlUql5rW0eXuwUdhqRrVaVbVabWtWAABawe+d3hJ5qDIMw3/ib2RkRJJk23bDq1Ke5zW8JdiMrStU//HvqvrJzzea3u9fDf6Wzhz+57ry4V81vd/vf7VPF49+I9CcAIDusbKyEvUI2EORh6rthoeH9eabb8q27YZXpRpdvWrVT36+ob/52cOmPz848JVA+wEAgN4S2dN/5XJZIyMjcl3X39bf3++/3v7U35btC9kBAADiJLJQZRiGTNOsufK0urqqkydP+q/Hx8drnvZzHKeuHBQAACAOIrv9Z5qmXnvtNeVyOb+fanx8XNPT0/5nFhcXNTc3J9d15XmeMpkMxZ8AACCWIl1T1cxXzmwPWQAAAHHFFyoDAACEgFAFAAAQAkIVAABACAhVAAAAISBUAQAAhIBQBQAAEAJCFQAAQAgIVQAAACEgVAEAAISAUAUAABACQhUAAEAICFUAAAAhIFQBAACEgFAFAAAQAkIVAABACAhVAAAAISBUAQAAhIBQBQAAEAJCFQAAQAgIVQAAACEgVAEAAISAUAUAABACQhUAAEAICFUAAAAheDbKk3uep+vXr0uSVldXlclkZNu2/77ruvI8T+l0Wp7nqVQqybKsqMYFAADYUaShan5+XrOzs5I2A9bIyIjy+bzS6bQkqVAoaH5+3v/81NQUoQoAAMRSZKHKdV3/SpRhGDIMQ5Zl6fr161pcXPQ/d+fOHf9qFQAAQFxFeqWqVCqpUqnIMAxJkmmaKpVKNZ8xTTOK0QAAAFoSWagyTVN3796t2eY4Tt3tvUKhIElaWVnRoUOHatZcAQAAxEWkV6q2K5fLqlQqmpqa8rel02kNDw/LMAzZtq0jR47INE1uBQIAgNiJTaXCxYsXlc/n/VuBkmRZVs3rVCqlXC4XxXgAAABPFItQNTc3pytXrtStn3o8QDVacwUAABAHkYeqXC6n48eP+7f0HMeRtPl0YDableu6/mcfPnyoAwcORDInAADAk0QaqrYClGEYcl1X5XJZ5XJZ0uZVqampqZqrV47j6OzZs5HMCgAA8CSRLVT3PE+nTp2q2759obpt2/4twHK5rNnZWRapAwCAWIosVBmGoXv37j3xM6ZpKpPJ7NFEAAAAwUW+pgoAAKAbEKoAAABCQKgCAAAIAaEKAAAgBIQqAACAEBCqAAAAQkCoAgAACAGhCgAAIASEKgAAgBAQqgAAAEJAqAIAAAgBoQoAACAEhCoAAIAQEKoAAABCQKgCAAAIAaEKAAAgBIQqAACAEBCqAAAAQvBs0B2LxaJmZmbkum7D99fW1gIPBQAA0GkCh6psNqtUKqXLly8rmUyGORMAAEDHCRyqKpWK3njjDfX19YU5DwAAQEcKvKbKtm2VSqUwZwEAAOhYga9UDQ4O6tSpU0qlUhoeHpZhGDXvX7hwoe3hAAAAOkXgULWysqLR0VFJqlusnkgk2psKAACgwwQOVQsLC22f3PM8Xb9+XZK0urqqTCYj27ZrPrO0tCTTNOW6rkzTrHsfAAAgDgKHqi3FYlGO40iSvvGNb+jo0aNN7zs/P6/Z2VlJmwFrZGRE+Xxe6XRa0uYThrZty7IsSdLk5KRM0/TfBwAAiIu2yj9Pnz6tc+fOqVwuq1wu67vf/a7+5E/+pKl9XdeV67ryPE+SZBiGLMvyr1xJUi6X8wOVJI2NjSmXy7UzMgAAwK4IfKVqZmZG/f39+vjjj2u2T05O6tKlS7p06dJTj1EqlVSpVPxF7qZp+k8UOo5Tt/h9azsAAEDcBL5SVSgU9MYbb9Rtv3r1qm7duvXU/U3T1N27d2Wapr/NcRwNDw9Lkn8Fa7tkMqlKpRJ0ZAAAgF0TOFQZhtEw4FQqFVWr1ZaPVy6XValUNDU15R/n8aZ2wzAahi0AAICoBQ5Vx44dUzab1cbGhr9tY2ND58+fVyaTafl4Fy9eVD6f92/5Nboq5Xlew1uCAAAAUQscqqanp/WVr3xFIyMjOnbsmI4dO6aRkRENDAy0XPw5NzenK1eu1NwKbHRVqtHVKwAAgDhoq1JhcXHRf/LP8zxZlqVUKtXSMXK5nI4fP+7XJDiOI8uyap7627J1DgAAgLhpu6cqnU4H7o3aepLPMAy/XqFcLvvBaXx83A9ZW5/fWnMFAAAQJ02HqpmZGdm27X81zbVr1574+afdAvQ8T6dOnarbvj00LS4uam5uzg9cmUyG4k8AABBLTYeqlZUVjY2N+a/L5fKOn23mu/8Mw9C9e/ee+rnp6enmBgQAAIhQ06Eqn8/XvL5x40bowwAAAHSqwE//ba9S2G59fV3r6+uBBwIAAOhEgUPVyMhIw+2u62pmZibwQAAAAJ0ocKjaqTV9eHjY//4+AACAXtFypcLRo0eVSCSUSCR07Nixuvdd1225qwoAAKDTtRyqZmdnVa1Wdfr06Ya1CaZpEqoAAEDPaTlUbfVUjY+Pa3x8PPSBAAAAOlHgNVULCwthzgEAANDR2vqamrW1Nf+rZrYzDEMnTpxo59AAAAAdJXCoun37ts6dO6dUKqW1tTUNDQ1JklZXVzU2NkaoAgAAPSVwqLp27Zpu3ryp0dFRTUxM6O2331ZfX58KhYKKxWKYMwIAAMRe4DVV9+/f9xetm6bpd1PZtq1CoRDOdAAAAB0icKjauu0nSZZlaWlpSZL07rvvyvO8cKYDAADoEIFv/509e1alUklDQ0PKZDLK5XIaGhpStVrVq6++GuaMAAAAsRc4VD3eUZXP5+W6rgYGBtTf39/2YAAAAJ0k8O2/YrGoDz74oGabaZp68803/duCAAAAvSJwqJqfn294RWp4eFjz8/NtDQUAANBpAoeqcrmsQ4cO1W23LKthISgAAEA3CxyqTNPU+vp63fYHDx7owIEDbQ0FAADQaQKHqpMnT2pyclI//elP/W3r6+s6f/68bNsOZTgAAIBOEfjpvzNnzsh1Xb300ktKJpOSJM/zdPLkSV24cCG0AQEAADpBW1+oPDs7qzNnzmh1dVXSZiGoaZqhDAYAANBJ2gpV0uYtv62vqNm3bx+hCgAA9KS2QtXp06dVKpU0PDwsScrlcrp+/bp+8IMfhDIcAABApwi8UH1mZkb9/f36+OOPdePGDd24cUMff/yxfvd3f1eXLl1q6hiu62pubk65XK7he+VyWdLmWi1qGgAAQJwFDlWFQkFvvPFG3farV6/q1q1bT93fcRyVy2Wtrq42/ALmQqGgiYkJHTx4UCMjI37AAgAAiKPAt/8Mw1ClUlFfX1/N9kqlomq1+tT9LcuSJL3//vs7fubOnTvyPE/pdDromAAAAHsi8JWqY8eOKZvNamNjw9+2sbGh8+fPK5PJhDKcaZoEKgAA0BECX6manp7W5OSkRkZG/Cf+XNeVZVmh9VQVCgVJ0srKig4dOkSpKAAAiK22nv5bXFxUuVxWuVyW53myLEupVCqUwdLptIaHh2UYhmzb1pEjR7hyBQAAYqvtnqp0Or0rQWdrzdWWVCqlXC6n2dnZ0M8FAADQrqZD1bVr11o6cLu3AHO5XM3aLNM0VSwW2zomAADAbmk6VLVSaZBIJAINs8V1XWWzWVmW5a/XevjwoQ4cONDWcQEAAHZL06Hqxo0buzlHDdM0NTU1VfOVN47jaGFhYc9mAAAAaEXba6qCchxHy8vLKhaLWl1d1YMHD/Tyyy/7Qcq2bb9pvVwua3Z2lkXqAAAgtiJbU2VZlizL0vT0dMP3TdMMre8KAABgt8VyTRUAAECnieWaKgAAgE4T+GtqAAAA8CttLVRfW1uT4zh12w3D0IkTJ9o5NAAAQEcJHKpu376tc+fOKZVKaW1tTUNDQ5Kk1dVVjY2NEaoAAEBPCRyqrl27pps3b2p0dFQTExN6++231dfXp0KhQPM5AADoOYHXVN2/f1+jo6OSNusPSqWSpM1+qUKhEM50AAAAHSJwqNq67Sdtdk4tLS1Jkt599115nhfOdAAAAB0i8O2/s2fPqlQqaWhoSJlMRrlcTkNDQ6pWqzpz5kyYMwIAAMRe4FA1Pj5e8zqfz8t1XQ0MDKi/v7/twQAAADpJqN/9t/0LkAEAAHoJ5Z8AAAAhIFQBAACEgFAFAAAQAkIVAABACHYlVNGoDgAAek3gUDUzM6O//uu/rtm2sbGhbDar06dPtz0YAABAJ2mrUf1b3/qWfvjDH0ravDr1R3/0R/rpT3+qDz74ILQBAQAAOkHgnqpMJiNJeuWVVzQ2NibHcXT58mWdPHkytOEAAAA6RVvln5lMRqZp6vTp05qeniZQAQCAntV0qHrrrbd2fM+yLM3Pz+v+/fsaHByUJH3nO99pfzoAAIAO0XSoWl5efuL7o6Ojcl1XrusqkUgQqgAAQE9pOlTduHFjN+cAAADoaJR/AgAAhCDwQvVisaiZmRm5rluzvVqtKpFIaG1t7anHcF1X77zzjgYHB/2nCbdbWlqSaZpyXVemacq27aDjAgAA7KrAoSqbzSqVSuny5ctKJpMt7+84jjzP0+rqqgYGBhoe37ZtWZYlSZqcnJRpmkqn00FHBgAA2DWBQ1WlUtEbb7yhvr6+QPtvhaX333+/4fu5XE6zs7P+67GxsbptAAAAcRF4TZVt2yqVSmHO4nMcR4ZhNNwOAAAQR4GvVA0ODurUqVNKpVIaHh6uC0EXLlwIPJTneXXbksmkKpVK4GMCAADspsChamVlRaOjo5JUt1g9kUi0NVSlUqlbp2UYRsOwBQAAEAeBQ9XCwkKYc9RodFXK87yGtwQBAADiIJY9VY2uSjW6egUAABAXbX2h8traWsPF44Zh6MSJE4GPu/Vk4Hae5zXcDgAAEAeBQ9Xt27d17tw5pVIpra2taWhoSJK0urqqsbGxtkKVJI2Pj8txHD9IOY6jqampto4JAACwWwKHqmvXrunmzZsaHR3VxMSE3n77bfX19alQKKhYLD51f8dxtLy8rGKxqNXVVT148EAvv/yyTNOUJC0uLmpubk6u68rzPGUyGYo/AQBAbAUOVffv3/ef/jNNU6VSSYcPH5Zt25qZmdHly5efuL9lWbIsS9PT0zt+5knvAQAAxEnghepbt/2kzYC0tLQkSXr33XepPgAAAD0n8JWqs2fPqlQqaWhoSJlMRrlcTkNDQ6pWqzpz5kyYMwIAAMRe4FA1Pj5e8zqfz8t1XQ0MDKi/v7/twQAAADpJ4Nt/jRajm6ZJoAIAAD0pcKjKZrP68MMPw5wFAACgYwUOVa+++qrm5ua0sbER5jwAAAAdKfCaqkQiob6+Pr300ksaHR31+6W2XLhwoe3hAAAAOkXgULW8vKyBgQENDAzI8zyVy2X/vUQiEcpwAAAAnSJwqFpYWAhzDgAAgI4WeE0VAAAAfoVQBQAAEAJCFQAAQAgIVQAAACEgVAEAAISAUAUAABACQhUAAEAICFUAAAAhIFQBAACEgFAFAAAQAkIVAABACAhVAAAAISBUAQAAhIBQBQAAEAJCFQAAQAhiHapc11W5XJYkeZ4nx3EinggAAKCxWIeqQqGgiYkJHTx4UCMjI37AAgAAiJtnox7gae7cuSPP85ROp6MeBQAAYEexD1WmaUY9AgAAwFPFPlQVCgVJ0srKig4dOiTbtiOeCAAAoF6sQ1U6ndbw8LAMw5Bt2zpy5IhM0+RWIAAAiJ1YL1S3LEuGYfivU6mUcrlchBMBAAA0FutQ9XiAMk1TpVIpomkAAAB2FttQ5bqustmsXNf1tz18+FAHDhyIcCoAAIDGYhuqTNPU1NRUzdN/juPo7NmzEU4FAADQWKwXqtu27d8CLJfLmp2dZZE6AACIpViHKtM0lclkoh4DAADgqWJ7+w8AAKCTEKoAAABCQKgCAAAIAaEKAAAgBIQqAACAEBCqAAAAQkCoAgAACAGhCgAAIASEKgAAgBAQqgAAAEJAqAIAAAgBoQoAACAEhCoAAIAQEKoAAABCQKgCAAAIAaEKAAAgBIQqAACAEBCqAAAAQkCoAgAACAGhCgAAIATPRj1Ap/j9r/btymcBAEB3IFQ9xcZnn+uLR1VdPPqNlvb74lFVz+xL7NJUAAAgbghVT/EPn36mZ/YlVPlf/lf9Yu2vm9pn/29/Tb/+p/9ulycDAABxQqhq0i/W/lr/5T/9p6Y++9wf/AGhCgCAHhP7ULW0tCTTNOW6rkzTlG3bUY8EAABQJ9ahKpvNyrZtWZYlSZqcnJRpmkqn0xFPBgAAUCvWlQq5XM4PVJI0NjamXC4X4UQAAACNxfZKleM4Mgyj4fYgPv/8c1WrVf3pP/s1/XLwuab3e3ZfQisrK6oe+SPpv//Xze20b59+vrISaE4AQPf47LPPlEjwJHiviG2o8jyvblsymVSlUgl0vK1/qQe+vD/Y/vv3S/uD7QsA6E2JRIJQ1UNiG6oqlYqSyWTNNsMwGoatZrz44othjAUAANBQbNdUNboq5Xlew1uCAAAAUYttqGp0VarR1SsAAIA4iG2o2v7U3xbP8xpuBwAAiFpsQ5UkjY+P1zzt5ziOMplMhBMBAAA0lqhWq9Woh3iSubk5DQ4OyvM8GtUBAEBsxT5UAQAAdIJY3/4DAADoFIQqAACAEBCqAAAAQkCoAgAACAGhCgAAIASEKgAAgBAQqgAAAELwbNQD7LW1tTX9xV/8hX74wx/6X9icTCY1Ojqqb3/723rhhRcinnDvFYtFjY6O+q/X1tZ0/fp1ua6r4eFhvfbaa/q93/u9CCcE4ou/U+rxdwp6VU+Vf96+fVsXL17U2bNnlUql/O2e52llZUXvvfeerl69qqNHj0Y45d4bGhrS2tqapM2/DE+dOqXx8XElk0mtr6+rWCzqP/yH/9Bzvxz4ZVmLX5T1+DulMf5OQc+q9pCjR49W79+/v+P79+/fr05MTOzhRPFw8OBB/88TExPVcrlc8/6tW7d67udSKBSqIyMj1e9///tVx3H8/xUKherc3Fz1m9/8ZvWDDz6Iesw99cILL/h/dhynevDgwerk5GT14sWL1VOnTlVfeOGF6traWoQT7j3+TmmMv1PQq3rq9t+DBw9kmuaO7w8MDGh9fX0PJ4qHRCLh/9nzPCWTyZr3bdtWNpvd67Eide3aNf3lX/5lw39fxsfHlclkdP78+Z66AlHddlF7fn5e+Xy+5upMoVDQd7/7Xf3lX/5lFONFgr9TGuPvFPSqnlqobtu2ZmZmtLGxUffe+vq6Ll68qMOHD0cwWbSq1ap++MMfqlgsKpVKaXJysuZntLa2pgMHDkQ44d7jl2W9Zn5Ruq6712NFir9TGuPvFPSqnrpSNTs7q8nJSf3Lf/kvNTg46G/f+kUwPj6uq1evRjVeZFKplL73ve/VbHvw4IH6+vr08OFDTU5Oanp6OqLporH1y3J6elp9fX01762vr2t+fr7nfllu/aKsVqv+L8o///M/938+vfiLkr9TGuPvFPSqnlqovl2xWJTneZIkwzA0PDys/v7+iKeKp9XV1ZrbPL1icnJSH3zwwY6/LK9cudJT/85MTEzUbVtcXNSBAwf08OFDTUxMaHp6WseOHYtguujxd0rzevXvFHS/ng1VQLP4ZdkcflEC6HWEKgBoQ7FYVC6Xk+u6evnll3XixAn/vYcPH+qVV17pqcX7W/i5oBf11EJ1IEzvvfde1CPETq/9TG7fvq1Tp07JMAz98R//sf7iL/5C58+fr/nM6upqNMNFqNHP5d//+39f85le/Lmg+3GlCgjo2LFj+uCDD6IeI1Z67Wdy7NgxXbhwQePj4/62bDarZDKpCxcu6OHDh/rmN7/pF2H2Cn4u6FU99fQf0Kxr16499TO9Vh/Az6Te/fv3a4KDtPlE4Pz8vN577z3Zth3RZNHi54JeRagCGkgmk8rlciy83oafSb1UKqUf/vCHdfUaU1NTOnfuXERTRY+fC3oVoQpoYHx8XI7jaGFhYcfP9FKbusTPpJGFhQWdPn1almXp8uXLde9NTk5GNFm0+LmgVxGqgAZM03xqkWUmk9mjaeKBn0k90zT14Ycf7rjoenFxsScXZPNzQa9ioToAAEAIqFQAAAAIAaEKAAAgBIQqAACAEBCqAAAAQkCoAjrQxMSEstls1GPUiONMALCXqFQAOtBrr70mwzCiHqNGHGcCgL1EpQLQBRzHUTab1Z07d7rqXADQSbj9BwAAEAJCFdCBTp06pbm5OUnS5OSkTp06Jdd1dfDgQR08eFCe5/mfzWazGhkZ0ZEjR5TL5fztk5OTyuVyyuVyOnLkiBzHkSQVCgVNTEzo4MGDOnLkiAqFQs0+O51r+0xb5ubmdOTIEY2MjNStt5qcnNTS0lLNfFszAEAnIlQBHW5xcVELCwsyTVP37t3TvXv3/LVNk5OTcl1XH330kRYWFjQ/P69yuSxJevjwoXK5nJaWljQ1NSXLsiRJlUpFV65c0b179zQ7O6tz5875+zzpXI+bnJzU6uqqbt68qY8++kgPHjzQqVOn/PcfPnyo+fl52batjz76SKlUioXuADoaC9WBLuW6rm7fvq27d+/KMAyl02lNTU3p/fffVzqd9j/z0Ucf1QSj7d/fZ1mW0um0HMfx92lGuVyuObe0GchGRkbkOI4f4NLptP/nl19+uSZ0AUCnIVQBXWrr6tJLL71Us314eNj/8+joaMMrTblcTsvLy1pfX1e5XNYf//Eft3TuUqkk0zTrjj08PKzl5WU/SG2fJZlMtnQOAIgbQhXQxdLptPL5/I7vm6ZZt21iYkLJZFJnzpzR8PCwXnnllZbPu31N15P09/e3fGwAiCvWVAFdKp1Oq1wuNx1wpM3bgeVyWTdv3pRlWYF7pyzLkuu6deculUo6dOhQoGMCQNwRqoAuYJqmH2Icx5HrujJNU5lMRufOnZPrupI2n+xbWlra8Thbt+C2nhIsFAr+bcQnnetxW2ulXnnlFf+zk5OTMk1Ttm2H9Y8NALFCqAK6QDqdVjqd1ksvvVQTmmZnZ5VKpTQxMaGRkRHlcjl/PVMjhmHo1Vdf9WsOthaVb79itdO5Hnfz5k2Njo5qYmJCL730kgYGBp54KxIAOh2N6gAAACHgShUAAEAICFUAAAAhIFQBAACEgFAFAAAQAkIVAABACAhVAAAAISBUAQAAhIBQBQAAEAJCFQAAQAgIVQAAACEgVAEAAITg/wc8W0tVXWYnVwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 500x800 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mf.plot_ra_evolution_col([rh_bot, rh_uni, rh_top],  dpi=100, figsize=(5, 8)) #filename=\"plots/ra_%s.pdf\"%filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc781a6d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.16"
  },
  "vscode": {
   "interpreter": {
    "hash": "cbb24b48ae2642bb0fe3c3a73dd180d2cf6b4e8df7de3f3850c06f8dd008ecd7"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}