{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](img/banner.png)\n",
    "%%HTML\n",
    "<script src=\"require.js\"></script>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bias-Variance Tradeoff\n",
    "\n",
    "**Mahmood Amintoosi, Fall 2024**\n",
    "\n",
    "Computer Science Dept, Ferdowsi University of Mashhad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* [Wiki: Bias variance tradeoff](https://en.wikipedia.org/wiki/Bias%E2%80%93variance_tradeoff)\n",
    "* [MLU-Explain bias-variance](https://mlu-explain.github.io/bias-variance/)\n",
    "* [MLU-Explain double-descent, part 1](https://mlu-explain.github.io/double-descent/)\n",
    "* [MLU-Explain double-descent, part 2](https://mlu-explain.github.io/double-descent2/)\n",
    "* [The Bias-Variance Tradeoff: A Newbie’s Guide, by a Newbie](https://medium.com/@DeepthiTabithaBennet/the-bias-variance-tradeoff-a-newbies-guide-by-a-newbie-95fb03dbebcb)\n",
    "* [bias-variance-trade-off](https://spotintelligence.com/2023/04/11/bias-variance-trade-off/)\n",
    "\n",
    "* Paper: [VC Theoretical Explanation of Double Descent](https://arxiv.org/abs/2205.15549)\n",
    "\n",
    "* Paper: [Reconciling modern machine-learning practice and the classical bias–variance trade-off](https://www.pnas.org/doi/10.1073/pnas.1903070116)\n",
    "  - [Double Descent](https://medium.com/mlearning-ai/double-descent-8f92dfdc442f), [Highlights](misc/medium-com_mlearning-ai_double-descent-highlightes.md)\n",
    "  - [Reproducing Deep Double Descent](https://hippocampus-garden.com/double_descent/), [Highlights](misc/hippocampus-garden-com_double_descent-highlightes.md)\n",
    "    + [deep_double_descent, colab](https://colab.research.google.com/drive/1lT2dUqal90NbLVQIGvseyAdKzH19MH2T?usp=sharing)\n",
    "* [Sec 22.3 of Zaki](https://fumdrive.um.ac.ir/index.php/f/4160875)\n",
    "\n",
    "* Paper: [Understanding the double descent curve in Machine Learning](https://arxiv.org/abs/2211.10322)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import polyfit\n",
    "from numpy import polyval\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import defaultdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return np.sin(x * np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def error_function(pred, actual):\n",
    "    return (pred - actual) ** 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(120)\n",
    "n_observations_per_dataset = 25\n",
    "n_datasets = 1000\n",
    "max_poly_degree = 12  # Maximum model complexity\n",
    "model_poly_degrees = range(1, max_poly_degree + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "NOISE_STD = .5\n",
    "percent_train = .8\n",
    "n_train = int(np.ceil(n_observations_per_dataset * percent_train))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1, 1, n_observations_per_dataset)\n",
    "x = np.random.permutation(x)\n",
    "x_train = x[:n_train]\n",
    "x_test = x[n_train:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_hat = defaultdict(list)\n",
    "\n",
    "pred_train = defaultdict(list)\n",
    "pred_test = defaultdict(list)\n",
    "\n",
    "train_errors = defaultdict(list)\n",
    "test_errors = defaultdict(list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "for dataset in range(n_datasets):\n",
    "\n",
    "    # Simulate training/testing targets\n",
    "    y_train = f(x_train) + NOISE_STD * np.random.randn(*x_train.shape)\n",
    "    y_test = f(x_test) + NOISE_STD * np.random.randn(*x_test.shape)\n",
    "\n",
    "    # Loop over model complexities\n",
    "    for degree in model_poly_degrees:\n",
    "        \n",
    "        # Train model\n",
    "        tmp_theta_hat = polyfit(x_train, y_train, degree)\n",
    "\n",
    "        # Make predictions on train set\n",
    "        tmp_pred_train = polyval(tmp_theta_hat, x_train)\n",
    "        pred_train[degree].append(tmp_pred_train)\n",
    "\n",
    "        # Test predictions\n",
    "        tmp_pred_test = polyval(tmp_theta_hat, x_test)\n",
    "        pred_test[degree].append(tmp_pred_test)\n",
    "\n",
    "        # Mean Squared Error for train and test sets\n",
    "        train_errors[degree].append(np.mean(error_function(tmp_pred_train, y_train)))\n",
    "        test_errors[degree].append(np.mean(error_function(tmp_pred_test, y_test)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_estimator_bias_squared(pred_test):\n",
    "    pred_test = np.array(pred_test)\n",
    "    average_model_prediction = pred_test.mean(0)\n",
    "\n",
    "    return np.mean((average_model_prediction - f(x_test)) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_estimator_variance(pred_test):\n",
    "    pred_test = np.array(pred_test)\n",
    "    average_model_prediction = pred_test.mean(0)\n",
    "    \n",
    "    return np.mean((pred_test - average_model_prediction) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "complexity_train_error = []\n",
    "complexity_test_error = []\n",
    "bias_squared = []\n",
    "variance = []\n",
    "for degree in model_poly_degrees:\n",
    "    complexity_train_error.append(np.mean(train_errors[degree]))\n",
    "    complexity_test_error.append(np.mean(test_errors[degree]))\n",
    "    bias_squared.append(calculate_estimator_bias_squared(pred_test[degree]))\n",
    "    variance.append(calculate_estimator_variance(pred_test[degree]))\n",
    "\n",
    "best_model_degree = model_poly_degrees[np.argmin(complexity_test_error)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Bias-Variance Tradeoff')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVgAAADgCAYAAABGgn7sAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAABROElEQVR4nO2dd3xVRfbAvychIZVUQknoJSIQQDoqoEgRpazIAgIJa2Fh7b3t2lZd22/tYtdgARWVRUWpFhCU3nsJktDSE5KQOr8/7n3P915e2kteGvP9fN4n986dmXvmlZO5Z845I0opNBqNRlPzeNS1ABqNRtNY0QpWo9Fo3IRWsBqNRuMmtILVaDQaN6EVrEaj0bgJrWA1Go3GTWgFq7FDRN4UkX/VtRzlISK7RWR4XcvhbkRkuIgk1lBfc0XktIicFZEwEblYRA6a5xNr4h6a0oj2gz2/EJEEoAVQDBQC64A5Sqnjbrrfm4CfUirWobwXsAFopZRKc8e9awMR+R641DxtCiigwDz/WCk1pxp9Dzf7iKqmjF5AFjBIKbXdLFsFLFFKvVydvjXlo2ew5yfjlFIBQCvgNPCqG+8VD1wjIv4O5TOBb6uiXEWkSY1KVgMopa5USgWY7+cnwHOWc1vlWseytwB8gN02Ze0czjVuQCvY8xil1DlgEXChpUxEPhSRJ83jEBH5VkSSRSTdPI6yqTtLRI6ISLaIHBWR6U7usR5IAibZtPMErgPmi0gnEVktIqkikiIin4hIsE3dBBG5X0R2ADki0sQsu8K8PkBE1otIhoicFJHXRMTbpr0SkTnm43CGiLwuImJz/SYR2WuOYY+IXGSWtxaRL82xHxWR26r6/pr3vllEDgIHzbKXReS4iGSJyGYRudSmvq/5/qeLyB6gv0N/ZcokIk1F5CUROWG+XjLLugL7zWoZ5nt9GOgIfGOaCJpWdWyayqEV7HmMiPgBU4DfyqjiAXyAMdtpC+QBr5lt/YFXgCuVUoHAEGBbGf3MB2xNBFcAXsBSQID/AK2BbkAb4DGH9tOAq4BgpVSRw7Vi4E4gHBgMjAD+4VDnagxlFQP8FRhtjmGyea9YoBkwHkgVEQ/gG2A7EGn2eYeIjC5jfOUxERjIn//ENgK9gVDgU+ALEfExrz0KdDJfo4E4SyeVkOlhYJDZdy9gAPBPpdQBoLtZJ1gpdblSqhPwB+aTjFIq34VxaSqDUkq/zqMXkACcBTIwbLAngJ421z8EniyjbW8g3Tz2N/uYBPhWcM+25r2izPNPgJfLqDsR2Oog7/VOxnBFGe3vAL62OVfAJTbnnwMPmMfLgNud9DEQ+MOh7EHggwrGaffemfe+vII26UAv8/gIMMbm2mwgsTIyAYeBsTbXRgMJ5nF7U5YmlXkP9avmXnoGe34yUSkVjGGXuwX4WURaOlYSET8ReUtEjolIFvALECwinkqpHIzZ7xzgpIh8JyIXmO12m4+eZ0XkUqXUH2bbGSISgKFE55t1W4jIQhFJMu/xMcZs1JYyF+BEpKtpujhltn/aSftTNse5QIB53AZDMTnSDmhtmhQyRCQDeAjDlllV7GQXkXtMk0Sm2W+QjbytHeofq4JMrR3qHzPLNHWIVrDnMUqpYqXUVxiP2Zc4qXI3EA0MVEo1A4aa5WK2X6aUGomxWLYPeMcs767+XOhZY7aJx1jYmgQcVUptNsufxphd9TTvMcPSv62o5QxjnnnvLmb7h5y0L4vjGI/jzsqPKqWCbV6BSqmxlezXFqvspr31PgwzRYj5Ty7TRt6TGErfQtsqyHQCQwnbtj3hgryaGkQr2PMYMZgAhAB7nVQJxLC7ZohIKIaN0NK2hYhMMG2x+Rhmh5Jybvclxo/+cQxla3uPs0CmiEQC91ZxGIEYLkhnzRn03Cq0fRe4R0T6mu9FZxFph+E+lm0urvmKiKeI9BCR/hX0VxlZi4BkoImIPIJh+7XwOfCgGIuLUcCtNtcqkmkB8E8RaS4i4cAjGE8DmjpEK9jzk29E5CyGYnoKiFNKOXPZeQnwBVIwFsJ+sLnmAdyFMUtKA4ZRjnIzTQpfAlEYNlgLjwMXYczkvgO+quJY7sHwSMjGmEF/VtmGSqkvMMb/qdl+MRCqlCrGWBjrDRzFGP+7GI/z1WEZxnt4AOMR/hz2JoHHzfKjwHLgIxtZK5LpSWATsAPYCWwxyzR1iA400Gg0GjehZ7AajUbjJtyqYEVkjIjsF5FDIvKAk+t3mc7dO0RklWn/slwrFpFt5muJO+XUaDQad+A2E4EY0ToHgJFAIoaD9TSl1B6bOpcBvyulckVkLjBcKTXFvHZWGeGHGo1G0yBx5wx2AHBIKXVEKVUALAQm2FZQSv2olMo1T3/DWADRaDSaRoE7FWwk9iukiWZZWdwAfG9z7iMim0TkN9Hp1DQaTQOkXmQnEpEZQD8MVx8L7ZRSSSLSEVgtIjuVUocd2s3GCCfE39+/7wUXXFBrMtckx4/DmTP2ZU2bQo8eVe9r9+7dnDt3DoDIyEhatiwVoFXv2LzZiDno27dvHUui0ZRm8+bNKUqp5q60daeCTcI+KiXKLLNDjKxIDwPDlE3SCaVUkvn3iIj8BPTBIaxRKfU28DZAv3791KZNm2p4CLXDgQMQHW1flp8PzzwDV1xRtb4++OADrr/+esDIM7Fu3Tq8vb0raFW3iJncqqF+fprGjYgcq7iWc9xpItgIdBGRDmKkj5sK2HkDiEgf4C1gvFLqjE15iCWFmhmVcjGwh0ZK164wYkTp8nnzqt7XddddZ521njhxgoULF1ZTOo1G4ypuU7DKSCt3C0b0yl7gc6XUbhF5QkTGm9Wex0i88YWDO1Y3YJOIbAd+BJ6x9T5ojMx1EgP1v/9BUqk5f/k0bdqU2277M3XpCy+8gA4m0WjqhkYTydWQTQQAhYXQrh2cPGlf/thj8OijTpuUSVpaGm3btiUnJweAZcuWMWrUqJoRVKM5zxCRzUqpfq60rReLXOcrhYWFJCYmWhelvvwSMjPt63h6wp49IJXND2WybNkysrOzAfD29mbvXme5XDQajQUfHx+ioqLw8vKqsT61gq1DEhMTCQwMpH379ogIBQWwY0fpeq1aQUhI1frOz89n586d1vN27drh5+dXTYk1msaJUorU1FQSExPp0KFDjfWrcxHUIefOnSMsLMy6iu7tDcHBpeslJ1e976ZNmxJio5VPnz7topTuZ8+ePezZ06hN7Jp6jogQFhZmfZqsKbSCrWPE4dk/IqJ0nawscOVzb9HizwT8aWlpFBQUlFO77sjNzSU3N7fiihqNG3H8LdYEWsHWMwIDjSADR1yZxQYEBBAQYKRzUEpxxjGaQaPRuBWtYOsZItDcScxISgqUlLdfQBnYRnIlJydTXFxcDek0Gk1V0Aq2HhIeXtproLgY0tKq3ldQUBA+Pj5mH8WkpKTUgIRls3jxYm666SamTJnC8uXL3Xovjaa+oxVsPaRJEwgNLV3uiplAROxssadPny4VePDTTz8xc+bMUm3XrVvHI488UqX7TZw4kXfeeYc333yTzz6r9O4tGk2jRLtp1VOaN4fUVPuynBzIzYWqeluFhYWRlJREUVERBQUFpKenE2qjwbdv306fPn1KtRsyZAhDhgxxRXyefPJJbr75ZpfaajSNBT2Draf4+ztXpK6sU3l4eBBh455w6tQpu1nstm3bSEpKYuDAgXTs2JGffvoJgMmTJ7NmjbHr9qJFixg0aBC9evXikksuIdmcTsfHx9O3b19iYmK45JJLUEpx//33c+WVV3LRRRdVSr7w8HDCw8OrPjCNpp6jZ7D1ADd4h5SioKA5J0+eRClFbm4uZ8+eJTAwEDBmsOPHj+f3339n+fLl/Otf/2LNmjXs2rWLmJgYAC677DKuvfZaAB5//HE+//xzYmNjefbZZ9m2bRve3t5kZGTw6quvsnLlSjIzMzl06BBz5sypULb27du7bdwaTV2iFex5gpeXF+Hh4daZ56lTpwgMDKSwsJCUlBQeeughAHr37k1KSgrnzp2joKCAoCBjV+gPP/yQzz77jPz8fE6dOsXTTz+Np6cneXl53H333cTFxdGvXz9uu+02u2QzGs35jDYRnEfYLnZlZmaSl5fHvn376Ny5szVn7JYtW+jVqxe7d+/mwgsvBGD+/Pls2LCB1atXs337dqKjo+nevTt+fn7s2rWLiy++mNmzZ/PGG2+4JFdOTo41MY1G05jQCvY8wsfHh2CbWNzTp0+zbds2jh49Sn5+PmfPnuXxxx/njjvuYOfOnVbzwM6dOxkyZAgBAQF8+eWXrFu3jp49e3Lw4EH8/f2ZOnUqV199tcthhnv37tXJaDSNEq1gzzNsZ7Gpqals3bqVa665hiFDhjBgwABuu+02Bg0aZKdgZ82axRtvvMGAAQPYunUrHTt2xN/fn6eeeoro6Gguuugijh49yj/+8Y+6GpZGUy/R+WDrkL1799KtW7cK66WlwZEj9mUiEBMDVc2sppRi37591kfyVq1aERlZ3l6U7sfyufXr51LKTY2mxnD2m6xOPlg9g20ABAcbwQe2KGWEz1YVx8ADHT6r0bgPrWAbAB4eRvisI8nJhqKtKiEhITQ1M8oUFRWR6hjRoNFoagStYBsIzhLAFBSU3gGhMoiIXeCBs/BZjUZTfbSCbSA0bQqmS6odruQnACN6ytPTEzB2P8jIyHBdOI1G4xStYBsQzpJxZ2ZCfn7V+/L09CwVPltXdOvWrVKLfRpNQ0Mr2AZEs2bGtjKOuDqLbd68uTWLe05ODmfPnq2GdK7j7++Pv79/ndxbo3EnWsE2IGo6Gbe3tzdhYWHW87qcxWo0jRGtYBsYzpJxFxVBerpr/dm6bGVkZFQpGmvs2LE1YrtNSEggISGh2v1oNPUNtypYERkjIvtF5JCIPODk+l0iskdEdojIKhFpZ3MtTkQOmq84d8rZkPDycr6Ft6tmAl9fX2tCF6jc7rNKKUpKSli6dKld6K2rpKSkuH2nBY2mLnBbNi0R8QReB0YCicBGEVmilLLdn3kr0E8plSsic4HngCkiEgo8CvQDFLDZbOviPK1+I4+7P1+herS0G9YDDzxAmzZtmDlzJpmZmbz99ts0adKEffv2kZGRQWFhIU8++SQTJkwgISGB0aNHM3DgQDZv3szSpUsZNmwYmzZtIjw8nIkTJ3L8+HHOnTvH7bffzuzZswFj48Xbb7+db7/9Fl9fX/73v//RokULTp8+zZw5czhy5Ah5eXncf//99OvXj48//phXXnmFgoICBg4cyBtvvGH1dtBoGhrunMEOAA4ppY4opQqAhcAE2wpKqR+VUpb9mn8Doszj0cAKpVSaqVRXAGPcKOt5yZQpU/j8888JDAzEz8+PlStXMnbsWObNm8eWLVv48ccfufvuu60+sgcPHuQf//gHu3fvpl27dnZ9vf/++2zevJlNmzbxyiuvWIMXcnJyGDRoENu3b2fo0KG88847ANx2220MGzaM7du389FHH9GpUyf27t3LZ599xq+//sq2bdvw9PTkk08+qd03RaOpQdyZDzYSOG5znggMLKf+DcD35bQtFTAvIrOB2QBt27atjqznJX369OHMmTOcPHmS5ORkAgMDCQ8P59FHH2XPnj14eHiQlJRkNRu0a9eOQYMGOe3rlVde4euvvwbg+PHjHDx4kLCwMLy9vbn66qsB6Nu3LytWrABg9erVzJ8/HzBcxgICAli1ahWbN2+mf//+AOTl5dm5kmk0DY16kXBbRGZgmAOGVaWdUupt4G0wkr24QbRGz+TJk1m0aBEnT55kzJgxfP/996SmpvLDDz/QunVr2rdvb134KsuV6qeffmLlypWsX78ePz8/hg8fbm3j5eVldQXz9PSkqKioTFmUUsTFxfGf//ynhkep0dQN7jQRJAFtbM6jzDI7ROQK4GFgvFIqvyptNdVnypQpLFy4kC+//JLrrruOs2fPEhoaSlpaGqtXr+bYsWMV9pGZmUlISAh+fn7s27eP3377rcI2I0aMYN68eYCxnfjZs2cZMWIEixYt4oy58VhaWlql7q/R1DTFxcV8//33TJ48uVr9uHMGuxHoIiIdMJTjVOA62woi0gd4CxijlLLdzm8Z8LSIWNbLRwEPulHWOsXZAlRlyM2FPXtKl19wAQQEVK6P7t27k52dTWRkJD169OCqq67ijjvuYOLEiQwcOJALLrigwj7GjBnDm2++Sbdu3YiOji7TjGDLyy+/zOzZs3nvvfcoLCzkkUceYfjw4Tz55JOMGjWKkpISvLy8eP3110vZezUad1FYWMjDDz9MfHw8SUnVn9O5NR+siIwFXgI8gfeVUk+JyBPAJqXUEhFZCfQETppN/lBKjTfbXg88ZJY/pZT6oLx7NeZ8sOX3YWznbUtoKHTs6Fp/x48ft9pcAwMDiY6OrpZ8Gk19p7i4mPT0dFJSUkhISODKK690rOJyPli32mCVUkuBpQ5lj9gcX1FO2/eB990nXeMgIgKOHrUvS0+HwsKqJ+MGI/DgzJkzKKXIzs4mJydHh7FqGh1KKXJyckhJSSEtLY0SJ6GQERERxMbG8sILL7h8Hx3J1cAJCXGejNvVFK/e3t6E2EQy6PBZTWOisLCQU6dOsXv3bvbt20dKSkop5Tp+/HgWL15MYmIizz//fLXuVy+8CDSuY0nG7agHk5OhRYvSYbWVoWXLlqSlpQGQnp5Ofn6+NUG3O9BbxmjciVKKzMxMUlJSyMzMdJr72MfHh/DwcHx8fPjf//5XY/fWCrYR4EzB5udDVpbzHLIV4efnR7NmzcjKygKM8FntZ6xpaOTl5ZGamkpqaiqFhYWlrnt4eBAaGkp4eDj+/v6ICOmuJvUoA61gGwE+PkYqQ1MfWklOdk3BgmGLtSjYlJQUWrduTRNHW4RGU88oLi4mLS2N1NTUMtNvBgQEEB4eTkhIiNvDsPUvppEQEVFawWZkGDNZV57umzVrhq+vL3l5eZSUlJCcnEyrVq1qRFaNpiZRSnH27FlSUlJIT093umDl5eVFWFiY1QxQWdYfX18t2fQiVyMhKMh5Mu7yklSlpqbSu3dvevfuTcuWLYmMjLSeFxYW2qUy/O9//1uphNzDhw/Hmbvc8OHDiY6OtvZ/7bXXVmpcGk1ZFBQUcPLkSXbt2sX+/ftJTU21U64iQnBwMJ07dyYmJoaoqKhKK9fDaYf56xd/Zcj7Q6olo57BNhJEDFvsiRP25Skp0KqVsRjmSFhYGNu2bQPgscceIyAggHvuucd6PTQ0lKSkJAoLC/nkk0+48cYbq+UX+8knn5S7kFVUVGRnhnA8r2w7TeOlpKTEbsHKGb6+voSHhxMaGopXFX0VS1QJdy27i9c2vEZhSWm7bVXR38pGRHg4nDxpv5V3YaFhKggNrVwfq1at4p577qGoqIj+/fvzyCOP8Nprr5GcnMy4ceOIjIzkxx9/ZO7cuWzcuJG8vDyuvfZaHn/8cZdknjVrFtnZ2ezfv59Ro0aRlpaGj48PW7du5eKLLyY2NpY5c+aQm5tLp06deP/99wkJCWH48OH07t2btWvXMm3aNO6++26X7q9pGOTl5ZGSkkJqaqrTfBaenp7WBSs/Pz9r/ovKUqJKOJNzhqSsJF787cWaElsr2PpCVb8QVWHfPlUpBXvu3DlmzZrFqlWr6Nq1K7GxsdYcBZ9++ilvvPGGNdPVU089RWhoKMXFxYwYMYIdO3YQExNTbv/Tp0/H19cXgJEjR1p9DLOzs1m5ciUtW7Zk1qxZJCYmsm7dOjw9PYmJieHVV19l2LBhPPLIIzz++OO89NJLgPGI2NCi9zSVp6ioyBphleMYrmhiyQAXHBzs0oKVUor0c+kkZSWRX5xPiXJh76Vy0Ar2PCA7G/LywNRtZVJcXEyHDh3o2rUrAHFxcbz++ut2CS9Onz5NUFAQn3/+OW+//TZFRUWcPHmSPXv2VKhgyzIRTJ8+nZYtW1rPJ0+ejKenJ5mZmWRkZDBs2DCrPLayTJkypcKxaxoWlghCy4KVM59Vy15y4eHh1fLPzs7PJjErkZxC58q7JtAK9jwhORlcdWW1XezKyspiz549vPDCC2zcuJGQkBBmzZpVpb28HHEMxa1saK4O4W08FBQUWE0A+U72obcsWIWHh9OsWbNqPfGdKzpHYlYiGecyyq3n5+XHvUPu5XFcM3+BVrD1hppMunP6NBw/bl+WmgqRkVDeU5SnpycJCQkcOnSIzp0789FHHzFs2DCaNm1KYGAgOTk5BAcHk5CQgL+/P0FBQZw+fZrvv/+e4cOHuyxvVlYWycnJNHfYMjcoKIiQkBDWrFnDpZdeapVH0zgoKSkhIyODlJQUq8+1I5YFq7CwsGovZBYWF3Ly7EmSc5JRlP178xAPru99PU9c9gStAltpBauxJywMkpLst/IuLoa0NOfbflvw8fHhgw8+YPLkydZFrjlz5gBw0003cdttt9G8eXPeeustevXqxQUXXECbNm24+OKLKyWXrQ02PDyclStXApCcnMyxY8dKKViA+Ph46yJXx44d+eCDcpOqaRoAubm51iQrZS1YWUwAfn5+1b6fZQHrZPZJilVxuXV9vXzZPmc7PSJ6VPu+4OZ0hbXJ+ZqusCwSEkr7wPr5wYUXut7nvn37rL6wLVu2JCoqqoIWlUPnImj8FBUVkZaWRkpKCrm5uU7rNGvWzLpg5eHMr7CKKKVIy0sjKTuJguKCcuv6NvGlTbM2JB1NKvWbFJH6ma5QU3c0b15awebmGrljXTVdtmzZkkOHDgFYI7v0jq+asqjsgpXFBFCTCYWy87M5nnWc3ELnytyCl4cXkc0iCfMNQ0RIquGNU7SCbaT4+xszVsfJQnKy6wo2KCgIHx8fzp07R3FxMSkpKXYLYBoNQH5+vnXBqqCg9MxRRAgJCSE8PJzAwMAadVHMK8wjKTupwgUsD/GgZUBLWvi3wNPDfZMErWAbMRERhqnAlrQ0iIoqnUO2MogILVq0sO6Tdfr0aZo3b14jj3Oahk1JSYnVZzU7O9tpHT8/P2uEVU1H3hUWF3Ii+wTJuckV1m3u15zWga3x8nQhI30V0Qq2ERMSYngTFNvY9UtKDI8CVyeeYWFhJCUlUVRUREFBAenp6YSFhdWMwJoGhVLKbsGquLj0AlKTJk0ICwsjLCysRhasHCkuKTYWsM6erDBIIKhpEFHNovD1qsAhvAbRCrYR4+lphM+aW2xZSU42ZreuPJl5eHgQERHBCTPpwenTpwkNDXVrJJqmflFQUGDNs1qW/3NQUBDh4eEEBQW55QlHKUVqXionsk9UuIDl5+VHVLMomjVtVuNyVIRWsI2c5s1LK9hz54zormYuft+aN2/OyZMnrTOY7OxsmrnaGdp7oCFQXFxMRkYGqampZfqsNm3a1Lpg5e0stVsNkZWfRWJWYoULWN6e3kQGRhLqW3cTAG08a+T4+EBgYOnyZNNUlZiYyIQJE+jSpQudOnXi9ttvd7owYUtOTg7ff/+99XzHjh01ln7wsccec7rJ3GOPPYaIWL0YAF566SVEpEr5CD788ENuueUWl+t8//339OvXjwsvvJA+ffq4PcnMrFmzWLRokUttx44dS0ZGBhkZGbzxxhtVbn/ixAnGjBnD0aNHee+992jXrh3jxo1j0qRJzJ49mzVr1uDh4UF4eDjR0dH06NGDVq1auU255hXmcTD1IAdSD5SrXD3Eg8jASNZ+uZbhA4fTp08fLrnkEvaYe9zv3LmTWbNmuUXGUrLUyl00dYqz4IL0dMjPV1xzzTVMnDiRgwcPcuDAAc6ePcvDDz9cbn8ZGRksWLDAeu7j48NHH31U02KXomfPnixcuNB6/sUXX9C9e3e339fCrl27uOWWW/j444/Zs2cPmzZtonPnzrV2/6qydOlSgoODq6xgz507R1JSEg899BAjR44kNTUVpRR9+vThk08+4csvv+TRRx/lxRdfJCUlhfbt21fLG0Ap5TRJtoXC4kKOZRxjd/JuMvOdpyi0EOEfQc+InrQKbMWMGTPYuXMn27Zt47777uOuu+4CjO9RYmIif/zxh0vyVgWtYM8DgoOdb+G9ePFqfHx8+Nvf/gYYETQvvvgi77//Prm5uXz44YdMmDCB4cOH06VLF2tKwgceeIAjR44wc+ZMXn75ZU6cOEHv3r0BY/Y3ceJERo4cSfv27Xnttdf473//S58+fRg0aJB1M8V33nmH/v3706tXL0aOHMnmzZsrHMfEiROtG9IdPnzYauezsGDBAnr27EmPHj24//77reUffPABXbt2ZcCAAfz666/W8uTkZCZNmkT//v3p37+/3TVnPPfcczz88MNccMEF1vdr7ty5ACQkJHD55ZcTExPDiBEjrD/eWbNmMXfuXAYNGkTHjh356aefuP766+nWrZvdLCogIIA777yT7t27M2LECJKTS6+Gb968mWHDhtG3b19Gjx7NyZMnyczMJDo6mv379wMwbdo03nnnHQDat29PSkoKDzzwAIcPH6Z3797ce++9xMbGsnjxYmu/06dP56uvviI5OZm9e/eya9cuTp48yYoVKxg8eLC1noeHB5GRkcTExDBhwgQee+wx5s2bV+57mZyczMiRI+nevTs33ngj7dq1IyUlhYSEBKKjo4mNjaVHjx4cP36c559/nv79+xMTE8Ojjz5KcUkxJ7JP8Oy8Zxl72ViuG3kdT9/3tNPFtGCfYHo070HboLZW7wBbs1VOTo7dP4Bx48bZ/bN2F25VsCIyRkT2i8ghEXnAyfWhIrJFRIpE5FqHa8Uiss18LXGnnPUFESnz9fbbb1vrvf322+XWdcSy86wjW7bs5qKL+tqVNWvWjLZt21ofxTds2MCXX37Jjh07+OKLL9i0aRPPPPMMnTp1YuPGjdx+++2AEalj2Vhu165dfPXVV2zcuJGHH34YPz8/tm7dyuDBg5k/fz4A11xzDRs3bmT79u20bdu2Ul/2Zs2a0aZNG3bt2sXChQvtsmmdOHGC+++/n9WrV7Nt2zY2btzI4sWLOXnyJI8++ii//vora9eutT4mAtx+++3ceeedbNy4kS+//JIbb7yx3Pvv2rWLvn37Or126623EhcXx44dO5g+fTq33Xab9Vp6ejrr16/nxRdfZPz48dx5553s3r3bOrsCQwH069eP3bt3M2zYsFL5dQsLC7n11ltZtGgRmzdv5vrrr+fhhx8mKCiI1157jVmzZrFw4ULS09O56aab7NpaPq9t27bx/PPPc8MNN/Dhhx+ilOLYsWP88ssvREZGcuzYMWtawKSkJJo1a4afnx8RERG0bduWwMBAOxPARRddxL59+8p9Lx9//HEuv/xydu/ezbXXXms3azx48CD/+Mc/2L17N/v37+fgwYNs2LCBrVu3sn7DeuL/F8+vW35l2f+W8d7i9/h0xad4eHrww1c/APDgnAeZMWoG1195Pddefi2D+g+id+/e1u8YwOuvv06nTp247777eOWVV6zl/fr1Y82aNeV+3jWB2xa5RMQTeB0YCSQCG0VkiVJqj021P4BZwD2leyBPKdXbXfKdbzRvbiTjtqW42NizqzxGjhxpdcO65pprWLt2LRMnTgSMbFa2Ga3OnDkDwGWXXUZgYCCBgYEEBQUxbtw4wHg027FjB2Aoq3/+859kZGSQlpbGoEGDKjWOqVOnsnDhQpYtW8aqVausuQk2btzI8OHDrfkMpk+fzi+//AJgVz5lyhQOHDgAwMqVK+0UblZWVqW2xXHG+vXr+eqrrwCYOXMm9913n/XauHHjEBF69uxJixYt6NmzJwDdu3cnISGB3r174+HhYf2HMWPGDK655hq7/vfv38+uXbsYOXIkYCw6WfZIGzlyJF988QU333wz27dvr1DW/v37s2fPHn766SeWL1/OsGHD7CLyRIT8/Hxat25NTEwMHh4eHDlypFQ/tpFZZb2Xa9eu5euvvwZgzJgxhISEWOu0a9fO+rkvX76c5cuXE9M7hsLiQnJycjh65Cj5e/LZt3MfsWNjAcg/l09oeCjent4sXLiwwgWsm2++mZtvvplPP/2UJ598kvj4eAA7Txh34k4vggHAIaXUEQARWQhMAKyfglIqwbxWs1luGyiVzQsxe/ZsZs+eXaW+vb0NU0FGxp9lHTpcyIcf2i+gZGVl8ccff9C5c2e2bNlS6strey4itGzZkoMHDwKGgi0pKbELefTw8LCee3h4WJN7zJo1i8WLF9OrVy8effTRSpkIAK6++mruvfde+vXrVy3PBTCc43/77bdK79PUvXt3Nm/eTK9evap0H9vxO743zpKdQOkE7Eopunfvzvr1pTfhKykpYe/evfj5+ZGenu40R4RSitOnT5OSkkJeXh6jR4/mm2++Yfny5TzyyCOAEQgQFhZGaGgoHh4eFBcXl+titXXrVmvcflXfS7BPN1lQXMD1t17P2Glj7ep89v5nXDX5Km550Fh09BRPWgW2IsI/gmlTp1lNI7bcddddxMbG2pVNnTrVas4Bw87sW1GC5BrAnSaCSMA2aV6iWVZZfERkk4j8JiITa1Sy8xTHxa4BA0aQk5PLe+8Zj1TFxcXcfffdzJo1y+oUvmLFCtLS0sjLy2Px4sVcfPHFBAYGWqN1goODrY+MxcXFlZ4BZmdn06pVKwoLC/nhhx8qPQY/Pz+effbZUgtxAwYM4OeffyYlJYXi4mIWLFjAsGHDGDhwID///DOpqakUFhbyxRdfWNuMGjWKV1991XpueVwvi3vvvZenn37aOgMuKSnhzTffBGDIkCFWM8cnn3zCpZdeWukxWfqyeAt8+umnXHLJJXbXo6OjSU5OtirYwsJCdu/eDcCLL75It27d+PTTT/nb3/5mNdWAYZ44c+YMaWlpHD9+nLy8PMD4R7VgwQJEhCFDhtC9e3cuvPBCWrRogZeXF127diXBMQzQhh07dvDvf/+bm2++GSj7vbz44ov5/PPPAWOWmp6ebtdPQXEBCRkJdB3Qlc8/+ZzcHMM74MzJM6SlpNH/kv6s/nY16SnpRPhH0LpJa/JT8/EQDz777DO2bdtW6mVRrpZ//ADfffcdXbp0sZ4fOHCAHj1qJmNWedRnP9h2SqkkEekIrBaRnUqpw7YVRGQ2MBugravZpM8jmjUztvC2mAVEhOee+5qXXvoHzzzzb0pKShg7dixPP/20tc2AAQOYNGkSiYmJzJgxw+qzevHFF9OjRw+uvPJKO1toVlZWpWbi//73vxk4cCDNmzenQ4cOZW4J4oypU6eWKmvVqhXPPPMMl112GUoprrrqKiZMmAAYLl6DBw8mODjYuhgH8Morr3DzzTcTExNDUVERQ4cOtSpMZ8TExPDSSy8xbdo0cnNzERGuvvpqAF599VX+9re/8fzzz9O8efMqp1X09/dnw4YNPPnkk0RERPDZZ5/ZXff29mbRokXcdtttZGZmUlRUxB133EGTJk1499132bBhA4GBgVx66aU88sgjzJkzh4KCAo4dO0ZwcDC9evViypQpDBkyhDvuuIMuXbrQvXt3Jk2a5HTG6+/vT6dOnay5gQHWrFlDnz59yM3NJSIigldeeYURI0aU+14++uijTJs2jY8++ojBgwfTsmVLAgMDyczKpKikiF1ndlGiShg0bBBHDx7l+vHXA8Y/0idefYKOXTtyzz/v4Z6Z96CUwsvLi9dff5127dpV+J6+9tprrFy5Ei8vL0JCQqzmAYAff/yRq666qkqfkSu4LV2hiAwGHlNKjTbPHwRQSv3HSd0PgW+VUk4d/iq6DjpdYWU5dQoSE+3LmjSBmJjSO89++OGHbNq0iddee63cPouLi9mxY4d1dbdTp052traK0OkKDS8CV+2/YCRYSUtLKze6KiAggLCwMEJCQigoKKBnz55s2bKFoKAgp/W//vprNm/ezJNPPlktuTw9PWnSpAnr169n7ty5rPh1BSeyT1S4a6u/lz9tmrUhoGmAy/cvS6Zhw4axdu3aUjkRnP0m62u6wo1AFxHpACQBU4HrKtNQREKAXKVUvoiEAxcDz7lN0vOI8HAjGbft/9WiIiMJjDNPg8rg6elJREQEJ81VtFOnThEcHFxpv8hwV298nlNcXEx6ejqpqallJlhp2rSpNReAxf67cuVKbrjhBu68884ylSvAX/7yF1JTU6sl4x9//MFf//pXSkpK8GziyQPPPsCxzGPltmnq2ZTIZpGE+IS4JQLrjz/+4JlnnqmVrd7dmnBbRMYCLwGewPtKqadE5Algk1JqiYj0B74GQoBzwCmlVHcRGQK8BZRg2IlfUkq9V9699Ay28hw9aiR8scXfH6ojSmFhITt27LCaBy644AICAmp25qH5M8dqamoq6enpTh30PTw8CA0NJSwsjICAgDrPE5FbmMvxzONkFzj/J2DBdgHLQ+rGRb8hzWBRSi0FljqUPWJzvBEoZQBSSq0DerpTtvOZ5s1LK9icHCN3rKsJj7y8vAgLCyPFzPJ96tSpeh3l1NA4d+6cNWtVWaHMzZo1IywszOUtrGuaguICkrKSSM0rfxYsCBH+EbQKbEUTj/q8LFR1GtdoGiBKqVqfYZSXjLsSawdl0qJFC6uCzcjI4Ny5c5Vy27EscOldYu2xbLOSmppa5iKgr6+v1bXKnQlWqkJxSTGnzp7iVM6pChc8Q3xCiGwWiU+Tyrt3uQt3PM1XqGBFxAMYZM4qNTWIj48PqamphIWF1aqSFTFmscccTGGWnWddNU35+voSFBREZqYRL3769OlKrfbu3bsXOL8XuSyUlJSQlZVFamoqGRkZTn/0TZo0sZoA/Pz86twEYEEpRXJuMieyT1BU4ty/10KAdwBRzaII8K4fZiSlFKmpqVXy460MFf6UlFIlIvI60KdG76whKiqKxMREp3Hn7saSeNvx97tli/PsW5XF8igLWBdfKnpctdS3KNrzkYKCAs6ePUtOTk6ZiU/8/Pzw9/fH29ubnJycKrm2uZvcwlwyzmVQWFy+Z0ATjyaE+IagvBTHk4+XW7e28fHxqbGNPC1Udq6ySkQmAV+pxrINbT3Ay8uLDh061Nn9582D11+3L7vwQti1y7Vk3GDMBPr168eWLVsAIxbdEilUFheaW92eb1+tkydP8sknnzB//nx27tzptE7//v2Ji4tj6tSp9W7niKKSItb+sZbHf36cnxJ+KrduqG8ojw57lDn95uDtWT9MGbVBpbwIRCQb8AeKgTxAAKWUqv0U4WXQEL0I6prdu8FZMMvPP8PQoa73u2DBAq67zvDIa968OceOHSs3LNHyiHs+KNi8vDyWLFlCfHw8y5YtczpbjYyMZObMmcTGxtaJl0l5ZOdns+zwMpbsX8LSg0srXMDy9vTm9oG389ClDxHsE1w7QtYw1fEicKubVm2iFaxrDB0KjkmFpk4Fm3SvVaawsJDOnTtbMye99dZb5eZOaOwKVinFunXriI+P5/PPP7faqG3x9fVl0qRJxMbGcvnll9cLLwALxzKO8c2Bb/jmwDf8ePTHCgMELEzrMY2nRzxN++D27hXQzdSKghWR8YBlXvOTUupbV27oLrSCdY0FC+A6h/APLy9js8Tq7Mj94osvWhMcR0dHs2fPnjIThzRWBZuQkMD8+fOZP38+hw8fdlpn+PDhxMbGcu211xJYHeN3DVKiSth8YjNL9i9hyYEl7Di9o0rtL217KS+MeoEBkQPcJGHt4nYFKyLPAP2BT8yiaRjBAg+6clN3oBWsa+TnQ9u2YGYatPL00/BgNT7drKws2rRpY92/acmSJda0hY40JgWblZXFokWLmD9/Pj///LPTOp07dyY2NpaZM2fSvn372hWwDHILc1l1ZJV1pnrq7Kkq99E1rCvPXfEc46PH1xvPhpqgNhTsDqC3Usa+uGau161KqRhXbuoOtIJ1nYcegv84ZIho1w4OHzZ2pnWV++67j+effx6AoUOHlqlwLKkKy0pmXd8pLi5m9erVxMfH89VXX1kzVtkSFBTElClTiIuLY/DgwfVCAZ06e4pvD3zLkv1LWHlkJXlFpeWuCD8vP0Z2HMmkbpOY2mOqdTeBxkRtKdjhSqk08zwUw0ygFWwjICEBOnYs7bL17bdQnYRDiYmJdOjQwZrz9Pfff2fAgMbx2AiGW1l8fDwff/wxSUlJpa57enoyevRo4uLiGDduXK3kHy0PpRQ7z+xkyf4lfHPgGzYkbXCpn9aBrRnXdRzjuo7j8g6X4+tVt+NyN7URKvs0sFVEfsTwIBgKlNoCRtMwad8exo6F776zL3/jjeop2KioKGuqOoD/+7//K5WGr6GRmprKggULmD9/Phs3bnRap2fPnsTFxTF9+nRatmxZyxLak1+Uz8/Hfuab/d+w5MAS/sh0baO/Pi37MD56POO6juOiVhfVixl4Q6DCGawZyXUtsAbDDguwQSlVdSONG9Ez2Orx3Xdgpja1ImKYCarjqrt9+3ZrDlYPDw8OHTpUyvfX4mFgu+9YfaKgoIClS5cyf/58vv32W7uE1haaN2/O9OnTiYuLo1evXnWqgFJzU1l6cClLDixh2aFlFSZZcYa3pzeXd7ic8V3Hc3XXq2kT1MYNkjYMasNEsMnVG9QWWsFWj+Ji6NSpdPjsAw+Uts9WlVGjRrFixQrA2BzvpZdesrteHxe5lFJs3ryZ+fPns2DBAmu0mS3e3t6MHz+euLg4Ro8ejZezrXtrif0p+62P/r8e/5USVfVdmML9wrm669WM6zqOkR1HEti0fng11DW15UWQAnwGWOPzLDbZ+oBWsNXnP/8xFrxsad7ccNmy2UqqyixbtowxY8YARkKX48eP2yXkrk8K9sSJE3z88cfEx8fbbeJny6BBg4iNjWXKlCmEhobWsoQGRSVF/PrHr3xz4BuW7F/CwbSDFTdywoXNL2Rc13GMjx7PwMiBeHrUH//b+kJtKNijToqVUqqjKzd1B1rBVp8zZyAqChyfgD/9FKZNc71fpRQxMTHs2rULgP/85z888MCfJvy6VrC5ubksXryY+Ph4Vq5c6TS6qk2bNtboqujo6DqQEjLPZdpFUaWfS6+4kQOe4snQdkOt9tROoZ3cIGnjwq0K1rTBTlZK1evVCa1ga4brrisdxXXppWDugO0y8fHxzJo1CzD2zzp69Kg1w35dKNiSkhLWrl1LfHw8X3zxhdMdAfz9/Zk0aRJxcXEMHz683B1W3UVCRoJ1geqnhJ8qzFLljKCmQYztMpZxXccxpvMYQnwrv52PRttgAa1ga4o1a5znIdi503negspSUFBAhw4drHvRf/DBB1aFW5sK9vDhw3z00UfMnz+fo0dLP5iJCJdddhlxcXFcc801tb4rQ4kqYUPSBr7Zbzj87zzjPAlMRXQM6cj4ruMZHz2eS9pe0ij9U2sLbYNFK9iaQino2dNIBGPLzTdDBXsfVsizzz5rNQ306NGDHTt2ICJuV7CZmZl88cUXxMfHs3btWqd1unbtSlxcHDNmzKj1HYpzCnJYeWQl3xz4hm8PfMvpnNNV7kMQBrcZzPiu4xkXPY5u4d20K1UNoW2waAVbk7z+Otxyi31ZYCCcOAHVmdBlZGTQpk0b6+6pP/zwA6NHj7ZGcFkiumqC4uJiVqxYQXx8PIsXL3a602pwcDBTp04lLi6OgQMH1qpCOpF9whpFteroKs4VOd8Jtjz8vfwZ1WkU46PHM7bLWCL8I9wgqUZn00Ir2JokKwtatzb26bLlzTfh73+vXt933HEHL7/8MgBXXHGF1X2rpti1axfz58/n448/tu5ya4unpydXXnklcXFxXH311TWewb4slFJsP73d6kq16YRr39XIwEjrAtVlHS6rF1utNHbcpmBF5D6l1HPm8WSl1Bc2155WSj1UZuNaRivYmuXvfwdHv/9evWDrVteTcYORYapTp07WlfqtW7daAxFcJTk5mQULFhAfH29N9O1I7969iYuLY9q0abSoTpqwKpBflM+PCT9a7anHs1zL4N+3VV+rK1Xvlr31o38t404Fu0UpdZHjsbPzukYr2Jpl2zbo42SToHXrYPDg6vU9depUa8jsjBkzrKG0VSE/P5/vvvuO+Ph4li5das13YEuLFi2YMWMGsbGxxMTUTtqM5JxkaxTV8sPLOVtwtsp9NPVsyoiOI6xRVJHNIt0gqaayuFPBblVK9XE8dnZe12gFW/MMGQLr19uXzZwJ8+dXr99NmzbRv78Rdd2kSROrcqzIXKWUYuPGjcTHx7Nw4ULS0kqvsTZt2pQJEyYQFxfHqFGjaOLqDo6VRCnFvpR91kf/dcfXoai62S3CP4Kru1zNuGgjisrfW++wW19wZ7IXVcaxs3NNI2Pu3NIK9vPP4cUXoTrbQ/Xr149hw4bx888/O515OpKYmGh1rdq3b5/TOkOGDCEuLo6//vWvBAcHuy5cJSgsLmTtH2utUVSH050n066IHhE9rI/+AyIH4CG172ercS8VKdheIpKFkUHL1zzGPK/Qui4iY4CXAU/gXaXUMw7XhwIvATHAVKXUIptrccA/zdMnlVLxFQ9HU5NMngx33mnsPmshPx8++ADuuad6fd9zzz1l5ocFyMnJ4auvvmL+/PmsWrXK6ey2Xbt2xMbGEhsbS+fOnasnUAVknMvgh0M/sGT/Er4/9D0Z5zKq3EcTjyYMazfMSPUXPY6OIfXGCUfjJtzmRWAm5T4AjAQSgY3ANKXUHps67YFmwD3AEouCNfPNbgL6YcyUNwN9lVJlxgZqE4F7uPdeeOEF+7JOneDAAahOYFNJSQkXXngh+/fvt5YVFxfzyy+/EB8fz6JFi6zuXLYEBAQwefJkYmNjGTp0qFujq46kH7E++v9y7BeXoqiCfYIZ22Us47uOZ0znMQT5BLlBUo07qY18sK4wADiklDoCICILgQmAVcEqpRLMa47B36OBFTYJvlcAY4BqbMWncYW//720gj18GFasgNGjXe/Xw8ODu+++224zxI4dO3LMMZ0XRnTViBEjiIuL4y9/+Qv+/u6xTxaXFPN70u/W0NQ9yc6TvVRE59DOVof/i9tcrKOozmPcqWAjAVu/lERgYDXa6qXUOqBzZxg1CpYvty+fN696ChZg5syZ/POf/+SMuSGYo3K94IILrNFVUVFR1btZGZwtOMuKwytYcmAJ3x34juTc5Cr34SEeDGkzxGpPjQ6L1q5UGsC9CtbtiMhsYDZQ6+GN5xNz55ZWsN98Y6QxbFONPMw+Pj7ceuut/Otf/7KWhYaGMm3aNOLi4ujXr59bFFViVqLVN3X10dXkF+dXuY8A7wBGdxptjaIK9wuvcTk1DR93KtgkwPbnF2WWVbbtcIe2PzlWUkq9DbwNhg3WFSE1FXP11UYaw8TEP8tKSuCdd+CJJ6rX9/3338+aNWvIysri3nvv5aqrrrJm2aoplFJsPbXV2IZ6/xK2ntrqUj9tmrWxRlENbz+cpk1qVk5N48Odi1xNMBa5RmAozI3AdUqp3U7qfgh867DItRmwBDJswVjkKjO5jF7kci///jc88oh9WatWxg4IdZjIv0zOFZ1j9dHV1plqUnZl/7fb0791f+ujf0yLGP3ofx5SLxe5lFJFInILsAzDTet9pdRuEXkC2KSUWiIi/YGvgRBgnIg8rpTqrpRKE5F/YyhlgCfqU+au85EbbzRmq7ZuqydPwv/+B9deW3dy2XIm5wzfHfiOJQeWsOLwCnIKcypu5IBPEx+u6HgF47uO56quV9E6sLUbJNWcL+hkL5pKM3kyLFpkX3b55bBqVfX6tWx2aOtRUBmUUuxJ3mN1pfot8TeXoqha+Lew+qZe0fEK/Lz8qtyHpvGis2mhFWxtsHo1jBhRunzvXrjgAtf7rUo+2MLiQn459os1iupohrNMmhXTM6Kn1Z7aP7K/jqLSlEm9NBFoGh+XXQbR0WATGwAYaQwdNoqtUdLz0vn+0Pcs2b+EHw79QGZ+ZpX78PLwYnj74daZavvg9jUvqEbjgFawmkojAnPmGOGztsTHw9NPg18NPlkfSjtkffRfc2wNxaq4yn2E+oZao6hGdx5Ns6bNak5AjaYSaAWrqRJxccbW3nl5f5ZlZMDChXD99dXr+9c/fjVcqQ4sYV+K86QuFdE1rKs1impImyE08dBfcU3dob99mioREgJTpxoJX2yZN69qCrZElbA3eS/rjq+zll3ywSVVlsdDPLik7SXGo3/XcUSH182W2hqNM7SC1VSZf/yjtILdtMl49StjKSArP4vfE39n3fF1rE9cz2+Jv7lkSwUI9A5kTOcxjI8ez5WdryTMrxq5EzUaN6IVrKbK9OtnvBydNubNg/feM7wBDqUdsirTdcfXsevMLpdcqCy0C2pnXfUf1n4Y3p7e1RyFRuN+tILVuMTcuXDDDeaJVy603sj8I+s5MX8dm06vJyU3pfKdPea8eEDkAMZ3Hc/46PH0iOiho6g0DQ6tYDVVQinFH5l/4NlrPd4T1lEQsR5abgPPIoqAH1xzSwXAt4kvIzuNZFzXcVzV5SpaBbaqKbE1mjpBK1hNueQX5bP11Fa7x/0T2SeMi9XckS3AO4CBkQMZ0mYIQ9oMYVi7Yfh6+VZfaI2mnqAVrMaOk9knWZ+4nvXH17MucR2bT2x2KZ2fMzqFdLIq08FRg+kR0QNPD0/69u3Ld3zH5s2ba+Q+Gk19QSvY85iikiJ2nN5hVabrjq8jISOhRvr2aeJD/9b9rcp0cJvBRPhHOK27ZcuWGrmnRlPf0Ar2PCI1N5XfEn+zPu7/nvQ7uYW5NdN5Zhs4PgSOD8bzxBD2retFuyi90q85v9EKtpFi68hvsZ3uT91fccNK4OXhRZ9WfRgSNYQBrYdwx6TBnDn055YuxcD8D8BmowKN5rxEK9hGgsWR36JMq+PI70iEf4RhO40awuA2g+nbqq/dYtTeaUZCblvefhsefBCa6G+Y5jxGf/0bIBZHfosyXZ+4np2nd1bLkd+Ch3gQ0yLGqkyHtBlCh+AO5fqg3nQTPPWUsY2MhcRE+O47mDCh2iJpNA0WrWAbALmFuWw6sYl1x9dZFWqVHPnLIcQnhMFtBjM4ylCm/Vv3J7BpYJX6aNMGxo0zdjew5Y03tILVnN9oBVvPUEpxPOu4oUjN1f1tp7ZRVFJUceNKcGHzC63KdHDUYKLDo2sk2fTcuaUV7PLlcOiQsfV3edx0003Vvr9GUx/RCraOsTjy27pKWR35q4mtI//gqMEMihpEiG9IjfTtyMiR0KkTHD5sX/7WW/D88+W3tWwZo9E0NrSCrWVOnT1lKNPj69zmyG+ZoVoc+WsDDw/4+9/hvvvsyz/4wFgA8/GpFTE0mnqFVrBupKikiJ2nd9q5Srm6h5QjFkd+6+N+OY78tcXf/ma4ZuXb/L9ITYUvvoCZM8tuZ4ng6tu3r5sl1GhqF61gaxBHR/4NSRtc2jraGW2atTFW9aOMUNNeLXvVu5R94eHw17/CRx/Zl8+bV76C7WcmkW0sG3BqNBa0gnURiyO/ZWbqLkd+ywp/m6A2NdK3u5k7t7SCXb8etm+HXr3qRiaNpq7QCraSZOVnsSFpg1WZusOR3/K47+jI35AYNMhQpNu325fPm2fsPqvRnE+4VcGKyBjgZcATeFcp9YzD9abAfKAvkApMUUoliEh7YC9gmRL+ppSa405ZbVFKcTj9sJ3faU078luUaWUc+RsSIsYsdo7Dp/Xxx/Dcc9BMb+yqOY9wm4IVEU/gdWAkkAhsFJElSqk9NtVuANKVUp1FZCrwLDDFvHZYKdXbXfLZYnHkt7hKrT++nuTc5BrpO8QnhEFRg6wz1AGRA6rsyN/QmD4d7r0XsrP/LMvJMUwHN99cd3JpNLWNO2ewA4BDSqkjACKyEJgA2CrYCfy5Ycgi4DVx81TO3Y783cK72T3u15Qjf0MiIMBY1HrjDfvyefOMDRMbyWRdo6kQdyrYSOC4zXkiMLCsOkqpIhHJBCxbhHYQka1AFvBPpdQaxxuIyGxgNkDbtm2dClFQXMDWk1utfqfrj68nKTupGsP6E4sjv0WZDowaSKhvaI303dCZO7e0gt29G9auhUsvrRuZNJrapr4ucp0E2iqlUkWkL7BYRLorpbJsKyml3gbeBujXr58Ce0f+9Ynr2XRiU4068tu6StWmI39Do0cPQ5Gucfi3OG9eaQW7yXF7Wo2mkeBOBZsE2PoWRZllzuokikgTIAhIVYZDZD6AUmqziBwGugJl/hKPZhyl48sda9SRv1/rflZlOihqEC0CWtRI3+cLc+eWVrCLFsFLL0GETUyEDjDQNFbcqWA3Al1EpAOGIp0KXOdQZwkQB6wHrgVWK6WUiDQH0pRSxSLSEegCHCnvZmm5aaRlpLksbFSzKDvbae+WveudI39D45proHlzSLZZLywshPffhwceqDu5NJrawm0K1rSp3gIsw3DTel8ptVtEngA2KaWWAO8BH4nIISANQwkDDAWeEJFCoASYo5RyXXs60MSjCRe1uqhBOvI3JJo2hRtugGeesS9/6y3Dy8DTtK7Mnj0b0ElfNI0PaSzhidJaFH93fi3CP8LO77QhO/I3NBISoGNHcPyaffstXHWVcWxxHGks30VN40JENiul+rnStr4ucrmMh3jQM6Kn3fbQHUM6NhpH/oZG+/Zw5ZWwdKl9+bx5fypYjaax0mhmsG26tVEffvfheeHI39D49ltjxwNbRODIEUMB6xmspj5TnRlso/GAb+HfghEdR2jlWg+58kpo186+TCljY0SNpjHTaBSspv7i6Wkk43bkvfegoKD25dFoagutYDW1wg03gJeXfdmZM/DVV3Ujj0ZTGzS6RS5N/SQiAiZNgoUL7cvnzYOLLrqoboTSlEl+PqSnG6+0tD+PnZVlZ0OLFkb0Xo8e0L27sT+bpw5ybDyLXP369VM65LJ+88svMGxY6fJdu4wfpaZmKSwsWymWVWY5z8ur3r19fKBbN+NztVW8bdsa+7c1JLSblqZBcOmlxo9s92778nnz4LXX6kam+k5REWRkuKYoc2pmtyKXOHcOtm41XrYEBPypdG2Vb8uWjTPLmp7BamqV116DW2+1LwsMhBMnjB9fY6S4GDIzK1aKzsqysiruvzEQEmI/07Uch4VV3NbdVGcGqxWsplbJzITWrSE317a0/vvBlpQYtsbKKEXHsszM0pFsmspha9u1KN/u3Wt3ZwxtItA0GIKCYMYM5z6wSrn3MVEpOHu28rNH27KMDEPJni94ehqzSsdXaGjpMj8/I2hk164/X2k1lDnk9GnjtWqVfXnbtqXtu926GbLUJ/QMVlPrbNsGffrYlhhadf16xaBB5bdVypj9VuUx21KWkWHYNM8XRCA42LlSdCxzPA8MdP2fnVKGUrQo2927/zw+e7ZGh2iHiOG94Kh4o6PBuxqJ8bSJAK1gGxqDB8Nvv1nOjF/yyJGKsWMrXuEuLKwzseuEoKCyZ4/lKcpmzerXir1ScPx4acW7Z4+xKOYumjSBrl1L23cr60qmFSxawTY05s+HuDjLmWWq1Di+i84IDCz/MbssRRkU1Pj9SYuLDROD7Ux3927Yt8+9TxxNmxpmBUfF6+hKphUsWsE2NM6dg8hIi62uYShYf/+K7ZHOFGVwsDGL0lSNggI4eLC04j10yL32cH9/e4V7111awWoF2wC55x74v/+D2lSwvr6VW7hxVlYdO56m5sjLM2a3top31y44dsxdd9ReBJoGyJ13wrvvQmbmW1Vq5+1deaXoWObj46bBaGoNX19jkdR+odRwo9uzp/TC2smTdSMn6Bmspo7ZtcuYxZ45YzxKV0ZR+vo2zqgfjXtISyttZti5syquZNpEoBWsRqOpNBZXMkfFu2uXMRO2R5sINA0Yy2aHls0PNRp3I2LkP2jZEkaM+LPc1pXMonDnz6/GffQMVlPX6C1jNPUZvWWMRqPR1EO0gtVoNBo34VYFKyJjRGS/iBwSkQecXG8qIp+Z138XkfY21x40y/eLyGh3yqnRaDTuwG0KVkQ8gdeBK4ELgWkicqFDtRuAdKVUZ+BF4Fmz7YXAVKA7MAZ4w+xPo9FoGgzunMEOAA4ppY4opQqAhcAEhzoTgHjzeBEwQowVjwnAQqVUvlLqKHDI7E+j0WgaDO5UsJHAcZvzRLPMaR2lVBGQCYRVsq1Go9HUaxq0H6yIzAYszpP5IrKrLuVxM+FASl0L4UbCRaRRj49G/vnReMcX7WpDdyrYJKCNzXmUWeasTqKINAGCgNRKtkUp9TbwNoCIbHLVV60hoMfXsNHja7iIiMsO9u40EWwEuohIBxHxxli0WuJQZwlgyQp6LbBaGd7mS4CpppdBB6ALsMGNsmo0Gk2N47YZrFKqSERuAZYBnsD7SqndIvIEsEkptQR4D/hIRA4BaRhKGLPe58AeoAi4WSlV7C5ZNRqNxh241QarlFoKLHUoe8Tm+BwwuYy2TwFPVeF2TrbRa1To8TVs9PgaLi6PrdHkItBoNJr6hg6V1Wg0GjfRoBSsiLQRkR9FZI+I7BaR253UERF5xQyz3SEiF9WFrK5QyfENF5FMEdlmvh5x1ld9RER8RGSDiGw3x/e4kzplhk/Xdyo5vlkikmzz+d1YF7K6ioh4ishWEfnWybUG+9lZqGB8Vf7sGpofbBFwt1Jqi4gEAptFZIVSao9NnSsxvA66AAOBeebfhkBlxgewRil1dR3IV13ygcuVUmdFxAtYKyLfK6V+s6ljDZ8WkakY4dNT6kJYF6jM+AA+U0rdUgfy1QS3A3uBZk6uNeTPzkJ544MqfnYNagarlDqplNpiHmdjvBGOEV4TgPnK4DcgWERa1bKoLlHJ8TVYzM/krHnqZb4cFwHKCp+u91RyfA0WEYkCrgLeLaNKg/3soFLjqzINSsHaYj5+9AF+d7jUKMJsyxkfwGDzMfR7Eeleu5JVD/MRbBtwBlihlCrz83MIn24QVGJ8AJNM89UiEWnj5Hp95SXgPqCsTbMb9GdHxeODKn52DVLBikgA8CVwh1Iqq67lqWkqGN8WoJ1SqhfwKrC4lsWrFkqpYqVUb4zovAEi0qOORapRKjG+b4D2SqkYYAV/zvjqNSJyNXBGKbW5rmVxB5UcX5U/uwanYE3b1pfAJ0qpr5xUqVSYbX2lovEppbIsj6Gmn7GXiITXspjVRimVAfyIkY7SFuvn5xA+3aAoa3xKqVSlVL55+i7Qt5ZFc5WLgfEikoCRGe9yEfnYoU5D/uwqHJ8rn12DUrCmPec9YK9S6r9lVFsCxJreBIOATKVUHe6MXnkqMz4RaWmxa4nIAIzPsEF8iUWkuYgEm8e+wEhgn0O1ssKn6z2VGZ/DesB4DDt7vUcp9aBSKkop1R4j4nK1UmqGQ7UG+9lVZnyufHYNzYvgYmAmsNO0cwE8BLQFUEq9iRE5NhYjh2wu8LfaF9NlKjO+a4G5IlIE5AFTG8qXGGgFxIuRPN0D+Fwp9a1UIny6gVCZ8d0mIuMxPEbSgFl1Jm0N0Ig+O6dU97PTkVwajUbjJhqUiUCj0WgaElrBajQajZvQClaj0WjchFawGo1G4ya0gtVoNBo3oRVsHSEiytaRWUSamJl6SmXxqaCfhIoCDcqqIyIBIvKWiBwWkc0i8pOIuC0xjoi0Fxc3phSR8SLygHk8UUQudKGPl0RkqHn8k4jsN0OOfxWRMje2q47c1cV23OXUmSUir5VRnmxmhzooIstEZIj7pC0fEfEWkV/MIITzAq1g644coIfpkA6GU3ptR5y9i+HP10Up1RfDZ7heRoUppZYopZ4xTycCVVKwIhIGDFJK/WJTPN0MOY4Hnq8RQWsYh3G7wmdKqT5KqS7AM8BXItKtunK5oiSVUgXAKhpehi2X0Qq2blmKkb0HYBqwwHJBREJFZLGZWOI3EYkxy8NEZLkY+UbfBcSmzQwx8pFuM2emnmXdWEQ6YaRx/KdSqgRAKXVUKfWdef0uEdllvu4wy9qLyD4R+VBEDojIJyJyhTkDPGhGliEij4nIRyKy3iy/ycn9PUXkeRHZaI7x72b5nSLyvnnc07y/n2WWZs7AxgPPm+PsJCJbbPrtYntuwyTghzLejl+AzmLwvHnPnSJSShGYM7DeNudrRaSXOeb3zZnxERG5zaZOdd5L6+xURMaJkWd1q4isFJEWZYzHKUqpHzG2P5lt9tdJRH4wn17WiMgFNuW/me/BkyJy1iwfbtZbAuwp6zM0695rU26bF3cxML0qcjdolFL6VQcv4CwQg5HWzQfYBgwHvjWvvwo8ah5fDmwzj18BHjGPr8JIhxcOdMNIRuFlXnsDiDWPE4Bwh/uPB74uQ7a+wE7AHwgAdmNk9mqPEcXSE+Of82bgfQwlPwFYbLZ/DNgO+JqyHQdam+13mXVmYyh3gKbAJqCD2e8vwF/MsovNOrOA18zjD4FrbeT9EehtHj8N3OpkTPHAOJvzn4B+5vG9wGcYSngFxiadLYA/MKKzbOWOA14yj7tiRPlYxrzOHEs4RviyVw28l7bjDuHP4KAbgf9zrOMw5lLlGLP/783jVRhPL2D8s11tHn8LTDOP5wBnzePhGE9eHSr4DEdhKHIxx/YtMNSs5wkk1/Xvr7Ze540tpD6ilNohRlrCaThsDglcgvGDRym12py5NgOGAteY5d+JSLpZfwTGj3mjGKkKfDFS5rnCJRjKNwdARL4CLsWINT+qlNpplu8GVimllIjsxFAaFv6nlMoD8kTkR2AAxj8RC6OAGBG51jwPwvixHxWRWcAO4C2l1K+VkPdd4G8ichfG4+cAJ3VaAckOZZ+ISB7GP6BbgbuABcrYwfi0iPwM9DdlsfAF8C8RuRe4HkPZW/hOGclA8kXkDIaSron30kIU8JkYMfHewNGK3hgnWPJYBABDgC/kz5StTc2/gzEUMcCnwAs27TcopSz3dfoZmuWjgK1meYBZ/otSqlhECkQkUBk5jxs1WsHWPUswvsDDqV7uTAHilVIPVrL+bqCXiHiqqm2Jnm9zXGJzXoL998kxBtvxXDBmmsuc3KMLxgy/dSVl+hJ4FFgNbFZKOUt+k4fxpGDLdKXUJqtAlcgNrZTKFZEVGLPMv2KfUcn2vSmm4t9XZd9LC68C/1VKLRGR4Riz5qrSByNJiQeQoYzUilUhx+bY6WcoIqOB/yil3iqjj6bAuSret0GibbB1z/vA45aZjA1rMG1V5o8pRRm5YX8BrjPLr8R4bATjce9aEYkwr4WKSLuybqqUOozxSPe4iDU7V3sRucq890TT9umP8bi+porjmiDGHlVhGP88NjpcX4aRtMbLvHdXEfEXkSAMM8hQIMxmdmRLNhBoM5ZzZn/zgA/KkGcv0LkCmdcAU0zbYnNThg1O6r1ryrhRKZXu5Lpjn9V9Ly0E8edCaFx5FZ0hIsMwHuvfMb9LR0VksnlNRKSXWfU3zKcnyk/Y4vQzNMuvN2fJiEikzfcyDOO7XFhV+RsiWsHWMUqpRKXUK04uPQb0FZEdGKu/lh/U48BQ85HyGgw7IcrYt+ufwHKzzQqMx+LyuBHjMfaQGG5IH2IkHd5iHm/A2FHhXaXU1rI6KYMdGLbR34B/K6VOOFx/F9gDbDHv/RbGrO1F4HWl1AGMPZ6esfw4bVgI3Gsu9nQyyz7BmPktL0Oe7zAUfXl8bcq9HWM2fJ9S6pRjJWUkZc6ibGVuW7cm3ksLj2E80m8GUirZZooYi4EHMDKzTVJKWdLsTQduEJHtGE80E8zyO4C7zO9RZ4ydCZzh9DNUSi3HMC2sN80di/jzH+JlGJ/FeYHOpqWpcUTkMYyFkRcqqluD97wHCFJK/aucOmuBq5WRDLs692qNsUh2gTI9MBoTIuIH5Jn24KkYC14TKmpXyb6/Ah4w/4E2erQNVtPgEZGvgU4Y3hblcTdGbt2MatwrFngKuKsxKleTvsBrpukoA2Mxr9qIiDeGd8R5oVxBz2A1Go3GbWgbrEaj0bgJrWA1Go3GTWgFq9FoNG5CK1iNRqNxE1rBajQajZvQClaj0WjcxP8DtXOkfQr2NNsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 360x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.figure(figsize=(5, 3))\n",
    "plt.plot(model_poly_degrees, bias_squared, color='blue', label='$bias^2$', linewidth=5)\n",
    "plt.plot(model_poly_degrees, variance, color='green', label='variance', linewidth=5)\n",
    "plt.plot(model_poly_degrees, np.array(bias_squared) + np.array(variance), linewidth=3, color='black', label='Total Error')\n",
    "plt.axvline(best_model_degree, color='black', linestyle='--', linewidth=2, label=f'Optimal Model Complexity (Degree={best_model_degree})')\n",
    "\n",
    "plt.xlabel('Model Complexity (Polynomial Degree)')\n",
    "plt.ylabel('Error')\n",
    "plt.ylim([0, .25])\n",
    "plt.xlim([2, 4.5])\n",
    "plt.legend()\n",
    "plt.title('Bias-Variance Tradeoff')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://spotintelligence.com/2023/04/11/bias-variance-trade-off/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7TUlEQVR4nO3deXhU5dn48e89W3YSEhJ2BBQQEAQJq7ZiqUAXgbpibV2qRa2K+rZVu6uvvrXb21+1CMVda91fWmxpcQNXlMUNAiIUsEQIhLBlT2bm/v1xJitZJmEmk2Tuz3Wda87ynDP3DOHc85znnOcRVcUYY0z8csU6AGOMMbFlicAYY+KcJQJjjIlzlgiMMSbOWSIwxpg454l1AG3Vq1cvHTx4cKzDMMaYLmXDhg0HVDW7qW1dLhEMHjyY9evXxzoMY4zpUkTks+a22aUhY4yJc5YIjDEmzlkiMMaYOGeJwBhj4pwlAmOMiXNd7q6htiothX/+E5KTISmp+deEBBCJdbTGGNPxun0i+PxzuOCC1suJtJwoIvnq81nSMcZ0Ht0+EZSXh1dOFcrKnCnavF7o0wf69m342nhd795OTaXbUIXqMig/VDdVHIGUbMgaBilZsY7QmLjU7RNBR5zY26q6GnbvdqbWZGaGlzQyMjqwlhHwOyfwisP1TuqHG53gm9kWrG7+uEk9IeskJyn0Cr1mnQSZQ8Gb2DGfzZg41O0TQbg1gs7q4EFn2ry55XIJCccmiKYSSO/ezqUp59d5eesn7mPWH4bKI60E0wMSMyApwzm554ysm0/qGdpWM98DSvbDgW1QtM153bEaPvpLvQMKZAyCXsNCyeHEuvke/ew6mzHHqdsngsxMOPdcp2ZQXt78a1VVrCM9PpWV8NlnztSarJTD9E0tIDPxACneMlJ8ZSR7y5x5b2g+QUlJySA5JZOUFCElzUNymo+UdB8pfZJIzkgmJSOZlJ5pJPXsgSulp3OyT0wHt7ftH2DY2Y0+UAkUbXemA9tC89vgszVQXVpXzpvsJIasYXXJoddJTk0iIa3tcRgTh6SrDVWZm5ur0ehryO93EkJLySISr2Vlznt1N0lJkJLiNIanpDScb2pd/XmXy6mgNDdBveWgouVH0JLC0HTAmUoL0bLDThkEVUF9PdCUbDSpF5ocmpKy0IQMVNzHvIfX69SsGk8+X9vWu91WSTHtU13d+nmkvBzOO6/t7YciskFVc5vaFtUagYjMBv4AuIEHVfWeRtt/D5wVWkwGclQ1I5oxNcfjgbQ0Z4q2khLYuxcKClp+LSysO0l1djV/oNEnQEZoGtYRb9hmIkqCL9hw8gZJ8AXweRsuO/M16wN1672B0PrQfL11zrIfnye03uMsJ3j9JHiD+Dz+0Hq/s97jx+cNICigoMHQH5Y28Uoz68N41aDzBbRn3+Zi6HsqTLwK3LG5eKHqnJxbOzG39Udhc9sCgfDi+tKXnEu9kRK1b1dE3MAi4GwgH1gnIstVtfZqt6reXK/8DcD4aMXTmaSmwrBhzoSqc/390E44tAsOhl4P7cJ/IJ/Cggr2FudQUJLD3uI+FJT1Z2/ViRRUDGJvSR8KjvZkb1Eq5RXuGH8qU0NVqKh0U1HZuf5NfO5KfO4qEtxVJHgq671W4nNXN1rnvNaWd1eS4KlyJrfz6vP4SfBU163zVpHgrsbnqa5b7wkd11MdKl+zrpoEbzVul4aqTwLiqjcvEAw4bUUfPglz7oN+44C6k3MkT8wt7RvuybkjRfpHVzTT7CRgu6ruABCRp4G5QHPNnhcDv4hiPNERDIK/3Gl4rZ3KwF/hvFY32uYvh9IDoZP9Tjj0GVQebXjM1N7QczCeoZPoO2EwfXsOgZ6DnSm1t3MtpR5VKC5uvYZRUODUMkx8qgokUBVIoCTWgdTjdrd0mc25DFh2oIjy27yUUUJ5IIWyMiEYjHXksRXpuyGjmQj6A/VvkMwHJjdVUEROAIYArzWzfQGwAGDQoEHti6ZkPxzeXXdy9tc7aVfXO2k3eQKvObE3UT5Q2fZY3Al1J/ZB0+rmew6GnieAL6VNhxOBHj2cafjwlstWV8P+/U5iOHrU+YMqLXWmpuZb215aChUVbf8KjAHn13bzz+/Uvwxo6utKNYK2mA88r6pNVsJUdSmwFJzG4na9w4dPwiu3t1xGXM5dKJ5E59Wb5Ny/7k127kBJyQmtqz+Fynkar09q4lhJdWVdsenmyeuF/v2dKVKCwbr/zK0ljaYSjIauDoQ7QdvKh3vM6mrn7qumpqqq8NfH+69V035ud10vBC31UNCjR2TfN5qJ4HNgYL3lAaF1TZkPXBfFWGDkHMgZ1egEXXMiD52s3db3Q3u4XE67R2pqrCPpHPz+tiWOlta3Z5+m1le38ByfaVnNybnxybg9Xcy0VsbbjjuvIyGaiWAdMExEhuAkgPnANxsXEpGTgZ7AmijGErrX/MSovoUx4NyB5vE4/7E7i2Cw+RpPLJNWe++K83iUZE8pSe5S50SalUlyqrddJ9/Wysbq5NyRopYIVNUvItcDK3FuH31YVfNE5E5gvaouDxWdDzytXe2BBmO6EJerrhG2s1B12ghaShzQ8Nd4zeT1CgQS4b1H4LW7weWBL/8Ccq+M2WXXrsweKDPGdG0Hd8Lfb3K6Jhk4BebcC9kjYh1Vp9PSA2WWOo0xXVvmEPj2X2HeYij8BJacAa//GvxdvN+YDmSJwBjT9YnAuG/C9etg5Dmw6m740xdh97pYR9YlWCIwxnQfqTlw/sNw8TPOg5oPnQ3/vBUqi2MdWadmicAY0/2MmA3XvQeTvgvv/QnunwrbXo51VJ2WJQJjTPeUkAZf/Q18Z6XznNCT58MLVzldvJgGLBEYY7q3QZPhmjdh+o8g76/wx4nw0dNdp2vfDmCJwBjT/XkSYPptTkLIOgmWXQ1/Ptfp/NFYIjDGxJGckc6loq/+FnavddoO1ixyuryOY5YIjDHxxeVyGpGvew8GfwFW/hge/DIUbIp1ZDFjicAYE5/SB8A3n4HzHoLD/4GlZ8KrdzrdzMcZSwTGmPglAmPOdx5EG3MhvPk7WHI67Hor1pF1KEsExhiTnAnfWAzfXgaBanj0a/DijVBxJNaRdQhLBMYYU+PEL8H31sDU6+H9x+GPk2DLi7GOKuosERhjTH2+FJh1N1z1KqRkwzPfcqbiglhHFjWWCIwxpin9T4MFq2DGL+DTl5zawYbHuuWDaJYIjDGmOW4vfOG/nMtFfcbAiwvhsXOg6N+xjiyiLBEYY0xrsk6Ey16Ec+6FvR87D6K9+b9Ow3I3YInAGGPC4XLBhMvg+rUwfBa8egcsPQs+fz/WkR23qCYCEZktIltFZLuI3NZMmQtFZLOI5InIX6IZjzHGHLe0PnDRE3DRn6G0EB6cASt/AlWlsY6s3aKWCETEDSwCvgKMAi4WkVGNygwDfgScrqqjgZuiFY8xxkTUyHOcbipOuwzW/NG5XPTv12IdVbtEs0YwCdiuqjtUtQp4GpjbqMx3gUWqeghAVfdHMR5jjImspAw45//B5SuchuUnvgHLroWyg7GOrE2imQj6A7vrLeeH1tU3HBguIm+LyLsiMrupA4nIAhFZLyLrCwsLoxSuMca00+DT4Zq34Qs/gI3PwqJJsPH5LnOraawbiz3AMGA6cDHwgIhkNC6kqktVNVdVc7Ozszs2QmOMCYc3EWb8DBa8DukD4YUr4an5cCQ/1pG1KpqJ4HNgYL3lAaF19eUDy1W1WlV3Ap/iJAZjjOma+pwCV70Cs34JO9+ARZNh7QMQDMY6smZFMxGsA4aJyBAR8QHzgeWNyvwVpzaAiPTCuVS0I4oxGWNM9LncMPV78L13YeAkWPEDeHgW7P8k1pE1KWqJQFX9wPXASmAL8Kyq5onInSIyJ1RsJVAkIpuBVcAPVbUoWjEZY0yH6nkCfOv/4BtLoWg7LDkDVt8D/spYR9aAaBdpzKiRm5ur69evj3UYxhjTNqUH4F+3wcbnIPtkmHOfU1voICKyQVVzm9oW68ZiY4yJDym94LwH4ZLnnYfPHpoJK34IlcWxjswSgTHGdKhhZzttB5OvdhqRF02GT1fGNCRLBMYY09ESUuErv4IrX4aEHvCXC+H570BJbJ6TskRgjDGxMnAiXP0GnPUTZyS0RRPhw790+INolgiMMSaWPD448xa45i3oNQL+ei08MQ8O7uywECwRGGNMZ5A9Aq74J3ztd5C/wenE7p37IOCP+ltbIjDGmM7C5YKJVzm9mg6dDi/91Onmeu/H0X3bqB7dGGNM26X3h4ufggsehaOfw9Lp8MrtUF0elbezRGCMMZ2RCIz+Bly3FsZdDG/9Ht77U1TeyhOVoxpjjImM5EyYuwhOvRj6T4jKW1giMMaYrmDwGVE7dIuXhsQxsKUyxhhjurYWE4E6PdKt6KBYjDHGxEA4jcXvi8jEqEdijDEmJsJpI5gMXCIinwGlgOBUFsZGNTJjjDEdIpxEMCvqURhjjImZVi8NqepnQAZwTmjKCK0zxhjTDbSaCETkRuBJICc0/VlEboh2YMYYYzpGOI3FVwKTVfXnqvpzYArw3XAOLiKzRWSriGwXkdua2H65iBSKyIeh6aq2hW+MMeZ4hdNGIECg3nIgtK7lnUTcwCLgbCAfWCciy1V1c6Oiz6jq9WHGa4wxJsLCSQSPAO+JyLLQ8jzgoTD2mwRsV9UdACLyNDAXaJwIjDHGxFBrTxa7gHeBK4CDoekKVf1/YRy7P7C73nJ+aF1j54nIxyLyvD3FbIwxHa/FGoGqBkVkkaqOB96Pwvu/CDylqpUicjXwGPClxoVEZAGwAGDQoEFRCMMYY+JXOI3Fr4rIeSLSartAI58D9X/hDwitq6WqRapaGVp8EGiyaz1VXaqquaqam52d3cYwjDHGtCScRHA18BxQKSJHRaRYRI6Gsd86YJiIDBERHzAfWF6/gIj0rbc4B9gSZtzGGGMipMVLQ6E2gtmq+nZbD6yqfhG5HlgJuIGHVTVPRO4E1qvqcmChiMwB/DjtD5e39X2MMcYcH3E6GG2hgMgHoTaCTiE3N1fXr18f6zCMMaZLEZENqprb1LZothEYY4zpAtrSRlDVxjYCY4wxXUCrD5SpalpHBGKMMSY2wul0TkTkWyLys9DyQBGZFP3QjDHGdIRwLg3dD0wFvhlaLsHpQ8gYY0w3ENYIZap6moh8AKCqh0LPBRhjjOkGwqkRVId6ElUAEckGglGNyhhjTIcJJxHcCywDckTkbuAt4H+iGpUxxpgOE85dQ0+KyAZgBs44BPNU1bqCMMaYbiKcNgJU9RPgkyjHElVb9h7l4/zDiAgCta8uFwiCCPW2Oetc4swT2u6qt93lFKo9lkvqHafm+KF5l6tuPxB6pfoYlJmMPaNnjOkMwkoE3cHrnxZyzz87Ty7r3SOByUOymDw0k8lDsjgxO8USgzEmJuImEVwyeRDnnNoPVUUVZ8KZD6qihNbVn0cJBuvK1d9HCe2nQO1x6vYPhgrVHKvuPZT8Q+W8t/Mg7+4oYvlHewDolZrA5KGZTBmSyeShWQzLSbXEYIzpEHGTCNISvaQlemMdRq1vTTkBVWVXURnv7Sji3R1FvLfzIP/4eC8AmSk+Jg/JdKahWYzonYbLZYnBGBN5zSYCESkmdMtoU1S1R1QiiiMiwpBeKQzplcL8SYNQVXYfLOfdnUW8t8OpMfxzUwEAGcleJg12ksLkIZmM7NsDtyUGY0wENJsIavoYEpH/BvYCT+C0fV4C9G1uP9N+IsKgrGQGZSVzYa4zuFv+obLapPDezoO8tHkfAD0SPUwaklnbzjCqbw887nDuBjbGmIbCGY/gI1U9tbV1HSXexyPYe6Sc93Yc5L2dRby74yA7D5QCkJrgIXdwT6aEagyn9E/Ha4nBGBPS0ngE4bQRlIrIJcDTOJeKLgZKIxifaYO+6UnMG9+feeP7A7DvaAXv7TxY286wemshAMk+N7mDnTaGKUMzGdM/A5/HEoMx5ljh1AgGA38ATsdJBG8DN6nqrmgH15R4rxG0prC4krU7nRrDezsOsnVfMQBJXjennZDBlCFZTB6axakD00nwuGMcrTGmo7RUI2g1EXQ2lgja5mBpFWtDl5He23mQTwqOogoJHhd90hNxhR6Gc17rHpxzuUIP0DXYTqPleuXl2PJuV8vbXS7neIkeN8k+N8kJbpK9bpITPM6yz02yz9PgNcXnIcnnttqNMW10XJeGRGQ4sBjoraqniMhYYI6q3hXGvrNxahNu4EFVvaeZcucBzwMTVdXO8hGUmeJj9il9mX2K075/uKyKtTsPsnbnQQpLKgnWPOOgzjMTQdXa5yFq5muelwjWWxcIBmu31ZQPBMPYN1izvW5dRXWQsio/wTb8JvG6hSSvm5QEJzHUJIiUBsnDTZLPQ4rP7WxL8DRKLg3nfR4XXrcLn9tlt+qauBLOpaHXgR8Cf6oZxF5ENqnqKa3s5wY+Bc4G8oF1wMWqurlRuTTgH4APuL61RGA1gu5JVan0BymrClBa6ae8OkBZVYCySr+zrspPeVWA0qoA5VX+0KtTtqy6rpwzNZyvDrS91utxCV63C69b8Hmc5OANJQqv2xVaJ7XLXreLBI9T3hsq6wuVq1lXe5zaSRquq91HGr1PXfnaMpasTBsdb2NxsqqubfSUqz+M/SYB21V1RyiIp4G5wOZG5f4b+BVOsjFxSkRI9LpJ9LrJTInscBdV/iDlVQHKqkMJorJhsqhJMlX+IFWBINWhqcofpDqgzjp/aF0gSJVfG5QprfRTFdBG+zXcP9CW6k6Y3C5xkkltwqlLHk7NRo5dF0oyvnqJpyYh1ZatXedusK0uGTbxnm4XXo80Oq7LnnXpIsJJBAdE5ETqxiM4H+e5gtb0B3bXW84HJtcvICKnAQNV9R8i0mwiEJEFwAKAQYMGhfHWxtTxeZyTVjqxe7I8EKxLFNUBrUsWxyQPbbgcKn/MOr9SFQg0ONYxyavePmXlgQbJrNofrE1eNWX8UUxWNcmiLnk0UUsK1bLqJ5hBmcnMHN2bEb3TrMuVKAonEVwHLAVOFpHPgZ04D5UdFxFxAf8LXN5aWVVdGoqB3NzcrtW6bQzOCdHtcmo8nVUwqFQHGyWqY5JPw/XONg0llkbJqEGZpmtZtTUwv1JeXn1MjeqFoxX878ufckJWMrNG92HW6N6MH9jTLotFWIuJIHSd/3uq+mURSQFcqloc5rE/BwbWWx4QWlcjDTgFWB3K9H2A5SIyxxqMjel4LpeQ4HKT4AESYh2NY39xBa9s3s/KvAIeeXsnS9/YQXZaAjNH9WbW6D5MGZpld5BFQDiNxe+q6pQ2H1jEg9NYPAMnAawDvqmqec2UXw38wBqLjTFNOVpRzapPnKSwemshZVUB0hI9zDg5h1mj+3DmiGySfXHTj2abHW9j8Qcishx4jnpPFKvq/7W0k6r6ReR6YCXO7aMPq2qeiNwJrFfV5WF/AmNM3OuR6GXuuP7MHdefiuoAb207wMq8Al7Zso+/friHBI+LLw7PZtboPnx5ZA4ZyZG96aA7C6dG8EgTq1VVvxOdkFpmNQJjTH3+QJC1uw7yUt4+VuYVsPdIBW6XMHlIJrNG92Hm6N70TU+KdZgxZ08WG2Pigqqy8fMjrMwrYGXePrbvLwHg1AHpzDqlD7NG9+HE7NQYRxkbx5UIRCQRuBIYDSTWrLcagTGms9u+v4SVeQW8lFfAR/lHADgpJ5VZo53G5jH90+PmttTjTQTP4Qxc/03gTpxbR7eo6o2RDjQclgiMMe2x90h57eWj93YeJBBU+qUnMjN0+WjS4MxuPabH8SaCD1R1vIh8rKpjRcQLvNmeO4kiwRKBMeZ4HSqt4pUt+1iZt483txVS6Q/SM9nLl0c6NYUzhvXq1M98tMfx3jVUHXo9LCKnAAVATqSCM8aYjtYzxccFuQO5IHcgZVV+Xt9ayMq8Av6VV8BzG/JJ9rmZPsK5A+msk3Po0YnGO4+GcBLBUhHpCfwMWA6kAj+PalTGGNNBkn0evjKmL18Z05cqf5B3dxTxr7wCXt68jxUbC/C6hWkn9mLW6D6cPao32Wmd5Gm7CLK7howxpgnBoPLB7kOsDLUrfFZUhghMGNQz1N1FHwZlJcc6zLAdbxtBk7/+VfXOCMTWZpYIjDEdTVXZuq+YlZucpLB571EATu6TxuzQbakn9+ncHeMdbyL4fr3FRODrOHcN2e2jxpi4tPtgWehZhQLWf3YIVRiUmVx7W+ppgzpfx3gRfaBMRBKAlao6PQKxtZklAmNMZ1JYXBm6A6mAt7cfoDqg9EpN4OxRvZk1ujfTTuzVKTrGi3Qi6AmsU9WTIhFcW1kiMMZ0VjUd472Ut49VW/fXdoz3pZqO8YZnk5IQm47xjnfM4o2EBqXB6TwuG+fBMmOMMfU07hjv7e01HePt52+hjvG+MKwXM0f34csje0d8NL72Cic1fb3evB/Yp6rhDFVpjDFxK9HrZsbI3swY2Rt/IMi6XYdqu7t4Zct+3C5h0uBMZo3uzczRfeiXEbuO8cJpLM5sabuqHoxoRK2wS0PGmK5MVdn0+dHaxuZtoY7xxg5Irx2F7aSctIi/7/HeNbQLZ6SxQ4AAGcB/QptVVYdGLNIwWCIwxnQnOwpLWJm3j3/lFfDR7sMAnJidUvuswtgBkekY73gTwQPAMlVdEVr+CjBPVa8+7sjawRKBMaa72nuknJc3O3cgvbvD6Rivb3pi7dCck4a0v2O8400EG1V1TGvrOoolAmNMPDhcVsWrW5yhOV//1OkY78dfPZkFXzyxXcc73k7n9ojIT4E/h5YvAfa0KxJjjDFhyUj2cd6EAZw3YQBlVX7e+LSQMQMyovJe4dQxLsa5ZXRZaMoOrWuViMwWka0isl1Ebmti+zUislFEPhSRt0RkVFuCN8aYeJDs8zD7lL70j9KdRa3WCEJ3Bd0IICJuIEVVj7a2X6jsIuBsIB9YJyLLVXVzvWJ/UdUlofJzgP8FZrf5UxhjjGm3VmsEIvIXEekhIinARmCziPwwjGNPArar6g5VrQKeBubWL9AooaRQ9+CaMcaYDhLOpaFRoRP2POCfwBDg22Hs1x/YXW85P7SuARG5TkT+DfwaWNjUgURkgYisF5H1hYWFYby1McaYcIWTCLyh4SnnActVtZoI/nJX1UWqeiJwK/DTZsosVdVcVc3Nzs6O1FsbY4whvETwJ2AXzqWbN0TkBKDVNgLgc5wH0WoMCK1rztM4ycYYY0wHajURqOq9qtpfVb+qzkMH/wHOCuPY64BhIjJERHzAfJyhLmuJyLB6i18DtoUfujHGmEhoc3+ooWTQaqdzquoXkeuBlTi9lj6sqnkiciewXlWXA9eLyJeBapwuLC5razzGGGOOT1Q7xg51S7Gi0bqf15u/MZrvb4wxpnWxHzbHGGNMTIVVIxCRacDg+uVV9fEoxWSMMaYDhTNC2RPAicCHQCC0WgFLBMYY0w2EUyPIxXmozJ76NcaYbiicNoJNQJ9oB2KMMSY2wqkR9MLpX2gtUFmzUlXnRC0qY4wxHSacRHB7tIMwxhgTO+F0Q/16RwRijDEmNsLphnqKiKwTkRIRqRKRgIiE09eQMcaYLiCcxuI/4oxItg1IAq7CGXDGGGNMNxDWk8Wquh1wq2pAVR/BRhEzxphuI5zG4rJQ76Efisivgb1Y1xTGGNNthHNC/3ao3PVAKc4YA+dFMyhjjDEdJ5y7hj4TkSSgr6re0QExGWOM6UDh9DV0DvBbwAcMEZFxwJ32QJkx8a26upr8/HwqKipiHYqpJzExkQEDBuD1esPeJ9wHyiYBqwFU9UMRGdKeAI0x3Ud+fj5paWkMHjwYEYl1OAZQVYqKisjPz2fIkPBP0+G0EVSr6pHG79em6Iwx3U5FRQVZWVmWBDoRESErK6vNtbRwagR5IvJNwB0aY3gh8E47YjTGdDOWBDqf9vybhFMjuAEYjdPh3FPAUeCmNr+TMcaYTqnVRKCqZar6E1WdqKq5ofmw6h0iMltEtorIdhG5rYnt/yUim0XkYxF5VUROaM+HMMbEn6KiIsaNG8e4cePo06cP/fv3r12uqqpqcd/169ezcOHCVt9j2rRpkQq3U2v20pCILG9px9buGhIRN05XFGcD+cA6EVmuqpvrFfsAyFXVMhG5Fvg1cFG4wRtj4ldWVhYffvghALfffjupqan84Ac/qN3u9/vxeJo+xeXm5pKbm9vqe7zzTsddBW8cb0vxt7Rfe7RUI5gKDADexLl99HeNptZMArar6g5VrQKeBubWL6Cqq1S1LLT4buj9jDGmXS6//HKuueYaJk+ezC233MLatWuZOnUq48ePZ9q0aWzduhWA1atX8/Wvfx1wksh3vvMdpk+fztChQ7n33ntrj5eamlpbfvr06Zx//vmcfPLJXHLJJdQM2rhixQpOPvlkJkyYwMKFC2uPW18gEOCHP/whEydOZOzYsfzpT3+qPe4XvvAF5syZw6hRo45Zrqio4IorrmDMmDGMHz+eVatWAfDoo48yZ84cvvSlLzFjxozj/t5aSiN9cH7NXwx8E/gH8JSq5oV57P7A7nrL+cDkFspfCfyzqQ0isgBYADBo0KAw394Y01HueDGPzXsi2ynxqH49+MU5o9u8X35+Pu+88w5ut5ujR4/y5ptv4vF4eOWVV/jxj3/MCy+8cMw+n3zyCatWraK4uJgRI0Zw7bXXHnMf/gcffEBeXh79+vXj9NNP5+233yY3N5err76aN954gyFDhnDxxRc3GdNDDz1Eeno669ato7KyktNPP52ZM2cC8P7777Np0yaGDBnC6tWrGyz/7ne/Q0TYuHEjn3zyCTNnzuTTTz+t3e/jjz8mMzOzzd9RY80mAlUNAP8C/iUiCTgJYbWI3KGqfzzud65HRL6FMzbymc3EshRYCpCbm2u3rhpjmnXBBRfgdrsBOHLkCJdddhnbtm1DRKiurm5yn6997WskJCSQkJBATk4O+/btY8CAhhcoJk2aVLtu3Lhx7Nq1i9TUVIYOHVp7z/7FF1/M0qVLjzn+Sy+9xMcff8zzzz9fG9e2bdvw+XxMmjSpwT3/9ZffeustbrjhBqoD1ST3TSarbxabNm8C4Oyzz45IEoBWbh8NJYCv4SSBwcC9wLIwj/05Tr9ENQaE1jV+jy8DPwHOVNXKxtuNMZ1fe365R0tKSkrt/M9+9jPOOussli1bxq5du5g+fXqT+yQkJNTOu91u/H5/u8o0R1W57777mDVrVoP1q1evbhBv4/iDGmR/yX427t9IUIMAFJYVHlPueDXbRiAijwNrgNOAO0J3Df23qh5zMm/GOmCYiAwJ9V46H2jQAC0i44E/AXNUdX+7PoExxjTjyJEj9O/fH3Cuq0faiBEj2LFjB7t27QLgmWeeabLcrFmzWLx4cW2N5NNPP6W0tLTZ4/qDfvKP5nPiuBN56qmnCGqQz/79GQWfF5DeP53qQNM1m/ZqqbH4W8Aw4EbgHRE5GpqKwxmhTFX9OD2WrgS2AM+qap6I3CkiNXcc/QZIBZ4TkQ9bu1PJGGPa4pZbbuFHP/oR48ePb9Mv+HAlJSVx//33M3v2bCZMmEBaWhrp6enHlLvqqqsYNWoUp512GqeccgpXX311k/EEggEq/BVs3LeRgpICzrv0PILBIPNnzOfH1/6YX/z+F/gSfByuOBzRzyE1Ld9dRW5urq5fvz7WYRgT97Zs2cLIkSNjHUbMlZSUkJqaiqpy3XXXMWzYMG6++eY2HSMQDLC/dD8FJQUENNBiWZe46J3Sm35p/Zp9iripfxsR2aCqTd4ze3w3nxpjTJx74IEHeOyxx6iqqmL8+PFcffXVYe8bCAYoLCukoKQAf7DlGotLXGQnZ9MntQ9ed/g9i4bDEoExxhyHm2++uc01gGAwSGFZIXtL9raaAAQhOyWbvql9I54AalgiMMaYDhLUIAfKDrC3eC/VwZYbfAWhV3Iv+qb1xef2RTUuSwTGGBNlQQ1SVFbE3pK9VAVa7gcJcBJAal8SPAmtlo0ESwTGGBMlqkpReRF7i/dSGWj9MamspCz6pvUl0ZPYAdHVsURgjDERpqocLD/InuI9YSWAzKRM+qX2I9HbsQmghiUCY0yXVFRUVNvhWkFBAW63m+zsbADWrl2Lz9fydfXVq1fj8/lqu5pesmQJycnJXHrppe2OSVU5VHGIPcV7qPC33lt/z8Se9EvrR5I3qd3vGQmWCIwxXVJr3VC3ZvXq1aSmptYmgmuuuabdsagqhysOs6d4D+X+8lbLZyRmkJOUQ4+kHrXrwu1OOhAI1PalFCnhjFBmjDFdwoYNGzjzzDOZMGECs2bNYu/evQDce++9jBo1irFjxzJ//nx27drFkiVL+P3vf8+4ceN48803uf322/ntb38LwPTp07n11luZNGkSw4cP58033wSgrKyMCy+8kFGjRvGNb3yDyZMns+rtVWw5sIV/H/p3bRLY8vEWFpy3gG/P/jY3fPMGDuw7AMB1F1zHI798hPkz5/PA/Q8wffp0brrpJnJzc/nDH/7Aq6++yvjx4xkzZgzf+c53qKx0LisNHjyYW2+9ldNOO43nnnsu4t+b1QiMMcfvn7dBwcbIHrPPGPjKPWEXV1VuuOEG/va3v5Gdnc0zzzzDT37yEx5++GHuuecedu7cSUJCAocPHyYjI4NrrrmmQS3i1VdfbXA8v9/P2rVrWbFiBXfccQevvPIK999/Pz179iQvL4/33n+PMyafwe6ju0mrTqvbr9rPb376G373yO/omdWTl/72Eg/89gEee+QxkrxJaECp6R3hxRdfpKqqivXr11NRUcGwYcN49dVXGT58OJdeeimLFy/mpptuApwa0Pvvv3+cX2rTLBEYY7qFyspKNm3axNlnnw04l1D69u0LwNixY7nkkkuYN28e8+bNC+t45557LgATJkyo7VTurbfe4rvXfpetRVvx9PVw0siTjtlv1793sWPrDq6bfx1ucePCRf9+/Un1OYPcXHRRw0EYa5a3bt3KkCFDGD58OACXXXYZixYtqk0EjfeLJEsExpjj14Zf7tGiqowePZo1a9Ycs+0f//gHb7zxBi+++CJ33303Gze2Xnup6Xa6psvpksoSSqpK2H10N72rercQCJw04iTeePsN0nxpx/QH1FK30y2JZLfTjVkbgTGmW0hISKCwsLA2EVRXV5OXl0cwGGT37t2cddZZ/OpXv+LIkSOUlJSQlpZGcXFxq8ctrSrFH/TzSdEnjJ4wmldefAWAHZ/uYPsn2xuUTfYmM2PiDEoOl5D3fl7tYDh5ea0P7DhixAh27drF9u3OMZ944gnOPLPJsboizhKBMaZbcLlcPP/889x6662ceuqpjBs3jnfeeYdAIMC3vvWt2nF/Fy5cSEZGBueccw7Lli2rbSxurKK6gm1F29h2cFvtoDDnX3Y+h4oOceH0C1n868UMHT6U1LRUkjxJnJR5EiN7jSS7R3aTcbQmMTGRRx55hAsuuIAxY8bgcrmO606mtrBuqI0x7dJdu6Euqy5jT/GeJvv8DwQC+Kv9JCQmkL8rn+vnX8+GjRvI6ZHTbJfQsWDdUBtjTDuUV5ezp3gPhyoONVumoryCay+4loA/gMflYemSpfROb6G9oIuwRGCMiWsV/gr2FO/hYPnBVstmpmfyxpo3yErK6lQ1gONlicAYE5cq/ZXsLdnLgbIDrZb1uX30Te1LVnIWLul+TatR/UQiMltEtorIdhG5rYntXxSR90XELyLnRzMWY4wBqPJX8dnhz9i0f1OrScDr8jIofRCn5JxCdkp2t0wCEMUagYi4gUXA2UA+sE5Elqvq5nrF/gNcDoTfQYgxxrRDVaCKgpICCksLUVq+Scbr8tIntU+3PvnXF81LQ5OA7aq6A0BEngbmArWJQFV3hbYFoxiHMSaOVQeqKSgpYH/Zflq7S9Lj8jgJIDkbtyuyHbt1ZtFMBP2B3fWW84HJ7TmQiCwAFgAMGjTo+CMzxnR7/oCfgtIC9pfur30OoDkel4feKb3JScmJqwRQo0s0FqvqUmApOM8RxDgcY0wn5g/62Veyj32l+1pNAG5x0zu1N71TesdlAqgRzUTwOTCw3vKA0DpjjIm4QDDAvtJ97CvZR0ADLZZ1iYveKb3pndobj6vuNPjXv/6Vf/zjHxw9epQrr7ySmTNnRjvsTiGarSDrgGEiMkREfMB8YHkU388YE4defe1Vzp9/Phv3b2RP8Z7aJPDRuo9Y8pslDcq6xEWf1D6MyRlD/x79GyQBgHnz5vHAAw+wZMkSnnnmmQ77DLEWtRqBqvpF5HpgJeAGHlbVPBG5E1ivqstFZCKwDOgJnCMid6jq6GjFZIzpPgLBAIVlhax8eyUDRwzEH/Q32H7qxFM5deKpAIgIOck59Entg9ftbfXYd911F9ddd11U4u6MonpflKquUNXhqnqiqt4dWvdzVV0eml+nqgNUNUVVsywJGNP1iER/qi+oQfaV7GPT/k3kH81ny6Yt7C/Yz+Vfv5y5U+ey4Z0NANy24DY+fO9DclJy2PbmNi6YdQG5p+VyxhlnUFhYCMBjjz3GhAkTGDt2LGeccQaqyq233spXvvIVTjvttI7+KmOmSzQWG2NMUIMcKDtAQUkBVYGq2vXbNm/jizO/yKN/f5R3X3+XJb9ZwoPLHmTXtl3M+eIcctJzSJmRwvwL5wNwxx138Oyzz3LppZfyq1/9ig8//BCfz8fhw4e57777eOWVVzhy5Ajbt2/vsN4/Y80SgTGm0yssLWRvyd4GCQCcYSEPHzzMFTdcAcDw0cMpPlzMST1OQv1KTlYOAI8++ijPPPMMlZWVFBQU8D//8z+43W7Ky8v5/ve/z2WXXUZubi4LFy5k4cKFHf75Yq37PzJnjOnyPjvy2TFJAGDX9l0MHDwQr89LZlImFbsryB2fy/at2xk1ahQAjz/+OGvXruW1117jo48+YsSIEYwePZrk5GQ2bdrE6aefzoIFC7j//vs7+mN1GlYjMMZ0WZ/mfUpBfgEn9TgJDx5+edcv+f3vf8/GjRsZO3YsABs3bmTatGmkpqbywgsv8M477zBmzBi2bdvGsGHDmD9/Pps3b6aioiLGnyZ2LBEYY9osEAwQ1CAV/gqKK/0EggH8QT/+oJ+A1ptvYn2kZCRmULSriAvPv5AZZ86gvLycn/3sZ0yZMoXnnnuOSZMmAXD55Zdz7rnn8uSTTzJz5kyGDh1KSkoKd999N2vWrCElJYXRo0fzwAMPRCy2rsZGKDMmjgWCAQ5XHOZg+UEOlh+kqLyodv5g+UGKyoo4WNFoufwghysOs2LmCnqd0KvDY05PSKdfWj9SfNEbzL2rsxHKjIlD/qC/4Qm9rOEJvcmTfHlRk8MxdlY9EnrQL60fqb7UWIfS7VgiMKYT8Qf9HCo/1OIJvKl1XemE3lZpvjT6pfUjLSEt1qF0W5YIjImCmhN6k5daapYrjv3lfqTySKxDjzqXuPC4PLWTW9x18666eY948Hl8+Ny+WIfc7VkiMB0qqEGqAlUNJkEanBjqT7EeF7Y6UM2hikNhX2qpmT9aeTSmcXcEEcHn9jV5Qm9wUpe6ebfLHRcDvXQ1lgji0JGKI3xQ8AGHyg9RFaiiMlB5zMm5panN5f115VvrFbKxxr8em5san4TaOgEcqTxyzC/34qriaPwTdCrpCelkJmXWTlnJWWQmNlquvz0pi4zEDLZ/up2RvUe2/gam07NE0M2pKluLtrJm9xrW5DtT3v68Vofq6yzq1yBMy9IT0ps8abe03DOp5zE9cJr4Y38B3UxxZTFrP19be9J/N/9dDpYfjHVYpg0yEjPCOonX/7WekZhhJ3TTbvaX04WpKtsPbndO+qFf/Bv3b2x1VCYTfYLUntAbXF4J45JLPI+U1ZG++tWv8pe//IWMjIxYhxJzlgi6kNKqUtbtWdfgMs+BsgOxDqvNfG5f7eR1OX3D1zx9Wn9qa3tCNAhCz6Sex/4Sb+XXejyd0OWO6Dfo6y8idylTVVFVVqxYEbFjdnWWCDopVWXn4Z0NTvofFXwUsZPjSZknMSp7FImeROek7PI1OEH73D4SPAnHrGvPlOCuO05b7gRS1QbdFTRIEsGm17c2tXS8Hgk9jvmlnp6QHjcn9K7mtttuY+DAgbUDyNx+++14PB5WrVrFoUOHqK6u5q677mLu3Lns2rWLWbNmMXnyZDZs2MCKFSs488wzWb9+Pb169WLevHns3r2biooKbrzxRhYsWABAamoqN954I3//+99JSkrib3/7G71792bfvn1cc8017NixA4DFixczbdo0/vznP3PvvfdSVVXF5MmTuf/++3G7O//fjyWCTqK8upz1e9bXnvTX7F7DvtJ9ETl2sjeZif0mMm3gNKYOmMqUAVPITsmOyLGjSUTwiMeufZsmXXTRRdx00021ieDZZ59l5cqVLFy4kB49enDgwAGmTJnCnDlzANi2bRuPPfYYU6ZMOeZYDz/8MJmZmZSXlzNx4kTOO+88srKyKC0tZcqUKdx9993ccsstPPDAA/z0pz9l4cKFnHnmmSxbtoxAIEBJSQlbtmzhmWee4e2338br9fK9732PJ598kksvvbRDv5f2sP9hMaCq/OfIfxpc2/+g4IOIdcg1JGMIUwdOZdqAaUwdOJWxvcfaydR0O+PHj2f//v3s2bOHwsJCevbsSZ8+fbj55pt54403cLlcfP755+zb5/ygOuGEE5pMAgD33nsvy5YtA2D37t1s27aNrKwsfD4fX//61wGYMGECL7/8MgCvvfYajz/+OABut5v09HSeeOIJNmzYwMSJEwEoLy8nJycnqt9BpET17CAis4E/4IxZ/KCq3tNoewLwODABKAIuUtVdkYxh56Gd/G7N75zrgjjXBoMarJ1Xjl1utky97U2tC/c4nxZ9yp7iPRH5fImeRHL75TJ1wFSmDZzGlAFT6JPaJyLHNqazu+CCC3j++ecpKCjgoosu4sknn6SwsJANGzbg9XoZPHhwbffSKSlNd1K3evVqXnnlFdasWUNycjLTp0+v3cfr9dZeynS73fj9zf9YU1Uuu+wyfvnLX0b4U0Zf1BKBiLiBRcDZQD6wTkSWq+rmesWuBA6p6kkiMh/4FXBRJOPYV7qPResWRfKQMTUofRBTB0ytPfGf2udUewTfxFQkG3Lb6qKLLuK73/0uBw4c4PXXX+fZZ58lJycHr9fLqlWr+Oyzz1o9xpEjR+jZsyfJycl88sknvPvuu63uM2PGDBYvXsxNN91Ue2loxowZzJ07l5tvvpmcnBwOHjxIcXExJ5xwQiQ+alRFs0YwCdiuqjsARORpYC5QPxHMBW4PzT8P/FFERCPYN3ZXfpzd5/Yxoe8E58Q/0Dn59+/RP9ZhGdNpjB49muLiYvr370/fvn255JJLOOeccxgzZgy5ubmcfPLJrR5j9uzZLFmyhJEjRzJixIhmLx/V94c//IEFCxbw0EMP4Xa7Wbx4MVOnTuWuu+5i5syZBINBvF4vixYt6hKJIGrjEYjI+cBsVb0qtPxtYLKqXl+vzKZQmfzQ8r9DZQ40OtYCYAHAoEGDJoST5Wus+3wdkx6cdLwfp0P0S+tX26A7beA0xvcZT4InIdZhGdOkpvq8N51DtxyPQFWXAkvBGZimLfvGutOy5nhcHk7re1rtZZ6pA6cysMfAThuvMab7imYi+BwYWG95QGhdU2XyRcQDpOM0GkfMoPRB3Dv7XkQEQXCJq3a+Letc4mqwPdx1TR07LSGNMTljSPImRfKjGmNMu0QzEawDhonIEJwT/nzgm43KLAcuA9YA5wOvRbJ9ACAnJYcbJt8QyUMaY0JU1WqxnUx7TqFRa0lVVT9wPbAS2AI8q6p5InKniMwJFXsIyBKR7cB/AbdFKx5jTGQlJiZSVFTUrhOPiQ5VpaioiMTExDbtZ4PXG2Papbq6mvz8/Np77k3nkJiYyIABA/B6vQ3Wd/nGYmNM5+P1ehkyZEiswzAR0HVvsjfGGBMRlgiMMSbOWSIwxpg41+Uai0WkEAj/0eLOqRfQ9UaUiR77PurYd9GQfR8NHc/3cYKqNtn/fJdLBN2BiKxvrvU+Htn3Uce+i4bs+2goWt+HXRoyxpg4Z4nAGGPinCWC2Fga6wA6Gfs+6th30ZB9Hw1F5fuwNgJjjIlzViMwxpg4Z4nAGGPinCWCDiQiA0VklYhsFpE8Ebkx1jHFmoi4ReQDEfl7rGOJNRHJEJHnReQTEdkiIlNjHVMsicjNof8nm0TkKRFpW5eaXZiIPCwi+0OjONasyxSRl0VkW+i1Z6TezxJBx/ID31fVUcAU4DoRGRXjmGLtRpxuyg38AfiXqp4MnEocfy8i0h9YCOSq6imAG2dMk3jxKDC70brbgFdVdRjwKhHstt8SQQdS1b2q+n5ovhjnP3rcjkYvIgOArwEPxjqWWBORdOCLOGN0oKpVqno4pkHFngdICo1emAzsiXE8HUZV3wAONlo9F3gsNP8YMC9S72eJIEZEZDAwHngvxqHE0v8DbgGCMY6jMxgCFAKPhC6VPSgiKbEOKlZU9XPgt8B/gL3AEVV9KbZRxVxvVd0bmi8AekfqwJYIYkBEUoEXgJtU9Wis44kFEfk6sF9VN8Q6lk7CA5wGLFbV8UApcTxiX+j691ycBNkPSBGRb8U2qs4jNKRvxO79t0TQwUTEi5MEnlTV/4t1PDF0OjBHRHYBTwNfEpE/xzakmMoH8lW1pob4PE5iiFdfBnaqaqGqVgP/B0yLcUyxtk9E+gKEXvdH6sCWCDqQOKN8PwRsUdX/jXU8saSqP1LVAao6GKcR8DVVjdtffKpaAOwWkRGhVTOAzTEMKdb+A0wRkeTQ/5sZxHHjechy4LLQ/GXA3yJ1YEsEHet04Ns4v34/DE1fjXVQptO4AXhSRD4GxgH/E9twYidUM3oeeB/YiHOuipvuJkTkKWANMEJE8kXkSuAe4GwR2YZTY7onYu9nXUwYY0x8sxqBMcbEOUsExhgT5ywRGGNMnLNEYIwxcc4SgTHGxDlLBCZmRCQQuoU2T0Q+EpHvi0in/psUkWwReS/UDcQXjvNYt4vIDyIVWxvf+8HWOjwUkdUiYgPHxwFPrAMwca1cVccBiEgO8BegB/CL4z2wiLhVNXC8x2nCDGCjql4VhWN3mK4ev4msTv3ry8QPVd0PLACuF4dbRH4jIutE5GMRuRpARFwicn+oz/6XRWSFiJwf2rZLRH4lIu8DF4jITBFZIyLvi8hzoT6eEJEJIvK6iGwQkZU1j+3XJyKDReS10Hu/KiKDRGQc8Gtgbqgmk9Ron10i8msR2Sgia0XkpOaO1Wi/E0Mx1ywPq1kOHfOO0GfYKCInh9ZnishfQ8d8V0TGhtbfLiKPicibIvKZiJxbL6Z/hbo4afBrX0QWi8j6UM3sjuP/1zRdjSUC02mo6g6cfudzgCtxepycCEwEvisiQ4BzgcHAKJyntBsP3lKkqqcBrwA/Bb4cWl4P/FfoRHgfcL6qTgAeBu5uIpz7gMdUdSzwJHCvqn4I/Bx4RlXHqWp5E/sdUdUxwB9xeldt8liNPve/gSOhRANwBfBIvSIHQp9hMVBzKekO4IPQMX8MPF6v/InAl4A5wJ+BVaGYynG6/W7sJ6qaC4wFzqxJKiZ+2KUh01nNBMbW/NoH0oFhwBnAc6oaBApEZFWj/Z4JvU7BSRZvO13V4CP0yD5wCvByaL0bp5vjxqbiJB2AJ3BqAuF4qt7r79twrAeBK0Tkv4CLgEn1ttV0Trih3nHOAM4DUNXXRCRLRHqEtv1TVatFZCPO5/tXaP1GnCTa2IUisgDnfNAX53v7uPWParoLSwSm0xCRoUAAp1dFAW5Q1ZWNyrTWN1NpTVHgZVW9uNH+Y4A8VY3WMJDazHxrXsBpG3kN2KCqRfW2VYZeA4T3f7YSQFWDIlKtdf3IBBvvH6pl/QCYqKqHRORRIG6GhDQOuzRkOgURyQaWAH8MnbhWAtfWu6Y9XJyBWt4Gzgu1FfQGpjdzyHeB0+tdp08RkeHAViBbQuMBi4hXREY3sf871A2NeAnwZpgf5aJ6r2vCPZaqVuB85sU0vCzUnDdDx0JEpuNcPmrP2BY9cJLnkdD3+ZV2HMN0cVYjMLGUJCIfAl6c8ZyfAGq6534Q5zLG++JcwynEGZrvBeq6aN6N0zvlkcYHVtVCEbkceEpEEkKrf6qqn4YuN90rzvCQHpxr+XmNDnEDzmhhPwy99xVhfqae4vQeWgnU1EbCPdaTwDeAcEbiuh14OPReZdR1T9wmqvqRiHwAfILzfb7dnuOYrs16HzVdjoikqmqJiGQBa4HTQ/35xzquXTiDrR9o5/4/ANJV9WcRDcyYVliNwHRFfxeRDJwG4P/uDEngeInIMuru9jGmQ1mNwBhj4pw1FhtjTJyzRGCMMXHOEoExxsQ5SwTGGBPnLBEYY0yc+/9S2oenmGPo1AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Generate some synthetic data with a non-linear relationship\n",
    "np.random.seed(0)\n",
    "x = np.linspace(-5, 5, num=100)\n",
    "# y = x ** 3 + np.random.normal(size=100)\n",
    "y = f(x) #np.sin(x * np.pi)\n",
    "\n",
    "# Split the data into training and testing sets\n",
    "x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=0)\n",
    "\n",
    "# Fit polynomial regression models with different degrees of polynomials\n",
    "degrees = [1, 2, 3, 4, 5, 8, 10]\n",
    "train_errors, test_errors = [], []\n",
    "bias_squared = []\n",
    "variance = []\n",
    "\n",
    "for degree in degrees:\n",
    "\n",
    "    # Transform the features to polynomial features\n",
    "    poly_features = PolynomialFeatures(degree=degree)\n",
    "    x_poly_train = poly_features.fit_transform(x_train.reshape(-1, 1))\n",
    "    x_poly_test = poly_features.transform(x_test.reshape(-1, 1))\n",
    "\n",
    "    # Fit the linear regression model to the polynomial features\n",
    "    model = LinearRegression()\n",
    "    model.fit(x_poly_train, y_train)\n",
    "\n",
    "    # Evaluate the model on the training and testing data\n",
    "    y_pred_train = model.predict(x_poly_train)\n",
    "    y_pred_test = model.predict(x_poly_test)\n",
    "    train_error = mean_squared_error(y_train, y_pred_train)\n",
    "    test_error = mean_squared_error(y_test, y_pred_test)\n",
    "    train_errors.append(train_error)\n",
    "    test_errors.append(test_error)\n",
    "    bias_squared.append(calculate_estimator_bias_squared(y_pred_test))\n",
    "    variance.append(calculate_estimator_variance(y_pred_test))\n",
    "\n",
    "\n",
    "# Plot the training and testing errors as a function of the degree of polynomial\n",
    "\n",
    "plt.plot(degrees, train_errors, label='Training error')\n",
    "plt.plot(degrees, test_errors, label='Testing error')\n",
    "plt.plot(degrees, bias_squared, color='blue', label='$bias^2$', linewidth=5)\n",
    "plt.plot(degrees, variance, color='green', label='variance', linewidth=5)\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel('Degree of polynomial')\n",
    "plt.ylabel('Mean squared error')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6939637643352092"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_error"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tf",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}