2020101089_Q4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "improving-pepper",
   "metadata": {},
   "source": [
    "# Assignment 2 - Question 4\n",
    "The objective of this assignment is to get you familiarize with  the  problem  of  `Linear Regression`.\n",
    "\n",
    "## Instructions\n",
    "- Write your code and analysis in the indicated cells.\n",
    "- Ensure that this notebook runs without errors when the cells are run in sequence.\n",
    "- Do not attempt to change the contents of other cells.\n",
    "- No inbuilt functions to be used until specified\n",
    "\n",
    "## Submission\n",
    "- Ensure that this notebook runs without errors when the cells are run in sequence.\n",
    "- Rename the notebook to `<roll_number>_Q4.ipynb`.\n",
    "- Fill the Name and Roll number in the below markdown"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f6269f1",
   "metadata": {},
   "source": [
    "Name: Siddik Ayyappa <br>\n",
    "Roll Number: 2020101089"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "final-transaction",
   "metadata": {},
   "source": [
    "## Background about the dataset\n",
    "\n",
    "TLDR: You have 4 independent variables (`float`) for each molecule. You can use a linear combination of these 4 independent variables to predict the bandgap (dependent variable) of each molecule.\n",
    "\n",
    "You can read more about the problem in [Li et al, Bandgap tuning strategy by cations and halide ions of lead halide perovskites learned from machine learning, RSC Adv., 2021,11, 15688-15694](https://doi.org/10.1039/D1RA03117A)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "lyric-olympus",
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "import random\n",
    "import numpy as np\n",
    "import pandas as pd # for visualisation\n",
    "import timeit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "hundred-receipt",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_molecules = list()\n",
    "\n",
    "with open('bg_data.txt', 'r') as infile:\n",
    "    input_rows = csv.DictReader(infile)\n",
    "    \n",
    "    for row in input_rows:\n",
    "        current_mol = ([float(row['Cs']), float(row['FA']), float(row['Cl']), float(row['Br'])], float(row['Bandgap']))\n",
    "        all_molecules.append(current_mol)\n",
    "\n",
    "random.shuffle(all_molecules)\n",
    "\n",
    "\n",
    "num_train = int(len(all_molecules) * 0.8)\n",
    "\n",
    "# each point in x_train has 4 values - 1 for each feature\n",
    "x_train = np.array([x[0] for x in all_molecules[:num_train]])\n",
    "# each point in y_train has 1 value - the bandgap of the molecule\n",
    "y_train = np.array([x[1] for x in all_molecules[:num_train]])\n",
    "\n",
    "x_test = np.array([x[0] for x in all_molecules[num_train:]])\n",
    "y_test = np.array([x[1] for x in all_molecules[num_train:]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "square-direction",
   "metadata": {},
   "source": [
    "### 4.1 Implement a Linear Regression model that minimizes the MSE **without using any libraries**. You may use NumPy to vectorize your code, but *do not use numpy.polyfit* or anything similar.\n",
    "\n",
    "4.1.1 Explain how you plan to implement Linear Regression in 5-10 lines.\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "frozen-forth",
   "metadata": {},
   "source": [
    "<!-- your answer to 1.1.1 -->\n",
    "We shall first look at the problem we are trying to solve. We have a four dimensional input (```x_train```), which is being mapped to a real number (```y```). We, plan to fit the data as follows. \n",
    " \n",
    " $\\theta.X + b = y $\n",
    "\n",
    "\n",
    " Where the paramter $\\theta$, is a 4 dimensional vector and $b$ is  The parameters we aim to predict are $\\theta, b$, which will result in the least MSE (Or rather the objective loss function). We shall achieve this with te help of Gradient Descent. <br>\n",
    " $\\theta^{i+1} = \\theta^i - \\eta*(\\nabla_{\\theta}J)$<br>\n",
    " $b^{i+1} = b^i - \\eta*(\\nabla_{b}J)$<br>\n",
    " $\\hat{y_i} = \\theta.X_i + b$<br>\n",
    " $J = \\frac{1}{M}\\sum_{i = 1}^m(y_i - \\hat{y_i})^2$<br>\n",
    "\n",
    " This is done for multiple number of epochs, and then the parameters are frozen, and tested on. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "addressed-winter",
   "metadata": {},
   "source": [
    "4.1.2 Implement Linear Regression using `x_train` and `y_train` as the train dataset.\n",
    "\n",
    "4.1.2.1 Choose the best learning rate and print the learning rate for which you achieved the best MSE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "angry-depression",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training Done:\n",
      " Time: 0.0840987409974332 s\n"
     ]
    }
   ],
   "source": [
    "class linear_regression():\n",
    "    def __init__(self, learning_rate):\n",
    "        self.parameters = None\n",
    "        self.learning_rate = learning_rate\n",
    "    def fit(self, x_train, y_train, epochs):\n",
    "        self.parameters, gradient = np.zeros(x_train[0].shape), np.zeros(x_train[0].shape)\n",
    "        bias_gradient = 0\n",
    "        self.bias = 0\n",
    "        for _ in range(epochs):\n",
    "            loss = 0\n",
    "            diff = 0\n",
    "            for i in range(len(x_train)):\n",
    "                y_temp = np.dot(self.parameters, x_train[i]) + self.bias\n",
    "                diff = y_temp - y_train[i]\n",
    "                loss += diff**2\n",
    "                gradient += 2*x_train[i]*diff\n",
    "                bias_gradient += diff\n",
    "            gradient /= len(x_train)\n",
    "            bias_gradient  /= len(x_train)\n",
    "            loss /= len(x_train)\n",
    "            self.parameters -= self.learning_rate*(gradient)\n",
    "            self.bias -= self.learning_rate*(bias_gradient)    \n",
    "    def predict(self, x):\n",
    "        y_pred = []\n",
    "        for i in range(len(x)):\n",
    "            y_pred += [np.dot(self.parameters, x[i]) + self.bias]\n",
    "        y_pred = np.array(y_pred)\n",
    "        return y_pred\n",
    "x = linear_regression(0.3)\n",
    "start = timeit.default_timer()\n",
    "x.fit(x_train, y_train, 100)\n",
    "end = timeit.default_timer()\n",
    "print(\"Training Done:\\n Time:\", end-start, \"s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "585f914c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean Square Error: \n",
      "0.08620944847648604\n"
     ]
    }
   ],
   "source": [
    "print(\"Mean Square Error: \")\n",
    "print(sum((x.predict(x_test) - y_test)**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "angry-tampa",
   "metadata": {},
   "source": [
    "4.1.3 Make a [Parity Plot](https://en.wikipedia.org/wiki/Parity_plot) of your model's bandgap predictions on the test set with the actual values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "foster-center",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcoAAAHZCAYAAAAG3nvpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABiQklEQVR4nO3dd1QU198G8GdBWJCyCIogEMCGBTtqEBUb1qDGJCZ2bD+jGFs0sSRBNIq9JLEkatRYorHF2LBi7wrGChYQo6jEAgKCwN73D142rrRd2AY8n3M4x52d8h1Yebgzd+6VCCEEiIiIKFdG+i6AiIjIkDEoiYiI8sGgJCIiygeDkoiIKB8MSiIionwwKImIiPLBoCQiIsoHg5KIiCgfDEoiIqJ8MChLMIlEgqlTp+q7DIWAgAC4ubnpuwwlSUlJGDJkCBwcHCCRSDBmzBh9l1TitGrVCq1atdJ3GVQER48ehUQiwdGjR/Vdil4wKHVszZo1kEgkii8zMzNUr14dI0eOxJMnT7R67NOnT2Pq1Kl4+fKlRvfbqlUrpXOytbVF48aN8euvv0Iul2vkGDNnzsSff/6pkX29u981a9Zg+PDhWLduHfr165fv+pmZmVi9ejVatWoFW1tbSKVSuLm5YeDAgbh48aLG6wOAR48eYerUqYiIiNDK/otCIpFg5MiR+i5Da2JiYpQ+20ZGRrC1tUWnTp1w5swZfZdHOlJG3wWUVtOmTYO7uztSU1Nx8uRJLFu2DHv37sW1a9dQtmxZjRzj9evXKFPmvx/x6dOnERwcjICAANjY2GjkGNmcnZ0REhICAIiPj8dvv/2GwYMHIyoqCrNmzSry/mfOnImPP/4Y3bt3L/K+3nbkyBG8//77CAoKKnDd169fo0ePHggNDUXLli0xefJk2NraIiYmBn/88QfWrl2L2NhYODs7a7TGR48eITg4GG5ubqhfv75G960LBw4c0HcJRdarVy907twZmZmZiIqKwtKlS9G6dWtcuHABderU0Xd5WteyZUu8fv0apqam+i5FLxiUetKpUyd4eXkBAIYMGQI7OzssWLAAO3fuRK9evQq9X7lcjjdv3sDMzAxmZmaaKrdAMpkMffv2VbweNmwYPDw88NNPP2H69OkwMTHRWS3qePr0KWrVqqXSuhMmTEBoaCgWLlyY4xJtUFAQFi5cqIUKiz9D/+WanJwMCwuLfNdp2LCh0ue7RYsW6NSpE5YtW4alS5dqu0QlqtSraUZGRjr9fWJoeOnVQLRp0wYAEB0dDQCYN28emjVrBjs7O5ibm6NRo0bYunVrju2yL31t2LABtWvXhlQqRWhoqOK97HuUU6dOxYQJEwAA7u7uiktJMTEx8PX1Rb169XKty8PDAx06dFD7fMqWLYv3338fycnJiI+Pz3O95ORkfPnll3BxcYFUKoWHhwfmzZuHtye1kUgkSE5Oxtq1axV1BwQE5Hv8p0+fYvDgwahYsSLMzMxQr149rF27VvF+9j2X6Oho7NmzR+n7kZt//vkHP//8M/z8/HK9j2lsbIzx48crWpN53Y+dOnUqJBKJ0rKDBw+iefPmsLGxgaWlJTw8PDB58mRFnY0bNwYADBw4UFHnmjVrFNtv2bIFjRo1grm5OcqXL4++ffvi4cOHSscICAiApaUlYmNj8cEHH8DS0hJOTk5YsmQJAODq1ato06YNLCws4Orqio0bN+b7/VXHu/cos7/3f/zxB2bMmAFnZ2eYmZmhbdu2uHPnTo7tz507h44dO0Imk6Fs2bLw9fXFqVOnlNa5f/8+RowYAQ8PD5ibm8POzg6ffPJJjp9n9q2PY8eOYcSIEbC3ty/UFYAWLVoAAO7evau0/OXLlxgzZozi81y1alXMnj07xy2IZ8+eoV+/frC2toaNjQ0GDBiAK1eu5PjZZv/c7t69i86dO8PKygp9+vQBkPVH8aJFi1C7dm2YmZmhYsWKGDZsGF68eKF0rIsXL6JDhw4oX748zM3N4e7ujkGDBimts2nTJjRq1AhWVlawtrZGnTp1sHjxYsX7ed2jVOez9/DhQ3Tv3h2WlpaoUKECxo8fj8zMTNW/6XrEFqWByP4PZ2dnBwBYvHgxunbtij59+uDNmzfYtGkTPvnkE+zevRtdunRR2vbIkSP4448/MHLkSJQvXz7XX9A9evRAVFQUfv/9dyxcuBDly5cHAFSoUAH9+vXD0KFDce3aNXh6eiq2uXDhAqKiovDNN98U6pzu3bsHY2PjPC/zCiHQtWtXhIWFYfDgwahfvz7279+PCRMm4OHDh4oW2rp16zBkyBA0adIE//vf/wAAVapUyfO4r1+/RqtWrXDnzh2MHDkS7u7u2LJlCwICAvDy5UuMHj0aNWvWxLp16zB27Fg4Ozvjyy+/VHw/crNv3z5kZGQUeA9TXdevX8cHH3yAunXrYtq0aZBKpbhz544iCGrWrIlp06bhu+++w//+9z/FL+hmzZoByPrFP3DgQDRu3BghISF48uQJFi9ejFOnTiE8PFzpe5+ZmYlOnTqhZcuWmDNnDjZs2ICRI0fCwsICU6ZMQZ8+fdCjRw8sX74c/fv3h7e3N9zd3TV6vm+bNWsWjIyMMH78eCQkJGDOnDno06cPzp07p1jnyJEj6NSpExo1aoSgoCAYGRlh9erVaNOmDU6cOIEmTZoAyPqsnj59Gp999hmcnZ0RExODZcuWoVWrVrhx40aO2xkjRoxAhQoV8N133yE5OVnt2rMDuFy5coplKSkp8PX1xcOHDzFs2DC89957OH36NCZNmoS4uDgsWrQIQFbA+fv74/z58xg+fDhq1KiBnTt3YsCAAbkeKyMjAx06dEDz5s0xb948xbkMGzZM8fMfNWoUoqOj8dNPPyE8PBynTp2CiYkJnj59ivbt26NChQqYOHEibGxsEBMTg+3btyv2f/DgQfTq1Qtt27bF7NmzAQA3b97EqVOnMHr06Dy/B+p+9jp06ICmTZti3rx5OHToEObPn48qVapg+PDhan//dU6QTq1evVoAEIcOHRLx8fHiwYMHYtOmTcLOzk6Ym5uLf/75RwghREpKitJ2b968EZ6enqJNmzZKywEIIyMjcf369RzHAiCCgoIUr+fOnSsAiOjoaKX1Xr58KczMzMTXX3+ttHzUqFHCwsJCJCUl5XtOvr6+okaNGiI+Pl7Ex8eLmzdvilGjRgkAwt/fX7HegAEDhKurq+L1n3/+KQCI77//Xml/H3/8sZBIJOLOnTuKZRYWFmLAgAH51pFt0aJFAoBYv369YtmbN2+Et7e3sLS0FImJiYrlrq6uokuXLgXuc+zYsQKACA8PV6mGd881W1BQkHj7v93ChQsFABEfH5/nvi5cuCAAiNWrVystf/PmjbC3txeenp7i9evXiuW7d+8WAMR3332nVA8AMXPmTMWyFy9eCHNzcyGRSMSmTZsUy2/dupXjs5MXACIwMDDfdXx9fYWvr6/idVhYmAAgatasKdLS0hTLFy9eLACIq1evCiGEkMvlolq1aqJDhw5CLpcr1ktJSRHu7u7Cz89Padm7zpw5IwCI3377TbEs+/9f8+bNRUZGRoHnFx0dLQCI4OBgER8fLx4/fixOnDghGjduLACILVu2KNadPn26sLCwEFFRUUr7mDhxojA2NhaxsbFCCCG2bdsmAIhFixYp1snMzBRt2rTJ8XPO/rlNnDhRaZ8nTpwQAMSGDRuUloeGhiot37FjhwAgLly4kOc5jh49WlhbW+f7/cj+mYWFhQkhCvfZmzZtmtI+GzRoIBo1apTnMQ0JL73qSbt27VChQgW4uLjgs88+g6WlJXbs2AEnJycAgLm5uWLdFy9eICEhAS1atMDly5dz7MvX11fl+2y5kclk6NatG37//XfFJc/MzExs3rwZ3bt3V+l+yK1bt1ChQgVUqFABNWvWxI8//oguXbrg119/zXObvXv3wtjYGKNGjVJa/uWXX0IIgX379hXqfPbu3QsHBwele70mJiYYNWoUkpKScOzYMbX3mZiYCACwsrIqVE15yf6re+fOnWr3EL548SKePn2KESNGKN0/6tKlC2rUqIE9e/bk2GbIkCFKx/bw8ICFhQV69uypWO7h4QEbGxvcu3dPzbNRz8CBA5XuX2a3lrOPGxERgdu3b6N379549uwZ/v33X/z7779ITk5G27Ztcfz4ccX37O3/L+np6Xj27BmqVq0KGxubXP/PDB06FMbGxirXGhQUhAoVKsDBwQEtWrTAzZs3MX/+fHz88ceKdbZs2YIWLVqgXLlyilr//fdftGvXDpmZmTh+/DgAIDQ0FCYmJhg6dKhiWyMjIwQGBuZ5/HdbXVu2bIFMJoOfn5/SsRo1agRLS0uEhYUB+O/ztXv3bqSnp+e6bxsbGyQnJ+PgwYMqfz8K89n7/PPPlV63aNFC658xTeGlVz1ZsmQJqlevjjJlyqBixYrw8PCAkdF/f7fs3r0b33//PSIiIpCWlqZY/u79LQAauTzWv39/bN68GSdOnEDLli1x6NAhPHnyROVLjW5ublixYoXikZdq1arB3t4+323u37+PSpUq5QifmjVrKt4vjPv376NatWpK38+i7tfa2hoA8OrVq0LVlJdPP/0UK1euxJAhQzBx4kS0bdsWPXr0wMcff5yj/ndln4eHh0eO92rUqIGTJ08qLTMzM8txaVkmk8HZ2TnH50omk+W416Vp7733ntLr7MuY2ce9ffs2AOR5SRIAEhISUK5cObx+/RohISFYvXo1Hj58qHSPOyEhIcd26v6f+d///odPPvkEqampOHLkCH744Ycc99du376Nv//+O8/L90+fPgWQ9XNzdHTMcTm4atWquW5XpkyZHPdRb9++jYSEhDz/j2Ufy9fXFx999BGCg4OxcOFCtGrVCt27d0fv3r0hlUoBZF2G/uOPP9CpUyc4OTmhffv26NmzJzp27Jjn90MTn71y5cpp/TOmKQxKPWnSpImi1+u7Tpw4ga5du6Jly5ZYunQpHB0dYWJigtWrV+fayeLtv6YLq0OHDqhYsSLWr1+Pli1bYv369XBwcEC7du1U2t7CwkLldYujGjVqAMjq9KLKIxq5/UEDIMcvV3Nzcxw/fhxhYWHYs2cPQkNDsXnzZrRp0wYHDhxQq9VTkLz2ldfyt8NGGwo6bnZrce7cuXl+zy0tLQEAX3zxBVavXo0xY8bA29sbMpkMEokEn332Wa4tdXX/z1SrVk3x+f7ggw9gbGyMiRMnonXr1or/x3K5HH5+fvjqq69y3Uf16tXVOmY2qVSa448muVwOe3t7bNiwIddtskNJIpFg69atOHv2LHbt2oX9+/dj0KBBmD9/Ps6ePQtLS0vY29sjIiIC+/fvx759+7Bv3z6sXr0a/fv3V+oAVxSa/BzrA4PSAG3btg1mZmbYv3+/4q8+AFi9enWR9pvXL28g64Pcu3dvrFmzBrNnz8aff/6p9uUpdbm6uuLQoUN49eqVUqvy1q1bivdVqT23/f7999+Qy+VKv2By26+qOnXqBGNjY6xfv16lVna5cuVyHdght9askZER2rZti7Zt22LBggWYOXMmpkyZgrCwMLRr1y7Pc88+j8jISEWv6WyRkZGFOk9Dkt1hy9rausA/wrZu3YoBAwZg/vz5imWpqakaH1wj25QpU7BixQp88803il7mVapUQVJSUoG1urq6IiwsDCkpKUqtytx6/OalSpUqOHToEHx8fFQK/ffffx/vv/8+ZsyYgY0bN6JPnz7YtGmT4lK8qakp/P394e/vD7lcjhEjRuDnn3/Gt99+m2tLt6R/9t7Fe5QGyNjYGBKJRKn1ERMTU+SRabLvNeb1y6Nfv3548eIFhg0bhqSkJKXnxrQh+wHun376SWn5woULIZFI0KlTJ8UyCwsLlX/pde7cGY8fP8bmzZsVyzIyMvDjjz/C0tISvr6+atfq4uKCoUOH4sCBA/jxxx9zvC+XyzF//nz8888/ALJ+kSUkJODvv/9WrBMXF4cdO3Yobff8+fMc+8puPWVfcs/r5+bl5QV7e3ssX75c6fL8vn37cPPmzRy9o4ubRo0aoUqVKpg3bx6SkpJyvP/2Y0fGxsY5WsA//vij1h4/sLGxwbBhw7B//37FiEk9e/bEmTNnsH///hzrv3z5EhkZGQCyrt6kp6djxYoVivflcrniUR1V9OzZE5mZmZg+fXqO9zIyMhSflRcvXuT4vrz7+Xr27JnS+0ZGRqhbt67SOu8q6Z+9d7FFaYC6dOmCBQsWoGPHjujduzeePn2KJUuWoGrVqkq/eNXVqFEjAFl/DX/22WcwMTGBv7+/4hdxgwYN4OnpiS1btqBmzZpo2LChRs4nL/7+/mjdujWmTJmCmJgY1KtXDwcOHMDOnTsxZswYpUdAGjVqhEOHDmHBggWoVKkS3N3d0bRp01z3+7///Q8///wzAgICcOnSJbi5uWHr1q04deoUFi1aVOgOOfPnz8fdu3cxatQobN++HR988AHKlSuH2NhYbNmyBbdu3cJnn30GAPjss8/w9ddf48MPP8SoUaOQkpKCZcuWoXr16kqdS6ZNm4bjx4+jS5cucHV1xdOnT7F06VI4OzujefPmALJC18bGBsuXL4eVlRUsLCzQtGlTuLu7Y/bs2Rg4cCB8fX3Rq1cvRRd9Nzc3jB07tlDnqa6LFy/i+++/z7G8VatWinMoDCMjI6xcuRKdOnVC7dq1MXDgQDg5OeHhw4cICwuDtbU1du3aBSDrcui6desgk8lQq1YtnDlzBocOHVI8bqUNo0ePxqJFizBr1ixs2rQJEyZMwF9//YUPPvgAAQEBaNSoEZKTk3H16lVs3boVMTExKF++PLp3744mTZrgyy+/xJ07d1CjRg389ddfij+aVLl64uvri2HDhiEkJAQRERFo3749TExMcPv2bWzZsgWLFy/Gxx9/jLVr12Lp0qX48MMPUaVKFbx69QorVqyAtbU1OnfuDCCrg9fz58/Rpk0bODs74/79+/jxxx9Rv359xX39d5mYmBjEZ09n9NjjtlTK7p6eX3dtIYRYtWqVqFatmpBKpaJGjRpi9erVOR4tECL/7vnIpYv/9OnThZOTkzAyMsr1UZE5c+bkeIygIL6+vqJ27doFrpfbIxOvXr0SY8eOFZUqVRImJiaiWrVqYu7cuUqPAwiR9chCy5Ythbm5uQBQ4KMiT548EQMHDhTly5cXpqamok6dOjkerxBC9cdDsmVkZIiVK1eKFi1aCJlMJkxMTISrq6sYOHBgjkdHDhw4IDw9PYWpqanw8PAQ69evz/EzPHz4sOjWrZuoVKmSMDU1FZUqVRK9evXK8YjBzp07Ra1atUSZMmVyPEKwefNm0aBBAyGVSoWtra3o06eP4jGjbAMGDBAWFhY5zievn52q3xcAeX5Nnz5dcYzcHg95+9EKIf57FOPdn1N4eLjo0aOHsLOzE1KpVLi6uoqePXuKw4cPK9Z58eKF4udtaWkpOnToIG7duiVcXV2VPiuq/v97t6a5c+fm+n5AQIAwNjZWPMr06tUrMWnSJFG1alVhamoqypcvL5o1aybmzZsn3rx5o9guPj5e9O7dW1hZWQmZTCYCAgLEqVOnBAClR3Xy+rll++WXX0SjRo2Eubm5sLKyEnXq1BFfffWVePTokRBCiMuXL4tevXqJ9957T0ilUmFvby8++OADcfHiRcU+tm7dKtq3by/s7e2FqampeO+998SwYcNEXFycYp13Hw/JVpTPXm6/zwyVRAgt37GnYmXx4sUYO3YsYmJicvRKJCLt+fPPP/Hhhx/i5MmT8PHx0Xc59BYGJSkIIVCvXj3Y2dkpnsMiIs17/fq1UieczMxMtG/fHhcvXsTjx4810pOdNIf3KAnJycn466+/EBYWhqtXr2Lnzp36LomoRPviiy/w+vVreHt7Iy0tDdu3b8fp06cxc+ZMhqQBYouSEBMTA3d3d9jY2GDEiBGYMWOGvksiKtE2btyI+fPn486dO0hNTUXVqlUxfPjwEj23Z3HGoCQiIsoHn6MkIiLKB4OSiIgoH6WuM49cLsejR49gZWWl1rBoRERUsggh8OrVK1SqVCnfSQhKXVA+evQILi4u+i6DiIgMxIMHD3LM0PK2UheU2cOXPXjwQDF1EhERlQ5JSUn49ttvERUVBTMzMxw6dKjAYS1LXVBmX261trZmUBIRlSJJSUmYPn067t27B1tbW0ybNg3169cv8DZcqQtKIiIqfZKSkjBhwgTcunUL1tbWWLhwIcqXL6/Stuz1SkREJVpuIVm5cmWVt2dQEhFRiVXUkAQYlEREVEJpIiQBBiUREZVAmgpJgEFJREQljCZDEmBQEhFRCaLpkAQYlEREVEJoIyQBBiUREZUA2gpJgEFJRETFnDZDEmBQEhFRMabtkAQYlEREVEzpIiQBjvVKREQGJlMucD76OZ6+SoW9lRmauNvC2Eh54HJdhSTAoCQiIgMSei0OwbtuIC4hVbHMUWaGIP9a6OjpCEC3IQnw0isRERmI0GtxGL7+slJIAsDjhFQMX38ZodfidB6SAFuURERkADLlAsG7bkDk8p4AIAEw9a/rqP9gGyJ1GJIAW5RERGQAzkc/z9GSfJsA8DgxDREPk3UakgCDkoiIDMDTV3mH5NtMZeV1GpIAg5KIiAyAvZWZSuuNGTZQpyEJMCiJiMgANHG3haPMDJK8VhACFSzKoFszT12WBYBBSUREBsDYSIIg/1q5duaBEIBEgukf1s3xPKUu6DUoQ0JC0LhxY1hZWcHe3h7du3dHZGRkgdstWrQIHh4eMDc3h4uLC8aOHYvUVNWubxMRUTEj0X04vk2vQXns2DEEBgbi7NmzOHjwINLT09G+fXskJyfnuc3GjRsxceJEBAUF4ebNm1i1ahU2b96MyZMn67ByIiLSpOzHQ/IiARC86wYy5bm2ObVKr89RhoaGKr1es2YN7O3tcenSJbRs2TLXbU6fPg0fHx/07t0bAODm5oZevXrh3LlzWq+XiIi0Q5XHQ+ISUnE++jm8q9jprjAY2D3KhIQEAICtrW2e6zRr1gyXLl3C+fPnAQD37t3D3r170blz51zXT0tLQ2JiotIXEREZlgfxCSqtp+pjJJpkMCPzyOVyjBkzBj4+PvD0zLtXU+/evfHvv/+iefPmEEIgIyMDn3/+eZ6XXkNCQhAcHKytsomIqIiSkpKwdf2vgKV3geuq+hiJJhlMizIwMBDXrl3Dpk2b8l3v6NGjmDlzJpYuXYrLly9j+/bt2LNnD6ZPn57r+pMmTUJCQoLi68GDB9oon4iICiF77Nan10/DJD0pz/UkyBocvYl73lcctcUgWpQjR47E7t27cfz4cTg7O+e77rfffot+/fphyJAhAIA6deogOTkZ//vf/zBlyhQYGSlnv1QqhVQq1VrtRERUOG8PcC6ztkZAZw8EH3wIAEqPiWT3eQ3yr6WXx0P0GpRCCHzxxRfYsWMHjh49Cnd39wK3SUlJyRGGxsbGiv0REZHhy2sWEIeKFXNMs+Xw1jRbqsxVqWl6DcrAwEBs3LgRO3fuhJWVFR4/fgwAkMlkMDc3BwD0798fTk5OCAkJAQD4+/tjwYIFaNCgAZo2bYo7d+7g22+/hb+/vyIwiYjIcOU3VVZHT0f41XLINQxVmatSGyRCj80wSR4Pka5evRoBAQEAgFatWsHNzQ1r1qwBAGRkZGDGjBlYt24dHj58iAoVKsDf3x8zZsyAjY1NgcdMTEyETCZDQkICrK2tNXQmRESkisLOJ5k9V+W7gZWdIsv6NlQ7LFXNA70GpT4wKImI9KOwIZkpF2g++0iez1lKkHV59uTXbdS6DKtqHhhMr1ciIiq5ChuSgHqDEWgDg5KIiLSqKCEJqD7IgLYGI2BQEhGR1hQ1JAHVBxnQ1mAEDEoiItIKTYQkUPBcldoejIBBSUREGqepkAT+m6sSQI6w1MVgBAxKIiLSKE2GZLaOno5Y1rchHGTKl1cdZGaFejREHQYxhB0REZUM2gjJbPkNRqBNDEoiItIIbYZkNmMjSemej5KIiIonXYSkvjAoiYioSEpySAIMSiIiKoKSHpIAg5KIiAqpNIQkwKAkIqJCKC0hCTAoiYhITaUpJAEGJRERqaG0hSTAoCQiIhWVxpAEGJRERKSC0hqSAIOSiIgKUJpDEmBQEhFRPkp7SAIc65WIiPJQ2JDMlAudD1yuTQxKIiLKobAhGXotDsG7biAuIVWxzFFmhiD/WlqdCkubeOmViIiUFCUkh6+/rBSSAPA4IRXD119G6LU4bZWsVQxKIiJSKMrl1uBdNyByeS97WfCuG8iU57aGYWNQEhERgKJ13Dkf/TxHS/JtAkBcQirORz/XULW6w6AkIqIi9259+irvkCzMeoaEQUlEVMpp4hEQeyszja5nSBiURESlmKaek2zibgtHmRnyeghEgqzer03cbYtUrz4wKImISilNDiZgbCRBkH8tAMgRltmvg/xrFcvnKRmURESlkDZG3Ono6YhlfRvCQaZ8edVBZoZlfRsW2+coOeAAEVEpo81h6Tp6OsKvlgNH5iEiouJJF2O3GhtJ4F3FTqP71CdeeiUiKiU4wHnhMCiJiEoBhmThMSiJiEo4hmTRMCiJiEowhmTRMSiJiEoohqRmMCiJiEoghqTmMCiJiEoYhqRmMSiJiEoQhqTmMSiJiEoIhqR2MCiJiEoAhqT26DUoQ0JC0LhxY1hZWcHe3h7du3dHZGRkgdu9fPkSgYGBcHR0hFQqRfXq1bF3714dVExEZHgYktql17Fejx07hsDAQDRu3BgZGRmYPHky2rdvjxs3bsDCwiLXbd68eQM/Pz/Y29tj69atcHJywv3792FjY6Pb4omIDABDUvv0GpShoaFKr9esWQN7e3tcunQJLVu2zHWbX3/9Fc+fP8fp06dhYmICAHBzc9N2qUREBochqRsGdY8yISEBAGBrm/cM2H/99Re8vb0RGBiIihUrwtPTEzNnzkRmZmau66elpSExMVHpi4iouGNI6o7BBKVcLseYMWPg4+MDT0/PPNe7d+8etm7diszMTOzduxfffvst5s+fj++//z7X9UNCQiCTyRRfLi4u2joFIiKdYEjqlkQIIfRdBAAMHz4c+/btw8mTJ+Hs7JznetWrV0dqaiqio6NhbGwMAFiwYAHmzp2LuLi4HOunpaUhLS1N8ToxMREuLi5ISEiAtbW15k+EiEiLGJKak5iYCJlMVmAeGMTEzSNHjsTu3btx/PjxfEMSABwdHWFiYqIISQCoWbMmHj9+jDdv3sDU1FRpfalUCqlUqpW6iYh0iSGpH3q99CqEwMiRI7Fjxw4cOXIE7u7uBW7j4+ODO3fuQC6XK5ZFRUXB0dExR0gSEZUUDEn90WtQBgYGYv369di4cSOsrKzw+PFjPH78GK9fv1as079/f0yaNEnxevjw4Xj+/DlGjx6NqKgo7NmzBzNnzkRgYKA+ToGISOsYkvql10uvy5YtAwC0atVKafnq1asREBAAAIiNjYWR0X957uLigv3792Ps2LGoW7cunJycMHr0aHz99de6KpuISGcYkvpnMJ15dEXVm7dERPrGkNQuVfPAYB4PISKi/zAkDQeDkojIwDAkDQuDkojIgDAkDQ+DkojIQDAkDRODkojIADAkDReDkohIzxiSho1BSUSkRwxJw8egJCLSE4Zk8cCgJCLSA4Zk8cGgJCLSMYZk8cKgJCLSIYZk8cOgJCLSEYZk8cSgJCLSAYZk8cWgJCLSMoZk8cagJCLSIoZk8cegJCLSEoZkycCgJCLSAoZkyaGRoHz58qUmdkNEVCIwJEsWtYNy9uzZ2Lx5s+J1z549YWdnBycnJ1y5ckWjxRERFTcMyZJH7aBcvnw5XFxcAAAHDx7EwYMHsW/fPnTq1AkTJkzQeIFERMUFQ7JkKqPuBo8fP1YE5e7du9GzZ0+0b98ebm5uaNq0qcYLJCIqDhiSJZfaLcpy5crhwYMHAIDQ0FC0a9cOACCEQGZmpmarIyIqBhiSJZvaLcoePXqgd+/eqFatGp49e4ZOnToBAMLDw1G1alWNF0hEZMgYkiWf2kG5cOFCuLm54cGDB5gzZw4sLS0BAHFxcRgxYoTGCyQiMlQMydJBIoQQ+i5ClxITEyGTyZCQkABra2t9l0NExRRDsvhTNQ8K9RzlunXr0Lx5c1SqVAn3798HACxatAg7d+4sXLVERMUIQ7J0UTsoly1bhnHjxqFTp054+fKlogOPjY0NFi1apOn6iIgMCkOy9FE7KH/88UesWLECU6ZMgbGxsWK5l5cXrl69qtHiiIgMCUOydFI7KKOjo9GgQYMcy6VSKZKTkzVSFBGRoWFIll5qB6W7uzsiIiJyLA8NDUXNmjU1URMRkUFhSJZuaj8eMm7cOAQGBiI1NRVCCJw/fx6///47QkJCsHLlSm3USESkNwxJUjsohwwZAnNzc3zzzTdISUlB7969UalSJSxevBifffaZNmokItILhiQBRXyOMiUlBUlJSbC3t9dkTVrF5yiJSBUMyZJP1TxQu0X5trJly6Js2bJF2QURkcFhSNLb1A5Kd3d3SCSSPN+/d+9ekQoiItInhiS9S+2gHDNmjNLr9PR0hIeHIzQ0lPNRElGxxpCk3KgdlKNHj851+ZIlS3Dx4sUiF0REpA8MScpLocZ6zU2nTp2wbds2Te2OiEhnGJKUH40F5datW2Fra6up3RER6QRDkgqi9qXXBg0aKHXmEULg8ePHiI+Px9KlSzVaHBGRNjEkSRVqB2X37t2VXhsZGaFChQpo1aoVatSooam6iIi0iiFJqtLrxM0hISHYvn07bt26BXNzczRr1gyzZ8+Gh4eHSttv2rQJvXr1Qrdu3fDnn3+qtA0HHCAihiQBGh5wIDExUeUDqxM+x44dQ2BgIBo3boyMjAxMnjwZ7du3x40bN2BhYZHvtjExMRg/fjxatGih8vGIiBiSpC6VWpRGRkb5DjIAZN2rlEgkiomcCyM+Ph729vY4duwYWrZsmed6mZmZaNmyJQYNGoQTJ07g5cuXbFESUYEYkvQ2jbYow8LCNFZYfhISEgCgwN6z06ZNg729PQYPHowTJ07ku25aWhrS0tIUr9VpHRNRycGQpMJSKSh9fX21XQfkcjnGjBkDHx8feHp65rneyZMnsWrVqlznxMxNSEgIgoODNVQlERVHDEkqikIPip6SkoLY2Fi8efNGaXndunULtb/AwEBcu3YNJ0+ezHOdV69eoV+/flixYgXKly+v0n4nTZqEcePGKV4nJibCxcWlUDUSUfHDkKSiUjso4+PjMXDgQOzbty/X9wtzj3LkyJHYvXs3jh8/Dmdn5zzXu3v3LmJiYuDv769YJpfLAQBlypRBZGQkqlSporSNVCqFVCpVuyYiKv4YkqQJao/MM2bMGLx8+RLnzp2Dubk5QkNDsXbtWlSrVg1//fWXWvsSQmDkyJHYsWMHjhw5And393zXr1GjBq5evYqIiAjFV9euXdG6dWtERESwpUhECgxJ0hS1W5RHjhzBzp074eXlBSMjI7i6usLPzw/W1tYICQlBly5dVN5XYGAgNm7ciJ07d8LKygqPHz8GAMhkMpibmwMA+vfvDycnJ4SEhMDMzCzH/UsbGxsAyPe+JhGVLgxJ0iS1W5TJycmwt7cHAJQrVw7x8fEAgDp16uDy5ctq7WvZsmVISEhAq1at4OjoqPjavHmzYp3Y2FjExcWpWyYRlVIMSdI0tVuUHh4eiIyMhJubG+rVq4eff/4Zbm5uWL58ORwdHdXalyqDAh09ejTf99esWaPWMYmo5GJIkjYUaj7K7BZeUFAQOnbsiA0bNsDU1JShRUR6w5AkbSnyWK8pKSm4desW3nvvPZUf2dAnjsxDVPIwJKkwVM0Dte9RvvucY9myZdGwYcNiEZJEVPIwJEnb1A7KNm3awN3dHZMnT8aNGze0URMRkUoYkqQLagflo0eP8OWXX+LYsWPw9PRE/fr1MXfuXPzzzz/aqI+IKFcMSdKVIt2jjI6OxsaNG/H777/j1q1baNmyJY4cOaLJ+jSO9yiJij+GJGmCqnlQ5M48mZmZ2LdvH7799lv8/fffRZpmSxcYlETFG0OSNEVrnXmynTp1CiNGjICjoyN69+4NT09P7Nmzp7C7IyIqEEOS9EHt5ygnTZqETZs24dGjR/Dz88PixYvRrVs3lC1bVhv1EREBYEiS/qgdlMePH8eECRPQs2dPPhJCRDrBkCR9UjsoT506pY06iIhyxZAkfSv0PUoiIm1jSJIhYFASkUFiSJKhYFASkcFhSJIhYVASkUFhSJKhUbszT7aLFy/i5s2bAICaNWvCy8tLY0URUenEkCRDpHZQ/vPPP+jVqxdOnToFGxsbAMDLly/RrFkzbNq0Cc7OzpqukYhKAYYkGSq1L70OGTIE6enpuHnzJp4/f47nz5/j5s2bkMvlGDJkiDZqJKISjiFJhkztsV7Nzc1x+vRpNGjQQGn5pUuX0KJFC6SkpGi0QE3jWK9EhoUhSfqitbFeXVxckJ6enmN5ZmYmKlWqpO7uiKgUY0hScaB2UM6dOxdffPEFLl68qFh28eJFjB49GvPmzdNocURUcjEkqbhQ+9JruXLlkJKSgoyMDJQpk9UXKPvfFhYWSus+f/5cc5VqCC+9EukfQ5IMgap5oHav10WLFhWlLiIq5RiSVNyoHZQDBgzQRh1EVAowJKk4KvSAAwCQmpqKN2/eKC3j5Uwiyg1DkoortTvzJCcnY+TIkbC3t4eFhQXKlSun9EVE9C6GJBVnagflV199hSNHjmDZsmWQSqVYuXIlgoODUalSJfz222/aqJGIijGGJBV3al963bVrF3777Te0atUKAwcORIsWLVC1alW4urpiw4YN6NOnjzbqJKJiiCFJJYHaLcrnz58rPujW1taKR0CaN2+O48ePa7Y6Iiq2GJJUUqgdlJUrV0Z0dDQAoEaNGvjjjz8AZLU0swdJJ6LSjSFJJYnaQTlw4EBcuXIFADBx4kQsWbIEZmZmGDt2LCZMmKDxAomoeGFIUkmj9sg877p//z4uXbqEqlWrom7dupqqS2s4Mg+R9jAkqTjR2sg873J1dYWrq2tRd0NExRxDkkoqtYPyhx9+yHW5RCKBmZkZqlatipYtW8LY2LjIxRFR8cCQpJJM7aBcuHAh4uPjkZKSohhg4MWLFyhbtiwsLS3x9OlTVK5cGWFhYXBxcdF4wURkWBiSVNKp3Zln5syZaNy4MW7fvo1nz57h2bNniIqKQtOmTbF48WLExsbCwcEBY8eO1Ua9RGRAGJJUGqjdmadKlSrYtm0b6tevr7Q8PDwcH330Ee7du4fTp0/jo48+QlxcnCZr1Qh25iHSDIYkFXeq5oHaLcq4uDhkZGTkWJ6RkYHHjx8DACpVqoRXr16pu2siKiYYklSaqB2UrVu3xrBhwxAeHq5YFh4ejuHDh6NNmzYAgKtXr8Ld3V1zVRKRwWBIUmmjdlCuWrUKtra2aNSoEaRSKaRSKby8vGBra4tVq1YBACwtLTF//nyNF0tE+sWQpNKo0AMO3Lp1C1FRUQAADw8PeHh4aLQwbeE9SqLCYUhSSaO1e5TZatSoga5du6Jr166FDsmQkBA0btwYVlZWsLe3R/fu3REZGZnvNitWrECLFi0U81+2a9cO58+fL9TxiUg1DEkqzVR6jnLcuHEq73DBggUqr3vs2DEEBgaicePGyMjIwOTJk9G+fXvcuHEDFhYWuW5z9OhR9OrVC82aNYOZmRlmz56N9u3b4/r163ByclL52ESkGoYklXYqXXpt3bq10uvLly8jIyND0ZKMioqCsbExGjVqhCNHjhS6mPj4eNjb2+PYsWNo2bKlSttkZmaiXLly+Omnn9C/f/8C1+elVyLVMSSpJNPoWK9hYWGKfy9YsABWVlZYu3at0sg82ZM4F0VCQgIAwNbWVuVtUlJSkJ6enuc2aWlpSEtLU7xOTEwsUo1EpQVDkiiL2p15nJyccODAAdSuXVtp+bVr19C+fXs8evSoUIXI5XJ07doVL1++xMmTJ1XebsSIEdi/fz+uX78OMzOzHO9PnToVwcHBOZazRUmUN4YklQZa68yTmJiI+Pj4HMvj4+OLNMhAYGAgrl27hk2bNqm8zaxZs7Bp0ybs2LEj15AEgEmTJiEhIUHx9eDBg0LXSFQaMCSJlKk9KPqHH36IgQMHYv78+WjSpAkA4Ny5c5gwYQJ69OhRqCJGjhyJ3bt34/jx43B2dlZpm3nz5mHWrFk4dOhQvvNgZj/rSUQFY0gS5aR2UC5fvhzjx49H7969kZ6enrWTMmUwePBgzJ07V619CSHwxRdfYMeOHTh69KjKo/nMmTMHM2bMwP79++Hl5aXuKRBRLhiSRLkr9IADycnJuHv3LoCsgdLzepwjPyNGjMDGjRuxc+dOpWcxZTIZzM3NAQD9+/eHk5MTQkJCAACzZ8/Gd999h40bN8LHx0exjaWlJSwtLQs8Jnu9EuXEkKTSSNU8KHRQaoJEIsl1+erVqxEQEAAAaNWqFdzc3LBmzRoAgJubG+7fv59jm6CgIEydOrXAYzIoiZQxJKm00ujjIW9LTk7GrFmzcPjwYTx9+hRyuVzp/Xv37qm8L1Uy+ujRo0qvY2JiVN4/EeWPIUlUMLWDcsiQITh27Bj69esHR0fHPFuFRGTYGJJEqlE7KPft24c9e/Yo3R8kouKFIUmkOrWfoyxXrpxaI+cQkWFhSBKpR+2gnD59Or777jukpKRoox4i0iKGJJH61L70On/+fNy9excVK1aEm5sbTExMlN6/fPmyxoojIs1hSBIVjtpB2b17dy2UQUTaxJAkKjy9PkepD3yOkkobhiRR7rQ2KDoRFR8MSaKiU/vSa2ZmJhYuXIg//vgDsbGxePPmjdL7z58/11hxRFR4DEkizVC7RRkcHIwFCxbg008/RUJCAsaNG4cePXrAyMhIpSHkiEj7GJJEmqN2UG7YsAErVqzAl19+iTJlyqBXr15YuXIlvvvuO5w9e1YbNRKRGhiSRJqldlA+fvwYderUAZA1Y0dCQgIA4IMPPsCePXs0Wx0RqYUhSaR5agels7Mz4uLiAGRNr3XgwAEAwIULFzhBMpEeMSSJtEPtoPzwww9x+PBhAMAXX3yBb7/9FtWqVUP//v0xaNAgjRdIRAVjSBJpT5Gfozxz5gzOnDmDatWqwd/fX1N1aQ2fo6SShiFJVDham4/yXd7e3vD29i7qboioEBiSRNpXqKCMjIzEjz/+iJs3bwIAatasiS+++AIeHh4aLY6I8saQJNINte9Rbtu2DZ6enrh06RLq1auHevXq4fLly/D09MS2bdu0USMRvYMhSaQ7at+jrFKlCvr06YNp06YpLQ8KCsL69etx9+5djRaoabxHSYWVKRc4H/0cT1+lwt7KDE3cbWFsJNF5HQxJIs1QNQ/UDsqyZcvi77//RtWqVZWW3759G/Xq1TP4eSoZlFQYodfiELzrBuISUhXLHGVmCPKvhY6ejjqrgyFJpDlaGxS9VatWOHHiRI7lJ0+eRIsWLdTdHZHBC70Wh+HrLyuFJAA8TkjF8PWXEXotTid1MCSJ9EOlzjx//fWX4t9du3bF119/jUuXLuH9998HAJw9exZbtmxBcHCwdqok0pNMuUDwrhvI7bKLACABELzrBvxqOWj1MixDkkh/VLr0amSkWsNTIpEgMzOzyEVpEy+9kjrO3H2GXisKHsP496Hvw7uKnVZqYEgSaYdGn6OUy+UaK4yoOHn6KrXgldRYT10MSSL948TNRPmwtzLT6HrqYEgSGQYGJVE+mrjbwlFmhrzuPkqQ1fu1ibutRo/LkCQyHAxKonwYG0kQ5F8LAHKEZfbrIP9aGu3Iw5AkMiwMSqICdPR0xLK+DeEgU7686iAzw7K+DTX6HCVDksjwFHlQdKLSoKOnI/xqOWh1ZB6GJJFhUikoExMTVd4hH7mgksrYSMJHQIhKIZWC0sbGBhKJan85G/pzlESGhiFJZNhUCsqwsDDFv2NiYjBx4kQEBAQo5qE8c+YM1q5di5CQEO1USVRCMSSJDJ/ag6K3bdsWQ4YMQa9evZSWb9y4Eb/88guOHj2qyfo0jiPzkKFgSBLpl9YGRT9z5gy8vLxyLPfy8sL58+fV3R1RqcSQJCo+1A5KFxcXrFixIsfylStXwsXFRSNFEZVkDEmi4kXtx0MWLlyIjz76CPv27UPTpk0BAOfPn8ft27exbds2jRdIVJIwJImKH7VblJ07d0ZUVBT8/f3x/PlzPH/+HP7+/oiKikLnzp21USNRicCQJCqe1O7MU9yxMw/pA0OSyPBorTMPAJw4cQJ9+/ZFs2bN8PDhQwDAunXrcPLkycJVS1SC5ReSmXKBM3efYWfEQ5y5+wyZ8lL1dytRsaD2Pcpt27ahX79+6NOnDy5fvoy0tDQAQEJCAmbOnIm9e/dqvEii4iq/kAy9FofgXTcQl/DfXJaOMjME+dfS6PixRFQ0arcov//+eyxfvhwrVqyAiYmJYrmPjw8uX76s1r5CQkLQuHFjWFlZwd7eHt27d0dkZGSB223ZsgU1atSAmZkZ6tSpw3Amg1RQSA5ff1kpJAHgcUIqhq+/jNBrcfoomYhyoXZQRkZGomXLljmWy2QyvHz5Uq19HTt2DIGBgTh79iwOHjyI9PR0tG/fHsnJyXluc/r0afTq1QuDBw9GeHg4unfvju7du+PatWvqngqR1hR0uTV41w3kdpE1e1nwrhu8DEtkINQOSgcHB9y5cyfH8pMnT6rdOSE0NBQBAQGoXbs26tWrhzVr1iA2NhaXLl3Kc5vFixejY8eOmDBhAmrWrInp06ejYcOG+Omnn9Q9FSKtKKjjzvno5zlakm8TAOISUnE++rkOqiWigqgdlEOHDsXo0aNx7tw5SCQSPHr0CBs2bMD48eMxfPjwIhWTkJAAALC1zXu2+DNnzqBdu3ZKyzp06IAzZ84U6dhEmqBK79anr/IOycKsR0TapXZnnokTJ0Iul6Nt27ZISUlBy5YtIZVKMX78eHzxxReFLkQul2PMmDHw8fGBp6dnnus9fvwYFStWVFpWsWJFPH78ONf109LSFB2OAPWmDCNSh6qPgNhbmeWydU6qrkdE2qV2UEokEkyZMgUTJkzAnTt3kJSUhFq1asHS0rJIhQQGBuLatWsaf8QkJCQEwcHBGt0n0bvUeU6yibstHGVmeJyQmut9SgkAB1nWxNBEpH9qX3odNGgQXr16BVNTU9SqVQtNmjSBpaUlkpOTMWjQoEIVMXLkSOzevRthYWFwdnbOd10HBwc8efJEadmTJ0/g4OCQ6/qTJk1CQkKC4uvBgweFqpEoL+oOJmBsJEGQfy0AWaH4tuzXQf61YGyk2hywRKRdagfl2rVr8fr16xzLX79+jd9++02tfQkhMHLkSOzYsQNHjhyBu7t7gdt4e3vj8OHDSssOHjyomBvzXVKpFNbW1kpfRJpS2BF3Ono6YlnfhnCQKV9edZCZYVnfhnyOksiAqHzpNTExEUIICCHw6tUrmJn99x88MzMTe/fuhb29vVoHDwwMxMaNG7Fz505YWVkp7jPKZDKYm5sDAPr37w8nJyfFpNCjR4+Gr68v5s+fjy5dumDTpk24ePEifvnlF7WOTVRURR2WrqOnI/xqOeB89HM8fZUKe6usy61sSRIZFpWD0sbGBhKJBBKJBNWrV8/xvkQiUfte4LJlywAArVq1Ulq+evVqBAQEAABiY2NhZPRfw7dZs2bYuHEjvvnmG0yePBnVqlXDn3/+mW8HICJN09TYrcZGEnhXsdNChUSkKSoPin7s2DEIIdCmTRts27ZN6REOU1NTuLq6olKlSlorVFM4KDoVFQc4JyoZVM0DlVuUvr6+AIDo6Gi89957kEh4eYhKH4YkUemjdmeeI0eOYOvWrTmWb9myBWvXrtVIUUSGiCFJVDqpHZQhISEoX758juX29vaYOXOmRooiMjQMSaLSS+2gjI2NzfUxDldXV8TGxmqkKCJDwpAkKt3UDkp7e3v8/fffOZZfuXIFdnbsvUclC0OSiNQOyl69emHUqFEICwtDZmYmMjMzceTIEYwePRqfffaZNmok0guGJBEBhRjrdfr06YiJiUHbtm1RpkzW5nK5HP379+c9SioxGJJElE3l5yjfFRUVhStXrsDc3Bx16tSBq6urpmvTCj5HSQVhSBKVDhp/jvJd1atXz3WEHqLijCFJRO9SKSjHjRuH6dOnw8LCAuPGjct33QULFmikMCJdY0gSUW5UCsrw8HCkp6cr/p0XjtZDxRVDkojyUuh7lMUV71HSuxiSRKWTqnmg9uMhRCUJQ5KICqLSpdcePXqovMPt27cXuhgibciUi1znfGRIEpEqVApKmUym+LcQAjt27IBMJoOXlxcA4NKlS3j58qVagUqkC6HX4hC86wbiElIVyxxlZvjKrzIOrJrDkCSiAqkUlKtXr1b8++uvv0bPnj2xfPlyGBsbAwAyMzMxYsQI3vOjPOXVqtOm0GtxGL7+Mt69Cf84IRVjt1yH/ZNMODIkiagAanfmqVChAk6ePAkPDw+l5ZGRkWjWrBmePXum0QI1jZ15dC+vVl2Qfy109HTUyjEz5QLNZx9ROqYSIWCSkYy9w71QrWoVrdRARIZNa515MjIycOvWrRzLb926Bblcru7uqITLbtW9G1iPE1IxfP1lhF6L08pxz0c/zzskAUAiQbqJJf6V2Gjl+ERUcqg9Ms/AgQMxePBg3L17F02aNAEAnDt3DrNmzcLAgQM1XiAVX5lygeBdN3Jc+gQAAUACIHjXDfjVctD4Zdinr/IJyUKsR0Sll9pBOW/ePDg4OGD+/PmIi8tqDTg6OmLChAn48ssvNV4gFV8FteoEgLiEVJyPfg7vKpqdos3eykyj6xFR6aV2UBoZGeGrr77CV199hcTERADgvT7KlT5bdU3cbeEoM8PjhNRcW7QSAA6yrE5FRET5KdSAAxkZGTh06BB+//13xbB1jx49QlJSkkaLo+JNn606YyMJvvKrDCEE8E5/teyLvEH+tbTe85aIij+1W5T3799Hx44dERsbi7S0NPj5+cHKygqzZ89GWloali9fro06qRjSZ6suKSkJB1bNgf2TTLxwb4t0E0vFew5a7nFLRCWL2kE5evRoeHl54cqVK7Cz++++0ocffoihQ4dqtDgq3oyNJAjyr4Xh6y9DAiiFpTZbdW+PuONobY0Nw73wr8RGp89wElHJoXZQnjhxAqdPn4apqanScjc3Nzx8+FBjhVHJ0NHTEcv6NszxHKW2WnV5DUtXTaNHIaLSRO2glMvlyMzMzLH8n3/+gZWVlUaKopKlo6cj/Go5aH1kHo7dSkTaoHZnnvbt22PRokWK1xJJ1uDSQUFB6Ny5syZroxLE2EgC7yp26FbfCd5V7BiSRFRsqD2E3YMHD9CxY0cIIXD79m14eXnh9u3bKF++PI4fPw57e3tt1aoRHMKu5GFIElFhqJoHhZq4OSMjA5s3b8aVK1eQlJSEhg0bok+fPjA3Ny9S0brAoCxZGJJEVFhaCcr09HTUqFEDu3fvRs2aNTVSqK4xKEsOhiQRFYVWBkU3MTFBairHxiT9Y0gSka6o3ZknMDAQs2fPRkZGhjbqISoQQ5KIdEntx0MuXLiAw4cP48CBA6hTpw4sLCyU3t++fbvGiiN6F0OSiHRN7aC0sbHBRx99pI1aiPLFkCQifVA7KFevXq2NOojyxZAkIn1R+R6lXC7H7Nmz4ePjg8aNG2PixIl4/fq1NmsjAsCQJCL9UjkoZ8yYgcmTJ8PS0hJOTk5YvHgxAgMDtVkbEUOSiPRO5aD87bffsHTpUuzfvx9//vkndu3ahQ0bNkAul2uzPirFGJJEZAhUDsrY2FilsVzbtWsHiUSCR48eaaUwKt0YkkRkKFQOyoyMDJiZKc9Eb2JigvT0dI0XRaUbQ5KIDInKvV6FEAgICIBUKlUsS01Nxeeff670LCWfo6RMuSj0lFoMSSIyNCoH5YABA3Is69u3b5EOfvz4ccydOxeXLl1CXFwcduzYge7du+e7zYYNGzBnzhzcvn0bMpkMnTp1wty5c2FnZ1ekWkgzQq/F5Zik2VHFSZoZkkRkiFQOSm08P5mcnIx69eph0KBB6NGjR4Hrnzp1Cv3798fChQvh7++Phw8f4vPPP8fQoUPZkjUAodfiMHz9Zbw7yv7jhFQMX38Zy/o2REdPx1xbnK9TkhmSRGSQ1B5wQJM6deqETp06qbz+mTNn4ObmhlGjRgEA3N3dMWzYMMyePVtbJZKKMuUCwbtu5AhJABAAJACCd92AXA5M36Pc4qxoZYpKj0/jJUOSiAyQ2oOi65O3tzcePHiAvXv3QgiBJ0+eYOvWrUq9cUk/zkc/Vwq/dwkAcQmpGLHxco71niSmIbxsI8id6jIkicjgFKug9PHxwYYNG/Dpp5/C1NQUDg4OkMlkWLJkSZ7bpKWlITExUemLNO/pqyJMvybJ6uiT4tEZrm7uGqqIiEgzilVQ3rhxA6NHj8Z3332HS5cuITQ0FDExMfj888/z3CYkJAQymUzx5eLiosOKSw97K7OCV8qPRIL45Aycj36umYKIiDREIoTI7baSzkkkkgJ7vfbr1w+pqanYsmWLYtnJkyfRokULPHr0CI6OOXtVpqWlIS0tTfE6MTERLi4uBc5oTerJlAs0n30EjxNSc71PqarFn9VHt/pOGquLiCgviYmJkMlkBeZBsWpRpqSkwMhIuWRjY2MAWc955kYqlcLa2lrpizTP2EiCIP9aALI67rxNtScosxS5ZUpEpGF6DcqkpCREREQgIiICABAdHY2IiAjExsYCACZNmoT+/fsr1vf398f27duxbNky3Lt3D6dOncKoUaPQpEkTVKpUSR+nQG/p6OmIZX0bwkGmHHYOMjMs7d0AFa1MgTz+oJEg63nLJu62OqiUiEh1en085OLFi2jdurXi9bhx4wBkDW6wZs0axMXFKUITAAICAvDq1Sv89NNP+PLLL2FjY4M2bdrw8RAD0tHTEX61HHJ9TrLS49N4UrZRVlhK/mtnZv8ryL+WyiP4EBHpisHco9QVVa9Jk2pUGa7u7RF35E51keLRGfHJGYr3VR25h4hIk1TNA722KKl4U2W4uhzD0k0bDVc390KPBUtEpGtsUVKh5DVcXbalvRugZWVrDktHRAaLLUrSmvyGq8s28vdw1E26iASGJBEVc8Xq8RAyDAUNVwcAcgFEcFg6IioBGJSkNnWGq+OwdERU3DEoSW0qDwrAYemIqARgUJLaGrmWgzqdVIs0YDoRkZ4xKEltl+6/gFyNvtIclo6IijMGJalNnRYih6UjouKOQUlqU6eFyGHpiKi4Y1CS2pq42+Y7wDkAGEmApb0bclg6Iir2GJSktuwBzgHkGZY/9WqAznUZkkRU/DEoSS3ZY7e+vHYUrg8PooKlidL7jjIzLO/bEJ3rctozIioZOIQdqYwDnBNRacSgJJXkCMm3hqXzrmKn5+qIiLSHl16pQPmFJBFRScegpHwxJImotGNQUp4YkkREDErKA0OSiCgLg5JyYEgSEf2HQUlKGJJERMoYlKSQkPgKgybOwsV4wMjBA/PmL2BIElGpx+coCQCw42I0Jv1xCamW3kBVIB5A/y0xCPIvy/FaiahUY4uyFMmUC5y5+ww7Ix7izN1nyPz/SSV3XIzG2C3XkSpRnhXkcUIqhq+/jNBrcfool4jIILBFWUqEXotD8K4biEv4by5JR5kZxretjClbLwESM0CiPPScACABELzrBvxqOXBoOiIqldiiLAVCr8Vh+PrLSiEJZLUYv9x2HalG5jlCMpsAEJeQivPRz3VQKRGR4WFQlnCZcoHgXTeQ22RYec8mmdPTV6kFr0REVAIxKEu489HPc7QkleTRknyXvZVZwSsREZVAvEdZwhW1JSgB4CDLmj6LiKg0YouyhFOnJfhu2zL7dZB/LXbkIaJSi0FZwjVxt4WjzCxHCGaTIKv369LeDeAgUw5VB5kZlvVtyOcoiahU46XXEs7YSIKv/Cpj7JbrWQveuif5douxo6cjOng64nz0czx9lQp7q6zLrWxJElFpx6As4ZKSknBg1RzYP8nEC/e2SDexVLznIDNThCSQFareVez0VSoRkUFiUJZgbw9w7mhtjQ3DvfCvxIYtRiIiNTAoS6i8ZgGppu/CiIiKGXbmKYE4VRYRkeYwKEsYhiQRkWbx0qsGZMqFQfQWZUgSEWkeg7KI8pqV4+3epLrAkCQi0g5eei2C/Gbl0OU8jgxJIiLtYVAWkiqzcgTvuqGYHFlbGJJERNql16A8fvw4/P39UalSJUgkEvz5558FbpOWloYpU6bA1dUVUqkUbm5u+PXXX7Vf7DsKmpVDF/M4GlJIZsoFztx9hp0RD3Hm7jOt/4FARKQrer1HmZycjHr16mHQoEHo0aOHStv07NkTT548wapVq1C1alXExcVBLpdrudKcVJ2VQ1vzOBpSSBrKfVoiIm3Qa1B26tQJnTp1Unn90NBQHDt2DPfu3YOtbda0T25ublqqLn+qzsqhjXkcDS0kh6+/nOMSdPZ9Wg6qTkTFXbG6R/nXX3/By8sLc+bMgZOTE6pXr47x48fj9evXOq9F1Vk5ND2PoyGFpKHcpyUi0qZi9XjIvXv3cPLkSZiZmWHHjh34999/MWLECDx79gyrV6/OdZu0tDSkpaUpXicmJmqkFmMjCYL8a2H4+suQAEphoa15HA0pJAH17tNysHUiKq6KVYtSLpdDIpFgw4YNaNKkCTp37owFCxZg7dq1ebYqQ0JCIJPJFF8uLi4aq6ejpyOW9W2ok3kck5KSMH7CBIQ/TEamcwMMmPA9XN3cNbb/wtD3fVoiIl0oVi1KR0dHODk5QSaTKZbVrFkTQgj8888/qFYt55DfkyZNwrhx4xSvExMTNR6WfrUctDoyT1JSEvpNmo+/zXyRWcsKADBudyzmnniq1w4z+rxPS0SkK8WqRenj44NHjx4hKSlJsSwqKgpGRkZwdnbOdRupVApra2ulL03LnsexW30neFex00pIhpdthExTS6X3dD2wwbv0dZ+WiEiX9BqUSUlJiIiIQEREBAAgOjoaERERiI2NBZDVGuzfv79i/d69e8POzg4DBw7EjRs3cPz4cUyYMAGDBg2Cubm5Pk5Bq7Ivt/5dxiNrgUQ5kvTdYSb7Pi2AHGGprfu0RES6ptegvHjxIho0aIAGDRoAAMaNG4cGDRrgu+++AwDExcUpQhMALC0tcfDgQbx8+RJeXl7o06cP/P398cMPP+ilfm3K7rgT8TAZmVKrHCGZTRcDG+RHl/dpiYj0QSKEKFV99xMTEyGTyZCQkKCVy7Ca8Hbv1kznBoh1alvgNos/q49u9Z10UF3uDGUGFSIiVamaB8WqM09pkJD4CkMnz0Z0PGDt4IHAoQEYv/dBgdvpu8NM9n1aIqKShkGpReq2snZcjMakPy4h1dIbqArEA5h38ilsypogISU91wf7Jci6zMkOM0RE2sGg1BJ1xz/dcTEaY7dcByTKLcMniWmKgNTVwAZERPSfYvV4SHGh7jyVCYmvMOmPS1kvcunZKgFgU9YEFa3ZYYaISNfYotSwgsY/lSDrcQ6/Wg4wNpIgKSkJQyfPzrrcmgcB4GVKOjYMbggjIwk7zBAR6RCDUsPUGf+0TkUpJkyYgOh4AFUL3ve/yWl67dlKRFQa8dKrhqk6rumD+IT/Bjg3Ue0JHX33bCUiKo0YlBqmaphtXf+rYhaQ5dMncCg4IiIDxaDUMFXGPzWTv8bT66cVU2VVq1qFQ8ERERkoBqWGFTT+qRAC1ncOQGZtpTSfJIeCIyIyTBzCTktye47STP4a1ncOwDHzSZ6TLnMoOCIi3VA1DxiUWpQdeg/iE7B1/a94ev10jpYkERHph6p5wEuvWmRsJEGdilKErZ2L+OunGJJERMUQg1KL3p4FJLvjDkOSiKh4YVBqCUOSiKhkYFBqAUOSiKjkYFBqGEOSiKhk4VivGlTYkOQjIUREhotBqSGFDUl1560kIiLd4qVXDcgrJDPlAmfuPsPOiIc4c/cZMuXKj6yqO28lERHpHluURZRXSBbUUlR33koiItIPtiiLIL+QLKilqM68lUREpD8MykLK73Jrfi1FIKul+DjhtUrHUXV+SyIi0g4GZSHk13FH1Zbi8+Q3Kh2LkzUTEekXg1JNBfVuVbUFaGsp5WTNRETFAINSDao8AqJqC9DB2oyTNRMRFQMMShWp+pxkE3dblVuKnKyZiMjw8fEQFagzmICxkQRB/rUwfP1lSAClTj25tRQ7ejrCr5YDR+YhIjJQnLi5ABxxh4ioZFI1D9iizEdRBjhnS5GIqGRgUOZBE7OAGBtJ4F3FTksVEhGRLrAzTy44VRYREWVjUL6DIUlERG9jUL6FIUlERO9iUP4/hiQREeWGQQmGJBER5a3UByVDkoiI8lOqg5IhSUREBSm1QcmQJCIiVZTaAQe+/fZb3Lt3jyFJRET5KrUtyqioKIYkEREVqNS1KLPHgDczM8O0adNQvnx5JCYm6rkqIiLStezf/QXNDVLqZg/5559/4OLiou8yiIjIQDx48ADOzs55vl/qglIul+PRo0ewsrKCRFL8ZvJITEyEi4sLHjx4oNI0YYaO52PYeD6GjedTNEIIvHr1CpUqVYKRUd53IkvdpVcjI6N8/3IoLqytrUvEf4xsPB/DxvMxbDyfwpPJZAWuU2o78xAREamCQUlERJQPBmUxI5VKERQUBKlUqu9SNILnY9h4PoaN56Mbpa4zDxERkTrYoiQiIsoHg5KIiCgfDEoiIqJ8MCiJiIjywaA0MMePH4e/vz8qVaoEiUSCP//8s8Bt0tLSMGXKFLi6ukIqlcLNzQ2//vqr9otVQWHOZ8OGDahXrx7Kli0LR0dHDBo0CM+ePdN+sQUICQlB48aNYWVlBXt7e3Tv3h2RkZEFbrdlyxbUqFEDZmZmqFOnDvbu3auDagtWmPNZsWIFWrRogXLlyqFcuXJo164dzp8/r6OK81fYn0+2TZs2QSKRoHv37torUg2FPZ+XL18iMDAQjo6OkEqlqF69ukF85gp7PosWLYKHhwfMzc3h4uKCsWPHIjU1VQcV/4dBaWCSk5NRr149LFmyROVtevbsicOHD2PVqlWIjIzE77//Dg8PDy1WqTp1z+fUqVPo378/Bg8ejOvXr2PLli04f/48hg4dquVKC3bs2DEEBgbi7NmzOHjwINLT09G+fXskJyfnuc3p06fRq1cvDB48GOHh4ejevTu6d++Oa9eu6bDy3BXmfI4ePYpevXohLCwMZ86cgYuLC9q3b4+HDx/qsPLcFeZ8ssXExGD8+PFo0aKFDipVTWHO582bN/Dz80NMTAy2bt2KyMhIrFixAk5OTjqsPHeFOZ+NGzdi4sSJCAoKws2bN7Fq1Sps3rwZkydP1mHlAAQZLABix44d+a6zb98+IZPJxLNnz3RTVBGocj5z584VlStXVlr2ww8/CCcnJy1WVjhPnz4VAMSxY8fyXKdnz56iS5cuSsuaNm0qhg0bpu3y1KbK+bwrIyNDWFlZibVr12qxssJR9XwyMjJEs2bNxMqVK8WAAQNEt27ddFOgmlQ5n2XLlonKlSuLN2/e6LCywlHlfAIDA0WbNm2Ulo0bN074+PhouzwlbFEWc3/99Re8vLwwZ84cODk5oXr16hg/fjxev36t79IKxdvbGw8ePMDevXshhMCTJ0+wdetWdO7cWd+l5ZCQkAAAsLW1zXOdM2fOoF27dkrLOnTogDNnzmi1tsJQ5XzelZKSgvT0dLW20RVVz2fatGmwt7fH4MGDdVFWoalyPn/99Re8vb0RGBiIihUrwtPTEzNnzkRmZqauylSZKufTrFkzXLp0SXF5/969e9i7d6/ufx/oNJZJLVChBdahQwchlUpFly5dxLlz58SePXuEq6urCAgI0E2RalDlfIQQ4o8//hCWlpaiTJkyAoDw9/c3uL+QMzMzRZcuXQr8y9bExERs3LhRadmSJUuEvb29NstTm6rn867hw4eLypUri9evX2upssJR9XxOnDghnJycRHx8vBBCGGyLUtXz8fDwEFKpVAwaNEhcvHhRbNq0Sdja2oqpU6fqqFLVqPN5W7x4sTAxMVH8Pvj88891UKEyBqUBUyVY/Pz8hJmZmXj58qVi2bZt24REIhEpKSlarlA9qpzP9evXhaOjo5gzZ464cuWKCA0NFXXq1BGDBg3STZEq+vzzz4Wrq6t48OBBvusVl6BU9XzeFhISIsqVKyeuXLmixcoKR5XzSUxMFG5ubmLv3r2KZYYalKr+fKpVqyZcXFxERkaGYtn8+fOFg4ODtktUi6rnExYWJipWrChWrFgh/v77b7F9+3bh4uIipk2bpqNKszAoDZgqwdK/f39RpUoVpWU3btwQAERUVJQWq1OfKufTt29f8fHHHystO3HihAAgHj16pMXqVBcYGCicnZ3FvXv3ClzXxcVFLFy4UGnZd999J+rWraul6tSnzvlkmzt3rpDJZOLChQtarKxwVD2f8PBwAUAYGxsrviQSiZBIJMLY2FjcuXNHRxXnT52fT8uWLUXbtm2Vlu3du1cAEGlpadoqUS3qnE/z5s3F+PHjlZatW7dOmJubi8zMTG2VmAPvURZzPj4+ePToEZKSkhTLoqKiiu28mykpKTkmUDU2NgaQNcmqPgkhMHLkSOzYsQNHjhyBu7t7gdt4e3vj8OHDSssOHjwIb29vbZWpssKcDwDMmTMH06dPR2hoKLy8vLRcperUPZ8aNWrg6tWriIiIUHx17doVrVu3RkREBFxcXHRUee4K8/Px8fHBnTt3IJfLFcuioqLg6OgIU1NTbZZboMKcj8H8PtBZJJNKXr16JcLDwxV/7S5YsECEh4eL+/fvCyGEmDhxoujXr5/S+s7OzuLjjz8W169fF8eOHRPVqlUTQ4YM0dcpKFH3fFavXi3KlCkjli5dKu7evStOnjwpvLy8RJMmTfR1CgrDhw8XMplMHD16VMTFxSm+3r7E3a9fPzFx4kTF61OnTokyZcqIefPmiZs3b4qgoCBhYmIirl69qo9TUFKY85k1a5YwNTUVW7duVdrm1atX+jgFJYU5n3cZ0qXXwpxPbGyssLKyEiNHjhSRkZFi9+7dwt7eXnz//ff6OAUlhTmfoKAgYWVlJX7//Xdx7949ceDAAVGlShXRs2dPndbOoDQwYWFhAkCOrwEDBgghsv4j+/r6Km1z8+ZN0a5dO2Fubi6cnZ3FuHHjDOb+ZGHO54cffhC1atUS5ubmwtHRUfTp00f8888/ui/+HbmdBwCxevVqxTq+vr6Kc8v2xx9/iOrVqwtTU1NRu3ZtsWfPHt0WnofCnI+rq2uu2wQFBem8/ncV9ufzNkMKysKez+nTp0XTpk2FVCoVlStXFjNmzFC6Z6kvhTmf9PR0MXXqVFGlShVhZmYmXFxcxIgRI8SLFy90Wjun2SIiIsoH71ESERHlg0FJRESUDwYlERFRPhiURERE+WBQEhER5YNBSURElA8GJRERUT4YlEQFmDp1KurXr6/vMtCqVSuMGTNGrW0CAgLQvXt3rdRTkDVr1sDGxkYvx9YGNzc3LFq0SN9lkB4wKElnHj9+jNGjR6Nq1aowMzNDxYoV4ePjg2XLliElJUXf5RXa0aNHIZFI8PLlS4PcHxEVTRl9F0Clw7179+Dj4wMbGxvMnDkTderUgVQqxdWrV/HLL7/AyckJXbt2zXXb9PR0mJiY6LhizXvz5o3eB6YmIvWxRUk6MWLECJQpUwYXL15Ez549UbNmTVSuXBndunXDnj174O/vr1hXIpFg2bJl6Nq1KywsLDBjxgwAwLJly1ClShWYmprCw8MD69atU2wTExMDiUSCiIgIxbKXL19CIpHg6NGjAP5rqR0+fBheXl4oW7YsmjVrhsjISKVaZ82ahYoVK8LKygqDBw9GampqnucVExOD1q1bAwDKlSsHiUSCgIAAAFmXSkeOHIkxY8agfPny6NChQ4F15rc/AJDL5fjqq69ga2sLBwcHTJ06VaXvf3BwMCpUqABra2t8/vnnePPmjeK90NBQNG/eHDY2NrCzs8MHH3yAu3fv5vjebt++Ha1bt0bZsmVRr149nDlzRukYa9aswXvvvYeyZcviww8/xLNnz3LU8f3338Pe3h5WVlYYMmQIJk6cqHRZ+8KFC/Dz80P58uUhk8ng6+uLy5cvK+0j+/PRqVMnmJubo3Llyti6dWue5/7LL7+gUqVKSjNqAEC3bt0waNAgAMDdu3fRrVs3VKxYEZaWlmjcuDEOHTqU5z5V+bwBwLVr19CpUydYWlqiYsWK6NevH/7991/F+1u3bkWdOnVgbm4OOzs7tGvXDsnJyXkel/REpyPLUqn077//ColEIkJCQlRaH4Cwt7cXv/76q7h79664f/++2L59uzAxMRFLliwRkZGRYv78+cLY2FgcOXJECCFEdHS0ACDCw8MV+3nx4oUAIMLCwoQQ/w3Q3rRpU3H06FFx/fp10aJFC9GsWTPFNps3bxZSqVSsXLlS3Lp1S0yZMkVYWVmJevXq5VprRkaG2LZtmwAgIiMjRVxcnGISbV9fX2FpaSkmTJggbt26JW7dulVgnQXtz9raWkydOlVERUWJtWvXColEIg4cOJDn93LAgAHC0tJSfPrpp+LatWti9+7dokKFCmLy5MmKdbZu3Sq2bdsmbt++LcLDw4W/v7+oU6eOYr6/7Jpr1Kghdu/eLSIjI8XHH38sXF1dRXp6uhBCiLNnzwojIyMxe/ZsERkZKRYvXixsbGyETCZTHGf9+vXCzMxM/PrrryIyMlIEBwcLa2trpe/t4cOHxbp168TNmzfFjRs3xODBg0XFihVFYmKi0ufDzs5OrFixQkRGRopvvvlGGBsbixs3buT6PXj+/LkwNTUVhw4dUix79uyZ0rKIiAixfPlycfXqVREVFSW++eYbYWZmppjlRoisAeGz5xZV5fP24sULUaFCBTFp0iRx8+ZNcfnyZeHn5ydat24thBDi0aNHokyZMmLBggUiOjpa/P3332LJkiUGMRMLKWNQktadPXtWABDbt29XWm5nZycsLCyEhYWF+OqrrxTLAYgxY8YordusWTMxdOhQpWWffPKJ6Ny5sxBCvaB8+xfmnj17BADx+vVrIYQQ3t7eYsSIEUrHadq0aZ5B+fZ+353RwNfXVzRo0EBpmTp15ra/5s2bKy1r3Lix+Prrr/OsbcCAAcLW1lYkJycrli1btkxYWlrmOfFtfHy8AKCYCiy75pUrVyrWuX79ugAgbt68KYQQolevXoqfRbZPP/1UKSibNm0qAgMDldbx8fHJ93ubmZkprKysxK5duxTLAIjPP/9cab2mTZuK4cOH57mfbt26iUGDBile//zzz6JSpUr5Tv5bu3Zt8eOPPypeqxuU06dPF+3bt1fa54MHDxR/BF26dEkAEDExMXnWQIaBl15Jb86fP4+IiAjUrl0baWlpSu+9OyHwzZs34ePjo7TMx8cHN2/eVPu4devWVfzb0dERAPD06VPFcZo2baq0flEmWW7UqFGht83N27UDWfVn156XevXqoWzZsorX3t7eSEpKwoMHDwAAt2/fRq9evVC5cmVYW1vDzc0NABAbG5vnsQvzfYuMjESTJk2Ulr37+smTJxg6dCiqVasGmUwGa2trJCUl5ajl3X17e3vn+1no06cPtm3bpvicbdiwAZ999pliUuCkpCSMHz8eNWvWhI2NDSwtLXHz5s0cx1XHlStXEBYWBktLS8VXjRo1AGRd6q1Xrx7atm2LOnXq4JNPPsGKFSvw4sWLQh+PtIedeUjrqlatColEkuNeYOXKlQEA5ubmObaxsLBQ6xjZv/DEW7PGpaen57ru2x2DJBIJAOS4f6Up756HOnXm5t1OTRKJpMi1+/v7w9XVFStWrFDcy/P09FS6j/nusbX1fRswYACePXuGxYsXw9XVFVKpFN7e3jlqUZe/vz+EENizZw8aN26MEydOYOHChYr3x48fj4MHD2LevHmoWrUqzM3N8fHHH+d5XFV+jklJSfD398fs2bNzbO/o6AhjY2McPHgQp0+fxoEDB/Djjz9iypQpOHfuHNzd3Yt0vqRZbFGS1tnZ2cHPzw8//fRToTsq1KxZE6dOnVJadurUKdSqVQsAUKFCBQBAXFyc4v23O1qoc5xz584pLTt79my+22T3ZM3MzCxw/6rUqc7+VHHlyhW8fv1a8frs2bOwtLSEi4sLnj17hsjISHzzzTdo27YtatasWahWjSrfNw8PD1y4cEFp2buvT506hVGjRqFz586oXbs2pFKpUueXvPZ99uxZ1KxZM8/6zMzM0KNHD2zYsAG///47PDw80LBhQ6XjBgQE4MMPP0SdOnXg4OCAmJiYPPenys+xYcOGuH79Otzc3FC1alWlr+w/oCQSCXx8fBAcHIzw8HCYmppix44deR6X9IMtStKJpUuXwsfHB15eXpg6dSrq1q0LIyMjXLhwAbdu3SrwEuWECRPQs2dPNGjQAO3atcOuXbuwfft2Rc9Ec3NzvP/++5g1axbc3d3x9OlTfPPNN2rXOXr0aAQEBMDLyws+Pj7YsGEDrl+/rmj95sbV1RUSiQS7d+9G586dYW5uDktLy1zXVaVOdfanijdv3mDw4MH45ptvEBMTg6CgIIwcORJGRkYoV64c7Ozs8Msvv8DR0RGxsbGYOHGi2scYNWoUfHx8MG/ePHTr1g379+9HaGio0jpffPEFhg4dCi8vLzRr1gybN2/G33//rfS9rVatGtatWwcvLy8kJiZiwoQJuV5x2LJlC7y8vNC8eXNs2LAB58+fx6pVq/KtsU+fPvjggw9w/fp19O3bV+m9atWqYfv27fD394dEIsG3336bb2tZlZ9jYGAgVqxYgV69eil6Kt+5cwebNm3CypUrcfHiRRw+fBjt27eHvb09zp07h/j4+HwDn/REv7dIqTR59OiRGDlypHB3dxcmJibC0tJSNGnSRMydO1epswkAsWPHjhzbL126VFSuXFmYmJiI6tWri99++03p/Rs3bghvb29hbm4u6tevLw4cOFBgJ5nw8HABQERHRyuWzZgxQ5QvX15YWlqKAQMGiK+++irfDidCCDFt2jTh4OAgJBKJGDBggBAiq/PN6NGjc6xbUJ3q7K9bt26K93MzYMAA0a1bN/Hdd98JOzs7YWlpKYYOHSpSU1MV6xw8eFDUrFlTSKVSUbduXXH06FGln4EqHVeEEGLVqlXC2dlZmJubC39/fzFv3jylzjzZ55X9vR00aJAYNWqUeP/99xXvX758WXh5eQkzMzNRrVo1sWXLFqVONEJkfT6WLFki/Pz8hFQqFW5ubmLz5s15fg+yZWZmCkdHRwFA3L17V+m96Oho0bp1a2Fubi5cXFzETz/9lOP7/W4dqvwco6KixIcffihsbGyEubm5qFGjhhgzZoyQy+Xixo0bokOHDqJChQpCKpWK6tWrK3UeIsMhEeKti+xERDrk5+cHBwcHpWdiCyKRSLBjxw69Dc1HpQ8vvRKRTqSkpGD58uXo0KEDjI2N8fvvv+PQoUM4ePCgvksjyheDkoh0QiKRYO/evZgxYwZSU1Ph4eGBbdu2oV27dvoujShfvPRKRESUDz4eQkRElA8GJRERUT4YlERERPlgUBIREeWDQUlERJQPBiUREVE+GJRERET5YFASERHlg0FJRESUj/8DatD4KHRwZbMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 500x1000 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get the predictions of x_test into `y_pred`\n",
    "\n",
    "\n",
    "# ...\n",
    "#\n",
    "\n",
    "y_pred = x.predict(x_test)\n",
    "import matplotlib.pyplot as plt\n",
    "fig, ax = plt.subplots(figsize=(5, 10))\n",
    "\n",
    "ax.scatter(y_test, y_pred)\n",
    "\n",
    "lims = [\n",
    "    np.min([ax.get_xlim(), ax.get_ylim()]),\n",
    "    np.max([ax.get_xlim(), ax.get_ylim()]),\n",
    "]\n",
    "ax.plot(lims, lims, 'k-', alpha=0.75, zorder=0)\n",
    "ax.set_aspect('equal')\n",
    "ax.set_xlim(lims)\n",
    "ax.set_ylim(lims)\n",
    "\n",
    "ax.set_title('Parity Plot of Custom Linear Regression')\n",
    "ax.set_xlabel('Ground truth bandgap values')\n",
    "ax.set_ylabel('Predicted bandgap values')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dominant-chaos",
   "metadata": {},
   "source": [
    "### 4.2 Implement Ridge regression\n",
    "4.2.1 Explain Ridge regression briefly in 1-2 lines."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "happy-cyprus",
   "metadata": {},
   "source": [
    "<!-- Your answer to 1.2.1 -->\n",
    "Ridge regression does the same thing that Vanilla regression does but with a variation in the loss function. It employs a regularization technique callled as L2 regularization. Regularization is used to create simpler models, and this avoid overfitting and feature selection. \n",
    "\n",
    "$J = \\frac{1}{M}\\sum_{i = 1}^m(y_i - \\hat{y_i})^2 + \\alpha \\sum \\theta_i^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tight-forwarding",
   "metadata": {},
   "source": [
    "4.2.2 Implement Ridge regression and make a table of different RMSE scores you achieved with different values of alpha. What does the parameter `alpha` do? How does it affect the results here? Explain in 5-10 lines in total. (You can use scikit-learn from this cell onwards)\n",
    "\n",
    "- The parameter alpha, determines how strongly the regularization is applied.\n",
    "- Must belong to [0, $\\infty$]\n",
    "- Smaller the value of alpha, smaller the regularization. Zero alpha gives the same results as a vanilla linear regression. (This is evident from the fact that Ridge and Lasso Regression techniques, give the same MSE, with $\\alpha = 0$)\n",
    "- If a sufficiently large value of $\\alpha$ would result in the paramters with L2 norm with small coefficients, but not zero. \n",
    "- Hence the alpha value must be chose in such a way that the MSE and the L2 Norm are minimized properly.\n",
    "- Too large of a alpha, would result in a simpler model, but maybe way too simpler that could be underfitting. ($\\uparrow$ Bias)\n",
    "- Too small of a alpha, would result in a complex model, but may have to risk the overfitting. ($\\uparrow$ Variance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "violent-northern",
   "metadata": {},
   "outputs": [],
   "source": [
    "# you should not have imported sklearn before this point\n",
    "import sklearn\n",
    "from sklearn.linear_model import Ridge\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "alpha_values = [0, 0.001, 0.01, 0.1, 1.0, 10]\n",
    "mse_dict = []\n",
    "for alpha in alpha_values:\n",
    "    ridge_regression = Ridge(alpha = alpha)\n",
    "    ridge_regression.fit(x_train, y_train)\n",
    "    ridge_y_pred = ridge_regression.predict(x_test)\n",
    "    mse_dict += [np.sqrt(sum((ridge_y_pred - y_test)**2))]\n",
    "    # implement Ridge regression and make a table where you explore the effect of different values of `alpha`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "69906b75",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Alpha</th>\n",
       "      <th>RMSE</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.283468</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.001</td>\n",
       "      <td>0.283459</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.010</td>\n",
       "      <td>0.283398</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.100</td>\n",
       "      <td>0.284385</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.000</td>\n",
       "      <td>0.384821</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>10.000</td>\n",
       "      <td>1.153260</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Alpha      RMSE\n",
       "0   0.000  0.283468\n",
       "1   0.001  0.283459\n",
       "2   0.010  0.283398\n",
       "3   0.100  0.284385\n",
       "4   1.000  0.384821\n",
       "5  10.000  1.153260"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame({\"Alpha\": alpha_values, \"RMSE\": mse_dict})\n",
    "df.set_index(\"Alpha\")\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adolescent-temperature",
   "metadata": {},
   "source": [
    "### 4.3 Implement Lasso regression\n",
    "4.3.1 Explain Lasso regression briefly in 1-2 lines."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "144cedf1",
   "metadata": {},
   "source": [
    "Lasso Regularization, just like Ridge regression, is a variant of the vanilla linear regression, which uses regularization techniques to produce simpler models, to address overfitting and feature selection. Instead of L2 regularization, like Ridge Regression, Lasso Regression uses L1 regularization technique. \n",
    "\n",
    "$J \\frac{1}{M}\\sum_{i = 1}^m(y_i - \\hat{y_i})^2 + \\alpha\\sum \\vert \\theta_i \\vert$ "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "popular-wonder",
   "metadata": {},
   "source": [
    "4.3.2 Implement Lasso regression and make a table of different RMSE scores you achieved with different values of alpha. What does the parameter `alpha` do? How does it affect the results here? Explain in 5-10 lines in total.\n",
    "\n",
    "- The parameter alpha, determines how strongly the regularization is applied.\n",
    "- Must belong to [0, $\\infty$]\n",
    "- Smaller the value of alpha, smaller the regularization. Zero alpha gives the same results as a vanilla linear regression. (This is evident from the fact that Ridge and Lasso Regression techniques, give the same MSE, with $\\alpha = 0$)\n",
    "- If a sufficiently large value of $\\alpha$ would result in the paramters with L2 norm with small coefficients, tend to make them zero. \n",
    "- Hence the alpha value must be chose in such a way that the MSE and the L2 Norm are minimized properly.\n",
    "- Too large of a alpha, would result in a simpler model, but maybe way too simpler that could be underfitting. ($\\uparrow$ Bias)\n",
    "- Too small of a alpha, would result in a complex model, but may have to risk the overfitting. ($\\uparrow$ Variance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "extra-brighton",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/13/_wyhp36n6jgdfcwhj2x6mm940000gn/T/ipykernel_28715/3429651262.py:8: UserWarning: With alpha=0, this algorithm does not converge well. You are advised to use the LinearRegression estimator\n",
      "  lasso_regression.fit(x_train, y_train)\n",
      "/usr/local/lib/python3.10/site-packages/sklearn/linear_model/_coordinate_descent.py:648: UserWarning: Coordinate descent with no regularization may lead to unexpected results and is discouraged.\n",
      "  model = cd_fast.enet_coordinate_descent(\n",
      "/usr/local/lib/python3.10/site-packages/sklearn/linear_model/_coordinate_descent.py:648: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations, check the scale of the features or consider increasing regularisation. Duality gap: 2.623e-01, tolerance: 1.784e-03 Linear regression models with null weight for the l1 regularization term are more efficiently fitted using one of the solvers implemented in sklearn.linear_model.Ridge/RidgeCV instead.\n",
      "  model = cd_fast.enet_coordinate_descent(\n"
     ]
    }
   ],
   "source": [
    "# implement Lasso regression and make a table where you explore the effect of different values of `alpha`\n",
    "from sklearn.linear_model import Lasso\n",
    "\n",
    "alpha_values = [0, 0.001, 0.01, 0.1, 1.0, 10]\n",
    "mse_dict = []\n",
    "for alpha in alpha_values:\n",
    "    lasso_regression = Lasso(alpha = alpha)\n",
    "    lasso_regression.fit(x_train, y_train)\n",
    "    ridge_y_pred = lasso_regression.predict(x_test)\n",
    "    mse_dict += [np.sqrt(sum((ridge_y_pred - y_test)**2))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "accompanied-worst",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Alpha</th>\n",
       "      <th>RMSE</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000</td>\n",
       "      <td>0.283468</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.001</td>\n",
       "      <td>0.282557</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.010</td>\n",
       "      <td>0.368660</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.100</td>\n",
       "      <td>2.048338</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.000</td>\n",
       "      <td>2.085298</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>10.000</td>\n",
       "      <td>2.085298</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Alpha      RMSE\n",
       "0   0.000  0.283468\n",
       "1   0.001  0.282557\n",
       "2   0.010  0.368660\n",
       "3   0.100  2.048338\n",
       "4   1.000  2.085298\n",
       "5  10.000  2.085298"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame({\"Alpha\": alpha_values, \"RMSE\": mse_dict})\n",
    "df.set_index(\"Alpha\")\n",
    "df"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.6 64-bit",
   "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.10.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}