activity1.answers.2023.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1a67f85d-6cac-47aa-ab01-3bb1654d7099",
   "metadata": {},
   "source": [
    "###  Activity 1: All Possible Regressions\n",
    "\n",
    "This is a quick assignment.  I would expect it to only take a couple of hours to complete. The start of this notebook develops the example.  You will then parallelize the exammple using `joblib`. \n",
    "\n",
    "Run this on a machine that has at least 4 cores. Typically this is your laptop. The ugrad machines are also OK. Google Colab is not adequate. Fill out the requested cells toward the bottom of the notebook.\n",
    "\n",
    "**Due date**: September 8, 2023, 5:00 pm EDT.\n",
    "\n",
    "**Instructions for Submission**: Submit via Gradescope.\n",
    "\n",
    "#### Preparing the Environment\n",
    "\n",
    "You will need a couple of packages installed to make this work.  \n",
    "If you are running the Docker file, it should just work. Otherwise,\n",
    "stop your instance of jupyter lab and then run\n",
    "\n",
    "```\n",
    "conda install pandas pandoc\n",
    "jupyter lab\n",
    "```\n",
    "At this point, you should be ready to go.\n",
    "\n",
    "#### Example code\n",
    "\n",
    "This is a Python reimplementation of the Section 3.4. in _Matloff_, Parallel Computing for Data Science. It is based on data from https://www.kaggle.com/divan0/multiple-linear-regression. The notebook asks the question what combination of variables best predict the price of a house."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a9525a93-0b4e-4a25-bf87-d4335132cffc",
   "metadata": {
    "tags": []
   },
   "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>id</th>\n",
       "      <th>date</th>\n",
       "      <th>price</th>\n",
       "      <th>bedrooms</th>\n",
       "      <th>bathrooms</th>\n",
       "      <th>sqft_living</th>\n",
       "      <th>sqft_lot</th>\n",
       "      <th>floors</th>\n",
       "      <th>waterfront</th>\n",
       "      <th>view</th>\n",
       "      <th>...</th>\n",
       "      <th>grade</th>\n",
       "      <th>sqft_above</th>\n",
       "      <th>sqft_basement</th>\n",
       "      <th>yr_built</th>\n",
       "      <th>yr_renovated</th>\n",
       "      <th>zipcode</th>\n",
       "      <th>lat</th>\n",
       "      <th>long</th>\n",
       "      <th>sqft_living15</th>\n",
       "      <th>sqft_lot15</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>7129300520</td>\n",
       "      <td>20141013T000000</td>\n",
       "      <td>221900.0</td>\n",
       "      <td>3</td>\n",
       "      <td>1.00</td>\n",
       "      <td>1180</td>\n",
       "      <td>5650</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>7</td>\n",
       "      <td>1180.0</td>\n",
       "      <td>0</td>\n",
       "      <td>1955</td>\n",
       "      <td>0</td>\n",
       "      <td>98178</td>\n",
       "      <td>47.5112</td>\n",
       "      <td>-122.257</td>\n",
       "      <td>1340</td>\n",
       "      <td>5650</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>6414100192</td>\n",
       "      <td>20141209T000000</td>\n",
       "      <td>538000.0</td>\n",
       "      <td>3</td>\n",
       "      <td>2.25</td>\n",
       "      <td>2570</td>\n",
       "      <td>7242</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>7</td>\n",
       "      <td>2170.0</td>\n",
       "      <td>400</td>\n",
       "      <td>1951</td>\n",
       "      <td>1991</td>\n",
       "      <td>98125</td>\n",
       "      <td>47.7210</td>\n",
       "      <td>-122.319</td>\n",
       "      <td>1690</td>\n",
       "      <td>7639</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5631500400</td>\n",
       "      <td>20150225T000000</td>\n",
       "      <td>180000.0</td>\n",
       "      <td>2</td>\n",
       "      <td>1.00</td>\n",
       "      <td>770</td>\n",
       "      <td>10000</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>6</td>\n",
       "      <td>770.0</td>\n",
       "      <td>0</td>\n",
       "      <td>1933</td>\n",
       "      <td>0</td>\n",
       "      <td>98028</td>\n",
       "      <td>47.7379</td>\n",
       "      <td>-122.233</td>\n",
       "      <td>2720</td>\n",
       "      <td>8062</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2487200875</td>\n",
       "      <td>20141209T000000</td>\n",
       "      <td>604000.0</td>\n",
       "      <td>4</td>\n",
       "      <td>3.00</td>\n",
       "      <td>1960</td>\n",
       "      <td>5000</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>7</td>\n",
       "      <td>1050.0</td>\n",
       "      <td>910</td>\n",
       "      <td>1965</td>\n",
       "      <td>0</td>\n",
       "      <td>98136</td>\n",
       "      <td>47.5208</td>\n",
       "      <td>-122.393</td>\n",
       "      <td>1360</td>\n",
       "      <td>5000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1954400510</td>\n",
       "      <td>20150218T000000</td>\n",
       "      <td>510000.0</td>\n",
       "      <td>3</td>\n",
       "      <td>2.00</td>\n",
       "      <td>1680</td>\n",
       "      <td>8080</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>8</td>\n",
       "      <td>1680.0</td>\n",
       "      <td>0</td>\n",
       "      <td>1987</td>\n",
       "      <td>0</td>\n",
       "      <td>98074</td>\n",
       "      <td>47.6168</td>\n",
       "      <td>-122.045</td>\n",
       "      <td>1800</td>\n",
       "      <td>7503</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 21 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           id             date     price  bedrooms  bathrooms  sqft_living  \\\n",
       "0  7129300520  20141013T000000  221900.0         3       1.00         1180   \n",
       "1  6414100192  20141209T000000  538000.0         3       2.25         2570   \n",
       "2  5631500400  20150225T000000  180000.0         2       1.00          770   \n",
       "3  2487200875  20141209T000000  604000.0         4       3.00         1960   \n",
       "4  1954400510  20150218T000000  510000.0         3       2.00         1680   \n",
       "\n",
       "   sqft_lot  floors  waterfront  view  ...  grade  sqft_above  sqft_basement  \\\n",
       "0      5650     1.0           0     0  ...      7      1180.0              0   \n",
       "1      7242     2.0           0     0  ...      7      2170.0            400   \n",
       "2     10000     1.0           0     0  ...      6       770.0              0   \n",
       "3      5000     1.0           0     0  ...      7      1050.0            910   \n",
       "4      8080     1.0           0     0  ...      8      1680.0              0   \n",
       "\n",
       "   yr_built  yr_renovated  zipcode      lat     long  sqft_living15  \\\n",
       "0      1955             0    98178  47.5112 -122.257           1340   \n",
       "1      1951          1991    98125  47.7210 -122.319           1690   \n",
       "2      1933             0    98028  47.7379 -122.233           2720   \n",
       "3      1965             0    98136  47.5208 -122.393           1360   \n",
       "4      1987             0    98074  47.6168 -122.045           1800   \n",
       "\n",
       "   sqft_lot15  \n",
       "0        5650  \n",
       "1        7639  \n",
       "2        8062  \n",
       "3        5000  \n",
       "4        7503  \n",
       "\n",
       "[5 rows x 21 columns]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "#importing dataset using panda\n",
    "dataset = pd.read_csv('../data/kc_house_data.csv')\n",
    "#to see what my dataset is comprised of\n",
    "dataset.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "06ad55ee-ba79-4673-b693-ebca04f1bb77",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_262/3445377822.py:5: FutureWarning: In a future version of pandas all arguments of DataFrame.any and Series.any will be keyword-only.\n",
      "  dataset = dataset[~dataset.isin([np.nan, np.inf, -np.inf]).any(1)]\n"
     ]
    }
   ],
   "source": [
    "#dropping the id and date column\n",
    "dataset = dataset.drop(['id','date'], axis = 1)\n",
    "\n",
    "# clean out NaN and inf values\n",
    "dataset = dataset[~dataset.isin([np.nan, np.inf, -np.inf]).any(1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f724288c-e17e-46a4-939b-c5ec0fc34a3b",
   "metadata": {},
   "source": [
    "Let's first do a simple regression. How does square footage predict price? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cbd8b494-d600-4bc9-9af7-15bef7796fb5",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "coefficient of determination: 0.4928817196006995\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "X = np.array(dataset.sqft_living)\n",
    "Y = np.array(dataset.price)\n",
    "\n",
    "# shape X into matrix of a single column\n",
    "X = X.reshape((X.shape[0],1))\n",
    "\n",
    "model = LinearRegression()\n",
    "model.fit(X,Y)\n",
    "r_sq = model.score(X,Y)\n",
    "print('coefficient of determination:', r_sq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "97519567-fe0a-439e-a7e7-d97e4d56dea7",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAGFCAYAAABg2vAPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9u0lEQVR4nO3deXxU1d3H8e9MVhJIIKBNVISIaI0RESuLKI8iVGR1oa1afFq1LihWqa07BYprrWtRXCraimgfV0A0LYh1QRQrUo2xihhcIBEhkIRAQpi5zx/jYJZZ7p2ZO3Nn5vN+vfyDZM7khGDO957ld1yGYRgCAABpy53oDgAAgMQiDAAAkOYIAwAApDnCAAAAaY4wAABAmiMMAACQ5ggDAACkuUwzL/J6vdq0aZO6desml8tld58AAEAMGIahxsZG7bfffnK7gz//mwoDmzZtUu/evWPWOQAAED9fffWVDjjggKCfNxUGunXrtvfNCgoKYtMzAABgq4aGBvXu3XvvOB6MqTDgXxooKCggDAAAkGTCLfGzgRAAgDRHGAAAIM0RBgAASHOEAQAA0hxhAACANEcYAAAgzREGAABIc4QBAADSHGEAAIA0Z6oCIQAATuLxGlpdXafNjc3at1uuBpcWKcPNRXqRIgwAAJJKRWWNZi+pUk19896PlRTmauaEMo0pL0lgz5IXywQAgKRRUVmjqQvWtAsCklRb36ypC9aoorImQT1LboQBAEBS8HgNzV5SJSPA5/wfm72kSh5voFcgFMIAACAprK6u6zQj0JYhqaa+Waur6+LXqRRBGAAAJIXNjcGDQCSvw/cIAwCApLBvt9yYvg7fIwwAAJLC4NIilRTmKtgBQpd8pwoGlxbFs1spgTAAAEgKGW6XZk4ok6ROgcD/55kTyqg3EAHCAAAgaYwpL9G8KYNUXNh+KaC4MFfzpgyizkCEKDoEAEgqY8pLNLqsmAqEMUQYAAAknQy3S8P69Ux0N1IGywQAAKQ5wgAAAGmOMAAAQJojDAAAkOYIAwAApDnCAAAAaY4wAABAmiMMAACQ5ggDAACkOcIAAABpjjAAAECaIwwAAJDmCAMAAKQ5wgAAAGmOMAAAQJojDAAAkOYyE90BAEBoHq+h1dV12tzYrH275WpwaZEy3K5EdwsphDAAAA5WUVmj2UuqVFPfvPdjJYW5mjmhTGPKSxLYM6QSlgkAwKEqKms0dcGadkFAkmrrmzV1wRpVVNYkqGdINYQBAHAgj9fQ7CVVMgJ8zv+x2Uuq5PEGegVgDWEAABxodXVdpxmBtgxJNfXNWl1dF79OIWURBgDAgTY3Bg8CkbwOCIUwAAAOtG+33Ji+DgiFMAAADjS4tEglhbkKdoDQJd+pgsGlRfHsFlIUYQAAHMTjNbRq/Va9+MEmnXlMbxlSp0Dg//PMCWXUG0BMUGcAABwiUE2B7nlZkqTtO1v3fqyYOgOIMcIAADiAv6ZAx4OC9d+FgOmjDlHfXnlUIIQtCAMAkGDhagq4JD317pd68+qRhADYgj0DAJBg1BRAohEGACDBqCmARCMMAECCUVMAiUYYAIAEo6YAEo0wAAAJluF2aeaEMknUFEBiEAYAwAHGlJdo3pRBKi5svxRQXJireVMGUVMAtuJoIQA4xJjyEo0uK9bq6jptbmympgDihjAAAA6S4XZpWL+eie4G0gzLBAAApDnCAAAAaY4wAABAmiMMAACQ5ggDAACkOcIAAABpjjAAAECaIwwAAJDmCAMAAKQ5wgAAAGmOMAAAQJojDAAAkOYIAwAApDnCAAAAaY4wAABAmiMMAACQ5ggDAACkOcIAAABpjjAAAECaIwwAAJDmCAMAAKQ5wgAAAGmOMAAAQJojDAAAkOYIAwAApDnCAAAAaY4wAABAmiMMAACQ5ggDAACkOcIAAABpjjAAAECaIwwAAJDmCAMAAKQ5wgAAAGmOMAAAQJojDAAAkOYIAwAApDnCAAAAaY4wAABAmiMMAACQ5ggDAACkOcIAAABpjjAAAECaIwwAAJDmCAMAAKQ5wgAAAGkuM9EdAIBk4/EaWl1dp82Nzdq3W64GlxYpw+1KdLeAiBEGAMCCisoazV5SpZr65r0fKynM1cwJZRpTXpLAngGRY5kAAEyqqKzR1AVr2gUBSaqtb9bUBWtUUVmToJ4B0SEMAIAJHq+h2UuqZAT4nP9js5dUyeMN9ArA2QgDAGDC6uq6TjMCbRmSauqbtbq6Ln6dAmKEMAAAJmxuDB4EInkd4CSEAQAwYd9uuTF9HeAkhAEAMGFwaZFKCnMV7AChS75TBYNLi+LZLSAmCAMAYEKG26WZE8okqVMg8P955oQy6g0gKREGADiKx2to1fqtWrR2o1at3+qo3fljyks0b8ogFRe2XwooLszVvCmDqDOApEXRIQCOkQwFfcaUl2h0WTEVCJFSXIZhhI3dDQ0NKiwsVH19vQoKCuLRLwBJIJZlef0FfTr+QvK/G0/egHVmx29mBgBEJJZP8eEK+rjkK+gzuqyYJ3DABuwZAGBZrMvyUtAHSCzCAABL7CjLS0EfILEIAwAsseMpnoI+QGIRBgBYYsdTPAV9gMQiDACwxI6neAr6AIlFGABgiV1P8RT0ARKHo4UALPE/xU9dsEYuqd1Gwmif4inoAyQGMwMALOMpHkgtzAwAiIgdT/HJUI4YSEWUIwbgCJQjBmLP7PjNMgGAhLOjkBEA8wgDABKOcsRAYhEGACQc5YiBxCIMAEg4yhEDiUUYAJBwlCMGEoswACDhKEcMJBZhAIAjUMgISByKDgFwDMoRA4lBGADgKBlul4b165nobgBphWUCAADSHGEAAIA0RxgAACDNEQYAAEhzhAEAANIcYQAAgDRHGAAAIM1RZwAA4sjjNSiqBMchDABAnFRU1mj2kirV1H9/FXNJYa5mTiij3DISimUCAIiDisoaTV2wpl0QkKTa+mZNXbBGFZU1CeoZQBgAANt5vIZmL6mSEeBz/o/NXlIljzfQKwD7EQYAwGarq+s6zQi0ZUiqqW/W6uq6+HUKaIMwAAA229wYPAhE8jog1ggDAGCzfbvlxvR1QKwRBgDAZoNLi1RSmKtgBwhd8p0qGFxaFM9uAXsRBoAk5fEaWrV+qxat3ahV67d22nwW7vOInwy3SzMnlElSp0Dg//PMCWXUG0DCUGcASELhzqtznt15xpSXaN6UQZ1+LsX8XOAALsMwwj4uNDQ0qLCwUPX19SooKIhHvwAE4T+v3vF/XP8z5YUjSvXQ69VBPz9vyiAGngSiAiHiyez4TRgAkojHa+i421YEPabmkuRyScFWBFzyPYm+efVIBiAgDZgdv9kzACQRM+fVQ20N4Dw7gEAIA0ASidU5dM6zA2iLMAAkkVidQ+c8O4C2CANAEjFzXj3UVgDOswMIhDAAJBEz59UvOL7Ut5EwyOc5zw6gI8IAkGT859WLC9tP9RcX5mrelEG6dmxZyM9zrBBARxwtBJJUuPPqnGcHYHb8pgIhkKQy3C4N69cz4s8DgB9hAEDUmIUAkhthAEBUuAcBSH5sIAQQMf89CR2rItbWN2vqgjWqqKxJUM8AWEEYABARj9fQ7CVVnS5EkrT3Y7OXVHF1MpAECAMAImLmngTuQQCSA2EAQETM3m/APQiA8xEGAETE7P0G3IMAOB9hAEBEzNyTwD0IQHIgDACIiJl7ErgHAUgOhAEAnXi8hlat36pFazdq1fqtQU8EhLsngToDQHKg6BCAdqwWERpTXqLRZcVUIASSGBcVAdjLX0So4y8F/7DO0z6QXMyO3ywTAJBEESEgnREGAEiiiBCQzggDACRRRAhIZ4QBAJIoIgSkM04TADbzeI2k2GnvLyJUW98ccN+AS74jgxQRAlIPYQCwkdVjeonkLyI0dcEauaR2gYAiQkBqY5kAiEKo4jz+Y3odN+XV1jdr6oI1qqisiXd3wxpTXqL7zj5KPfKz2n2cIkJAamNmAIhQqKf+0WXFIY/pueQ7pje6rNhRT9oVlTWas/Rj1TW17v1YUX62Zoxz3kwGgNhhZgCIQLin/rkrPku6Y3rBvqdtTbt16UJnzmQAiA3CAGCRmeI8j75Vbeq9nHJMj4JDQHojDAAWmSnOs31na9DPt+WUY3oUHALSG2EAsMjs03xedkanq339XPLtL3DKMT0KDgEJsHGjdMIJkssl9e8vXXmltHt3QrpCGAAsMvs0v3O3Z+9mwbaceEyPgkNAnOzaJV16qS8AHHCA9Nprvo9/9pl0553S4sUJ6RanCZDyQhX9iaQgULjiPG31yMtSTqZbtQ0tez9W7MA6AxQcAmxkGNK8eb4QEM6339rfnwAIA0hpoY7/SYqoIFCG26WJR5bowdfDbxLctrNVT/xqiNwul6MrEFJwCLDBq69Kp5witbSEf60klZdLkyfb26cgXIZhhN0ebPY+ZMBJ/EflOv4D7zjYdfycpJAFdioqa3TxgjWm+3HPmQM1aeD+pl+fSMlUMRFwpM8/l376U+m996y1+8lPpPnzpa5dY9ods+M3MwNISWaOygXS9hhdoIJAHq+hWYs/stSXZFpnH1NeotFlxUlxlwLgGI2N0rRp0t/+Zq1dWZn07LPSD39oT78sIAwgJYU7KheO/xjdsH49O71v2/X/cJx0YsCsDLer0/cNoAOvV7rjDumqq6y3XbpUGjs29n2KAmEAKSkWR+BqGzq/h9X3ZZ0dSDEvvSSNG2e93Z/+JE2fLrmdeYiPMICUFIup+bodnWcArLzv9FH9WWcHUsHHH0unny7997/W2p17rnTvvTHfB2AHZ0YUIEqDS4tU1OHmPauK8rMDvm9xQU7YtiWFuZo2sn9UXx9AAm3b5tvZ73L51vbNBoFjjvFtIjQMWzYE2oUwgJSU4XbptCh38BcXdgn4vrMmHh62LcsDQBLas0eaNcsXAIqKfJv7zOjSxXeM0DCk1aul0lJbu2kHwgBS1qiy4ojbhtr4N6a8RA9MGaTueZ1nHnrkZemBEMcSATjQ00/7AkBWljR7tvl28+b5NhLu3OkrK5zE2DOAlGWlUmBbLoV/sh9dVqxuOVlauf5bbdrerP17dNGx/Xpp6EE9mREAksH770sTJ0pff22t3bRp0h//6JsNSCGEAThSJGWCA7U/pbxY81duCFhVz5DUPS+r3Q2DbQvsBOtDsMI8R+xfSBAAnGzzZukXv5AqKqy1O+EEacECaf/kKB4WCcIAHCfaKniB2rtcvuU8P//9AMEK7ATrw8QjS/TQ69WdZhpq65s1dcGakJULASTA7t3Sddf5agJYsc8+0pIl0pAh9vTLYShHDEcJVUJYCl0mOFR7v/OH99WosuKQMw3h3iMY/2U+b149khkCIJEMQ/rrX31H+6x6/HFpypTY9ylBzI7fbCCEY5gpITx7SZU83sDDdKj2km+wfqmyNmQQCPceoRj6vnIhgAR4+22pZ09fYR8rQeDqq32XCRlGSgUBK1gmgGOEKyHcdrANVC432vZm3sOMWFQ/9It274STWPleUun7hs2+/lo6+2zpjTestRs/3lcHYJ997OlXkiEMwDGWVdWael2wwdbsIBzqdbEYyGN1MZFdNwgmYqC18r1wcyLC2rVL+s1vpAcesNaub19p0SJpwABbupXMCANwhIrKGs1fucHUa4MNtmYH4VCvi2Yg9+8ZiMXFRMH2LUS7UTERA62V78Wu7xspwDCk++/3He2z6plnpDPOiH2fUgh7BpBw/nV6M0IVA/LXFQj2jOsK097Me4QTi8qD0e6dCMY/0HZcBvEPtBWVNZF1OAQr34td3zeS3IoVUk6Obx+AlSDwhz/4KgoaBkHABMIAEs7KOn2owTbD7dLMCWWSFHAwNyTNGHdYyME63HsEU5SfFbOnVit7H8xK1EBr5Xux4/tGklq/Xho0yHcm+KSTfMcDzfjZz3x3ChiGNGOGlJFhbz9TCMsESDiz6/TnD+/babDtuP49uqxY86YM6jQV7jdn6cdyu10hB+0x5SUh36OjnvnZWnXtScrOjE22jsXeh47MDrSPraxWr2456tU1RzKkLU0tUe0rsON7ieUGTThIY6N06aW+o31WHH647w6BQw+1p19pgjCAhDO7Tr9f9y7yeI29g1Ko9e8Z48p0ycI1nd7D7NrzmPKSvQWJllXVBq1iKEk3nVYesyAgxWbvQ0dmB9A5Sz8O+PFI9xXY8b3EaoMmHMDrlW6/XbrmGuttX35ZGjMm9n1KUywTIOHMrtPPWfqxjrtthSoqa0Kuf1+8YI2ue+HDgO9hZUo8w+3SsH499fsJh+uBKYNUXNh+ECouzLVlQ9vg0qKAlyC11SMvy9JGxWgH0Ej3FVjZxxGLPR9IEi++6FsCyMiwFgTuvFPyeHzLAASBmGJmAAnnX6efumBNp6fvjvyDUn5OZsj177b3DQR6Tbh6Ax21nSlwwtl3qyv7kV7a1PbrueQLUaPLik1/36F+tv53aLsPxMprkWSqqqTTTpM+/dRau/POk+65R+ra1Z5+QRIzA4gTj9fQqvVbtWjtRq1av7XTU7l/nb7j03dHxnf/7WjZE3Wf2k6dh+uf9P1MwaSB+2tYP/tuJ1xdXRcyzEi+sGNlI12kGyPbinQDX7CfbaCZFSuvRRKoq/Pt5He5fGv7ZoPAkCFSdbVvBuCRRwgCccDMAGxXUVmjWYurVNvw/eBbXJCrWRO/X4P2eA0VdsnWVScfqjVfbtPjb39pe7/8U+dOK3Jjx6Y7yfrGyFh9Xf/XNjuz4rRZGFi0Z4/vWN+cOdba5edLL70kjRhhT78QEmEAtqqorNHFCwJs5Gvwre0/MGWQJEU9QFnRtjiQE4vc2LHpzq/jQPvGp1v0zBpr97lHuv/AP7MS69fCIf7v/3xH+6x68EHpggt8swdIGMIAbOPxGrrmucAb+fyu/L//aOduT0Tr2KHkZ2do526PpOBrz5JCnr2PZI08FsKt70db6dA/0Hq8hm55KfDpgUBiWWERKWLNGmniRGnjRmvtLrtM+uMfpVxOhjgFewZgm7c/3xp27bvJhiAgSVmZbt13due158K8LF0x6pC9T8dWityY2VcQC6HW961upAvV59XVdaptaDHVJzbwYa9vvpFOPtn3JH/00eaDwMiRvtcahnTvvQQBh2FmALZZtX5rwr729p2t6pGfrTevHqm5K9bp0ZUbtH1Xq7bvbNVdyz/VU+9+qbHlxabea3Njc9z3FQRb3y+28DXD9dnK2n+or8sNg2mgpUW67jrf0T4r9t1XWrJEGjzYnn4hZggDsFFia8gvq6pV/a7dunv5uoB7Ah4xeTHShi1NQd/Dzn0F/vX9tz/f+l2wMjTsoF4aamIt3cxeCLNr/zPGHaZfDi8NOMA7bfMlYsgwpEcflc4/33rbxx+XpkyJfZ9gG5dhGGF/Yzc0NKiwsFD19fUqKCiIR7+QAlau26KfP/JO1O+T12b934qi/GxlZ7hMT4UH4nZJ+3TN1jeNgWuj+9fR37x6pC1Pw5EMth6voeNuWxF0CcTf59d+d6L+5/ZXw+5NCPa9BQsc/ldyFDBJvfWWNHasVF9vrd2110qzZ0tZoQtmIb7Mjt/sGYBt6neF3i8g+ZYdQw2hLimiICBJdU27owoCkuQ1FDQISPZenhPpLYNm90K898W2iPcmcMNgivnqK+m443z/Qw4fbj4ITJggffutbxbh5psJAkmMMICY8m9Ye37N17phUWXY1//quL6SggeCZBlKQq2/R7LxMJrB1kqdgkiL/CTrDYPx2gSaFHbulC6+2BcADjxQWrnSXLuDDpI++MAXABYvlnr1srefiAv2DCBmAk1phzPyh8U6uk9Rp3Zul++pPFkEW3/3FVz6qN0MRXFBjmZNPDzkFLqVwbbjeXyrdQoiKfJjV2EkO7G/Qb4BfO5c6de/tt72ued85YSRkggDiIlg68fhbG5s1qSB+7cbjLY0tgS9Pc+KLllu7Wr1Rv0+4XQPcmlQ8IJLLXsLLgUbhKIZbCOpU2C1yI+dhZHs4MTiUnH1yiu+i332WCzjfeONvouEMjLs6Rccg2UCRC3UlHY4/sGibd3/Xt1yYtKveE0sBHp+NlNw6ZrnPgw6TR3NYBvLOgXBJNMNg2m7v2H9eumoo3zLAKNGmQ8CZ54pbd/um0W4/nqCQJogDCBq4aa0Awk1WMTqabK51auifPs3NG0LcGnQ2+vDF1zavrNVbwepxRDtYGv3hT/xCByxkqz7GyLS0OA70udySQcfLK1da67dEUdIn3ziCwBPPikVFtraTTgPywSImtV14XCDRbTX7bZ12sD9NX/lBttnCTr+Haz6fIupdqs+36Lh/TtvwLJ69W8gdl/4E4vCSPGQjPsbLPF4pNtv9x3ts8Ltll5+Wfrxj+3pF5IKYQBRs/ok33GwCFTBLthAaNWosmIdU9p5g2Ksdf47MDvghh7MAw22hXlZOvfYUo0uC19B0e4Lf5LhhsFk299g2uLF0qRJ1tvddZdvA6GbiWF8jzCAqJl5ki/Kz9KM8YeruKD9YBFqh3e01+26XdLRfXooO9PdbsBa980OzX31s4jes6Ngl/cM69fT1NcIN1D7B9tAJZWfXP2Fzhp8oPr2yk/oIOz0Gwbtvvgprj76SDr1VOkzi/9+f/Ur6e67fdcEAwEQBhC1UFPafr8Y1lcTj9xPkvYOyhu2NOmu5es6vba23ne98XnD++pPk4/UuxvqdPcrnV8XjteQ3vtim4b169luwFq1fmtMwkCo6fqhB/VU97yskPsGeuRlaehB4QfRZVW1gcshN7S0+/tLxmNy8bjXIBZLLgm1datvMH/hBWvthg6VnnpK6tPHlm4htVCOGDETrs5A9zzfZr5wG+s6imap4J4zB2rSwP3bfcxfrtfsnoRgNQ/CDb7Bjhb6hTpa2LGvZmZHkq0McLzP/SdVnYE9e3ylfW+80Vq7bt2kpUul44+3p19IOmbHb8IAYsrjNTR3xWe6a/mnie6KJOnJC4YGnML2nzuXAgeN6aP6q2FXa8jLjO4/e5DGDgg9iPiKDlWptiGyAWjV+q066+G3w77Oz+67EmIlUfcaOP6Gxaeeks46y3q7hx7yzR64HPS9wBHMjt8sEyDmnnr3y0R3QZJU2CVTXq8hj9fo9As/2OY8/0A9uqxYx922Iuh7uyTNWVqlk8uLbd3Rb3WHe6jKhE4R7ty/S75z/6PLQv/dRsKR+xvWrPHV+N+0yVq7yy+Xbr1Vyk2yjY9wJMIAYiqSmgN2qd+1Rz9/5J2gT+KhBupV67eaPps+uLQo5GAfzQAU6Q53Jx+Ti6bUcsqorZX+93+lZcustTvpJN/1wCUOW9ZA0iMMICwrU6tOHIRqQpScDTZQm/0+llfV6jf/t9a2dehIay44+Zhcyp/7D6alxVfa9+67rbUrLpaWLJF+9CNbugVIhAF00HHg39a0W3OWmt905dRByJA0a/FHpqeezX4fgfYUxLLevZmTGm0lwzG5lD33H4hhSPPn+9bzrXriCenss2PfJyAAwgD2MnvrYKjBLpbVA2OttqFFc1d8pstH9Q/72m1NLWFfE2xwblvvPhbr3sH2NwTqj+TwY3JKsXP/waxcKY0d6ysPbMV110mzZklZ9pfRBtqiBBUkfb+728x6f7DLXTxeQ2+v36qBvQtNBwHXd/+NH1Ci/Gz7L0S5a/mnqqis2fvnQPfbe7yGqVsTw32Psax3P6a8RG9ePVJPXjBU95w5UNNHHaLigtjcOxDo78BOyXSvgSVffikde6xvR/9xx5kPApMmSd9+65tFuOkmggASgpkBRHTrYMdNXhWVNbrmuQ9D1hDo3iVTe7zSjpbvb08rLszVxCNL9NDr1XGbSfA/sS+rqg14muDMY3rHbBNkzfZdMXkfqfP+hmkjD476mFyizt4ny70GYTU1SdOnSw8/bK1dv37S88/7LggCHIA6A7B8lr2t84b31eDSopDFdSTfk/+/N2xrd96+KD9Lf5hwuG56+b9xP4EwfVT/gFX9or0LoaNuuRm6ffKRjhnc2u4JCVYBMp7Fixx/7j8Qw5DuvVe64grrbV94IbL7BIAIUXQIpi1au1GXP7U24vbdu2Rq+y6Td6W3EeuB14ruXbK0fZe1SoiRcin2A2vHQfToPj303hfbQg6qZveE+PucDMWL4mrZMumUU3y3BFpx883SVVdJGfYvgwEdUXQIpkW7azuSICAlLghIilsQ8ItlEZ1Ag3rHkskdp/qDVfwLJi3O+puxbp00ebL0wQfW2p19tnT//VJhoT39AmKMDYTYu7s7XfjvSIiXtgNrtIJt9Oy4589/4qOisiaiPSF+8TjrH+8NjGE1NEg//7lvI+Ahh5gPAgMGSJ9+6ltGeOIJggCSCjMDaSTY+qx/d3e4df9UMbhvD/2zanPcv+7mxuao1sg9XkOzFn9kalBvW9a3W05WxHsy7D7r75jLgzwe6bbbpOuvt9YuI0N6+WVp9Gh7+gXECWEgTYT7pTu6rFiTBx2gZ9Z8ncBexkcigoAkbdiys9MNhFYGvrkrPlNtQ/j6B37+GYlVn2+x3Nd4nPUPtnQRy6JNYS1aJJ16qvV299wjXXYZFwMhZbBMkAaCTS37f+ne8lKVjrttRVoEgURwybc0cffyT4P+DNrWPgikorImipsgrQ1Y8TjrH+6yIqlzHYuYqayUDj7YN5BbCQIXXCDt2OFbBvj1rwkCSCmEgRQX7peuIenB16sdc7mQkxTlZ1scRjvzn5ho9XgjHvj8P8NIDevXUyWFuaa/l0iLF1lh5bKimNi61Xekz+Xyne1fv95cu2OPlb74whcAHnpIys+PTX8Ah2GZIMU56RbBZNI9L0s3TirXpQvN3QkQTJfsDLldrnaFljoKt3M/0p+hf6p/6EE9g95v4A8IV4w6RH175cXtrH9cLitqbfWV9r35ZmvtCgqkpUt9VQSBNEEYSHEpd/NbnGzf2Sq3W7pwRKkefqNa4atxBLZzt/kz6cF+VpH8DDtO9Tut4p+tlxUtXOg7DWDVX/4inXce0/+IK6cU3iIMpLiUuPktQa74+1o1t3rj9vWC/awi+RkGGuT9G0Wd8Isn5pcV/fvf0oQJUm2ttY5Mny7dcouUk2OtHRADjjlNI8JAynPyLYJ2y8pwqdUT+XcdryAQbuAz8zMsLsjRHT8dqC07WkIO8h3vN7AqVk8xoa5mNr2BsaZGOucc6ZVXrH3x0aOlv/1NKi622GsgdhxxmqYNyhGnAf8/OimxVf8Q3PRR/dW3V37IUsKBfobxvEfAjqcYy+/Z3Cxdc43vaJ8VJSXSkiXS0UdH1E8gljxeo9Mx47ZiWQ6cuwnQjpW69IifTLdLXXMz2932GGwwTOSUYrCnmFiEkbCzDYbhW8+/8ELrb75woXTWWRH1C7CL2cvhnrxgaNTlwLmbAO3414sfW1mtOUs/TnR3UsbkQfvrmTUbI26/x2t0uvY52DRhotb8wx1P9Vc6jPTuhaBLF2+8IY0d6zvbb8UNN0i//72UFd+y04BZcTlNYxF1BtJIhtuljdt3JbobKcEl31P5bZOP1ANTBnW626GkMFcXjSiNqE5BqNoD/oFz0sD9Naxfz7hs/otrTYAvvpCGDfPt6B8xwnwQOPVUacsW3yzCnDkEATiaradpIsTMQBqpqKzR/JUbEt2NlGBIOqXc95Q+uqw46BP71h27I5o5cNKtgbY/xTQ1SZdfLj3yiLV2/ftLzz0nlZdH9nWBBIn5aZoYIAykiWir2KGz+Ss3aP7KDSHX7fNyovtfzAl1Imx5ijEM3ybA6dOtd2jRImniROvtAIeIyWmaGGOZIE1QidA+oe4X6FOUF9V7O6FOhP8pJtivJf+SiamnmH/+U3K7ff9ZCQK33CLt2eMLEQQBpAB/IbDiDkuM8SgHHggzA2nCCU+YoWS5pUiO9RflZ6muqTX8C2PULpBQm+jOGdZXN730sazet5OIacJgon6K+fRT6YwzfBcEWTFlinTffb7ywEAKclIhMGYG0kS8njAzLPwbnnbiwbrrp0dqxrjDdMvpAyx/rVPKizVj/OGW2hR2ydQT5w+x3C6cYJvosjPduuD4UkvvlahpwlAsP8XU1/uO9Llc0qGHmg8CAwdK69b5ZgAef5wggJSXiE3BgTAzkITMVIHr+Jqj+/SISyVCjyFluqQ9Jr5IVoZbf/zHJxEvX6yurtOUIX0stbntjAEa3r+XVq3fGtHXDGdzY3Onv/urxhwmSXr4jep2MwRul3TSYfuqcmND1PcFxKO+edinGI9HuvVW39E+KzIzpZdflkaNiml/AZhHGEgyZgrPBHpNUX6WBh3YPS77BsxWAL5r+adRfZ2tTbt19/JPTN8qeMVJ/TW6zFeC1q4yzRu27OxUWcz/8/nvnB/q8VUb9EXdTvUpytM5w/oqO9Md9UAez2JEAWsCvPCCdNpp1t/s3nuladO4GAhwACoQJgmP19DcFet01/J1nT7XtgqcpICV4qIRbY1/J2k7SAarqhcJl6TCvCzV72yNaZW+cEHBzsqAIX34oe9s/+efW2t30UXSnXdKedFtrARgDuWIU0hFZY1mLa5SbUPwp3r/hjPDMFTb0BK/ziWpB74bJGNZprl7XlanaoJ+kdQaD/fEH8/65h6voTXvrdMBv52mkteXWWs8fLj05JNS795R9QGAdYSBFBHLp1d8r3telt67YbQy3K69T9+1Dc3a3NCse1esU1OLx/R7lRTm6sxjegectenIbK1xM0/8hV2y7a9v3tqqz6f+Rgc9MtdSs/qcfH388BMaes6kdh93yt3tQLrgboIUEKomPKKzfWer5q74TJeP6q8Mt0uDS4u0urpOn21uNBUECnIzNfnoAzS6rFiDS4v04gebTH1dM0c8zd4FcNWYH8bsa3byxBO+o32SDrLQ7HenXK6njxgll8ul4q9z9abX2DvYO+nudgDtEQYcjEJB9nr0rWpNG3mwllXVWl4qmD2pXKcdtf/eP5s9urnum0atWr815BOx2bsA6naYWw4yfaz03Xel8eOlzZvNvf47Dx9zqm4f8Qvtzvz+PoCO5ZSddnc7gPYIAw7m9EJByc4/O3D38k8tz74UF7QfYI/u00NF+dmqa9odst3cV9dr7qvrQz4Rm/25F+Vnh9ynIEk98rJCFy7atMk3A/Dqq6a+pt9rpYN05bjp2pLfI+Tr/Ect7bz1EED0CAMO5oRStKnu0ZXVloNAxwHWP/0dLgi01fGJuO1a+pZGc0/8W3bs1u49ocs2tnq8WvyfTSouaLM+39wsXXWV9Oc/m+6vJNV07alfnTFDHxUfbLrNvt1yLd16mOhLmYB0RRhwMLvOwuN723dZL0nc9mcR6QbPtk/EXq+hOUs/bjdgul0KWcLYJemmlz4O+3V2tHg0/e9rJcPQxZ++omteuNtiT6VpE6/Si4eNsNSmbTnlWO6nAGAPwoCDhaoJj8TZvrNVf35lnY7pW6Srn/0w4p+L/4n4koXvd/pcuLsMzH7NwV9V6tGnZym/1eJAO2OGPDfM0HF3vK7aMHspwt1X4MS72wG0RxhwoI7Hr+47e5DmLI3NWXjExt2vhD9GGAvhZggCOaD+G81ddJsG1lis8Hj66dLDD0tFviWQDCnsBUUXjijV4v/UhCyn7MS72wG0R50Bhwl2/Gr8gGI9/MaGxHUMCTNj3GEqys/W7xd/pMbmPQFf02V3s2Ytf1A/+9BiQaBDDpGef14qKwv6EjPFj8LVDvAvp0iBQwWnCQB7UHQoCVFgCIGcP7yvRpUVdyow5DK8Ou/dRZrx6iOW3/Ptux/V0Mt/afr1sSgWRJ0BIP4oOpRkKDCUmlwu32280Xh+7UaVH9B9759HfP6e/vb0TMvvc8sJv9RDg0+X4XLryfFDLbUNeEGRCR1DxGu/O1HvfbGNCoSAwxAGHIICQ6kp2iAgSXVNrfJUfax/PHKJDt3ypaW2zx5+omaOnqodOb6LgVzyPY3HY30+1EzApIH7h2gJIN4IAw7BsSp0VNC8Qzf/4z6N/+8bltrVH3aEJh7/a33ZoyTkLn87UXEQSC7uRHcAPhyrgiS5vR5dtvJJbbhtvD6450zTQaAlI1MLb3pEMgwVVn2gay8bp+LC9v+migtz4zIIh6s4KPnqK3isHpMAYBtmBhzCf/yKpYL0dPInb+nBF2623G7G6Iv10nGnac6pR+hn5cVatX5rwtfnqTgIJB/CQAJ13Fw1fkCJHn6jOtHdQpz8cHO1Hn7uRvWu/8ZSO+9FF+ndaderdo9bY7vlalZpkZZV1eq421Y4Yn3e7JIXS2OAcxAGEiTQ5ioXm6pTXtHOet3+0t06af27ltq93/cIDXhtqTIO7C23pCFtPue09XkqDgLJhzCQAMF+ecdi5zmcJ8vTqt+88YSmvvOMpXbbc7vqvDNmas0Bh0mSnmzN07AOr3HijYBUHASSD2HAJsGKtFBPIH2c+tGruvvFOyy3++3YK/RM+UmdpooCTas7cX0+1J0a8TzRAMA8woANQp2vLuySzSbBFHbkpk80/5nZ6rmrwVK7Bwefrj+NOEetGVlBXxNoWj2a9flYVBUMZkx5ieZNGdTp/4OO9xYAcAbCQIyFW78dU16ckH7BPvs2btU9L/5Jw7780FK7f5Uerd+Ou0Jb8nuEfW2XLHfAaXWz6+4btjS1+3M8SgOPKS/R6LJi2wIHgNghDMSQmfPVL1fWxrNLsElOa4uu+9d8/WLNUkvtNnbbRxeecYM++kE/S+3+55B9Ag6ig0uLVFyQq9qG0DMET67+UtNG9leG2xXXDYeRljEGEF+EgRiipHCKMwz9fO3Luumf91tueunEq7X0sOMj/tLnDO0b8OMZbpfOGnyg7loe+rri2oYWra6u0+DSIsdtOASQeISBGFpexVN/Khry5Yd67OlZ6rKnxVK7e449S/cOP1Med0ZUX797XpaGhni67tsrz9T7bG5sduSGQwCJRxiIEY/X0PNrNya6G4iRA7bX6r5Ft+nI2nWW2r10yLG6dsxlqu/SLWZ9ufX0I0I+pVs5109BIACBEAZiZHV1neqaWhPdDUQhb/cuzV72oH5SudxSu8+KDtBFp12v9b16x7Q/3btk6tYzBoRdv7dyrn91dZ2pr01BICC9EAZihCWC5OQyvDr/3Rd0w6vzLbc9d/JMvdrvGBt6JQ0rLdJPBx+owi7Z8niNTjMDHY8FzhhXpksXhj/XT0EgAIEQBiLU9pdxr645eu79rxPdJVhwwvp/67FnZllud/MJ5+rhwafJcNl74eeq6jqt+u4pvuORv2DHAi8cUapFa2vanSzoeK6fgkAAAiEMRCDQL2M4X7+tX+mB529W/61fWWr3TPlJmjnqIjXlmNuoF2s19c26eMEaPTBlkCQFPRb44OvVKuzS/n9pI0CNawoCAejIZQT6bdFBQ0ODCgsLVV9fr4KCgnj0y7GCndGGMxU079DNFXM1/pM3LbX7oPhgXTrpGn3V3TlForrnZSk3063aBmunGiTp/rOP0tgB+7X7mJ0VCAE4g9nxm5kBC7hXIDlkeD2a9tbfNX3lQkvtmjOzde7kWVrVZ4BNPYvO9p2Rb1Cd9uT7miuXxg74/qmfgkAA/AgDJnm8hua/+TlLAw52yn/f1LxFt1pud8OPL9GCgaek9B3SXkO6ZOEaPeCO73XGAJIDYcCEisoaXfPch1E9mcEeh23+XA8/O0cHNHxrqd3fjhqnm048Ty1ZOTb1zJmoLgggEMJAGBWVNbp4wZpEdwNtFO2s15+W3qWRn//bUrt3epfr8vG/VW1BL5t65nxUFwQQCGEgBI/X0KzFHyW6G5CU5WnVla8/rotXP2epXV2XAp1/xu/1/v4/tKlnycdfXTCaDYRsPgRSC2EghNXVdRHt3EaMGIZO/2iF7lx6l+WmV46drmfLR6b0PoC2euRlaZvJZax9u+VGdYVxoLZF+dm6cVJ5uw2KAJIHYSAE6rMnxsBNn2j+M7NVtKvBUrsHBp+uO0aco9aMLJt65kwzxh2mXw4v1T8qazTtyfflDXLcxV9dcFtTiy5d+H5EVxgHO1pb17Rblyxco4u+LtW1Y8ui/ZYAxBlhIATqs8fPDxq36N4lf9KQryottVtx0I/0u7FXaGt+d3s6lgR6dctRhtulsQP201y5dMnCzntc/PMjM8aVac7SyK4wNnO09sHXq3XkAd071TQA4GyEge90XAM9uk8PeQ1DmW5pjzfRvUtNOa0tuv7V+frf95daavd1wT668PQZqvrBQTb1LLm0Da1jB5ToAXfw6oKFXbIjvsI43PXHfjcsqtTJ5SXsIQCSCGFAgddA3S4FnW5FFAxDU95/STcum2e56dRJ1+jlHx5nQ6eSU7BLhcaUl2h0WXHADX6LTF6zHWiJzOyyWV1TKycWgCST9mEg2BooQSC2hn3xgR57eqZyPNZqNdw1/GzNPfZn8rgzbOpZfORnZ6hpt8dyu+55War/bmOglUuFglUXNLv0Feh1VpbN2G8DJJe0DgOUF7ZX7+21um/RrRpQ+5mldi8eepyuGzNNDbldbepZ/J15TG/NX7lBkkz/e3NJuvX0IyQpZpcKRXOF8eDSIhXlZ6uuaXfYr8N+GyC5pHUYMLsGCvPyW3bqD8sf0BmVKyy1W9ezty4+7Tqt79nbpp4lVkGX7IA3BQbTIy9Lt5x+xN7BPti0v1XRXGGc4XbpxknlATcotlUSJEwAcK60vbVw9x6vpvzlba3esC3RXUl6LsOrC1Y/r+v+9ajltr+cPFP/6neMDb1ylq45mVozY7Qy3K69g/rn3zbpb6s2tKsP0DUnQ8f376UpQ/pqaL+etm3Ci6bOwC0vVenB16sDfs4lhTyaCCC+zI7faRUG/CcGHnp9vV79xFote3R2wvp39dgzsy23u/HE8/TIMafKcLlt6JVz5WS6dc+ZA9sNlP5/k8uravX82o2qa/o+GJgdnCMVropgqM+/9MEm3bCoMq79BWAdYaANj9fQ3BXr9OjKDdq+i8uGotFv61d68LmbdHDd15baPV0+SjNHX6Sd2V1s6lnyeKDDk3OwTaz+YTkRT9pmZg4oSQw4H2HgO9w4GL2C5h269eV7NfbTtyy1W1vSX9MmXq2vuxfb1LPkVFKYqzevHqkMt0ser6HjblsRdB+Bf0Of//Xx4MRwAiAyZsfvlN5AGOyXGsLL8Hp02VtP6YqVT1pqtyszR+f+ZKbePnCATT1Lfm2L+oTbxBqqCJAdQp2wCVehEEDyStkwwLHByIz975u6f9Gtlttd/+NL9MTAU9LmYqBo+c/hmz2PH69z+04LJwDiI2XDAMcGzTv8m/V6+Nk52q9xi6V2jw0ar1tOOFctWTk29Sx1bWlskcdrRFUEyA7RhBP2EADJK2XDABXQQuvZtF13LL1LJ1S/Z6ndqgOP0BXjr9Q33XrZ1LP0MGfpx/rLm9WaMa4s4iJAdog0nERzVBFA4qVsGHi5clOiu+A4WZ5W/e61v+nCd5+31G5rlwKdP3mm1u53qE09S0+19c26dOEaXTiiVA+9Xm25CJAdIqlQGGxvjpkrkQE4Q0qGgRuXfKSKys2J7oYzGIZ+8uFy3f7yPZabTh/3Gz1fPtKGTqU+M7dd+jfkLf5Pje47+yjNWfpxTEoOR8NqhUI2HAKpIeXCwJwXq/TIdzXg09mgjR9r/jOz1b15h6V284ZM1p3H/1ytGVk29czZpp14sIaUFuniBf9W0+7I767+63lD9G51ne5+ZV3I1/k35PXIz9GbV490xJr7mPKSgKWTA4UTNhwCqSGlwsCcFz/SI29uSHQ3Eqa4YYvuWXK7hnz9kaV2r/Q7Rr8be4Xq8gpt6lny6P+DrsrMcEcVBCRpc2OL/v7vryy8vjnoTYOJEOoa5LacdhoCQGRSJgzctDQ9g0BOa4tmrPiLpqx92VK7rwp/oAtPv14f73uQTT1LTlsaW/TeF9HfV1G3o8XSaRYn3vJnJpw47TQEgMikRBh46YMaPfzGhkR3I34MQ+e8v1Rzlj1guenFp16rikOH29Cp5Od2+Xb5R8O/wa4oP9t0m2S+5S+aK5EBOEfShwGP19D1L3yY6G7ExbAv/qPHnp6pHM8eS+3uOO7num/YT+V1Z9jUs9TgjbJCVdsNdoVdzIeBeJ4WiLVorkQG4BxJHwZWV9e1uwI21Ry4rUb3L7pV5d+st9TuxR8er+tOvlQNuV1t6lnqcLmk8Dd0hNc9L0u3nH6ExpSXyOM1Qj4xS76ZiLlnJf+xOysbDgE4U1KGgbaVzioqaxLdnZjLb9mpPyybpzM+etVSu096HahLTr1W63v2tqlnqSkWQUBSu1Aa6onZb+5ZR2nsgNQYKM1uOATgTEl3a2GgSmepwGV4ddE7z+ma1x6z3PaXk2fpX/1+FPtOwbKSDjcMUpkPQCKl5K2FqXgL4cjPVmv+s3+w3G7Oiedr/jGTZLjcNvQKkep4pp4nZgDJIGnCgMdr6JrnPkyJIHDwli/14PM3qV/dRkvt/n7EaM0edaF2ZnexqWfoyL8bfsa4Ml33/Ifaviv8/pRlVbXtjuQ5qX4AAASSNGFg7op12p7EGwULdzXq1oo/65RP37LUbm3JIZo26Wp9XfgDm3qGYNruhh9TXqLCLln6+SPvhG23aO0mXT+OHfQAkkdShAGP19CDr3+e6G5YluH16Ncrn9Tlbz1lqV1TVq7OmzxT7xx4hE09QyD5ORlqavHs/XPH3fBD+/VUUX6W6ppCh9KtTbspvwsgqTg6DPhPDbyxbrN27vaEb+AQ4z5+Q/ctvs1yu2tPnqYnjzzZd9YNcecPAt27ZOnc4aWaNvLgdk/3GW6XThu4v6m7Lyi/CyCZODYMVFTWaNbiKtU2JMcv1cO/Wa+/PPMHlezYaqndo0dP0K0nnKuWTPNFamCv+l2tunv5pzq0uGunHf+jyopNhQHK7wJIJo4MAxWVNbp4wZpEdyOsXk3bdOeLd2rEhvcttVvZZ4Cmj7tSm7sxjexEoa7epfwugFTkuDDgPzXgVNl7WvW71/+qC959wVK7b/O661dnzNB/9jvUno4hpoJdvUv5XQCpyHFh4O3Ptzrv1IBh6CcfLtPtL99ruekV46/UC4efaEOn0FFupkvNe2J7+DTQ2j/ldwGkGkeFAY/X0DMW7oC326CvP9Zjz8xSQUuTpXb3D52sO4+boj0ZjvrrTXl3/vQozVlaFfI+gKL8LM0Yf7jqdrSYuqEw2No/xYQApBLHjFZOKTNc0vCt7l18u47ZWGWp3bKDB+vqUy5XXV6hTT1DMPnZGbrjp0dqTHmJ3G6FnMK/+bTvLxL6y5vVUa39U0wIQKpwRBh46YNNumShtU14sZTb2qwZK/6in6+tsNTuy8If6MLTb9B/9y21qWcIJT/LrV+NOEi/PumQvU/kZqfwWfsHgO8l/KKilz6o0aULE3DfgGHoF2te1OzlD1puetFp1+kfhxxrQ6dg1vRRh3SqA9BW25stQ03hc5EQgFRmdvxOaBhIxBHCYzes1WNPz1K2d4+ldrcff47mDZ0srzvDpp7Bb0T/XpKkdzfUaVert93nuudl6dbTj4jpQG02OABAsnH8rYUer6FZiz+Ky9fqs22Tfvv64ypu3GppL8Diw0bo+pMvVWNOvo29Sz8/6JatU4/aX4v/UxPyidzjNfT251u1av1WSYaGHdRLQ/v1jPlAzdo/gHSXsDCwurpOtQ0ttr1/15admvPP+3Va1b8stft4n7665NRrVV20vz0dS1O5mW6dPeRAjS4r3vvkfdWYw0I+kWe4XRp+cC8NP7hXAnsOAKkvYWHAjtrtLsOrqW8/o6te/5uldl659MufzNLrBx0d8z6lO5ekcQNKdM+ZR3V6oueJHACcIWFhIJa120/67B098uwcy+3mjPyV5v9oogyXO2Z9SWclhbm6/pTD9E1js76o26k+RXk6Z1hfZWfy9wsATpawMLCtaXdU7ft/+4Ueev5GlW6rsdTuwtOu16v9fqTWjKyovj6kLllunVJeouMP2UfFBWy8A4BklZAwUFHpO05oVeGuRv3x5Xt08rq3LbVbs9+humzi1dpYuK/lr4n2euZna9LA/dqt/QMAklvcw4DHa2j2kirTdQUyvB5d/uZC/XrV3y19ncbsLjpv8ky927vceifTSIZLGnXYvvrfY0t1TN8ivffFNm1ubFavrjmS4dvbUde0W0Vdc3j6B4AUFfcwsLq6zlTJ4QlVr+nPS263/P5Xj7lMfx/wY8nFgBVIj7wsDT2op/rtkx/wqB4b+gAg/cQ9DIQ6RXB47Wd65Nk/qHhHnaX3nH/0RN12wi/VkpkdbfeSSnaGS8f1K9JdZx6tqpqGvU/0Xq+hd6rrJBka0ren3BkubdnRQkEdAEBAcQ8DHU8R9GrapjtfvFMjNli7m+DNPkdq+vgr9W3X4BfJpIKcTJcu/p9+uvh/DtbCd74Iuku/4xP98YfsE++uAgCSVNzDwODSIpUU5uqb7Tt11Wt/1cXvPGu67eb8Hjr/jN/rw5L+NvYwPrpkSL265Sg/J0vbd7WqW06mBvYpVM+8XH3T0Kz9unfRsf3aT+Off/xBCe41ACAVxT0M+G+Le/73c00HgcvHX6lFh59oc89ir0uWS4Yheb2GigtzdXz/fXTD+MPVJZv7DQAAzpGQo4VjykvU59DQRYfmDvup7h5+tvZkOOKW5Xb2ycvQvoV56tUtW+u+2aHmVo+6dcnSkNIiZWdmqLRnPsV2AABJI2Ej7WGXXyDjuQVyffrp3o8tO3iIrjrl19qWV5iobik7w6VTyn+gbxt36/MtTeqWk6mj+/bQMX2KtF+PPDbgAQBSTuIeu4uL5Xr/fWndOqlHD1U0ZGn2kiptM3HsMBYKcjP0oz49dPaQvlr71XbZeSseAABO5jIMI2z9H7P3IUer7b3yvfJzJJe0adtO/bPqG329bZc2bG1qd799XpZUlJ+rhuZWte7xqmtOhvKyM9UjP1utHq9yMjN0YM98nXbU/sp0u7WlieN1AID0YXb8dtSCfLBb7H5yzIGS2ocFBnUAAGLDUWEgHK68BQAg9tjuDgBAmiMMAACQ5ggDAACkOcIAAABpjjAAAECaIwwAAJDmCAMAAKQ5wgAAAGmOMAAAQJojDAAAkOZMlSP232XU0NBga2cAAEDs+MftcHcSmgoDjY2NkqTevXtH2S0AABBvjY2NKiwsDPp5U1cYe71ebdq0Sd26dZPLxS2BAAAkA8Mw1NjYqP32209ud/CdAabCAAAASF1sIAQAIM0RBgAASHOEAQAA0hxhAACANEcYAAAgzREGAABIc4QBAADS3P8D0Ik/ejF4kB0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# scatter data\n",
    "plt.plot(X,Y,'o')\n",
    "# best fit line\n",
    "y_pred = model.intercept_ + model.coef_ * X\n",
    "plt.plot(X, y_pred, color='red', linewidth=3)\n",
    "\n",
    "plt.xticks(())\n",
    "plt.yticks(())\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c639cbb5-6bdd-4032-812a-5bf8e0174833",
   "metadata": {},
   "source": [
    "We see that there is a strong correlation between square footage and house price.  The coefficient of determination measures the strength of the correlation and varies between 0 (no correlation) and 1.0 (perfectly correlated)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e2d9af2-3882-4c1d-9f58-038567ffd5ac",
   "metadata": {},
   "source": [
    "### Multiple Linear Regression\n",
    "\n",
    "Adding more variables often improves the score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fccb721f-db9f-4c58-ad3b-d7ed629791a2",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "coefficient of determination: 0.5249639858277426\n"
     ]
    }
   ],
   "source": [
    "X = np.array(dataset[['sqft_living','condition','yr_built']])\n",
    "Y = np.array(dataset.price)\n",
    "\n",
    "# shape X into matrix of a single column\n",
    "X = X.reshape((X.shape[0],3))\n",
    "\n",
    "model = LinearRegression()\n",
    "model.fit(X,Y)\n",
    "r_sq = model.score(X,Y)\n",
    "print('coefficient of determination:', r_sq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e900ab88-873c-4290-9d99-723cef811c6d",
   "metadata": {},
   "source": [
    "but some variables are confounding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ff60bd53-6ebc-4f13-8c78-5e81992aebd4",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "coefficient of determination: 0.1358874997770878\n"
     ]
    }
   ],
   "source": [
    "X = np.array(dataset[['zipcode','floors','waterfront']])\n",
    "Y = np.array(dataset.price)\n",
    "\n",
    "# shape X into matrix of a single column\n",
    "X = X.reshape((X.shape[0],3))\n",
    "\n",
    "model = LinearRegression()\n",
    "model.fit(X,Y)\n",
    "r_sq = model.score(X,Y)\n",
    "print('coefficient of determination:', r_sq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d948a45-a861-423a-b2ce-8fcf1bf976c3",
   "metadata": {},
   "source": [
    "This leads to a first parallel program.  What are the right set of variables? A brute force approach called _All Possible Regressions_ examines all combinations. So, let's build a big matrix that and we will regress on subsets.  We will look at all combinations of 1, 2, or 3 variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "615005bd-f622-43a0-a844-c32f11e2676d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from itertools import combinations, chain\n",
    "\n",
    "Y = np.array(dataset.price)\n",
    "X = np.array(dataset.drop(['price'], axis=1))\n",
    "\n",
    "## Let's choose all combinations of 1, 2, and 3 columns.\n",
    "col_idxs = np.array(range(X.shape[1]))\n",
    "combos = list(chain(combinations(col_idxs, 1), combinations(col_idxs, 2), combinations(col_idxs, 3)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bc919ff9-24aa-4a9b-a4ee-efc84f8b51d9",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6095149101819037 (2, 6, 14)\n"
     ]
    }
   ],
   "source": [
    "model = LinearRegression()\n",
    "\n",
    "# do in a for loop (dumbest way)\n",
    "r_sq_best = 0.0\n",
    "for combo in combos:\n",
    "    Xp = X[:,combo]\n",
    "    model = model.fit(Xp,Y)\n",
    "    r_sq = model.score(Xp,Y)\n",
    "    if r_sq > r_sq_best:\n",
    "        r_sq_best = r_sq\n",
    "        combo_best = combo\n",
    "    \n",
    "print(r_sq_best, combo_best)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "603e45b6-80e4-4c6c-9a2f-1e8a728b1e5b",
   "metadata": {},
   "source": [
    "The outcome is kind of crazy. The fields are `sqft-living`, `view`, and `latitude`. Latitude is probably a somewhat accurate proxy for wealth in this area, e.g. N of town richer than south of town. But, this is the kind of outcome that would not translate to other regions, i.e. is likely specific to this data.  _Neat_.\n",
    "\n",
    "Back to performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0e8377d0-63a9-439c-be04-92d19f25e71c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.02 s ± 30 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "model = LinearRegression()\n",
    "\n",
    "r_sq_best = 0.0\n",
    "for combo in combos:\n",
    "    Xp = X[:,combo]\n",
    "    model = model.fit(Xp,Y)\n",
    "    r_sq = model.score(Xp,Y)\n",
    "    if r_sq > r_sq_best:\n",
    "        r_sq_best = r_sq\n",
    "        combo_best = combo"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f8f2fb1-6df7-4902-9e5a-ed4ccfc231a1",
   "metadata": {},
   "source": [
    "Let's see if flattening the loop matters. Replace for loop with a list comprehension of all combination of variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7993f4bd-6331-42f5-b909-31c6f1b659cf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def r_sq_regression (combo):\n",
    "    Xp = X[:,combo]\n",
    "    model = LinearRegression()\n",
    "    model = model.fit(Xp,Y)\n",
    "    return model.score(Xp,Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "96184a66-012a-4068-a462-677fcafb2ba7",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(2, 6, 14) 0.6095149101819037]\n"
     ]
    }
   ],
   "source": [
    "r_sq_list = [ (combo, r_sq_regression(combo)) for combo in combos ]\n",
    "r_sq_arr = np.array(r_sq_list, dtype=object)\n",
    "r_sq_idx = np.argmax(r_sq_arr[:,1])\n",
    "print(r_sq_arr[r_sq_idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "efbf14c0-c47c-4b51-aed3-cde5f1b9af7e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "992 ms ± 17.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "r_sq_list = [ (combo, r_sq_regression(combo)) for combo in combos ]\n",
    "r_sq_arr = np.array(r_sq_list, dtype=object)\n",
    "r_sq_idx = np.argmax(r_sq_arr[:,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3298cb7-ab53-4a34-8442-883faa04a8ab",
   "metadata": {},
   "source": [
    "It didn't seem to help. But, this is a step toward parallelization.\n",
    "\n",
    "### Exercise (This is the assignment)\n",
    "\n",
    "1. Use `joblib.Parallel` and `joblib.delayed` to parallelize the computation of the calls to `r_sq_regression`. \n",
    "    1. In one cell, print the answer to verify that your parallel program is correct.\n",
    "    2. In another cell time the computation. What is the speedup at `n_jobs=4`?\n",
    "    3. Based on this single number, estimate the Amdahl number for this computation. Show your work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d1194559-1102-4978-bc80-a4911dddf37f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6095149101819037 (2, 6, 14)\n"
     ]
    }
   ],
   "source": [
    "# TODO code for 1A\n",
    "from joblib import Parallel, delayed\n",
    "r_sq_arr = Parallel(n_jobs=4)(delayed(r_sq_regression)(combo) for combo in combos)\n",
    "r_sq_idx = np.argmax(r_sq_arr)\n",
    "print(r_sq_arr[r_sq_idx], combos[r_sq_idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "3643a3bb-383c-4da9-aa94-e2f35e9b67ad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "586 ms ± 4.66 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit \n",
    "# code for 1B \n",
    "r_sq_arr = Parallel(n_jobs=4)(delayed(r_sq_regression)(combo) for combo in combos)\n",
    "r_sq_idx = np.argmax(r_sq_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8ab0daba-7e41-4d74-8bb1-df5cc13181d8",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7889546351084813"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# estimate for 1C\n",
    "#. we have a 586 ms with 4 cores.  so speedup of 1.69\n",
    "#1.69 = 1 / ((1-p)+p/4)\n",
    "1/1.69*4/3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec4d94a2-fb3d-4d04-a315-6c913e677fad",
   "metadata": {},
   "source": [
    "2. Use the batch size parameter to vary the number of tasks in each batch from 1,2,...128 @ n_jobs=4. You will need to look at the joblib documentation to read up about batch size.\n",
    "   1. Plot your results\n",
    "   2. Model the problem as having two performance components: a fixed startup cost per batch ($C_B$) and perfect parallelism.  Estimate the batch startup cost.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "0ebbd6e8-16ae-4a89-bed8-05f10bda3f17",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.07 s ± 23.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "747 ms ± 3.53 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "497 ms ± 4.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "454 ms ± 4.89 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "434 ms ± 5.71 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "440 ms ± 22.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "421 ms ± 2.28 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "424 ms ± 12.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "[(1, 1.0721810302929953), (2, 0.7469500536930614), (4, 0.49679515513292116), (8, 0.45430990514744607), (16, 0.43380886944110636), (32, 0.4400174761373949), (64, 0.4206638155737892), (128, 0.42407120855724706)]\n"
     ]
    }
   ],
   "source": [
    "# code for 2A\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "timings = list()\n",
    "for i in [1,2,4,8,16,32,64,128]: \n",
    "  result = %timeit -o Parallel(n_jobs=4, batch_size=i)(delayed(r_sq_regression)(combo) for combo in combos)\n",
    "  timings.append((i,result.average))\n",
    "    \n",
    "print(timings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "f3b5cf1b-4914-4294-8077-844c16e8aabc",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAIPCAYAAABg71YVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABDKUlEQVR4nO3de3hU1b3/8c/ck5ALQiDcQkCLiqKooVJAVKon/uKt9lhLtYpa6CmnKmKqRYpWpdp4qYjVggVF6tFWngpY26JtahEvtCIhqQioqGACJMYg5sIlk8zs3x/AyMzkMpPZM3smvF/PMw/OytqZlYWPfviy9nfbDMMwBAAAAKQgu9ULAAAAALqLMAsAAICURZgFAABAyiLMAgAAIGURZgEAAJCyCLMAAABIWYRZAAAApCzCLAAAAFIWYRYAAAApizALAACAlOW08sNff/11PfTQQyovL1dNTY1Wrlypyy67rMP5NTU1+slPfqLy8nJt3bpVM2bM0Pz586P6TL/fr127dikrK0s2my22HwAAAACmMwxDTU1NGjRokOz2zmuvlobZvXv3avTo0br++ut1+eWXdzm/paVF/fr105w5c/TII4906zN37dql/Pz8bl0LAACAxKmurtaQIUM6nWNpmC0uLlZxcXHE84cNG6ZHH31UkrRkyZJufWZWVpakg5uTnZ3dre8BAACA+GlsbFR+fn4gt3XG0jCbCC0tLWppaQm8b2pqkiRlZ2cTZgEAAJJYJEdCe/wNYKWlpcrJyQm8OGIAAADQc/T4MDt79mw1NDQEXtXV1VYvCQAAACbp8ccMPB6PPB6P1csAAABAHPT4MAsAAIDEMgxDbW1t8vl8Hc5xuVxyOBwxf5alYba5uVkfffRR4P22bdtUWVmpPn36aOjQoZo9e7Z27typZ555JjCnsrIycO3nn3+uyspKud1unXTSSYlePgAAAEJ4vV7V1NRo3759nc6z2WwaMmSIMjMzY/o8m2EYRkzfIQavvfaaJk2aFDZ+7bXXaunSpbruuuu0fft2vfbaa4GvtXdXW0FBgbZv3x7RZzY2NionJ0cNDQ10MwAAADCR3+/X1q1b5XA41K9fP7nd7nazm2EY+vzzz7Vv3z6NGDEirEIbTV6ztDJ77rnnqrMsvXTp0rAxC7M3AAAAOuH1euX3+5Wfn6+MjIxO5/br10/bt29Xa2trTMcNenw3AwAAACRWV4+glSLrIRvRZ5nyXQAAAAALEGYBAACQsgizAAAASFmEWQAAAKQswiwAAABMFUn3KbM6VBFmAQAAYAqXyyVJXT4wQTrYxktSzE8B43G2AAAAMIXD4VDv3r1VV1cnScrIyGi3BZff79fnn3+ujIwMOZ2xxVHCbBy9X9uorZ81q9XnV6vPr/w+GRp/XK7VywIAAIibAQMGSFIg0HbEbrdr6NChMfebJczG0Z8qd2nhax8H3n/rtEGEWQAA0KPZbDYNHDhQ/fv3V2tra4fz3G53RA9X6AphNo5c9uA/abT5eBQvAAA4OjgcjpjPw0aCG8DiyOkI3l6vz2/RSgAAAHomwmwcuULCbBthFgAAwFSE2ThyOUKOGfg5ZgAAAGAmwmwchVZmvW1UZgEAAMxEmI0jJ5VZAACAuCLMxlFoZbaVM7MAAACmIszGUeiZ2VZacwEAAJiKMBtHTjvdDAAAAOKJMBtHHDMAAACIL8JsHHHMAAAAIL4Is3FEZRYAACC+CLNxRGsuAACA+CLMxpE7tDLLQxMAAABMRZiNI2domPUTZgEAAMxEmI0jpz3kmAE3gAEAAJiKMBtHbmdIn1m/IcMg0AIAAJiFMBtHoZVZifZcAAAAZiLMxlFoay6J9lwAAABmIszGUXthlnOzAAAA5iHMxlFon1mJjgYAAABmIszGEccMAAAA4oswG0eudiqzHDMAAAAwD2E2jtqrzHqpzAIAAJiGMBtH7bXmojILAABgHsJsHNlstrCjBpyZBQAAMA9hNs6c9uAtJswCAACYhzAbZ6Htudr8HDMAAAAwC2E2ztwhN4G1tlGZBQAAMAthNs5CK7OtVGYBAABMQ5iNs9D2XFRmAQAAzEOYjbPQMNvG42wBAABMQ5iNs/DWXBwzAAAAMAthNs5ozQUAABA/hNk4C63M8gQwAAAA8xBm4yz0zKyXyiwAAIBpCLNxFvbQBMIsAACAaQizcRbWmotjBgAAAKYhzMZZWJilNRcAAIBpCLNx5rRzAxgAAEC8EGbjzOWkNRcAAEC8EGbjzGXnoQkAAADxQpiNs/AbwKjMAgAAmMXSMPv666/rkksu0aBBg2Sz2fTiiy92ec2aNWtUWFiotLQ0HXvssXriiSfiv9AYOEPCLK25AAAAzGNpmN27d69Gjx6txx9/PKL527Zt04UXXqiJEyeqoqJCP/vZzzRjxgwtX748zivtPndIn1kvxwwAAABM47Tyw4uLi1VcXBzx/CeeeEJDhw7V/PnzJUkjR47U+vXr9atf/UqXX355nFYZGyqzAAAA8ZNSZ2b/9a9/qaioKGjsggsu0Pr169Xa2truNS0tLWpsbAx6JVLYE8D8VGYBAADMklJhtra2Vnl5eUFjeXl5amtrU319fbvXlJaWKicnJ/DKz89PxFID3CGVWS+VWQAAANOkVJiVJJstuNJpGEa744fNnj1bDQ0NgVd1dXXc13gkp51jBgAAAPFi6ZnZaA0YMEC1tbVBY3V1dXI6nerbt2+713g8Hnk8nkQsr10uJ31mAQAA4iWlKrPjxo1TWVlZ0Njf//53jRkzRi6Xy6JVdc5lp88sAABAvFgaZpubm1VZWanKykpJB1tvVVZWqqqqStLBIwJTpkwJzJ8+fbo+/fRTlZSUaMuWLVqyZImeeuop3XrrrVYsPyKu0BvAqMwCAACYxtJjBuvXr9ekSZMC70tKSiRJ1157rZYuXaqamppAsJWk4cOHa9WqVbrlllv0m9/8RoMGDdKvf/3rpG3LJYW35qIyCwAAYB5Lw+y5554buIGrPUuXLg0bO+ecc7Rhw4Y4rspcoZXZVlpzAQAAmCalzsymIldoZbaNyiwAAIBZCLNxFvYEMD9hFgAAwCyE2Thzhx4z4AYwAAAA0xBm4yz0oQncAAYAAGAewmycOWnNBQAAEDeE2Thz05oLAAAgbgizcUafWQAAgPghzMZZWJ9ZjhkAAACYhjAbZ6F9ZmnNBQAAYB7CbJyFPTTBZ3T61DMAAABEjjAbZ6HdDCSpjUfaAgAAmIIwG2cue/gW054LAADAHITZOHM5wyuzXjoaAAAAmIIwG2ehTwCTpDbCLAAAgCkIs3EW+tAEifZcAAAAZiHMxll7N4Dx4AQAAABzEGbjLLQ1l0Q3AwAAALMQZuMs9AlgEpVZAAAAsxBm48xms8lhD32kLWEWAADADITZBAitznIDGAAAgDkIswkQ+uAEWnMBAACYgzCbAC5n8DZTmQUAADAHYTYBnJyZBQAAiAvCbAKEtudq8xNmAQAAzECYTYDQG8C8bRwzAAAAMANhNgGcVGYBAADigjCbAKHHDDgzCwAAYA7CbALQZxYAACA+CLMJEHYDGGEWAADAFITZBKA1FwAAQHwQZhOAM7MAAADxQZhNAM7MAgAAxAdhNgHCWnNRmQUAADAFYTYB3KHHDPxUZgEAAMxAmE0AZ9gxAyqzAAAAZiDMJkB4ay7CLAAAgBkIswnADWAAAADxQZhNAKed1lwAAADxQJhNAPrMAgAAxAdhNgFCjxnwOFsAAABzEGYTIKwyS2suAAAAUxBmEyCsNVcbxwwAAADMQJhNgLDWXH7CLAAAgBkIswkQembWy5lZAAAAUxBmEyC0NRcPTQAAADAHYTYBXM7QMEtlFgAAwAyE2QRw2UOPGVCZBQAAMANhNgG4AQwAACA+CLMJEN6ai2MGAAAAZiDMJkD4QxOozAIAAJiBMJsAYWGWM7MAAACmIMwmQOgxA7oZAAAAmIMwmwDusMosYRYAAMAMhNkEcIa05uKYAQAAgDksD7MLFizQ8OHDlZaWpsLCQr3xxhudzv/Nb36jkSNHKj09XSeccIKeeeaZBK20+8IfmkCYBQAAMIPTyg9ftmyZZs6cqQULFmjChAn67W9/q+LiYm3evFlDhw4Nm79w4ULNnj1bixcv1te//nWtW7dOP/zhD3XMMcfokksuseAniIzLzjEDAACAeLAZhmFZsho7dqzOOOMMLVy4MDA2cuRIXXbZZSotLQ2bP378eE2YMEEPPfRQYGzmzJlav3693nzzzYg+s7GxUTk5OWpoaFB2dnbsP0QEttQ0qvjRryrONpu0rfSihHw2AABAqokmr1l2zMDr9aq8vFxFRUVB40VFRVq7dm2717S0tCgtLS1oLD09XevWrVNra2uH1zQ2Nga9Ei20NZdhSD4/1VkAAIBYWRZm6+vr5fP5lJeXFzSel5en2tradq+54IIL9OSTT6q8vFyGYWj9+vVasmSJWltbVV9f3+41paWlysnJCbzy8/NN/1m64gppzSVxExgAAIAZLL8BzGYLDnqGYYSNHXbnnXequLhY3/jGN+RyufStb31L1113nSTJ4XC0e83s2bPV0NAQeFVXV5u6/kiEVmYlwiwAAIAZLAuzubm5cjgcYVXYurq6sGrtYenp6VqyZIn27dun7du3q6qqSsOGDVNWVpZyc3Pbvcbj8Sg7OzvolWihD02QuAkMAADADJaFWbfbrcLCQpWVlQWNl5WVafz48Z1e63K5NGTIEDkcDj3//PO6+OKLZbdbXmTuUGg3A4n2XAAAAGawtDVXSUmJrrnmGo0ZM0bjxo3TokWLVFVVpenTp0s6eERg586dgV6yH374odatW6exY8dqz549mjdvnt577z397ne/s/LH6FJon1lJ8hJmAQAAYmZpmJ08ebJ2796tuXPnqqamRqNGjdKqVatUUFAgSaqpqVFVVVVgvs/n08MPP6wPPvhALpdLkyZN0tq1azVs2DCLfoLIhD4BTJLaOGYAAAAQM0v7zFrBij6zPr+h4362KmjsHyVn62v9sxLy+QAAAKkkJfrMHk0cdptCi7PetqPqzxAAAABxQZhNkND2XG1+zswCAADEijCbIKFhlj6zAAAAsSPMJkhor1n6zAIAAMSOMJsgVGYBAADMR5hNEFfIHWC05gIAAIgdYTZBQh+cQGUWAAAgdoTZBAl9cAJnZgEAAGJHmE0QWnMBAACYjzCbIKFh1ttGmAUAAIgVYTZBQltztfk5ZgAAABArwmyChB0z4AYwAACAmBFmE8QVUpn1cgMYAABAzAizCUJlFgAAwHyE2QRx2ukzCwAAYDbCbIKEHjOgzywAAEDsCLMJEnrMgMosAABA7AizCUJrLgAAAPMRZhPETWUWAADAdITZBAmtzBJmAQAAYkeYTZDw1lwcMwAAAIgVYTZBQsOsl8osAABAzAizCeK0h9wARmUWAAAgZoTZBAk7ZuCnMgsAABArwmyChD40wdtGZRYAACBWhNkEoTILAABgPsJsgjjpMwsAAGA6wmyChB4zaOUGMAAAgJgRZhMk9JgBlVkAAIDYEWYThNZcAAAA5iPMJojbSWUWAADAbITZBHHaCbMAAABmI8wmSOgNYG1+jhkAAADEijCbIGE3gLVRmQUAAIgVYTZBnKGtuajMAgAAxIwwmyBhTwDjzCwAAEDMCLMJwkMTAAAAzEeYTRAemgAAAGA+wmyC0JoLAADAfITZBHE7g48Z+A3Jx01gAAAAMSHMJkhoZVaiOgsAABArwmyChLbmknhwAgAAQKwIswnidoRvNe25AAAAYkOYTRBnO2HWS5gFAACICWE2QUL7zEpSG71mAQAAYkKYTZDQPrMSN4ABAADEijCbIE57eGWWp4ABAADEhjCbIA67TbaQPNvmpzILAAAQC8JsgthsNrlCnwLWRmUWAAAgFoTZBAq9CayVyiwAAEBMCLMJFNqeq7WNMAsAABALwmwChXY04AlgAAAAsSHMJlDoMQMemgAAABAby8PsggULNHz4cKWlpamwsFBvvPFGp/Ofe+45jR49WhkZGRo4cKCuv/567d69O0GrjY0zJMzy0AQAAIDYWBpmly1bppkzZ2rOnDmqqKjQxIkTVVxcrKqqqnbnv/nmm5oyZYqmTp2qTZs26Y9//KPeeecdTZs2LcEr756wYwZUZgEAAGJiaZidN2+epk6dqmnTpmnkyJGaP3++8vPztXDhwnbn//vf/9awYcM0Y8YMDR8+XGeddZZ+9KMfaf369QleefeEtubimAEAAEBsLAuzXq9X5eXlKioqChovKirS2rVr271m/Pjx2rFjh1atWiXDMPTZZ5/phRde0EUXXdTh57S0tKixsTHoZRWXk2MGAAAAZrIszNbX18vn8ykvLy9oPC8vT7W1te1eM378eD333HOaPHmy3G63BgwYoN69e+uxxx7r8HNKS0uVk5MTeOXn55v6c0TDGfrQBCqzAAAAMbH8BjBbyDNeDcMIGzts8+bNmjFjhn7+85+rvLxcr7zyirZt26bp06d3+P1nz56thoaGwKu6utrU9UfDHdpnltZcAAAAMXFa9cG5ublyOBxhVdi6urqwau1hpaWlmjBhgm677TZJ0qmnnqpevXpp4sSJuvfeezVw4MCwazwejzwej/k/QDeEdzOgMgsAABALyyqzbrdbhYWFKisrCxovKyvT+PHj271m3759sof8Vb3D4ZB0sKKb7MKeAEaYBQAAiImlxwxKSkr05JNPasmSJdqyZYtuueUWVVVVBY4NzJ49W1OmTAnMv+SSS7RixQotXLhQn3zyid566y3NmDFDZ555pgYNGmTVjxExd0hltpUbwAAAAGJi2TEDSZo8ebJ2796tuXPnqqamRqNGjdKqVatUUFAgSaqpqQnqOXvdddepqalJjz/+uH7yk5+od+/e+uY3v6kHHnjAqh8hKtwABgAAYC6bkQp/P2+ixsZG5eTkqKGhQdnZ2Qn97Jv+UKE//2dX4P2Nk76mWy84IaFrAAAASHbR5DXLuxkcTVz20GMGVGYBAABiQZhNoNBuBpyZBQAAiA1hNoFcId0M2vxUZgEAAGJBmE2g0DDLMQMAAIDYEGYTyMUxAwAAAFMRZhOIhyYAAACYizCbQGFnZqnMAgAAxIQwm0C05gIAADAXYTaBOGYAAABgLsJsAoXeANbm55gBAABALAizCRR6ZtbbRmUWAAAgFoTZBAp/aAKVWQAAgFgQZhMo/HG2VGYBAABiQZhNIHfYDWBUZgEAAGJBmE2g0MpsG5VZAACAmHQrzH788ce64447dOWVV6qurk6S9Morr2jTpk2mLq6ncdppzQUAAGCmqMPsmjVrdMopp+jtt9/WihUr1NzcLEl69913ddddd5m+wJ7E7Qw9M8sxAwAAgFhEHWZvv/123XvvvSorK5Pb7Q6MT5o0Sf/6179MXVxPQ2UWAADAXFGH2Y0bN+rb3/522Hi/fv20e/duUxbVU9GaCwAAwFxRh9nevXurpqYmbLyiokKDBw82ZVE9VegTwFp5aAIAAEBMog6zV111lWbNmqXa2lrZbDb5/X699dZbuvXWWzVlypR4rLHHcIa25vITZgEAAGIRdZi97777NHToUA0ePFjNzc066aSTdPbZZ2v8+PG644474rHGHiO0MtvGDWAAAAAxcUZ7gcvl0nPPPae5c+eqoqJCfr9fp59+ukaMGBGP9fUo7Z2ZNQxDNputgysAAADQmajD7GHHHXecjjvuODPX0uOFhlnpYHuu0JZdAAAAiEzUYdYwDL3wwgtavXq16urq5A8597lixQrTFtfTOO3hobXV55fbyYPYAAAAuiPqMHvzzTdr0aJFmjRpkvLy8vgr8ii0F1o5NwsAANB9UYfZZ599VitWrNCFF14Yj/X0aO1WZuloAAAA0G1R//12Tk6Ojj322HispcdrrzLbQq9ZAACAbos6zN5999265557tH///nisp0fLcIcXwve2tFmwEgAAgJ4h6mMGV1xxhf7whz+of//+GjZsmFwuV9DXN2zYYNriehqH3aZebof2en2BsaYDhFkAAIDuijrMXnfddSovL9fVV1/NDWDdkJnmDAmzrRauBgAAILVFHWb/+te/6m9/+5vOOuuseKynx8v0OPWZWgLvmzlmAAAA0G1Rn5nNz89XdnZ2PNZyVMhMCz6W0cwxAwAAgG6LOsw+/PDD+ulPf6rt27fHYTk9X3ZacDGcyiwAAED3RX3M4Oqrr9a+fft03HHHKSMjI+wGsC+++MK0xfVEmZ7gLecGMAAAgO6LOszOnz8/Dss4eoSGWSqzAAAA3Rd1mL322mvjsY6jRmboMQMqswAAAN0WUZhtbGwM3PTV2NjY6VxuDutcFpVZAAAA00QUZo855hjV1NSof//+6t27d7u9ZQ3DkM1mk8/na+c74LDQymwTYRYAAKDbIgqz//znP9WnTx9J0tNPP638/Hw5HI6gOX6/X1VVVeavsIfJ9IS25uKhCQAAAN0VUZg955xzAv/8gx/8IFClPdLu3bt1/vnnc6a2C2FnZqnMAgAAdFvUfWYPHycI1dzcrLS0NFMW1ZOFnpmlNRcAAED3RdzNoKSkRJJks9l05513KiMjI/A1n8+nt99+W6eddprpC+xp6GYAAABgnojDbEVFhaSDldmNGzfK7XYHvuZ2uzV69Gjdeuut5q+whwnrM+ttk99vyG4Pr3YDAACgcxGH2dWrV0uSrr/+ej366KO04OqmrJDKrGFI+1p9YSEXAAAAXYs6QT399NPxWMdRIyukm4F08KgBYRYAACB6Ud8Ahtj08jjCxppbaM8FAADQHYTZBHM67Ep3BQdaOhoAAAB0D2HWAvSaBQAAMAdh1gKhvWZpzwUAANA9hFkLhFZmm6jMAgAAdAth1gJhvWapzAIAAHQLYdYCoWGWG8AAAAC6hzBrgfAbwGjNBQAA0B2Wh9kFCxZo+PDhSktLU2Fhod54440O51533XWy2Wxhr5NPPjmBK45d2A1gnJkFAADoFkvD7LJlyzRz5kzNmTNHFRUVmjhxooqLi1VVVdXu/EcffVQ1NTWBV3V1tfr06aMrrrgiwSuPTVZa8FPAOGYAAADQPZaG2Xnz5mnq1KmaNm2aRo4cqfnz5ys/P18LFy5sd35OTo4GDBgQeK1fv1579uzR9ddfn+CVx4Y+swAAAOawLMx6vV6Vl5erqKgoaLyoqEhr166N6Hs89dRTOv/881VQUNDhnJaWFjU2Nga9rEY3AwAAAHNYFmbr6+vl8/mUl5cXNJ6Xl6fa2tour6+pqdHLL7+sadOmdTqvtLRUOTk5gVd+fn5M6zZDFpVZAAAAU1h+A5jNZgt6bxhG2Fh7li5dqt69e+uyyy7rdN7s2bPV0NAQeFVXV8eyXFPQmgsAAMAczq6nxEdubq4cDkdYFbauri6sWhvKMAwtWbJE11xzjdxud6dzPR6PPB5PzOs1U9gxAyqzAAAA3WJZZdbtdquwsFBlZWVB42VlZRo/fnyn165Zs0YfffSRpk6dGs8lxk17N4AZhmHRagAAAFKXZZVZSSopKdE111yjMWPGaNy4cVq0aJGqqqo0ffp0SQePCOzcuVPPPPNM0HVPPfWUxo4dq1GjRlmx7JhleYJbc/n8hg60+pXudli0IgAAgNRkaZidPHmydu/erblz56qmpkajRo3SqlWrAt0JampqwnrONjQ0aPny5Xr00UetWLIpQiuzktR0oJUwCwAAECWbcZT9/XZjY6NycnLU0NCg7OxsS9bgbfPr+DteDhp79Sfn6Lh+mZasBwAAIJlEk9cs72ZwNHI77fI4g7eeXrMAAADRI8xahF6zAAAAsSPMWoReswAAALEjzFqkvfZcAAAAiA5h1iJhD0440GrRSgAAAFIXYdYimSG9ZqnMAgAARI8wa5HQG8CaCLMAAABRI8xaJPyYAWEWAAAgWoRZi3ADGAAAQOwIsxahNRcAAEDsCLMWCXtoAmEWAAAgaoRZi3ADGAAAQOwIsxYJb81Fn1kAAIBoEWYtQjcDAACA2BFmLRJ2ZralTYZhWLQaAACA1ESYtUhoZbbVZ6ilzW/RagAAAFITYdYioX1mJXrNAgAARIswa5HQyqzEuVkAAIBoEWYt4nHa5XLYgsaozAIAAESHMGsRm83GU8AAAABiRJi1UOi52aYD9JoFAACIBmHWQuEPTqAyCwAAEA3CrIXa6zULAACAyBFmLZTFmVkAAICYEGYtFHpmlsosAABAdAizFgrtZkCfWQAAgOgQZi1EZRYAACA2hFkLcWYWAAAgNoRZC4UdM2ihzywAAEA0CLMWykyjzywAAEAsCLMW4gYwAACA2BBmLRT60ATOzAIAAESHMGuh0MpsE8cMAAAAokKYtVBoZdbb5ldLm8+i1QAAAKQewqyFQvvMStLeFsIsAABApAizFsryuMLGuAkMAAAgcoRZC6W57HLYbUFjTfSaBQAAiBhh1kI2m432XAAAADEgzFos/ClghFkAAIBIEWYtFtrRgDALAAAQOcKsxcJ6zXLMAAAAIGKEWYuFtueiMgsAABA5wqzFwiuzdDMAAACIFGHWYmFnZjlmAAAAEDHCrMWy0oIfnNDEMQMAAICIEWYtRp9ZAACA7iPMWow+swAAAN1HmLUY3QwAAAC6jzBrsSyOGQAAAHQbYdZioZVZbgADAACIHGHWYtwABgAA0H2EWYuF9pnd3+pTm89v0WoAAABSC2HWYpkeV9jY3hafBSsBAABIPZaH2QULFmj48OFKS0tTYWGh3njjjU7nt7S0aM6cOSooKJDH49Fxxx2nJUuWJGi15gs9MytJjTzSFgAAICLhSSqBli1bppkzZ2rBggWaMGGCfvvb36q4uFibN2/W0KFD273mu9/9rj777DM99dRT+trXvqa6ujq1taXuOdMMl0M2m2QYX43RngsAACAylobZefPmaerUqZo2bZokaf78+frb3/6mhQsXqrS0NGz+K6+8ojVr1uiTTz5Rnz59JEnDhg1L5JJNZ7fblOlxqumIG78IswAAAJGx7JiB1+tVeXm5ioqKgsaLioq0du3adq956aWXNGbMGD344IMaPHiwjj/+eN16663av39/h5/T0tKixsbGoFeyodcsAABA91hWma2vr5fP51NeXl7QeF5enmpra9u95pNPPtGbb76ptLQ0rVy5UvX19frxj3+sL774osNzs6WlpbrnnntMX7+ZMtOcUsNX7+k1CwAAEBnLbwCz2WxB7w3DCBs7zO/3y2az6bnnntOZZ56pCy+8UPPmzdPSpUs7rM7Onj1bDQ0NgVd1dbXpP0Os6DULAADQPZZVZnNzc+VwOMKqsHV1dWHV2sMGDhyowYMHKycnJzA2cuRIGYahHTt2aMSIEWHXeDweeTwecxdvssy04PZczS10MwAAAIiEZZVZt9utwsJClZWVBY2XlZVp/Pjx7V4zYcIE7dq1S83NzYGxDz/8UHa7XUOGDInreuOJM7MAAADdY+kxg5KSEj355JNasmSJtmzZoltuuUVVVVWaPn26pINHBKZMmRKYf9VVV6lv3766/vrrtXnzZr3++uu67bbb9IMf/EDp6elW/RgxCz1mwJlZAACAyFjammvy5MnavXu35s6dq5qaGo0aNUqrVq1SQUGBJKmmpkZVVVWB+ZmZmSorK9NNN92kMWPGqG/fvvrud7+re++916ofwRShD06gMgsAABAZm2Ec2a6/52tsbFROTo4aGhqUnZ1t9XIkSY+UfahHX90aeF88aoAWXl1o4YoAAACsE01es7ybAaSskMpsE5VZAACAiBBmkwBnZgEAALqHMJsEskJbcx2gNRcAAEAkCLNJIOwGMCqzAAAAESHMJgGeAAYAANA9hNkkEHoD2F6vTz7/UdVkAgAAoFsIs0kgtDIrSXu9VGcBAAC6QphNAqFnZiWOGgAAAESCMJsEernbCbPcBAYAANAlwmwScNht6uV2BI3x4AQAAICuEWaTBO25AAAAokeYTRJhTwHjwQkAAABdIswmicywp4BRmQUAAOgKYTZJZHPMAAAAIGqE2SQRfsyAMAsAANAVwmySCHukLZVZAACALhFmk0RYNwMqswAAAF0izCaJLCqzAAAAUSPMJonQymwTYRYAAKBLhNkkkekJbc1Fn1kAAICuEGaTBE8AAwAAiB5hNkmEnZnlBjAAAIAuEWaTRNiZWcIsAABAlwizSSIr9JiBt01+v2HRagAAAFIDYTZJhD40wTCkfa0+i1YDAACQGgizSSIrpJuBxLlZAACArhBmk0QvjyNsrLmF9lwAAACdIcwmCafDrnRXcKDlJjAAAIDOEWaTCL1mAQAAokOYTSL0mgUAAIgOYTaJhPWapTILAADQKcJsEgltz0VlFgAAoHOE2SQSFmapzAIAAHSKMJtEwh9pS2suAACAzhBmk0h2WvCDE6jMAgAAdI4wm0RCjxnQZxYAAKBzhNkkQp9ZAACA6BBmkwjdDAAAAKJDmE0iWVRmAQAAokKYTSKcmQUAAIgOYTaJ0GcWAAAgOoTZJNLeDWCGYVi0GgAAgORHmE0iWZ7gPrM+v6EDrX6LVgMAAJD8CLNJJLQyK0lNLTwFDAAAoCOE2SQSemZW4iYwAACAzhBmk4jbaZfHGfxbQq9ZAACAjhFmkwy9ZgEAACJHmE0y9JoFAACIHGE2ybTXngsAAADtI8wmmbAHJxygmwEAAEBHCLNJJjOk1yyVWQAAgI4RZpNM6A1gTYRZAACADhFmk0z4MQPCLAAAQEcIs0mGG8AAAAAiZ3mYXbBggYYPH660tDQVFhbqjTfe6HDua6+9JpvNFvZ6//33E7ji+KIyCwAAEDlLw+yyZcs0c+ZMzZkzRxUVFZo4caKKi4tVVVXV6XUffPCBampqAq8RI0YkaMXxF3ZmljALAADQIUvD7Lx58zR16lRNmzZNI0eO1Pz585Wfn6+FCxd2el3//v01YMCAwMvhcCRoxfHHDWAAAACRsyzMer1elZeXq6ioKGi8qKhIa9eu7fTa008/XQMHDtR5552n1atXdzq3paVFjY2NQa9kFt6aiz6zAAAAHbEszNbX18vn8ykvLy9oPC8vT7W1te1eM3DgQC1atEjLly/XihUrdMIJJ+i8887T66+/3uHnlJaWKicnJ/DKz8839ecwG4+zBQAAiJyz6ynxZbPZgt4bhhE2dtgJJ5ygE044IfB+3Lhxqq6u1q9+9SudffbZ7V4ze/ZslZSUBN43NjYmdaDNTg/+LflyX6tWv1+nSSf2t2hFAAAAycuyymxubq4cDkdYFbauri6sWtuZb3zjG9q6dWuHX/d4PMrOzg56JbPj87LUp5c7aGzOyo1q4rG2AAAAYSwLs263W4WFhSorKwsaLysr0/jx4yP+PhUVFRo4cKDZy7OMy2HX7f/vxKCxXQ0H9OArH1i0IgAAgORl6TGDkpISXXPNNRozZozGjRunRYsWqaqqStOnT5d08IjAzp079cwzz0iS5s+fr2HDhunkk0+W1+vVs88+q+XLl2v58uVW/himu2LMEP3pPzv11ke7A2P/9+9PdcnoQTpzeB8LVwYAAJBcLA2zkydP1u7duzV37lzV1NRo1KhRWrVqlQoKCiRJNTU1QT1nvV6vbr31Vu3cuVPp6ek6+eST9de//lUXXnihVT9CXNhsNpV++1RdMP917W/1BcZnLX9XL988UWmuntOKDAAAIBY2wzAMqxeRSI2NjcrJyVFDQ0PSn5996s1t+sVfNgeN/e+5x2lWyDEEAACAniSavGb542zRsevGD9Np+b2Dxha9/one29lgzYIAAACSDGE2iTnsNj34nVPlcnzVqsznN/TTF95Vq89v4coAAACSA2E2yR2fl6UbJ40IGttc06hFr39i0YoAAACSB2E2BfzvucfphLysoLFHX92qj+qaLVoRAABAciDMpgC3064HvnOq7Ec8GM3b5tfty9+V339U3b8HAAAQhDCbIk7L762pZw0PGlv/6R49+/anFq0IAADAeoTZFFLyXydoaJ+MoLEHXn5fO7/cb9GKAAAArEWYTSHpbofu/+9Tgsb2en2as3KjjrJ2wQAAAJIIsyln/Ndy9b2v5weNvfbB53qxcqdFKwIAALAOYTYFzb5wpPpneYLG7vnzZtU3t1i0IgAAAGsQZlNQTrpL9142Kmjsy32tuvulTRatCAAAwBqE2RRVdPIAXXTqwKCxv7xbo79vqrVoRQAAAIlHmE1hd19ysnpnuILG7vzTe2rY32rRigAAABKLMJvC+mV59POLTwoa+6yxRfe/vMWiFQEAACQWYTbFffv0wTrn+H5BY39YV621H9VbtCIAAIDEIcymOJvNpvu+PUq93I6g8dtXbNR+r8+iVQEAACQGYbYHGHJMhmYVnxg0VvXFPs0r+8CiFQEAACQGYbaHuHpsgcYUHBM09tSb21RZ/aU1CwIAAEgAwmwPYbfb9MB3TpXb+dVvqd+QZr3wrrxtfgtXBgAAED+E2R7kuH6Zuvm8EUFjH3zWpCfWfGzRigAAAOKLMNvD/M/Zx+qkgdlBY4/9c6vmlX2oT3fvtWhVAAAA8WEzDMOwehGJ1NjYqJycHDU0NCg7O7vrC1LQezsb9K3fvCWfP/y39uvDjtF3CofowlMGKivN1c7VAAAA1oomrxFme6j7X36/0+MFaS67Ljh5gL5TOETjj8uVw25L4OoAAAA6Fk1ecyZoTUiwmeePUE3Dfv2pcle7Xz/Q6tefKnfpT5W7NCA7Td8+Y7AuP2OIvtY/M8ErBQAA6D4qsz2YYRh6b2ejlm/YoT9V7tSefa1dXnNafm9dXjhEl546SDkZHEMAAACJxzGDThxNYfZI3ja//vl+nZZv2KHV79eprZ3ztEdyO+z6r5PydHnhYJ09op+cDu4VBAAAiUGY7cTRGmaPtLu5RS/9Z5deKN+hTbsau5yfm+nRt08fpMsLh+jEAUfnngEAgMQhzHaCMBvs/dpGLS/foZUVu1Tf3NLl/JMHZevyM4boW6cNUt9MTwJWCAAAjjaE2U4QZtvX5vPr9a2fa3n5TpVt/kxeX+dPDXPabZp0Yn9dfsYQffPE/kFPHgMAAIgFYbYThNmufbnPqz+/W6Pl5TtUWf1ll/OPyXDpW6cd7IYwanC2bLajt81Xq8+v5gNtajrQpqaW1sA/N7e0qelAq5paDr0/cPB986H3TQfaZEjq28ut3Ey3+mZ6lJvpUd9Mt/od+vXwe4/TYfWPCQBAXBFmO0GYjc5Hdc1avmGHVm7YqdrGA13OPyEvS5cXDtZlpw1W/+y0BKzQHD6/ob3e4KDZ1NJ2RBhtDYTOI983H5rTeGjsQGvnFW0zZKU5lZvpORh6e3mUm3Xo10Dg9QQCcXaa86j+wwUAIDURZjtBmO0en9/Q2o/rtbx8h17ZVNtlaLPbpLOP76fvFA7R+SPzlOaKTzXRMAztb/UdETS/qnY2H2g7VAkNrpA2hsxpbjn46oncDntQVbfdau+hQNwnw03XCgBAUiDMdoIwG7umA61atbFGy8t3at32L7qcn53m1MWjB+k7hUN0en7vQKWwpc13RCX04F/Lh/0VfCd/LX84hLb32F5Ez2aTjslwHzrq8FXYDa345h46ApHu5rgDACA+CLOdIMya69Pde7V8w06t2LBDO/bs73J+bqZHfsNQ84G2Lm8yS3V2m5TpcSorzaWsNOehf3YqM82lTI9T2SFjkvTF3hbVN3tV33zw193NLapvbtHuZm+XvYETLcPtCAu9h487HD7zezj85qS7ZOeRyQCACBFmO0GYjQ+/39Db277Q8g07tGpjjfZ5fVYvKSa93A5lpbmUmeYMCqJZnuCx7ENzAl9POxheMz1OZbgdpp1XNQxDDftbA0F3d+DXFn1+ZOjd61V9U4v2Jtn+O+029enlDqnutn+jW99eHrpjAMBRjjDbCcJs/O3ztumV92r1QvkO/euT3Urkv2FpLrsyPa6DVc8jK58e1xFh8+D7w6E061D19PD7Xm6nHCleRdzv9QWF292dVHy/2OdN6O9RJLLTnMrN8ig36Aa39o4+uJXp4SY3AOhpCLOdIMwm1s4v92vlhh16oXyHtu/e1+E8l8MWqGi2V+XMOhROA8EzaOxgUO3lcVLR64Y2n1979rUGVXyPDL27935VDf68uUXetuQ6HuJx2oMqvX17uZWbdfDXflmeoI4PfXq5U/4PKgBwNCDMdoIwaw3DMLRpV6N2fbn/YGANCaoep53qWgowDEPNLW1B1d3wow9fBeLGA8nVJcJmk/pkhHd36OhGt3h14QAAdI4w2wnCLJA4LW0+fbHXq/omr+r3thw68uD96tcjwvAXe71J15ki0+M8ooWZ+9DRh8OV36+qwf0yPcpO57gDAJglmrzmTNCaAByFPE6HBuaka2BOepdz/X5DX+5vPXRTW3uV3kPv97aovsmr/a3xv8ntcPu3Tzs5InOYy2FT317tV3wPHnX46uhDn15uuejpCyQdwzDkNyS/YchvGDIC/3zwV8P/1df8RgTzj/y6XyFzovgefkU03+hqDYHrj/z6wf/+RjP/qrFDdXxeltW/XQGEWQBJwX6o40GfXm6NiOA/kvu8bYFzvEd2d6g/4pjD4fE9+1rjvv5Wn6HaxgMRPSlPknpnuAI9fdvr7nDk0QczO2Ogfe2HgChCh7+bIeWIIHHkfEPRB5kOg4i/vZ8pfkHHCL2+wz054vqIgl4H6zvie0Q1v521IjITR+QSZgEgVhlupzL6OJXfJ6PLuW0+/8HjDiHV3fpDvx7s9vBV+G31xf//al/ua9WX+1r18ed7u5yb5rIHzvP2C3mMcXa6K27VGEPRB5muK0rmBJ1IgmckcwgyQPSS7EQYYRZAz+d02NU/O039s9O6nGsYhhoPtHXez/eIIxBNCXgU8oFWv3bs2R/Rg0kAIN6S7XYrwiwAHMFmsykn3aWcdJeO69f1/AOtPu3eG0l3B6++2NuSdBUNAPFns0l2m01228H/xtgD721BXzv4/sivH5pvb3++TUe8t3dyfaSfZz88/8ivh88f1Lvr+yASiTALADFIczk0uHe6BkfwH3ef39CX+7yB1mbt3uh2qNtDfXOLWpKspy+6p/MgcfC8eHvBw6bOg0xUwcXUoNP5/INjnX1eB0HJHs38dtbfaRDren9s6mSOvRvf84j5nHmPL8IsACSIw247+GCHTI+kzm+eMAxDe72+Tiu+nze3aJ+3LWHVmCNDSuTB6sigEsX8iINY1yEj9NeOgl7U4clmk80eHlxC5wOIL8IsACQhm80WeCJeQd9eVi8HAJIWjQ4BAACQsgizAAAASFmEWQAAAKQswiwAAABSFmEWAAAAKYswCwAAgJRleZhdsGCBhg8frrS0NBUWFuqNN96I6Lq33npLTqdTp512WnwXCAAAgKRlaZhdtmyZZs6cqTlz5qiiokITJ05UcXGxqqqqOr2uoaFBU6ZM0XnnnZeglQIAACAZ2QzDsOxJ4WPHjtUZZ5yhhQsXBsZGjhypyy67TKWlpR1e973vfU8jRoyQw+HQiy++qMrKyog/s7GxUTk5OWpoaFB2dnYsywcAAEAcRJPXLKvMer1elZeXq6ioKGi8qKhIa9eu7fC6p59+Wh9//LHuuuuuiD6npaVFjY2NQS8AAAD0DJaF2fr6evl8PuXl5QWN5+Xlqba2tt1rtm7dqttvv13PPfecnM7InsRbWlqqnJycwCs/Pz/mtQMAACA5WH4DmM1mC3pvGEbYmCT5fD5dddVVuueee3T88cdH/P1nz56thoaGwKu6ujrmNQMAACA5RFbejIPc3Fw5HI6wKmxdXV1YtVaSmpqatH79elVUVOjGG2+UJPn9fhmGIafTqb///e/65je/GXadx+ORx+OJzw8BAAAAS1lWmXW73SosLFRZWVnQeFlZmcaPHx82Pzs7Wxs3blRlZWXgNX36dJ1wwgmqrKzU2LFjE7V0AAAAJAnLKrOSVFJSomuuuUZjxozRuHHjtGjRIlVVVWn69OmSDh4R2Llzp5555hnZ7XaNGjUq6Pr+/fsrLS0tbBwAAABHB0vD7OTJk7V7927NnTtXNTU1GjVqlFatWqWCggJJUk1NTZc9ZwEAAHD0srTPrBXoMwsAAJDcoslrllZmrXA4u9NvFgAAIDkdzmmR1FyPujDb1NQkSfSbBQAASHJNTU3KycnpdM5Rd8zA7/dr165dysrKarefbawaGxuVn5+v6upqjjHECXscf+xxfLG/8ccexx97HH9H8x4bhqGmpiYNGjRIdnvnzbeOusqs3W7XkCFD4v452dnZR92/eInGHscfexxf7G/8scfxxx7H39G6x11VZA+z/AlgAAAAQHcRZgEAAJCyCLMm83g8uuuuu3iEbhyxx/HHHscX+xt/7HH8scfxxx5H5qi7AQwAAAA9B5VZAAAApCzCLAAAAFIWYRYAAAApizALAACAlEWYNdmCBQs0fPhwpaWlqbCwUG+88YbVS0pJpaWl+vrXv66srCz1799fl112mT744IOgOYZh6O6779agQYOUnp6uc889V5s2bbJoxamttLRUNptNM2fODIyxv+bYuXOnrr76avXt21cZGRk67bTTVF5eHvg6+9x9bW1tuuOOOzR8+HClp6fr2GOP1dy5c+X3+wNz2N/ovP7667rkkks0aNAg2Ww2vfjii0Ffj2Q/W1padNNNNyk3N1e9evXSpZdeqh07diTwp0hune1xa2urZs2apVNOOUW9evXSoEGDNGXKFO3atSvoe7DHwQizJlq2bJlmzpypOXPmqKKiQhMnTlRxcbGqqqqsXlrKWbNmjW644Qb9+9//VllZmdra2lRUVKS9e/cG5jz44IOaN2+eHn/8cb3zzjsaMGCA/uu//ktNTU0Wrjz1vPPOO1q0aJFOPfXUoHH2N3Z79uzRhAkT5HK59PLLL2vz5s16+OGH1bt378Ac9rn7HnjgAT3xxBN6/PHHtWXLFj344IN66KGH9NhjjwXmsL/R2bt3r0aPHq3HH3+83a9Hsp8zZ87UypUr9fzzz+vNN99Uc3OzLr74Yvl8vkT9GEmtsz3et2+fNmzYoDvvvFMbNmzQihUr9OGHH+rSSy8NmscehzBgmjPPPNOYPn160NiJJ55o3H777RatqOeoq6szJBlr1qwxDMMw/H6/MWDAAOP+++8PzDlw4ICRk5NjPPHEE1YtM+U0NTUZI0aMMMrKyoxzzjnHuPnmmw3DYH/NMmvWLOOss87q8Ovsc2wuuugi4wc/+EHQ2H//938bV199tWEY7G+sJBkrV64MvI9kP7/88kvD5XIZzz//fGDOzp07DbvdbrzyyisJW3uqCN3j9qxbt86QZHz66aeGYbDH7aEyaxKv16vy8nIVFRUFjRcVFWnt2rUWrarnaGhokCT16dNHkrRt2zbV1tYG7bfH49E555zDfkfhhhtu0EUXXaTzzz8/aJz9NcdLL72kMWPG6IorrlD//v11+umna/HixYGvs8+xOeuss/Tqq6/qww8/lCT95z//0ZtvvqkLL7xQEvtrtkj2s7y8XK2trUFzBg0apFGjRrHn3dTQ0CCbzRb4Gx32OJzT6gX0FPX19fL5fMrLywsaz8vLU21trUWr6hkMw1BJSYnOOussjRo1SpICe9refn/66acJX2Mqev7557Vhwwa98847YV9jf83xySefaOHChSopKdHPfvYzrVu3TjNmzJDH49GUKVPY5xjNmjVLDQ0NOvHEE+VwOOTz+XTffffpyiuvlMS/x2aLZD9ra2vldrt1zDHHhM3h/4XRO3DggG6//XZdddVVys7OlsQet4cwazKbzRb03jCMsDFE58Ybb9S7776rN998M+xr7Hf3VFdX6+abb9bf//53paWldTiP/Y2N3+/XmDFj9Mtf/lKSdPrpp2vTpk1auHChpkyZEpjHPnfPsmXL9Oyzz+r3v/+9Tj75ZFVWVmrmzJkaNGiQrr322sA89tdc3dlP9jx6ra2t+t73vie/368FCxZ0Of9o3mOOGZgkNzdXDocj7E9FdXV1YX+KReRuuukmvfTSS1q9erWGDBkSGB8wYIAksd/dVF5errq6OhUWFsrpdMrpdGrNmjX69a9/LafTGdhD9jc2AwcO1EknnRQ0NnLkyMBNofx7HJvbbrtNt99+u773ve/plFNO0TXXXKNbbrlFpaWlkthfs0WynwMGDJDX69WePXs6nIOutba26rvf/a62bdumsrKyQFVWYo/bQ5g1idvtVmFhocrKyoLGy8rKNH78eItWlboMw9CNN96oFStW6J///KeGDx8e9PXhw4drwIABQfvt9Xq1Zs0a9jsC5513njZu3KjKysrAa8yYMfr+97+vyspKHXvsseyvCSZMmBDWUu7DDz9UQUGBJP49jtW+fftktwf/b8zhcARac7G/5opkPwsLC+VyuYLm1NTU6L333mPPI3Q4yG7dulX/+Mc/1Ldv36Cvs8ftsOrOs57o+eefN1wul/HUU08ZmzdvNmbOnGn06tXL2L59u9VLSzn/+7//a+Tk5BivvfaaUVNTE3jt27cvMOf+++83cnJyjBUrVhgbN240rrzySmPgwIFGY2OjhStPXUd2MzAM9tcM69atM5xOp3HfffcZW7duNZ577jkjIyPDePbZZwNz2Ofuu/baa43Bgwcbf/nLX4xt27YZK1asMHJzc42f/vSngTnsb3SampqMiooKo6KiwpBkzJs3z6ioqAjcSR/Jfk6fPt0YMmSI8Y9//MPYsGGD8c1vftMYPXq00dbWZtWPlVQ62+PW1lbj0ksvNYYMGWJUVlYG/f+vpaUl8D3Y42CEWZP95je/MQoKCgy3222cccYZgVZSiI6kdl9PP/10YI7f7zfuuusuY8CAAYbH4zHOPvtsY+PGjdYtOsWFhln21xx//vOfjVGjRhkej8c48cQTjUWLFgV9nX3uvsbGRuPmm282hg4daqSlpRnHHnusMWfOnKD/6bO/0Vm9enW7/+299tprDcOIbD/3799v3HjjjUafPn2M9PR04+KLLzaqqqos+GmSU2d7vG3btg7//7d69erA92CPg9kMwzASVwcGAAAAzMOZWQAAAKQswiwAAABSFmEWAAAAKYswCwAAgJRFmAUAAEDKIswCAAAgZRFmAQAAkLIIswAAAEhZhFkAiMG5556rmTNnJvxzt2/fLpvNpsrKStO+57BhwzR//nzTvh8AJILT6gUAwNHutdde06RJk7Rnzx717t3bsnW888476tWrl2WfDwDdQZgFAEiS+vXrZ/USACBqHDMAgBi1tbXpxhtvVO/evdW3b1/dcccdMgwj8PVnn31WY8aMUVZWlgYMGKCrrrpKdXV1kg4eF5g0aZIk6ZhjjpHNZtN1110nSfL7/XrggQf0ta99TR6PR0OHDtV9990X9NmffPKJJk2apIyMDI0ePVr/+te/Ol3r3XffraFDh8rj8WjQoEGaMWNG4GtHHjNYunSpbDZb2Ovuu+8OzH/66ac1cuRIpaWl6cQTT9SCBQu6u4UA0G2EWQCI0e9+9zs5nU69/fbb+vWvf61HHnlETz75ZODrXq9Xv/jFL/Sf//xHL774orZt2xYIrPn5+Vq+fLkk6YMPPlBNTY0effRRSdLs2bP1wAMP6M4779TmzZv1+9//Xnl5eUGfPWfOHN16662qrKzU8ccfryuvvFJtbW3trvOFF17QI488ot/+9rfaunWrXnzxRZ1yyintzp08ebJqamoCrz/84Q9yOp2aMGGCJGnx4sWaM2eO7rvvPm3ZskW//OUvdeedd+p3v/tdTHsJAFEzAADdds455xgjR440/H5/YGzWrFnGyJEjO7xm3bp1hiSjqanJMAzDWL16tSHJ2LNnT2BOY2Oj4fF4jMWLF7f7PbZt22ZIMp588snA2KZNmwxJxpYtW9q95uGHHzaOP/54w+v1tvv1goIC45FHHgkb/+ijj4y+ffsaDz74YGAsPz/f+P3vfx807xe/+IUxbty4dr83AMQLlVkAiNE3vvEN2Wy2wPtx48Zp69at8vl8kqSKigp961vfUkFBgbKysnTuuedKkqqqqjr8nlu2bFFLS4vOO++8Tj/71FNPDfzzwIEDJSlwhCHUFVdcof379+vYY4/VD3/4Q61cubLDKu5hDQ0Nuvjii1VcXKzbbrtNkvT555+rurpaU6dOVWZmZuB177336uOPP+70+wGA2bgBDADiaO/evSoqKlJRUZGeffZZ9evXT1VVVbrgggvk9Xo7vC49PT2i7+9yuQL/fDhQ+/3+dufm5+frgw8+UFlZmf7xj3/oxz/+sR566CGtWbMm6Psc5vP5NHnyZGVnZ2vx4sWB8cPff/HixRo7dmzQNQ6HI6J1A4BZCLMAEKN///vfYe9HjBghh8Oh999/X/X19br//vuVn58vSVq/fn3QfLfbLUmBSq4kjRgxQunp6Xr11Vc1bdo009aanp6uSy+9VJdeeqluuOEGnXjiidq4caPOOOOMsLm33HKLNm7cqHfeeUdpaWmB8by8PA0ePFiffPKJvv/975u2NgDoDsIsAMSourpaJSUl+tGPfqQNGzboscce08MPPyxJGjp0qNxutx577DFNnz5d7733nn7xi18EXV9QUCCbzaa//OUvuvDCC5Wenq7MzEzNmjVLP/3pT+V2uzVhwgR9/vnn2rRpk6ZOndqtdS5dulQ+n09jx45VRkaG/u///k/p6ekqKCgIm/v0009rwYIFWrlypex2u2prayUpcKTg7rvv1owZM5Sdna3i4mK1tLRo/fr12rNnj0pKSrq1PgDoDs7MAkCMpkyZov379+vMM8/UDTfcoJtuukn/8z//I+lg79alS5fqj3/8o0466STdf//9+tWvfhV0/eDBg3XPPffo9ttvV15enm688UZJ0p133qmf/OQn+vnPf66RI0dq8uTJHZ6HjUTv3r21ePFiTZgwQaeeeqpeffVV/fnPf1bfvn3D5q5Zs0Y+n0+XXnqpBg4cGHgdXvu0adP05JNPaunSpTrllFN0zjnnaOnSpRo+fHi31wcA3WEzjCOaIQIAAAAphMosAAAAUhZhFgAAACmLMAsAAICURZgFAABAyiLMAgAAIGURZgEAAJCyCLMAAABIWYRZAAAApCzCLAAAAFIWYRYAAAApizALAACAlPX/AYpIpL2JPuN7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "df = pd.DataFrame(data=timings, columns=['batch_size', 'time'])\n",
    "\n",
    "plt.figure(figsize=(8,6))\n",
    "plt.plot(df.batch_size, df.time, linewidth=3)\n",
    "plt.xlabel('batch size')\n",
    "plt.ylabel('time')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b762ada-44d1-41e5-b2b0-142ff9402977",
   "metadata": {},
   "source": [
    "# code for 2B\n",
    "# code for 2B\n",
    "import numpy as np\n",
    "import math\n",
    "\n",
    "# total number of jobs\n",
    "print(len(r_sq_list))\n",
    "\n",
    "# for batch of size 1, there are 987 job\n",
    "#987*b + rw = 1.08\n",
    "# for a batch of size 128, there are 8 jobs\n",
    "#8*b + rw = .42\n",
    "\n",
    "\n",
    "      \n",
    "ar =np.array(timings)\n",
    "\n",
    "\n",
    "\n",
    "#result = (ar[:,1]-ar[-1,1])\n",
    "#result2 = np.ceil(987/ar[:,0])-4\n",
    "#print(result2)\n",
    "#result/result2\n",
    "ar"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d662db6-2c99-45c2-a71b-a4482c454371",
   "metadata": {},
   "source": [
    "**TODO Answer for 2B**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57b3cb77-8d34-4eab-bc1a-09fd8266a07f",
   "metadata": {},
   "source": [
    "3. Run the job with `prefer='threads'` and `prefer='processes'`. Given that `LinearRegression.fit()`. You do not need to vary batch size for this part. Which is more efficient?  Why? Consider our discussion of parallel threads in python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "192d3aed-d2c7-45c8-83bd-83462fdc24bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "585 ms ± 21.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit \n",
    "# code for 1B \n",
    "r_sq_arr = Parallel(n_jobs=4, prefer='processes')(delayed(r_sq_regression)(combo) for combo in combos)\n",
    "r_sq_idx = np.argmax(r_sq_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "39bccac6-7435-4034-b320-ff15b329b10a",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.17 s ± 34 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit \n",
    "# code for 1B \n",
    "r_sq_arr = Parallel(n_jobs=4, prefer='threads')(delayed(r_sq_regression,)(combo) for combo in combos)\n",
    "r_sq_idx = np.argmax(r_sq_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b831431-437f-48e8-b24c-3952a8ab2ebf",
   "metadata": {},
   "source": [
    "**TODO Answer to 3**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23e2c81e-aa2d-487e-bbce-f5c48ab45b0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "Processes is faster, because threads holds on to the GIL."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}