quantization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6dc9f1eb",
   "metadata": {},
   "source": [
    "# The task\n",
    "\n",
    "From taks 3, we are familiar with the model to predict default probability and potential losses. The personal loans and risk team is using the model as a guide to loss provisions for the upcoming year. The team now wants to look at their mortgage book. They suspect that FICO scores will provide a good indication of how likely a customer is to default on their mortgage. The team wants to build a machine learning model that will predict the probability of default. The architecture they are using requires categorical data. As FICO ratings can take integer values in a large range, they will need to be mapped into buckets. They ask if we can find the best way of doing this to analyze the data. The process of doing this is known as quantization. We could consider many ways of solving the problem by optimizing different properties of the resulting buckets, such as the mean squared error (MSE) or log-likelihood. \n",
    "\n",
    "The data is provided here (https://cdn.theforage.com/vinternships/companyassets/Sj7temL583QAYpHXD/JiwEkbBq8pFwMRYLc/1684245611456/Task%203%20and%204_Loan_Data.csv).\n",
    "\n",
    "# The Lloyd-Max algorithm for quantization\n",
    "\n",
    "The task is essentially a one-dimensional quantization problem. The goal is to find the boundaries of different buckets and the representative value of each bucket. We can use sklearn.cluster.KMeans to solve it. However, KMeans cannot guarantee the centroids, i.e., the representative value of each bucket, are integers; and KMeans doesn't give boundries, but only tells which bucket (cluster) each data point is in. The extra steps are: (1) to approximate those centroids as integers by rounding; and (2) using each cluster's maximum and minimum values to find the boundaries. The other option is to implement the Lloyd-Max algorithm using dynamic programming to solve it. The details of the algorithm can be found here: https://en.wikipedia.org/wiki/Lloyd%27s_algorithm and https://ocw.mit.edu/courses/6-450-principles-of-digital-communications-i-fall-2006/926689aaa62a0315473fa9b982de1b07_book_3.pdf. I will first implement the Lloyd-Max algorithm without limiting the representative values as integers and use KMeans to validate the implementation, and then limit the representative values to integers only to see whether it makes big differences to the MSE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8576393d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import bisect\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "968bf546",
   "metadata": {},
   "outputs": [],
   "source": [
    "borrowers = pd.read_csv('./Task 3 and 4_Loan_Data.csv')\n",
    "borrowers.set_index('customer_id', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a2e2fbfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# implement Lloyd-Max algorithm\n",
    "def dp_lloyd_quantization(data, n, tol=0.00001, max_iter=500, is_integer=False):\n",
    "    data = sorted(data)\n",
    "    l = len(data)\n",
    "    mse = float('inf')\n",
    "    centroids = [data[i] for i in range(0, l//n*n, l//n)]\n",
    "    while max_iter > 0:\n",
    "        interval = [-float('inf')] + [0.5*(centroids[i]+centroids[i+1]) for i in range(len(centroids)-1)] \\\n",
    "                    + [float('inf')]\n",
    "        se = 0\n",
    "        counter = []\n",
    "        for j in range(len(interval) - 1):\n",
    "            nums = data[bisect.bisect_right(data, interval[j]):bisect.bisect_right(data, interval[j+1])]\n",
    "            centroids[j] = np.mean(nums)\n",
    "            if is_integer:\n",
    "                centroids[j] = int(round(centroids[j], 0))\n",
    "            se += np.var(nums) * len(nums)\n",
    "            counter.append(len(nums))\n",
    "        if mse - se/l > tol:\n",
    "            mse = se/l\n",
    "            max_iter -= 1\n",
    "        else:\n",
    "            break\n",
    "    return centroids, mse, interval, counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "025fb3fa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The representative values of each bucket are: [498.64057971014495, 550.6336206896551, 586.3137931034482, 616.5715915573303, 645.2651708312086, 676.5701086956522, 712.6266129032258, 761.1635610766045]\n",
      "\n",
      "The MSE is: 125.66896199112011\n",
      "\n",
      "The boundaries of buckets are: [-inf, 524.6371001999, 568.4737068965517, 601.4426923303893, 630.9183811942694, 660.9176397634303, 694.598360799439, 736.8950869899152, inf]\n",
      "\n",
      "The number of data points in each bucket are: [345, 928, 1450, 1753, 1961, 1840, 1240, 483]\n"
     ]
    }
   ],
   "source": [
    "# when the number of buckets is 8 and without limiting the representative values to be integers\n",
    "centroids, mse, interval, counter = dp_lloyd_quantization(data=borrowers.fico_score.values, n=8, \n",
    "                                                          tol=0.00001, max_iter=500)\n",
    "print(f\"The representative values of each bucket are: {centroids}\\n\")\n",
    "print(f\"The MSE is: {mse}\\n\")\n",
    "print(f\"The boundaries of buckets are: {interval}\\n\")\n",
    "print(f\"The number of data points in each bucket are: {counter}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6729c9a8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The representative values of each bucket are: [499, 552, 587, 618, 647, 678, 714, 763]\n",
      "\n",
      "The MSE is: 124.99724103265251\n",
      "\n",
      "The boundaries of buckets are: [-inf, 525.5, 569.5, 602.5, 632.5, 662.5, 696.0, 738.5, inf]\n",
      "\n",
      "The number of data points in each bucket are: [356, 948, 1469, 1822, 1979, 1769, 1209, 448]\n"
     ]
    }
   ],
   "source": [
    "# when the number of buckets is 8 and the the representative values are limited to integers\n",
    "centroids, mse, interval, counter = dp_lloyd_quantization(data=borrowers.fico_score.values, n=8, \n",
    "                                                          tol=0.00001, max_iter=500, is_integer=True)\n",
    "print(f\"The representative values of each bucket are: {centroids}\\n\")\n",
    "print(f\"The MSE is: {mse}\\n\")\n",
    "print(f\"The boundaries of buckets are: {interval}\\n\")\n",
    "print(f\"The number of data points in each bucket are: {counter}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "18796bae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1IAAAINCAYAAAA0iU6RAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABpGUlEQVR4nO3dd3wUdeL/8fdm0zvpCV16lI5AEBAFCchxICgWDkE8vQtBaaLyOxGQO9opKipgBf2iINwJCh4gohQxEgQj1YgYCCUhkZCEYtru/P6IrK60LEJmA6/n47EPszOzO+8JMcubz8xnLIZhGAIAAAAAVJiH2QEAAAAAoKqhSAEAAACAiyhSAAAAAOAiihQAAAAAuIgiBQAAAAAuokgBAAAAgIsoUgAAAADgIooUAAAAALjI0+wA7sBut+vIkSMKCgqSxWIxOw4AAAAAkxiGoRMnTiguLk4eHucfd6JISTpy5Ihq1qxpdgwAAAAAbuLgwYOqUaPGeddTpCQFBQVJKv9mBQcHm5wGAAAAgFkKCwtVs2ZNR0c4H4qU5DidLzg4mCIFAAAA4KKX/DDZBAAAAAC4iCIFAAAAAC6iSAEAAACAi7hGCgAAAFcFwzBUVlYmm81mdhS4MavVKk9Pzz982yOKFAAAAKq8kpISZWVl6fTp02ZHQRXg7++v2NhYeXt7X/J7UKQAAABQpdntdmVkZMhqtSouLk7e3t5/eLQBVyfDMFRSUqLc3FxlZGSoQYMGF7zp7oVQpAAAAFCllZSUyG63q2bNmvL39zc7Dtycn5+fvLy8dODAAZWUlMjX1/eS3ofJJgAAAHBVuNSRBVx7LsfPCj9tAAAAAOAiihQAAAAAuIhrpAAAAABJNruh1Iw85ZwoUlSQr9rWDZPV49qZtGLIkCHKz8/XsmXLzI5SJVCkAAAAcM1btTNLk5bvVlZBkWNZbIivJvSOV48bYk1MBnfFqX0AAAC4pq3amaWkBducSpQkZRcUKWnBNq3amWVSMrgzipQbsdkNpew7pg/TDitl3zHZ7IbZkQAAAK5qNruhSct361x/6zqzbNLy3Vfs72VdunTRI488opEjR6patWqKjo7W66+/rlOnTumBBx5QUFCQ6tevr5UrV0qSjh8/roEDByoyMlJ+fn5q0KCB5s2b53i/gwcPasCAAQoNDVVYWJj69Omj/fv3X1K2LVu2KDIyUtOnT5ckTZw4US1atNBbb72lWrVqKTAwUMOGDZPNZtOMGTMUExOjqKgo/etf/3J6n/z8fP31r39VZGSkgoODdeutt+rbb791rN+3b5/69Omj6OhoBQYG6sYbb9Snn37q9B516tTRlClTNHToUAUFBalWrVp67bXXHOtLSko0fPhwxcbGytfXV7Vr19bUqVMv6bgriiLlJlbtzFLH6Z/p+fmztP3z7np+/ix1nP4Z/wICAABwBaVm5J01EvVbhqSsgiKlZuRdsQxvv/22IiIilJqaqkceeURJSUm666671KFDB23btk3du3fXoEGDdPr0aY0fP167d+/WypUrtWfPHs2ZM0cRERGSpNLSUiUmJiooKEgbN27Upk2bFBgYqB49eqikpMSlTJ999pluu+02/etf/9ITTzzhWL5v3z6tXLlSq1at0sKFC/Xmm2+qV69eOnTokNavX6/p06frqaee0ubNmx2vueuuu5STk6OVK1dq69atatWqlbp27aq8vPLv6cmTJ3X77bdr7dq1+uabb9SjRw/17t1bmZmZTpmee+45tWnTRt98842GDRumpKQkpaenS5JmzZqljz76SIsXL1Z6erreffdd1alT51L+OCqMIuUGfh1O/ln3+XyoxWEeus/nQ2UX/MxwMgAAwBWUc+L8JepStrsUzZs311NPPaUGDRpo3Lhx8vX1VUREhB566CE1aNBATz/9tI4dO6bt27crMzNTLVu2VJs2bVSnTh1169ZNvXv3liS9//77stvteuONN9S0aVM1adJE8+bNU2ZmptatW1fhPEuXLlWfPn306quv6uGHH3ZaZ7fb9dZbbyk+Pl69e/fWLbfcovT0dL3wwgtq1KiRHnjgATVq1Eiff/65JOmLL75QamqqlixZojZt2qhBgwZ69tlnFRoaqv/85z+O4//b3/6mG264QQ0aNNDkyZNVr149ffTRR077vv322zVs2DDVr19fTzzxhCIiIhz7yczMVIMGDdSxY0fVrl1bHTt21L333nupfyQVQpEy2W+Hkzt7bFdDj0OSpIYeh9TJY7ukKzucDAAAcC2LCvK9rNtdimbNmjm+tlqtCg8PV9OmTR3LoqOjJUk5OTlKSkrSokWL1KJFCz3++OP68ssvHdt9++23+uGHHxQUFKTAwEAFBgYqLCxMRUVF2rdvX4WybN68WXfddZf+7//+T3ffffdZ6+vUqaOgoCCnbPHx8U43uI2OjlZOTo4j08mTJxUeHu7IFBgYqIyMDEemkydP6rHHHlOTJk0UGhqqwMBA7dmz56wRqd9+nywWi2JiYhz7GTJkiNLS0tSoUSM9+uij+uSTTyp0vH+EqUWqTp06slgsZz2Sk5MlSUVFRUpOTnZ84/v376+jR486vUdmZqZ69eolf39/RUVFaezYsSorKzPjcC5Jakaesk/lyMP3kPr5L9ZObx9J0k5vH/XzXyyL7yFln8q5osPJAAAA16q2dcMUG+Kr801yblH57H1t64ZdsQxeXl7O+7RYnJZZLOXp7Ha7evbsqQMHDmjUqFE6cuSIunbtqscee0xSeSFp3bq10tLSnB7ff/+97rvvvgplqVevnho3bqy33npLpaWlLmc9s8xutzsyxcbGnpUpPT1dY8eOlSQ99thjWrp0qaZMmaKNGzcqLS1NTZs2Pet0xAvtp1WrVsrIyNDkyZP1888/a8CAAbrzzjsrdMyXytTpz7ds2SKbzeZ4vnPnTt1222266667JEmjRo3Sxx9/rCVLligkJETDhw9Xv379tGnTJkmSzWZTr169FBMToy+//FJZWVm6//775eXlpSlTpphyTK7KOVEkr9DN8olcq/GSpPL/SSdHhkmyKUAvqzi3q3JOdDYxJQAAwNXJ6mHRhN7xSlqwTRbJadKJM+VqQu94t7qfVGRkpAYPHqzBgwerU6dOGjt2rJ599lm1atVK77//vqKiohQcHHxJ7x0REaEPPvhAXbp00YABA7R48eKzCowrWrVqpezsbHl6ep73mqVNmzZpyJAhuuOOOySVl69LmSAjODhYd999t+6++27deeed6tGjh/Ly8hQWdmVKsKkjUpGRkYqJiXE8VqxYoXr16unmm29WQUGB3nzzTc2cOVO33nqrWrdurXnz5unLL7/UV199JUn65JNPtHv3bi1YsEAtWrRQz549NXnyZL3yyisuX1BnlqggX5Xmt9Xkg1a9d+ioJuYekyRNzD2m9w4d1eSDVpXmt72iw8kAAADXsh43xGrOX1opJsT571sxIb6a85dWbnUfqaeffloffvihfvjhB+3atUsrVqxQkyZNJEkDBw5URESE+vTpo40bNyojI0Pr1q3To48+qkOHDlV4H1FRUfrss8/03Xff6d577/1DZ3t169ZNCQkJ6tu3rz755BPt379fX375pf7xj3/o66+/liQ1aNBAH3zwgdLS0vTtt9/qvvvuc4w0VdTMmTO1cOFCfffdd/r++++1ZMkSxcTEKDQ09JKzX4zb3JC3pKRECxYs0OjRo2WxWLR161aVlpaqW7dujm0aN26sWrVqKSUlRe3bt1dKSoqaNm3qOG9UkhITE5WUlKRdu3apZcuW59xXcXGxiouLHc8LCwuv3IFdRNu6Yerjd0R9SzMkSVZL+Q9Nk5ISxZeWqqkytM7vyBUdTgYAALjW9bghVrfFxyg1I085J4oUFVR+Op87jURJkre3t8aNG6f9+/fLz89PnTp10qJFiyRJ/v7+2rBhg5544gn169dPJ06cUPXq1dW1a1eXR6hiYmL02WefqUuXLho4cKDee++9S8prsVj0v//9T//4xz/0wAMPKDc3VzExMercubPj7/AzZ87U0KFD1aFDB0VEROiJJ55w+e/nQUFBmjFjhvbu3Sur1aobb7xR//vf/5yu3brcLIZhuMUsBosXL9Z9992nzMxMxcXF6b333tMDDzzgVHgkqW3btrrllls0ffp0Pfzwwzpw4IBWr17tWH/69GkFBATof//7n3r27HnOfU2cOFGTJk06a3lBQcElD4NeMsNQ/qyOCs7bJQ+Lod3eXrq7eqzeP5yl+JJS2Q2LCsOuV+ijX0gW9/ofGQAAwB0UFRUpIyNDdevWla8vZ/Hg4i70M1NYWKiQkJCLdgO3mbXvzTffVM+ePRUXF3fF9zVu3DgVFBQ4HgcPHrzi+zwvW4lCS3LkYSnvs5E2m5KOFyjyl2vHPCyGQktzJFvVOFURAAAAuBa4xal9Bw4c0KeffqoPPvjAsSwmJkYlJSXKz893Orfx6NGjiomJcWyTmprq9F5nZvU7s825+Pj4yMfH5zIewR/g6SM9/Ll06ifZDENZhwsVf7pEWf7eCqseLKvFIgVElm8HAAAAXKLAwMDzrlu5cqU6depUiWmqPrcoUvPmzVNUVJR69erlWNa6dWt5eXlp7dq16t+/vyQpPT1dmZmZSkhIkCQlJCToX//6l3JychQVFSVJWrNmjYKDgxUfH1/5B3KpQmpIITVkldSsutlhAAAAcDVKS0s777rq1flLqKtML1J2u13z5s3T4MGD5en5a5yQkBA9+OCDGj16tMLCwhQcHKxHHnlECQkJat++vSSpe/fuio+P16BBgzRjxgxlZ2frqaeeUnJysvuMOAEAAABuoH79+mZHuKqYXqQ+/fRTZWZmaujQoWete/755+Xh4aH+/furuLhYiYmJmj17tmO91WrVihUrlJSUpISEBAUEBGjw4MF65plnKvMQAAAAAFxj3GbWPjNVdGYOAAAAuB9m7YOrrqpZ+wAAAACgqqBIAQAAAICLKFIAAAAA4CKKFAAAAIBzGjJkiPr27Wt2DLdEkQIAAADO2Pe59HLb8v8CF0CRAgAAACTJMKS1k6Sf0sv/y+TWuACKFAAAACBJ+9ZKR74p//rIN+XPr7AuXbrokUce0ciRI1WtWjVFR0fr9ddf16lTp/TAAw8oKChI9evX18qVKyVJx48f18CBAxUZGSk/Pz81aNBA8+bNc7zfwYMHNWDAAIWGhiosLEx9+vTR/v37K5TFZrNp9OjRCg0NVXh4uB5//HH9/k5JXbp00fDhwzV8+HCFhIQoIiJC48ePP2u7awFFCgAAADAM6bN/ShZr+XOLtfx5JRSEt99+WxEREUpNTdUjjzyipKQk3XXXXerQoYO2bdum7t27a9CgQTp9+rTGjx+v3bt3a+XKldqzZ4/mzJmjiIgISVJpaakSExMVFBSkjRs3atOmTQoMDFSPHj1UUlJy0RzPPfec5s+fr7feektffPGF8vLytHTp0nPm9fT0VGpqql588UXNnDlTb7zxxmX/vrg7bsgrbsgLAABQlV2WG/L+8Km0oP/Zy//yX6l+tz8W8AK6dOkim82mjRs3SiofFQoJCVG/fv30zjvvSJKys7MVGxurlJQUTZkyRREREXrrrbfOeq8FCxbon//8p/bs2SOLxSJJKikpUWhoqJYtW6bu3btfMEtcXJxGjRqlsWPHSpLKyspUt25dtW7dWsuWLXPkzcnJ0a5duxz7ePLJJ/XRRx9p9+7dl+V7Uhm4IS8AAADwR/1+NOqMShqVatasmeNrq9Wq8PBwNW3a1LEsOjpakpSTk6OkpCQtWrRILVq00OOPP64vv/zSsd23336rH374QUFBQQoMDFRgYKDCwsJUVFSkffv2XTBDQUGBsrKy1K5dO8cyT09PtWnT5qxt27dv7yhRkpSQkKC9e/fKZrO5fvBVmKfZAQAAAABT/fbaqN8ybL9eK3UFR6W8vLycnlssFqdlZ0qL3W5Xz549deDAAf3vf//TmjVr1LVrVyUnJ+vZZ5/VyZMn1bp1a7377rtn7SMyMvKK5b9WMSIFAACAa9eZ0ajz/rXYo9KulaqoyMhIDR48WAsWLNALL7yg1157TZLUqlUr7d27V1FRUapfv77TIyQk5ILvGRISotjYWG3evNmxrKysTFu3bj1r299uI0lfffWVGjRoIKvVeta2VzOKFAAAAK5dthKp4LAk+3k2sEuFh8u3cwNPP/20PvzwQ/3www/atWuXVqxYoSZNmkiSBg4cqIiICPXp00cbN25URkaG1q1bp0cffVSHDh266HuPGDFC06ZN07Jly/Tdd99p2LBhys/PP2u7zMxMjR49Wunp6Vq4cKFeeukljRgx4nIfqtvj1D4AAABcuzx9pIc/l079dP5tAiLLt3MD3t7eGjdunPbv3y8/Pz916tRJixYtkiT5+/trw4YNeuKJJ9SvXz+dOHFC1atXV9euXSs0odqYMWOUlZWlwYMHy8PDQ0OHDtUdd9yhgoICp+3uv/9+/fzzz2rbtq2sVqtGjBihhx9++Iocrztj1j4xax8AAEBVdllm7UOFdOnSRS1atNALL7xgdpQ/hFn7AAAAAMAEFCkAAADgGnBmSvRzPc7cxwoVxzVSAAAAwDUgLS3tvOuqV69eofdYt27d5QlzFaBIAQAAANeA+vXrmx3hqsKpfQAAAADgIooUAAAAALiIIgUAAAAALqJIAQAAAICLKFIAAAAA4CKKFAAAAIBzGjJkiPr27XvZ3q9OnTp64YUXLtv7mYnpzwEAAABUii1btiggIKDC269bt0633HKLjh8/rtDQ0CsX7BIwIgUAAAD8Ivd0rmanzVbu6Vyzo1yVIiMj5e/vb3aMy4IiBQAAAPwi9+dczfl2jnJ/rpwi1aVLFz3yyCMaOXKkqlWrpujoaL3++us6deqUHnjgAQUFBal+/fpauXKlJOn48eMaOHCgIiMj5efnpwYNGmjevHmO9zt48KAGDBig0NBQhYWFqU+fPtq/f3+FsthsNo0ePVqhoaEKDw/X448/LsMwzso7fPhwDR8+XCEhIYqIiND48ePP2u58fn9qn8Vi0RtvvKE77rhD/v7+atCggT766CNJ0v79+3XLLbdIkqpVqyaLxaIhQ4ZUaD+VgSIFAAAAmOjtt99WRESEUlNT9cgjjygpKUl33XWXOnTooG3btql79+4aNGiQTp8+rfHjx2v37t1auXKl9uzZozlz5igiIkKSVFpaqsTERAUFBWnjxo3atGmTAgMD1aNHD5WUlFw0x3PPPaf58+frrbfe0hdffKG8vDwtXbr0nHk9PT2VmpqqF198UTNnztQbb7xxycc/adIkDRgwQNu3b9ftt9+ugQMHKi8vTzVr1tR///tfSVJ6erqysrL04osvXvJ+LjeukQIAAMA1Lfd0rmMEas+xPU7/laRIv0hF+kdesf03b95cTz31lCRp3LhxmjZtmiIiIvTQQw9Jkp5++mnNmTNH27dvV2Zmplq2bKk2bdpIKh/hOeP999+X3W7XG2+8IYvFIkmaN2+eQkNDtW7dOnXv3v2COV544QWNGzdO/fr1kyTNnTtXq1evPmu7mjVr6vnnn5fFYlGjRo20Y8cOPf/88468rhoyZIjuvfdeSdKUKVM0a9YspaamqkePHgoLC5MkRUVFud01UhQpAAAAXNOWfL9Ec76d47RsYspEx9dJzZM0rMWwK7b/Zs2aOb62Wq0KDw9X06ZNHcuio6MlSTk5OUpKSlL//v0dI1V9+/ZVhw4dJEnffvutfvjhBwUFBTm9f1FRkfbt23fBDAUFBcrKylK7du0cyzw9PdWmTZuzTttr3769o6hJUkJCgp577jnZbDZZrVYXj975+AMCAhQcHKycnByX36eyUaQAAABwTbur4V3qUrOLpPKRqIkpEzUxYaKahDeRVD4idSV5eXk5PbdYLE7LzpQWu92unj176sCBA/rf//6nNWvWqGvXrkpOTtazzz6rkydPqnXr1nr33XfP2kdk5JU9hj/iXMdvt9tNSlNxFCkAAABc0yL9zz51r0l4E8WHx5uU6MIiIyM1ePBgDR48WJ06ddLYsWP17LPPqlWrVnr//fcVFRWl4OBgl94zJCREsbGx2rx5szp37ixJKisr09atW9WqVSunbTdv3uz0/KuvvlKDBg0uaTTqYry9vSWVT4ThbphsAgAAAKginn76aX344Yf64YcftGvXLq1YsUJNmpSPnA0cOFARERHq06ePNm7cqIyMDK1bt06PPvqoDh06dNH3HjFihKZNm6Zly5bpu+++07Bhw5Sfn3/WdpmZmRo9erTS09O1cOFCvfTSSxoxYsTlPlRJUu3atWWxWLRixQrl5ubq5MmTV2Q/l4IiBQAAAPwi0i9SSc2TrvjpfJfK29tb48aNU7NmzdS5c2dZrVYtWrRIkuTv768NGzaoVq1a6tevn5o0aaIHH3xQRUVFFRqhGjNmjAYNGqTBgwcrISFBQUFBuuOOO87a7v7779fPP/+stm3bKjk5WSNGjNDDDz982Y9VkqpXr65JkybpySefVHR0tIYPH35F9nMpLEZFJ32/ihUWFiokJEQFBQUuD4MCAADAXEVFRcrIyFDdunXl6+trdpyrWpcuXdSiRQune0FVRRf6maloN2BECgAAAABcRJECAAAArgGBgYHnfWzcuPEPv//GjRsvuI+rDbP2AQAAANeAtLS0866rXr16hd5j3bp1513Xpk2bC+7jakORAgAAAK4B9evXv6Lv7+fnd8X34U44tQ8AAABXBeZQQ0Vdjp8VihQAAACqNC8vL0nS6dOnTU6CquLMz8qZn51Lwal9AAAAqNKsVqtCQ0OVk5Mjqfx+ShaLxeRUcEeGYej06dPKyclRaGiorFbrJb8XRQoAAABVXkxMjCQ5yhRwIaGhoY6fmUtFkQIAAECVZ7FYFBsbq6ioKJWWlpodB27My8vrD41EnUGRAgAAwFXDarVelr8kAxfDZBMAAAAA4CKKFAAAAAC4iCIFAAAAAC6iSAEAAACAiyhSAAAAAOAi04vU4cOH9Ze//EXh4eHy8/NT06ZN9fXXXzvWG4ahp59+WrGxsfLz81O3bt20d+9ep/fIy8vTwIEDFRwcrNDQUD344IM6efJkZR8KAAAAgGuEqUXq+PHjuummm+Tl5aWVK1dq9+7deu6551StWjXHNjNmzNCsWbM0d+5cbd68WQEBAUpMTFRRUZFjm4EDB2rXrl1as2aNVqxYoQ0bNujhhx8245AAAAAAXAMshmEYZu38ySef1KZNm7Rx48ZzrjcMQ3FxcRozZowee+wxSVJBQYGio6M1f/583XPPPdqzZ4/i4+O1ZcsWtWnTRpK0atUq3X777Tp06JDi4uIumqOwsFAhISEqKChQcHDw5TtAAAAAAFVKRbuBqSNSH330kdq0aaO77rpLUVFRatmypV5//XXH+oyMDGVnZ6tbt26OZSEhIWrXrp1SUlIkSSkpKQoNDXWUKEnq1q2bPDw8tHnz5nPut7i4WIWFhU4PAAAAAKgoU4vUjz/+qDlz5qhBgwZavXq1kpKS9Oijj+rtt9+WJGVnZ0uSoqOjnV4XHR3tWJedna2oqCin9Z6engoLC3Ns83tTp05VSEiI41GzZs3LfWgAAAAArmKmFim73a5WrVppypQpatmypR5++GE99NBDmjt37hXd77hx41RQUOB4HDx48IruDwAAAMDVxdQiFRsbq/j4eKdlTZo0UWZmpiQpJiZGknT06FGnbY4ePepYFxMTo5ycHKf1ZWVlysvLc2zzez4+PgoODnZ6AAAAAEBFmVqkbrrpJqWnpzst+/7771W7dm1JUt26dRUTE6O1a9c61hcWFmrz5s1KSEiQJCUkJCg/P19bt251bPPZZ5/JbrerXbt2lXAUAAAAAK41nmbufNSoUerQoYOmTJmiAQMGKDU1Va+99ppee+01SZLFYtHIkSP1z3/+Uw0aNFDdunU1fvx4xcXFqW/fvpLKR7B69OjhOCWwtLRUw4cP1z333FOhGfsAAAAAwFWmTn8uSStWrNC4ceO0d+9e1a1bV6NHj9ZDDz3kWG8YhiZMmKDXXntN+fn56tixo2bPnq2GDRs6tsnLy9Pw4cO1fPlyeXh4qH///po1a5YCAwMrlIHpzwEAAABIFe8Gphcpd0CRAgAAACBVkftIAQAAAEBVRJECAAAAABdRpAAAAADARRQpAAAAAHARRQoAAAAAXESRAgAAAAAXUaQAAAAAwEUUKQAAAABwEUUKAAAAAFxEkQIAAAAAF1GkAAAAAMBFFCkAAAAAcBFFCgAAAABcRJECAAAAABdRpAAAAADARRQpAAAAAHARRQoAAAAAXESRAgAAAAAXUaQAAAAAwEUUKQAAAABwEUUKAAAAAFxEkQIAAAAAF1GkAAAAAMBFFCkAAAAAcBFFCgAAAABcRJECAAAAABdRpAAAAADARRQpAAAAAHARRQoAAAAAXESRAgAAAAAXUaQAAAAAwEUUKQAAAABwEUUKAAAAAFxEkQIAAAAAF1GkAAAAAMBFFCkAAAAAcBFFCgAAAABcRJECAAAAABdRpAAAAADARRQpAAAAAHARRQoAAAAAXESRAgAAAAAXUaQAAAAAwEUUKQAAAABwEUUKAAAAAFxEkQIAAAAAF1GkAAAAAMBFFCkAAAAAcBFFCgAAAABcRJECAAAAABeZWqQmTpwoi8Xi9GjcuLFjfVFRkZKTkxUeHq7AwED1799fR48edXqPzMxM9erVS/7+/oqKitLYsWNVVlZW2YcCAAAA4BriaXaA66+/Xp9++qnjuafnr5FGjRqljz/+WEuWLFFISIiGDx+ufv36adOmTZIkm82mXr16KSYmRl9++aWysrJ0//33y8vLS1OmTKn0YwEAAABwbTC9SHl6eiomJuas5QUFBXrzzTf13nvv6dZbb5UkzZs3T02aNNFXX32l9u3b65NPPtHu3bv16aefKjo6Wi1atNDkyZP1xBNPaOLEifL29q7swwEAAABwDTD9Gqm9e/cqLi5O1113nQYOHKjMzExJ0tatW1VaWqpu3bo5tm3cuLFq1aqllJQUSVJKSoqaNm2q6OhoxzaJiYkqLCzUrl27zrvP4uJiFRYWOj0AAAAAoKJMLVLt2rXT/PnztWrVKs2ZM0cZGRnq1KmTTpw4oezsbHl7eys0NNTpNdHR0crOzpYkZWdnO5WoM+vPrDufqVOnKiQkxPGoWbPm5T0wAAAAAFc1U0/t69mzp+PrZs2aqV27dqpdu7YWL14sPz+/K7bfcePGafTo0Y7nhYWFlCkAAAAAFWb6qX2/FRoaqoYNG+qHH35QTEyMSkpKlJ+f77TN0aNHHddUxcTEnDWL35nn57ru6gwfHx8FBwc7PQAAAACgotyqSJ08eVL79u1TbGysWrduLS8vL61du9axPj09XZmZmUpISJAkJSQkaMeOHcrJyXFss2bNGgUHBys+Pr7S8wMAAAC4Nph6at9jjz2m3r17q3bt2jpy5IgmTJggq9Wqe++9VyEhIXrwwQc1evRohYWFKTg4WI888ogSEhLUvn17SVL37t0VHx+vQYMGacaMGcrOztZTTz2l5ORk+fj4mHloAAAAAK5iphapQ4cO6d5779WxY8cUGRmpjh076quvvlJkZKQk6fnnn5eHh4f69++v4uJiJSYmavbs2Y7XW61WrVixQklJSUpISFBAQIAGDx6sZ555xqxDAgAAAHANsBiGYZgdwmyFhYUKCQlRQUEB10sBAAAA17CKdgO3ukYKAAAAAKoCihQAAAAAuIgiBQAAAAAuokgBAAAAgIsoUgAAAADgIooUAAAAALiIIgUAAAAALqJIAQAAAICLKFIAAAAA4CKKFAAAAAC4iCIFAAAAAC6iSAEAAACAiyhSAAAAAOAiihQAAAAAuIgiBQAAAAAuokgBAAAAgIsoUgAAAADgIooUAAAAALiIIgUAAAAALqJIAQAAAICLKFIAAAAA4CKKFAAAAAC4iCIFAAAAAC6iSAEAAACAiyhSAAAAAOAiihQAAAAAuMjT7ACoWmx2Q6kZeco5UaSoIF+1rRsmq4fF7FgAAABApaJIocJW7czSpOW7VfP0BjUN+4925N2pg/6dNaF3vHrcEGt2PAAAAKDScGofKmTVziwlLdimrIKfdZ/Ph1oc5qH7fD5UdsHPSlqwTat2ZpkdEQAAAKg0FClclM1uaNLy3TIkdfbYroYehyRJDT0OqZPHdknSpOW7ZbMbJqYEAAAAKg9FCheVmpGn7FM58vA9pH7+i7XT20eStNPbR/38F8vie0jZp3KUmpFnclIAAACgcnCNFC4q50SRvEI3yydyrcZLksIkSZMjwyTZFKCXVZzbVTknOpuYEgAAAKg8FClcVFSQr0rz22pG0TrVsxzR9z6emhgZrom5x9SwuEz7jDiNKmqrqCBfs6MCAAAAlYIihYtqWzdMffyOqG9phiTJarFLkpqUlCi+tFRNlaF1fkfUtm6YmTEBAACASsM1Urgoq0WaGLhUduPc94uyGxZNDFwqK7eTAgAAwDWCIoWLs5UotCRHHpbyWfkibTYlHS9QpM0mSfKwGAotzZFsJWamBAAAACoNp/bh4jx9pIc/l079JJthKOtwoeJPlyjL31th1YNltVikgMjy7QAAAIBrAEUKFRNSQwqpIaukZtXNDgMAAACYi1P7AAAAAMBFFCkAAAAAcBFFCgAAAABcRJECAAAAABdRpAAAAADARRQpAAAAAHARRQoAAAAAXESRAgAAAAAXUaQAAAAAwEUUKQAAAABwEUUKAAAAAFxEkQIAAAAAF1GkAAAAAMBFblOkpk2bJovFopEjRzqWFRUVKTk5WeHh4QoMDFT//v119OhRp9dlZmaqV69e8vf3V1RUlMaOHauysrJKTg8AAADgWuIWRWrLli169dVX1axZM6flo0aN0vLly7VkyRKtX79eR44cUb9+/RzrbTabevXqpZKSEn355Zd6++23NX/+fD399NOVfQgAAAAAriGmF6mTJ09q4MCBev3111WtWjXH8oKCAr355puaOXOmbr31VrVu3Vrz5s3Tl19+qa+++kqS9Mknn2j37t1asGCBWrRooZ49e2ry5Ml65ZVXVFJSYtYhAQAAALjKmV6kkpOT1atXL3Xr1s1p+datW1VaWuq0vHHjxqpVq5ZSUlIkSSkpKWratKmio6Md2yQmJqqwsFC7du067z6Li4tVWFjo9AAAAACAivI0c+eLFi3Stm3btGXLlrPWZWdny9vbW6GhoU7Lo6OjlZ2d7djmtyXqzPoz685n6tSpmjRp0h9MDwAAAOBaZdqI1MGDBzVixAi9++678vX1rdR9jxs3TgUFBY7HwYMHK3X/AAAAAKo204rU1q1blZOTo1atWsnT01Oenp5av369Zs2aJU9PT0VHR6ukpET5+flOrzt69KhiYmIkSTExMWfN4nfm+ZltzsXHx0fBwcFODwAAAACoKNOKVNeuXbVjxw6lpaU5Hm3atNHAgQMdX3t5eWnt2rWO16SnpyszM1MJCQmSpISEBO3YsUM5OTmObdasWaPg4GDFx8dX+jEBAAAAuDaYdo1UUFCQbrjhBqdlAQEBCg8Pdyx/8MEHNXr0aIWFhSk4OFiPPPKIEhIS1L59e0lS9+7dFR8fr0GDBmnGjBnKzs7WU089peTkZPn4+FT6MQEAAAC4Npg62cTFPP/88/Lw8FD//v1VXFysxMREzZ4927HearVqxYoVSkpKUkJCggICAjR48GA988wzJqYGAAAAcLWzGIZhmB3CbIWFhQoJCVFBQQHXSwEAAADXsIp2A9PvIwUAAAAAVQ1FCgAAAABcRJECAAAAABdRpAAAAADARS4VqRkzZujnn392PN+0aZOKi4sdz0+cOKFhw4ZdvnQAAAAA4IZcmrXParUqKytLUVFRkqTg4GClpaXpuuuukyQdPXpUcXFxstlsVybtFcKsfQAAAACkKzRr3+87FzOnAwAAALgWcY0UAAAAALiIIgUAAAAALvJ09QVvvPGGAgMDJUllZWWaP3++IiIiJJVPNgEAAAAAVzuXJpuoU6eOLBbLRbfLyMj4Q6EqG5NNAAAAAJAq3g1cGpHav3//H80FAAAAAFUe10gBAAAAgItcKlIpKSlasWKF07J33nlHdevWVVRUlB5++GGnG/QCAAAAwNXIpSL1zDPPaNeuXY7nO3bs0IMPPqhu3brpySef1PLlyzV16tTLHhIAAAAA3IlLRSotLU1du3Z1PF+0aJHatWun119/XaNHj9asWbO0ePHiyx4SAAAAANyJS0Xq+PHjio6Odjxfv369evbs6Xh+44036uDBg5cvHQAAAAC4IZeKVHR0tGNq85KSEm3btk3t27d3rD9x4oS8vLwub0IAAAAAcDMuFanbb79dTz75pDZu3Khx48bJ399fnTp1cqzfvn276tWrd9lDAgAAAIA7cek+UpMnT1a/fv108803KzAwUPPnz5e3t7dj/VtvvaXu3btf9pAAAAAA4E4shmEYrr6ooKBAgYGBslqtTsvz8vIUFBRU5U7vq+jdiwEAAABc3SraDVwakRo6dGiFtnvrrbdceVsAAAAAqFJcKlLz589X7dq11bJlS13CQBYAAAAAXBVcKlJJSUlauHChMjIy9MADD+gvf/mLwsLCrlQ2AAAAAHBLLs3a98orrygrK0uPP/64li9frpo1a2rAgAFavXo1I1QAAAAArhmXNNnEGQcOHND8+fP1zjvvqKysTLt27VJgYODlzFcpmGwCAAAAgFTxbuDSiNRZL/bwkMVikWEYstlsf+StAAAAAKDKcLlIFRcXa+HChbrtttvUsGFD7dixQy+//LIyMzOr5GgUAAAAALjKpckmhg0bpkWLFqlmzZoaOnSoFi5cqIiIiCuVDQAAAADckkvXSHl4eKhWrVpq2bKlLBbLebf74IMPLku4ysI1UgAAAACkK3RD3vvvv/+CBQoAAAAArgUu35AXAAAAAK51f2jWPgAAAAC4FlGkAAAAAMBFFCkAAAAAcBFFCgAAAABcRJECAAAAABdRpAAAAADARRQpAAAAAHARRQoAAAAAXESRAgAAAAAXUaQAAAAAwEUUKQAAAABwEUUKAAAAAFxEkQIAAAAAF1GkAAAAAMBFFCkAAAAAcBFFCgAAAABcRJECAAAAABdRpAAAAADARaYWqTlz5qhZs2YKDg5WcHCwEhIStHLlSsf6oqIiJScnKzw8XIGBgerfv7+OHj3q9B6ZmZnq1auX/P39FRUVpbFjx6qsrKyyDwVuxmY3lLLvmD5MO6yUfcdksxtmRwIAAMBVxNPMndeoUUPTpk1TgwYNZBiG3n77bfXp00fffPONrr/+eo0aNUoff/yxlixZopCQEA0fPlz9+vXTpk2bJEk2m029evVSTEyMvvzyS2VlZen++++Xl5eXpkyZYuahwUSrdmZp0vLduu7EFk30fEcTy+7Xj0E3akLvePW4IdbseAAAALgKWAzDcKt/qg8LC9O///1v3XnnnYqMjNR7772nO++8U5L03XffqUmTJkpJSVH79u21cuVK/elPf9KRI0cUHR0tSZo7d66eeOIJ5ebmytvbu0L7LCwsVEhIiAoKChQcHHzFjg1X3qqdWUpasE2GDM33e0q7Q/IUXxCmB37+pySL5vylFWUKAAAA51XRbuA210jZbDYtWrRIp06dUkJCgrZu3arS0lJ169bNsU3jxo1Vq1YtpaSkSJJSUlLUtGlTR4mSpMTERBUWFmrXrl3n3VdxcbEKCwudHqj6bHZDk5bvliGps8d2RXkd0pxqIYryOqROHtslSZOW7+Y0PwAAAPxhphepHTt2KDAwUD4+Pvr73/+upUuXKj4+XtnZ2fL29lZoaKjT9tHR0crOzpYkZWdnO5WoM+vPrDufqVOnKiQkxPGoWbPm5T0omCI1I09ZBUWSDI3xXCKbUf7jbTM8NMZziQwZyiooUmpGnrlBAQAAUOWZXqQaNWqktLQ0bd68WUlJSRo8eLB27959Rfc5btw4FRQUOB4HDx68ovtD5cg5USSLZ6Ha+K2Tl+9Bfe9Tfgng9z6e8vI9qDZ+62TxLFTOiSKTkwIAAKCqM3WyCUny9vZW/fr1JUmtW7fWli1b9OKLL+ruu+9WSUmJ8vPznUaljh49qpiYGElSTEyMUlNTnd7vzKx+Z7Y5Fx8fH/n4+FzmI4HZooJ85RW6WemRa3W3fr0OamJk+C9frZZXbpmigrqbExAAAABXDdNHpH7PbreruLhYrVu3lpeXl9auXetYl56erszMTCUkJEiSEhIStGPHDuXk5Di2WbNmjYKDgxUfH1/p2WGutnXD1LU0TO8fztL7h7M0MfeYJGli7jHHsq6lYWpbN8zkpAAAAKjqTB2RGjdunHr27KlatWrpxIkTeu+997Ru3TqtXr1aISEhevDBBzV69GiFhYUpODhYjzzyiBISEtS+fXtJUvfu3RUfH69BgwZpxowZys7O1lNPPaXk5GRGnK5BVos01X+1gvPK5GH5dUKJJiUlii8pld2waGrYalkt401MCQAAgKuBqUUqJydH999/v7KyshQSEqJmzZpp9erVuu222yRJzz//vDw8PNS/f38VFxcrMTFRs2fPdrzearVqxYoVSkpKUkJCggICAjR48GA988wzZh0SzGQrUWhJjmQ596x8HhZDoaU5kq1E8qRoAwAA4NK53X2kzMB9pK4iBYekUz/JZhj6Yv9BfZy7Xr0ib1bHOjVltVikgEgppLrZKQEAAOCmKtoNKFKiSAEAAAAoV+VuyAsAAAAAVQVFCgAAAABcRJECAAAAABdRpAAAAADARRQpAAAAAHARRQoAAAAAXESRAgAAAAAXUaQAAAAAwEUUKQAAAABwEUUKAAAAAFxEkQIAAAAAF1GkAAAAAMBFFCkAAAAAcBFFCgAAAABcRJECAAAAABdRpAAAAADARRQpAAAAAHARRQoAAAAAXESRAgAAAAAXUaQAAAAAwEUUKQAAAABwEUUKAAAAAFxEkQIAAAAAF1GkAAAAAMBFFCkAAAAAcBFFCgAAAABcRJECAAAAABdRpAAAAADARRQpAAAAAHARRQoAAAAAXESRAgAAAAAXUaQAAAAAwEUUKQAAAABwEUUKAAAAAFxEkQIAAAAAF1GkAAAAAMBFFCkAAAAAcBFFCgAAAABcRJECAAAAABdRpAAAAADARRQpwCQ2u6FN6xbo+dnNtGndAtnshtmRAAAAUEEUKcAEq3ZmqeO0tcr94lm9FWAo94tn1XHaWq3amWV2NAAAAFQARQqoZKt2ZilpwTY1OJmqhh6HJEkNPQ6p4clUJS3YRpkCAACoAihSQCWy2Q1N+DhFFt9D6ue/WDu9fSRJO719dIf/Ynn4HtKEj1M4zQ8AAMDNeZodALiWpGbk6bh1gwLqrtV4SVKYJGlyZJgkm/z1so7ndlVqxs1KqBduYlIAAABcCEUKqEQ5J4pUmt9WM4rWqZ7liL738dTEyHBNzD2mhsVl2mfEaVRRW+WcKDI7KgAAAC6AIgVUoqggX3Wy71ffsgxJktVilyQ1KSlRfGmpmipDH9j3Kyoo0cyYAAAAuAiKFFCJ2tapphCf/8put8jDcvZ1UHbDoid9/qtGdZ40IR0AAAAqiskmgEpkNUpVz/u4o0RF2mxKOl6gSJtNkuRhMVTPO19Wo9TMmAAAALgIRqSAyuTpI59hG/Tl9u/06oYf9dPJEum09LGkiEBv/a3zderQrInk6WN2UgAAAFyAqSNSU6dO1Y033qigoCBFRUWpb9++Sk9Pd9qmqKhIycnJCg8PV2BgoPr376+jR486bZOZmalevXrJ399fUVFRGjt2rMrKyirzUICKC6mhDp266a1xD+mpv96rh+++Q0/99V69Ne4hdejUTQqpbnZCAAAAXISpRWr9+vVKTk7WV199pTVr1qi0tFTdu3fXqVOnHNuMGjVKy5cv15IlS7R+/XodOXJE/fr1c6y32Wzq1auXSkpK9OWXX+rtt9/W/Pnz9fTTT5txSECFWT0sSqgXrj4tqiuhXrisHhazIwEAAKCCLIZhuM2dP3NzcxUVFaX169erc+fOKigoUGRkpN577z3deeedkqTvvvtOTZo0UUpKitq3b6+VK1fqT3/6k44cOaLo6GhJ0ty5c/XEE08oNzdX3t7eF91vYWGhQkJCVFBQoODg4Ct6jAAAAADcV0W7gVtNNlFQUCBJCgsrv0np1q1bVVpaqm7dujm2ady4sWrVqqWUlBRJUkpKipo2beooUZKUmJiowsJC7dq165z7KS4uVmFhodMDAAAAACrKbYqU3W7XyJEjddNNN+mGG26QJGVnZ8vb21uhoaFO20ZHRys7O9uxzW9L1Jn1Z9ady9SpUxUSEuJ41KxZ8zIfDQAAAICrmdsUqeTkZO3cuVOLFi264vsaN26cCgoKHI+DBw9e8X0CAAAAuHq4xfTnw4cP14oVK7RhwwbVqFHDsTwmJkYlJSXKz893GpU6evSoYmJiHNukpqY6vd+ZWf3ObPN7Pj4+8vFhemkAAAAAl8bUESnDMDR8+HAtXbpUn332merWreu0vnXr1vLy8tLatWsdy9LT05WZmamEhARJUkJCgnbs2KGcnBzHNmvWrFFwcLDi4+Mr50AAAAAAXFNMHZFKTk7We++9pw8//FBBQUGOa5pCQkLk5+enkJAQPfjggxo9erTCwsIUHBysRx55RAkJCWrfvr0kqXv37oqPj9egQYM0Y8YMZWdn66mnnlJycjKjTgAAAACuCFOnP7dYzn3fnHnz5mnIkCGSym/IO2bMGC1cuFDFxcVKTEzU7NmznU7bO3DggJKSkrRu3ToFBARo8ODBmjZtmjw9K9YTmf4cAAAAgFTxbuBW95EyC0UKAAAAgFRF7yMFAAAAAFUBRQoAAAAAXESRAgAAAAAXUaQAAAAAwEUUKQAAAABwEUUKAAAAAFxEkQIAAAAAF1GkAAAAAMBFFCkAAAAAcBFFCgAAAABcRJECAAAAABdRpABckM1uaPcXH6rw2Zba/cWHstkNsyMBAACYjiIF4LxW7cxSx2lrlbN2ghZ4HlPO2gnqOG2tVu3MMjsaAACAqShSAM5p1c4sJS3YpgYnUxXldUhzqoUoyuuQGp5MVdKCbZQpAABwTaNIATiLzW5o0vLdMmRojOcS2YzyXxU2w0OjPZdIKl/PaX4AAOBaRZECcJbUjDxln8pRG7918vI9qO99PCVJ3/t4ysv3oFr7rVP2qRylZuSZnBQAAMAcnmYHAOB+ck4UySt0s9Ij1+puxTqWT4wM/+Wr1fLKLVPOic7mBAQAADAZRQrAWaKCfNW6IEjjS8qvg9rj7a2JkeGamHtMTUpKJEmTfw5SVJCvmTEBAABMQ5ECcJa2dappvOf/1Li4TB6WX6+DalJSoviSUtkNi8Z7/k+N6kwyMSUAAIB5uEYKwFmsRqnqeR93KlG/5WExVM87X1ajtJKTAQAAuAdGpACczdNHPsM26Mvt3+nVDT/q6OnjCv5pj0afbqJov2r6W+fr1KFZE8nTx+ykAAAApqBIATi3kBrq0KmG2t1kKDUjTzknihQV5Ku2dcNk9bCYnQ4AAMBUFCkAF2T1sCihXvjFNwQAALiGcI0UAAAAALiIIgUAAAAALqJIAQAAAICLKFIAAAAA4CKKFAAAAAC4iCIFAAAAAC6iSAEAAACAiyhSAAAAAOAiihQAAAAAuIgiBQAAAAAu8jQ7AAD8UTa7odSMPOWcKFJUkK/a1g2T1cNidiwAAHAVo0gBqNJW7czSpOW7lVVQ5FgWG+KrCb3j1eOGWBOTAQCAqxmn9gGoslbtzFLSgm3KKihSW69UPRj9uNp6pSq7oEhJC7Zp1c4ssyMCAICrFEUKQJVksxuatHy3DEmSoft8PtTiMA/d5/OhjF+WTlq+Wza7YWZMAABwlaJIAaiSUjPyHKfzdfbYroYehyRJDT0OqbPHdhmSsgqKlJqRZ2JKAABwtaJIAaiSck4UyeJZKA/fQ+rnv1g7vX0kSTu9fdTPf7E8fA/J4lmonBNFF3knAAAA1zHZBIAqKSrIV16hm+UTuVbjJUlhkqTJkWGSbArQyyrO7aqooO4mpgQAAFcrihSAKqlt3TBVK+ukcQfXqZ7liL738dTEyHBNzD2mhsVl2mfEaapnJ7WtG2Z2VAAAcBWiSAGokqweFr3SxqY2X2SUP7fYJUlNSkoUX1qqpspQjfY27icFAACuCK6RAlA1GYba/Dhbhs5dlAxZ1ObH2ZLBrH0AAODyo0gBqJpsJVLBYVl+meo80mZT0vECRdpsklS+vPBw+XYAAACXGaf2AaiaPH2khz+XTv0kSYqUNOz32wRElm8HAABwmVGkAFRdITXKHwAAAJWMIgUAlcRmN5SakaecE0WKCvJV27phTIYBAEAVRZECgEqwameWJi3frayCX28QHBviqwm949XjhlgTkwEAgEvBZBMAcIWt2pmlpAXblFVQpJs8dmiN91jd5LFD2QVFSlqwTat2ZpkdEQAAuIgiBQBXkM1uaNLy3b/MLWjocc/31cDjsB73fF/GL0snLd8tm51p2gEAqEpMLVIbNmxQ7969FRcXJ4vFomXLljmtNwxDTz/9tGJjY+Xn56du3bpp7969Ttvk5eVp4MCBCg4OVmhoqB588EGdPHmyEo8CAM4vNSPPcTpfZ4/tivPar9mhIYrz2q/OHttlSMoqKFJqRp65QQEAgEtMLVKnTp1S8+bN9corr5xz/YwZMzRr1izNnTtXmzdvVkBAgBITE1VU9Os1BgMHDtSuXbu0Zs0arVixQhs2bNDDDz9cWYcAABeUc+LM7ytDYzyXKNvDS3OqhSjbw0tjPJdIv4xK/bodAACoCkydbKJnz57q2bPnOdcZhqEXXnhBTz31lPr06SNJeueddxQdHa1ly5bpnnvu0Z49e7Rq1Spt2bJFbdq0kSS99NJLuv322/Xss88qLi6u0o4FAM4lKshXUvloVHOPH7Xb4iVJslrsaurxozp7bNcGe3PHdgAAoGpw22ukMjIylJ2drW7dujmWhYSEqF27dkpJSZEkpaSkKDQ01FGiJKlbt27y8PDQ5s2bz/vexcXFKiwsdHoAwJXQtm6YokOL1M9/sXZ4+WiPt7ckaY+3t3Z4+aif/2JFhxapbd0wk5MCAABXuO3059nZ2ZKk6Ohop+XR0dGOddnZ2YqKinJa7+npqbCwMMc25zJ16lRNmjTpMicGgLNZPSy6reFnGn/CJunX32cTI8N/+cqmPkGfyerR35R8AADg0rhtkbqSxo0bp9GjRzueFxYWqmbNmiYmAnDVMgyN+GmH7v0pWxYZ2uPtrYmR4ZqYe0xNSkpkyKKoiB2SYUgWbs4LAEBV4bZFKiYmRpJ09OhRxcb+erPKo0ePqkWLFo5tcnJynF5XVlamvLw8x+vPxcfHRz4+Ppc/NAD8nq1EkQVZiiwpcVrcpKRE8SWl5U8KsyRbieTJ7yUAAKoKty1SdevWVUxMjNauXesoToWFhdq8ebOSkpIkSQkJCcrPz9fWrVvVunVrSdJnn30mu92udu3amRUdAH7l6SM9/Ll06qfy54X7pc3jpX5vSMF1ypcFRFKiAACoYkwtUidPntQPP/zgeJ6RkaG0tDSFhYWpVq1aGjlypP75z3+qQYMGqlu3rsaPH6+4uDj17dtXktSkSRP16NFDDz30kObOnavS0lINHz5c99xzDzP2AXAfITXKH5IiQ6srqShJkbVukvwjTQ4GAAAulcUwDMOsna9bt0633HLLWcsHDx6s+fPnyzAMTZgwQa+99pry8/PVsWNHzZ49Ww0bNnRsm5eXp+HDh2v58uXy8PBQ//79NWvWLAUGBlY4R2FhoUJCQlRQUKDg4ODLcmwAUFXlns7Vku+X6K6GdymSsgcAuMZUtBuYWqTcBUUKAMrZ7IZWfj5X4w7N1tQaw9Tzlr/L6sEkGACAa0dFu4Hb3kcKAFC5Vu3MUsdpa1WW8qokqSzlVXWctlardmaZnAwAAPdDkQIAaPE3O5X8n+WqXrpKNt/y2VBtvjmqUbpKyf9ZrsXf7DQ5IQAA7oVT+8SpfQCubTa7oRtnj1Vp0OrzbuN1IlFbhv2b0/wAAFc9Tu0DAFRIakaejme3VqP9iXr/cJYm5h6TJE3MPab3D2ep0f5EHc9urdSMPJOTAgDgPihSAHCNyzlRJKMsSONta9Ww2KYmv9w8uElJiRoW2zTetlZGWZByThSZnBQAAPdBkQKAa1xUkK86e2xXc48f5WmxO63ztNjV3ONHdfbYrqggX5MSAgDgfky9IS8AwHxt61RTiM9/Zbdb5GExFGmzKel4gSJtNkmS3bDoSZ//qlGdJ01OCgCA+6BIAcA1zmqUqp73cXkUl889FGmza1h+gWO9h8VQPe98WY1SST4mpQQAwL1QpADgWufpI59hG/Tl9u/06oYf9dPJEseqiEBv/a3zderQrInk6R4lymY39En6Xq088IF61u6n7o0aMJsgAKDSMf25mP4cAM6w2Q2lZuQp50SRooJ81bZumFuVlFU7szRp+W7FlKzWD3VWqv7+nsr2TtSE3vHqcUOs2fEAAFcBpj8HALjM6mFRQr1w9WlRXQn1wt2uRCUt2Kasgp812PqJJGmw9RNlF/yspAXbtGpnlskJAQDXEooUAMDt2eyGJnycIovvYbXxWyebb075ct8ctfZbJw/fw5rwcYps9mv+JAsAQCXhGikAgNtLzcjTcesGBdRdq3RJExUuSZoYGS5Frpa/Vut4blelZtyshHrh5oYFAFwTGJECALi9nBNFKs1vp0b7E/X+4SxNzD0mSZqYe0zvH85So/2JKs1v51Y3DbbZDaXsO6YP0w4rZd8xRssA4CrDiBQAwO1FBfnKKAvSeI+1ami3SSqfWbBJSYkaFts03lirPmVd3OamwWcmxah5eoOahv1HO/Lu1EH/zkyKAQBXEUakAABur23dMPUN+k7NPX6Up8XutM7TYldzjx/VN+g7ta0bZlLCX/12Uoz7fD7U4jAP3efzIZNiAMBVhiIFAHB7Vos0MXCp7Eb5LIKRNpuSjhco0maTJNkNiyYGLpXV5EkGbXZDk5bvliGps8d2NfQ4JElq6HFInTy2S5ImLd/NaX4AcBWgSAEA3J+tRKElOfKwlBeQSJtdw/ILFGkrH53ysBgKLc2RbCUXepcrLjUjT9mncuThe0j9/Bdrp3f5TYx3evuon/9iWXwPKftUjlIz8kzNeQbXcQHApeMaKQCA+/P0kR7+XDr1k2yGoV2HC5V3ukRh/t66vnqwrBaLFBBZvp2Jck4UySt0s3wi12q8JKn8VMPJkWGSbArQyyrO7aqcE51NTFnuzHVcWQW/TtARG+LLdVwAUEEUKQBA1RBSQwqpIaukZtXNDnNuUUG+Ks1vqxlF61TPckTf+3hqYmS4JuYeU8PiMu0z4jSqqK3pk2KcuY7r9+NP2QVFSlqwTXP+0ooyBQAXwal9AABcJm3rhqmP3xH1LctQ09JiNSn5dXbBpqXF6luWoT5+R0ydFOO313FJUluvVD0Y/bjaeqU6lnEdFwBcHEUKAIDL5PeTYvyeO0yKkZqR95vT+QynmQUlQ4akrIIit7mOS+JaLgDuiVP7AAC4XH6ZFEOOSTGcZxd0mhTDpOu5fnvT4l9nFoxVQ49D6uyxXRvszc/azkxcywXAXVGkAAC4XH43KUbW4ULFny5Rlr+3wtxkUoyoIF9ZPAtl8SxUP8+zZxb8oixcRlmw6ddxSc7XcrX1SnXc3HhLQVuu5QJgOooUAACXk5tPitG2bpiqxWxVadDq884s6HUiUW3r3mNiyt9fy1V+CuL4MA9NPvWhUktvlEUWTVq+W7fFx8jqYfINxABck7hGCgCAa4jVw6InOwzW5INWvXfoqCbmHpMkTcw9pvcOHdXkg1Y92WGw6eXkt9dy/f7mxp09trvdtVw2u6FN6xbo+dnNtGndAq7jAq4BFCkAAK4xA4KPXnBmwQHBR01OWH6NlsWz8Lw3N/bwPSSLZ6FbXMu1ameWOk5bq9wvntVbAYZyv3hWHaet1aqdWWZHO4vNbmjlnu/16KppWrnnewof8Adwah8AANcSw5A++6fK/y3Vfo4NPMrX1+sqWcwblYoK8q3QzY2jgrqbllH69TquTh7fqmHgrxN3NDyZqqQFxW51HdeZiTtiSlbrhzordXBznp7xTmTiDuASMSIFAMC1xFYiFRzWmRL1+5kFJbtUeLh8OxO1rRumamWdLngKYrWyTqbfk2vCxymynGPU7I5fRs0mfJziFqM+ZwpfVsHPGmz9RJI02PqJsgt+VtKCbW45ega4O0akAAC4lvxmZkFJipQ07PfbmDyzoFR+LdcrbWxq80VG+XNLefFrUlKi+NJSNVWGarS3mXotV2pGno5bNyig7rlHzfz1so7ndlVqxs1KqBduWs5fC99Pam35XjbfHEnhsvnmqLVtnb4xGmrCx8W6Lf4O06+NOyN3zzIt2ThJd3WaoMgmfc2OA5wTRQoAgGvNLzMLujXDUJsfZ8uQRRadPaJjyKI2P86WjLtMOwUx50SRSvPbakbROtWzHNH3Pp6aGBmuibnH1LC4TPuMOI0qamv6dVy/LXzpkiaqvNRNjAyXIlfLX6vdovCdYbPZlfnZdM0JLFO7z6YrrOGfZbVyEhXcDz+VAADA/fxyCuKZEvX7UxAtMkw/BTEqyFed7PsvOHFHJ/t+0+/JVV742qnR/kS9fzjL6TTJ9w9nqdH+RJXmtzO98EnlpyCOmTZTfgV7JUl+BXs1ZtpMtz31MHfPMs1+raVy9ywzOwpMwIgUAABwP1XgFMS2daopxOe/stst8rCcPWpmNyx60ue/alTnSRPS/SoqyFdGWZDGe6xVQ7tN0q+Fr2GxTeONtepT1sX0wrf4m50a9+EmPee5QDt9fr3erEvRAiX/J0hTS2/SgJY3mJrxt6rKyBmnSV457venDQAAIJWffhjX4vyPEHPveGw1SlXP+7ijRP1+1MzDYqied76sRqmZMdW2bpj6Bn2n5h4/ytPiPFOjp8Wu5h4/qm/Qd6ZP3DFt03z5131Z42vafrnOrPx6s/E1bfKv+7KmbZrvFhN3SFVo5Mwon45/jk+Zcr94tnzWTjdUVUf2GJECAAC4FJ4+8hm2QV9u/06vbvhRP50skU5LH0uKCPTW3zpfpw7Nmpg/cYdFmhi4VPa88pGz3xc+u2HRxMClslrGmJYxNSNPx7Nb6fn8Ty9wvVkrpWbkmX4dV1UaObPtXSvlfidVj5Vyv5Nt71pZG3YzO5az35S9Ll88q8jGfUy99YIrKFIAAACXKqSGOnSqoXY3GUrNyFPOiSJFBfmqbd0wt5kBT7YShZbkSI6RM7uG5Rc4VntYDIWW5pRfb2ZS6cs5UeS43kw69yyNH9j3m34d168jZ6vPO1PjtE171b/5v0398889nasPd+xWzGdjdfo30/L/+N+xyr71BfVpGq9I/0jT8jnZ51z2tG+tVN/Nyt55UKQAAAD+IKuHxfSRkvP6zfVmNsPQrsOFyjtdojB/b11fPVhWi8X0682iAn00xnOJ7Mb5rzcb47lEpwMfNiHdr6rKyNnUL97SmqwFUpzkXPbKpO+Ga3fBXzSz+xOm5ZPKy17u6RwZn0/Unl9G9nb7+Mj4fKIsobGK9I9yn7J3HhQpAACAq90vU95bJTUz99Kyc2pbK1D5HsfkcZ5ZGj0shmpY8hRaK9DMmFVi5MxmN7RpW31NLrOet+xNzaovWzfD1FGzJd8v0Zxv50jekiLKy96kiDBJBdLH9yipeZKGtThrihm3QpECAACAqazevtpx+zI9u/RLSSqvU79cb3bmr/pj7rhJt3ibO7NgVRg5S83IU6PCdPX1vkDZO51u+qhZhL2TJh987bxlr6RpJ9OyVRRFCgAAAKbr0raVivxjNWn5bmUV/DqiExviqwm943XLDbEmpitXFUbOcgp/rlDZ2194vwnpytnshr5cuVovXGBkb+TK1erf4gb3udbwHChSAAAAcAs9bojVbfExbjtxR1UYOYsO8FCc5dgFp+WP1TH9HGDeXZBSfzymB4rfld1y/rL3QPG7Sv1xiBLqR5iQsGIoUgAAAHAbbj1xh9x/5OzG+rHq5/Nv2U7kylFRfil7Unnh8wiK1NL65uX8qeCE2leg7G0uOCGJIgUAAABcFdx55MzqYVHSnzsracE2SdJvx3vOpJvz51amZo0IDdafi/+pMEvhrwt/U/Yk6ZgRrOdDgys9mysoUgAAAICL3HnkrMcNsZrzl1ZnjZrF/DJq1sPkUbO2dcOkkOraXRCus0/sKy98MSHl5dSdUaQAAACAq4y7j5pN6B2vpAXbZNG5R80m9I53i6wXYjEM41xF8JpSWFiokJAQFRQUKDjYvYcQAQAAgKvBqp1Z573WzMxRs4p2A0akAAAAAFQ6dx41qwiKFAAAAABTuPO1Zhdj3gTyAAAAAFBFUaQAAAAAwEUUKQAAAABw0VVTpF555RXVqVNHvr6+ateunVJTU82OBAAAAOAqdVUUqffff1+jR4/WhAkTtG3bNjVv3lyJiYnKyckxOxoAAACAq9BVUaRmzpyphx56SA888IDi4+M1d+5c+fv766233jI7GgAAAICrUJUvUiUlJdq6dau6devmWObh4aFu3bopJSXlnK8pLi5WYWGh0wMAAAAAKqrKF6mffvpJNptN0dHRTsujo6OVnZ19ztdMnTpVISEhjkfNmjUrIyoAAACAq0SVL1KXYty4cSooKHA8Dh48aHYkAAAAAFWIp9kB/qiIiAhZrVYdPXrUafnRo0cVExNzztf4+PjIx8enMuIBAAAAuApV+REpb29vtW7dWmvXrnUss9vtWrt2rRISEkxMBgAAAOBqVeVHpCRp9OjRGjx4sNq0aaO2bdvqhRde0KlTp/TAAw+YHQ0AAADAVeiqKFJ33323cnNz9fTTTys7O1stWrTQqlWrzpqAAgAAAAAuB4thGIbZIcxWUFCg0NBQHTx4UMHBwWbHAQAAAGCSwsJC1axZU/n5+QoJCTnvdlfFiNQfdeLECUliGnQAAAAAkso7woWKFCNSKp+c4siRIwoKCpLFYjE1y5kG7M6jY1Uho1Q1claFjFLVyFkVMkpVIycZL5+qkLMqZJSqRs6qkFGqGjmrQkapauQko+sMw9CJEycUFxcnD4/zz83HiJQkDw8P1ahRw+wYToKDg93iB+lCqkJGqWrkrAoZpaqRsypklKpGTjJePlUhZ1XIKFWNnFUho1Q1claFjFLVyElG11xoJOqMKj/9OQAAAABUNooUAAAAALiIIuVmfHx8NGHCBPn4+Jgd5byqQkapauSsChmlqpGzKmSUqkZOMl4+VSFnVcgoVY2cVSGjVDVyVoWMUtXIScYrh8kmAAAAAMBFjEgBAAAAgIsoUgAAAADgIooUAAAAALiIIgUAAAAALqJIuYk5c+aoWbNmjhuRJSQkaOXKlWbHuqBp06bJYrFo5MiRZkdxmDhxoiwWi9OjcePGZsc6p8OHD+svf/mLwsPD5efnp6ZNm+rrr782O5ZDnTp1zvpeWiwWJScnmx3Nic1m0/jx41W3bl35+fmpXr16mjx5stxtHp0TJ05o5MiRql27tvz8/NShQwdt2bLF1EwbNmxQ7969FRcXJ4vFomXLljmtNwxDTz/9tGJjY+Xn56du3bpp7969bpXxgw8+UPfu3RUeHi6LxaK0tLRKzVeRnKWlpXriiSfUtGlTBQQEKC4uTvfff7+OHDniNhml8t+fjRs3VkBAgKpVq6Zu3bpp8+bNlZqxIjl/6+9//7ssFoteeOGFSssnXTzjkCFDzvrd2aNHD7fKKEl79uzRn//8Z4WEhCggIEA33nijMjMz3SrnuT6HLBaL/v3vf7tNxpMnT2r48OGqUaOG/Pz8FB8fr7lz51ZavormPHr0qIYMGaK4uDj5+/urR48elf47ferUqbrxxhsVFBSkqKgo9e3bV+np6U7bFBUVKTk5WeHh4QoMDFT//v119OjRSs1ZURQpN1GjRg1NmzZNW7du1ddff61bb71Vffr00a5du8yOdk5btmzRq6++qmbNmpkd5SzXX3+9srKyHI8vvvjC7EhnOX78uG666SZ5eXlp5cqV2r17t5577jlVq1bN7GgOW7Zscfo+rlmzRpJ01113mZzM2fTp0zVnzhy9/PLL2rNnj6ZPn64ZM2bopZdeMjuak7/+9a9as2aN/u///k87duxQ9+7d1a1bNx0+fNi0TKdOnVLz5s31yiuvnHP9jBkzNGvWLM2dO1ebN29WQECAEhMTVVRU5DYZT506pY4dO2r69OmVlul8Oc6X8/Tp09q2bZvGjx+vbdu26YMPPlB6err+/Oc/u01GSWrYsKFefvll7dixQ1988YXq1Kmj7t27Kzc3161ynrF06VJ99dVXiouLq6Rkv6pIxh49ejj9Dl24cGElJrx4xn379qljx45q3Lix1q1bp+3bt2v8+PHy9fV1q5y//R5mZWXprbfeksViUf/+/d0m4+jRo7Vq1SotWLBAe/bs0ciRIzV8+HB99NFHlZbxYjkNw1Dfvn31448/6sMPP9Q333yj2rVrq1u3bjp16lSlZVy/fr2Sk5P11Vdfac2aNSotLVX37t2dMowaNUrLly/XkiVLtH79eh05ckT9+vWrtIwuMeC2qlWrZrzxxhtmxzjLiRMnjAYNGhhr1qwxbr75ZmPEiBFmR3KYMGGC0bx5c7NjXNQTTzxhdOzY0ewYLhkxYoRRr149w263mx3FSa9evYyhQ4c6LevXr58xcOBAkxKd7fTp04bVajVWrFjhtLxVq1bGP/7xD5NSOZNkLF261PHcbrcbMTExxr///W/Hsvz8fMPHx8dYuHChCQnPzvhbGRkZhiTjm2++qdRM53KhnGekpqYakowDBw5UTqjfqUjGgoICQ5Lx6aefVk6oczhfzkOHDhnVq1c3du7cadSuXdt4/vnnKz3bGefKOHjwYKNPnz6m5DmXc2W8++67jb/85S/mBDqPivxc9unTx7j11lsrJ9A5nCvj9ddfbzzzzDNOy8z+/f77nOnp6YYkY+fOnY5lNpvNiIyMNF5//XUTEpbLyckxJBnr1683DKP8c8bLy8tYsmSJY5s9e/YYkoyUlBSzYp4XI1JuyGazadGiRTp16pQSEhLMjnOW5ORk9erVS926dTM7yjnt3btXcXFxuu666zRw4MBKP02hIj766CO1adNGd911l6KiotSyZUu9/vrrZsc6r5KSEi1YsEBDhw6VxWIxO46TDh06aO3atfr+++8lSd9++62++OIL9ezZ0+RkvyorK5PNZjvrX3r9/PzccsRUkjIyMpSdne30/3lISIjatWunlJQUE5NdHQoKCmSxWBQaGmp2lHMqKSnRa6+9ppCQEDVv3tzsOE7sdrsGDRqksWPH6vrrrzc7znmtW7dOUVFRatSokZKSknTs2DGzIznY7XZ9/PHHatiwoRITExUVFaV27dpd8DRKd3D06FF9/PHHevDBB82O4qRDhw766KOPdPjwYRmGoc8//1zff/+9unfvbnY0h+LiYkly+hzy8PCQj4+PqZ9DBQUFkqSwsDBJ0tatW1VaWur02dO4cWPVqlXLLT97KFJuZMeOHQoMDJSPj4/+/ve/a+nSpYqPjzc7lpNFixZp27Ztmjp1qtlRzqldu3aaP3++Vq1apTlz5igjI0OdOnXSiRMnzI7m5Mcff9ScOXPUoEEDrV69WklJSXr00Uf19ttvmx3tnJYtW6b8/HwNGTLE7ChnefLJJ3XPPfeocePG8vLyUsuWLTVy5EgNHDjQ7GgOQUFBSkhI0OTJk3XkyBHZbDYtWLBAKSkpysrKMjveOWVnZ0uSoqOjnZZHR0c71uHSFBUV6YknntC9996r4OBgs+M4WbFihQIDA+Xr66vnn39ea9asUUREhNmxnEyfPl2enp569NFHzY5yXj169NA777yjtWvXavr06Vq/fr169uwpm81mdjRJUk5Ojk6ePKlp06apR48e+uSTT3THHXeoX79+Wr9+vdnxzuvtt99WUFCQ253m9dJLLyk+Pl41atSQt7e3evTooVdeeUWdO3c2O5rDmTIybtw4HT9+XCUlJZo+fboOHTpk2ueQ3W7XyJEjddNNN+mGG26QVP7Z4+3tfdY/MrnrZ4+n2QHwq0aNGiktLU0FBQX6z3/+o8GDB2v9+vVuU6YOHjyoESNGaM2aNZV+DnVF/XYUolmzZmrXrp1q166txYsXu9W/YNntdrVp00ZTpkyRJLVs2VI7d+7U3LlzNXjwYJPTne3NN99Uz549TbkW4WIWL16sd999V++9956uv/56paWlaeTIkYqLi3Or7+X//d//aejQoapevbqsVqtatWqle++9V1u3bjU7GipRaWmpBgwYIMMwNGfOHLPjnOWWW25RWlqafvrpJ73++usaMGCANm/erKioKLOjSSr/1+oXX3xR27Ztc7vR8d+65557HF83bdpUzZo1U7169bRu3Tp17drVxGTl7Ha7JKlPnz4aNWqUJKlFixb68ssvNXfuXN18881mxjuvt956SwMHDnS7v4O89NJL+uqrr/TRRx+pdu3a2rBhg5KTkxUXF+c2Z+94eXnpgw8+0IMPPqiwsDBZrVZ169ZNPXv2NG1ypuTkZO3cudNtz8yoCEak3Ii3t7fq16+v1q1ba+rUqWrevLlefPFFs2M5bN26VTk5OWrVqpU8PT3l6emp9evXa9asWfL09HSbf2n7rdDQUDVs2FA//PCD2VGcxMbGnlWQmzRp4panIR44cECffvqp/vrXv5od5ZzGjh3rGJVq2rSpBg0apFGjRrndqGm9evW0fv16nTx5UgcPHlRqaqpKS0t13XXXmR3tnGJiYiTprJmSjh496lgH15wpUQcOHNCaNWvcbjRKkgICAlS/fn21b99eb775pjw9PfXmm2+aHcth48aNysnJUa1atRyfQwcOHNCYMWNUp04ds+Od13XXXaeIiAi3+SyKiIiQp6dnlfkcksr/7NPT093us+jnn3/W//t//08zZ85U79691axZMw0fPlx33323nn32WbPjOWndurXS0tKUn5+vrKwsrVq1SseOHTPlc2j48OFasWKFPv/8c9WoUcOxPCYmRiUlJcrPz3fa3l0/eyhSbsxutzvOaXUHXbt21Y4dO5SWluZ4tGnTRgMHDlRaWpqsVqvZEc9y8uRJ7du3T7GxsWZHcXLTTTedNd3n999/r9q1a5uU6PzmzZunqKgo9erVy+wo53T69Gl5eDj/KrNarY5/cXU3AQEBio2N1fHjx7V69Wr16dPH7EjnVLduXcXExGjt2rWOZYWFhdq8ebNbXrvp7s6UqL179+rTTz9VeHi42ZEqxN0+hwYNGqTt27c7fQ7FxcVp7NixWr16tdnxzuvQoUM6duyY23wWeXt768Ybb6wyn0NS+ZkRrVu3drtr9kpLS1VaWlqlPodCQkIUGRmpvXv36uuvv67UzyHDMDR8+HAtXbpUn332merWreu0vnXr1vLy8nL67ElPT1dmZqZbfvZwap+bGDdunHr27KlatWrpxIkTeu+997Ru3Tq3+mAICgpynMN6RkBAgMLDw89abpbHHntMvXv3Vu3atXXkyBFNmDBBVqtV9957r9nRnIwaNUodOnTQlClTNGDAAKWmpuq1117Ta6+9ZnY0J3a7XfPmzdPgwYPl6emevy569+6tf/3rX6pVq5auv/56ffPNN5o5c6aGDh1qdjQnq1evlmEYatSokX744QeNHTtWjRs31gMPPGBappMnTzr9C3lGRobS0tIUFhamWrVqaeTIkfrnP/+pBg0aqG7duho/frzi4uLUt29ft8mYl5enzMxMxz2ZzvzFMCYmplL/9fJCOWNjY3XnnXdq27ZtWrFihWw2m+Nc/7CwMHl7e5ueMTw8XP/617/05z//WbGxsfrpp5/0yiuv6PDhw5V+y4OL/Zn/voR6eXkpJiZGjRo1couMYWFhmjRpkvr376+YmBjt27dPjz/+uOrXr6/ExES3yFirVi2NHTtWd999tzp37qxbbrlFq1at0vLly7Vu3bpKy1iRnFL5P+IsWbJEzz33XKVmq2jGm2++WWPHjpWfn59q166t9evX65133tHMmTPdKueSJUsUGRmpWrVqaceOHRoxYoT69u1bqZNiJCcn67333tOHH36ooKAgx+/CkJAQ+fn5KSQkRA8++KBGjx6tsLAwBQcH65FHHlFCQoLat29faTkrzNQ5A+EwdOhQo3bt2oa3t7cRGRlpdO3a1fjkk0/MjnVR7jb9+d13323ExsYa3t7eRvXq1Y27777b+OGHH8yOdU7Lly83brjhBsPHx8do3Lix8dprr5kd6SyrV682JBnp6elmRzmvwsJCY8SIEUatWrUMX19f47rrrjP+8Y9/GMXFxWZHc/L+++8b1113neHt7W3ExMQYycnJRn5+vqmZPv/8c0PSWY/BgwcbhlE+Bfr48eON6Ohow8fHx+jatWul/yxcLOO8efPOuX7ChAluk/PM1Oznenz++edukfHnn3827rjjDiMuLs7w9vY2YmNjjT//+c9GampqpeWrSM5zMWP68wtlPH36tNG9e3cjMjLS8PLyMmrXrm089NBDRnZ2tttkPOPNN9806tevb/j6+hrNmzc3li1bVqkZK5rz1VdfNfz8/Ez7nXmxjFlZWcaQIUOMuLg4w9fX12jUqJHx3HPPVfrtQi6W88UXXzRq1KhheHl5GbVq1TKeeuqpSv+sPN/vwnnz5jm2+fnnn41hw4YZ1apVM/z9/Y077rjDyMrKqtScFWUxDJOuMAMAAACAKoprpAAAAADARRQpAAAAAHARRQoAAAAAXESRAgAAAAAXUaQAAAAAwEUUKQAAAABwEUUKAAAAAFxEkQIAuIX9+/fLYrEoLS3N7CgO3333ndq3by9fX1+1aNHinNt06dJFI0eOvKI5Jk6ceN79AwDMQZECAEiShgwZIovFomnTpjktX7ZsmSwWi0mpzDVhwgQFBAQoPT1da9euNTvOJZs/f75CQ0PNjgEAVxWKFADAwdfXV9OnT9fx48fNjnLZlJSUXPJr9+3bp44dO6p27doKDw+/jKkAAFUdRQoA4NCtWzfFxMRo6tSp593mXKeZvfDCC6pTp47j+ZAhQ9S3b19NmTJF0dHRCg0N1TPPPKOysjKNHTtWYWFhqlGjhubNm3fW+3/33Xfq0KGDfH19dcMNN2j9+vVO63fu3KmePXsqMDBQ0dHRGjRokH766SfH+i5dumj48OEaOXKkIiIilJiYeM7jsNvteuaZZ1SjRg35+PioRYsWWrVqlWO9xWLR1q1b9cwzz8hisWjixInn/Z6UlZVp+PDhCgkJUUREhMaPHy/DMJzea9myZU6vCQ0N1fz58x3PDx06pHvvvVdhYWEKCAhQmzZttHnz5nPub9++fbruuus0fPhwGYah4uJiPfbYY6pevboCAgLUrl07rVu3TpK0bt06PfDAAyooKJDFYnE6ltmzZ6tBgwby9fVVdHS07rzzzvMeIwDAGUUKAOBgtVo1ZcoUvfTSSzp06NAfeq/PPvtMR44c0YYNGzRz5kxNmDBBf/rTn1StWjVt3rxZf//73/W3v/3trP2MHTtWY8aM0TfffKOEhAT17t1bx44dkyTl5+fr1ltvVcuWLfX1119r1apVOnr0qAYMGOD0Hm+//ba8vb21adMmzZ0795z5XnzxRT333HN69tlntX37diUmJurPf/6z9u7dK0nKysrS9ddfrzFjxigrK0uPPfbYeY/17bfflqenp1JTU/Xiiy9q5syZeuONNyr8vTp58qRuvvlmHT58WB999JG+/fZbPf7447Lb7Wdtu337dnXs2FH33XefXn75ZVksFg0fPlwpKSlatGiRtm/frrvuuks9evTQ3r171aFDB73wwgsKDg5WVlaW41i+/vprPfroo3rmmWeUnp6uVatWqXPnzhXODADXPAMAAMMwBg8ebPTp08cwDMNo3769MXToUMMwDGPp0qXGbz8uJkyYYDRv3tzptc8//7xRu3Ztp/eqXbu2YbPZHMsaNWpkdOrUyfG8rKzMCAgIMBYuXGgYhmFkZGQYkoxp06Y5tiktLTVq1KhhTJ8+3TAMw5g8ebLRvXt3p30fPHjQkGSkp6cbhmEYN998s9GyZcuLHm9cXJzxr3/9y2nZjTfeaAwbNszxvHnz5saECRMu+D4333yz0aRJE8NutzuWPfHEE0aTJk0czyUZS5cudXpdSEiIMW/ePMMwDOPVV181goKCjGPHjp1zH2e+55s2bTKqVatmPPvss451Bw4cMKxWq3H48GGn13Tt2tUYN26cYRiGMW/ePCMkJMRp/X//+18jODjYKCwsvODxAQDOjREpAMBZpk+frrffflt79uy55Pe4/vrr5eHx68dMdHS0mjZt6nhutVoVHh6unJwcp9clJCQ4vvb09FSbNm0cOb799lt9/vnnCgwMdDwaN24sqfx0tzNat259wWyFhYU6cuSIbrrpJqflN9100yUdc/v27Z0m5EhISNDevXtls9kq9Pq0tDS1bNlSYWFh590mMzNTt912m55++mmNGTPGsXzHjh2y2Wxq2LCh0/dl/fr1Tt+T37vttttUu3ZtXXfddRo0aJDeffddnT59ukJ5AQCSp9kBAADup3PnzkpMTNS4ceM0ZMgQp3UeHh5O1/9IUmlp6Vnv4eXl5fTcYrGcc9m5Tl87n5MnT6p3796aPn36WetiY2MdXwcEBFT4PSuDxWK54PfMz8/vou8RGRmpuLg4LVy4UEOHDlVwcLCk8u+J1WrV1q1bZbVanV4TGBh43vcLCgrStm3btG7dOn3yySd6+umnNXHiRG3ZsoUZ/gCgAhiRAgCc07Rp07R8+XKlpKQ4LY+MjFR2drZTMbic93766quvHF+XlZVp69atatKkiSSpVatW2rVrl+rUqaP69es7PVwpT8HBwYqLi9OmTZuclm/atEnx8fEuZ/79pBBfffWVGjRo4Cg2kZGRysrKcqzfu3ev0+hPs2bNlJaWpry8vPPuw8/PTytWrJCvr68SExN14sQJSVLLli1ls9mUk5Nz1vckJiZGkuTt7X3O0TFPT09169ZNM2bM0Pbt27V//3599tlnLh8/AFyLKFIAgHNq2rSpBg4cqFmzZjkt79Kli3JzczVjxgzt27dPr7zyilauXHnZ9vvKK69o6dKl+u6775ScnKzjx49r6NChkqTk5GTl5eXp3nvv1ZYtW7Rv3z6tXr1aDzzwQIVPoztj7Nixmj59ut5//32lp6frySefVFpamkaMGOFy5szMTI0ePVrp6elauHChXnrpJaf3ufXWW/Xyyy/rm2++0ddff62///3vTqNz9957r2JiYtS3b19t2rRJP/74o/773/+eVWIDAgL08ccfy9PTUz179tTJkyfVsGFDDRw4UPfff78++OADZWRkKDU1VVOnTtXHH38sSapTp45OnjyptWvX6qefftLp06e1YsUKzZo1S2lpaTpw4IDeeecd2e12NWrUyOXjB4BrEUUKAHBezzzzzFmn3jVp0kSzZ8/WK6+8oubNmys1NfWCM9q5atq0aZo2bZqaN2+uL774Qh999JEiIiIkyTGKZLPZ1L17dzVt2lQjR45UaGio0/VYFfHoo49q9OjRGjNmjJo2bapVq1bpo48+UoMGDVzOfP/99+vnn39W27ZtlZycrBEjRujhhx92rH/uuedUs2ZNderUSffdd58ee+wx+fv7O9Z7e3vrk08+UVRUlG6//XY1bdpU06ZNO+tUPan8dL2VK1fKMAz16tVLp06d0rx583T//fdrzJgxatSokfr27astW7aoVq1akqQOHTro73//u+6++25FRkZqxowZCg0N1QcffKBbb71VTZo00dy5c7Vw4UJdf/31Lh8/AFyLLMbvT9oGAAAAAFwQI1IAAAAA4CKKFAAAAAC4iCIFAAAAAC6iSAEAAACAiyhSAAAAAOAiihQAAAAAuIgiBQAAAAAuokgBAAAAgIsoUgAAAADgIooUAAAAALiIIgUAAAAALqJIAQAAAICL/j8doHu+bLXwBQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# compare my implementation results with the KMeans results\n",
    "# it is obvious that the more buckets we have, the smaller the MSE tends to be. So we need to compare the MSE\n",
    "# at the same number of buckets (or clusters).\n",
    "from sklearn.cluster import KMeans\n",
    "fico = borrowers.fico_score.values.reshape(-1, 1)\n",
    "mse_kmeans = []\n",
    "mse_dp = []\n",
    "mse_dp_int = []\n",
    "for i in range(3, 21):\n",
    "    kmeans = KMeans(n_clusters=i, n_init='auto', max_iter=500, tol=0.00001, verbose=0, \n",
    "                    random_state=100, copy_x=True, algorithm='lloyd')\n",
    "    kmeans.fit(fico)\n",
    "    mse_kmeans.append(kmeans.inertia_ / len(fico))\n",
    "    centroids, mse, interval, counter = dp_lloyd_quantization(data=borrowers.fico_score.values, \n",
    "                                                                  n=i, tol=0.00001, max_iter=500)\n",
    "    mse_dp.append(mse)\n",
    "    centroids, mse, interval, counter = dp_lloyd_quantization(data=borrowers.fico_score.values, \n",
    "                                                                n=i, tol=0.00001, max_iter=500, \n",
    "                                                                  is_integer=True)\n",
    "    mse_dp_int.append(mse)\n",
    "df = pd.DataFrame({'n_buckets': range(3,21), 'mse_kmeans': mse_kmeans, 'mse_dp': mse_dp, \n",
    "                   'mse_dp_int': mse_dp_int})\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.plot(df['n_buckets'], df['mse_kmeans'], 'o', label='mse_kmeans')\n",
    "plt.plot(df['n_buckets'], df['mse_dp'], '^', label='mse_dp')\n",
    "plt.plot(df['n_buckets'], df['mse_dp_int'], '+', label='mse_dp_int')\n",
    "plt.xlabel('Number of buckets')\n",
    "plt.ylabel('MSE')\n",
    "plt.xticks(df['n_buckets'])\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa3d52d7",
   "metadata": {},
   "source": [
    "# Comparison with KMeans\n",
    "\n",
    "From the comparison above, we can see that my implementation with representative values of each bucket being limited to integers have very close MSE to the KMeans results. We can use the implementation to quantize the FICO scores."
   ]
  }
 ],
 "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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}