{ "cells": [ { "cell_type": "markdown", "id": "fd45f376", "metadata": {}, "source": [ "# Choosing Features and Avoiding Overfitting" ] }, { "cell_type": "code", "execution_count": 1, "id": "194d625a", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:29.573484Z", "iopub.status.busy": "2026-10-24T05:29:29.573247Z", "iopub.status.idle": "2026-10-24T05:29:30.442743Z", "shell.execute_reply": "2026-10-24T05:29:30.442280Z" }, "tags": [ "hide-input", "thebe-init", "runnable" ] }, "outputs": [], "source": [ "path_data = '../../data/'\n", "\n", "import numpy as np\n", "import pandas as pd\n", "from itertools import combinations\n", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "\n", "from sklearn.linear_model import LinearRegression\n", "from sklearn.model_selection import train_test_split, cross_val_score\n", "\n", "def rmse(actual, predicted):\n", " \"Root mean squared error, in the units of the target.\"\n", " return np.sqrt(np.mean((actual - predicted) ** 2))" ] }, { "cell_type": "markdown", "id": "6ee50f61", "metadata": {}, "source": [ "```{contents}\n", ":local:\n", ":depth: 2\n", "```" ] }, { "cell_type": "markdown", "id": "9d69f820", "metadata": {}, "source": [ "```{index} feature selection\n", "```\n", "\n", "A real estate analyst at a brokerage wants a model that prices homes before they are listed. In {doc}`1302-multiple-regression` you saw that adding predictors can help: a home's price depends on more than its size. That raises an obvious question. If three predictors are better than one, why not use thirty, or all the data the county records?\n", "\n", "This section shows why \"more is better\" is the wrong rule. You will see that a model's error on the data it was fitted to always falls as features are added, but its error on new homes can rise. You will learn to spot this problem, called **overfitting**, to choose among models fairly with a **validation set**, and to weigh accuracy against how easily a model can be explained.\n", "\n", "The data are the same Ames, Iowa home sales used in {doc}`1301-linear-regression`: real sales from 2006 to 2010, recorded by the Ames Assessor's Office and published by Dean De Cock in the *Journal of Statistics Education* (2011). The file keeps normal sales only (foreclosures, family sales, and similar are removed) and drops the few homes over 4,000 square feet that De Cock recommends removing. That leaves 2,412 homes. Prices are in dollars; we report errors in thousands of dollars." ] }, { "cell_type": "code", "execution_count": 2, "id": "9cdc0149", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.444074Z", "iopub.status.busy": "2026-10-24T05:29:30.443960Z", "iopub.status.idle": "2026-10-24T05:29:30.449655Z", "shell.execute_reply": "2026-10-24T05:29:30.449273Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2412 homes\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sale_priceliving_areaoverall_qualityyear_builtgarage_carstotal_roomsfull_bathslot_areaneighborhoodyear_sold
021500016566196027131770NAmes2010
11050008965196115111622NAmes2010
217200013296195816114267NAmes2010
324400021107196828211160NAmes2010
418990016295199726213830Gilbert2010
\n", "
" ], "text/plain": [ " sale_price living_area overall_quality year_built garage_cars \\\n", "0 215000 1656 6 1960 2 \n", "1 105000 896 5 1961 1 \n", "2 172000 1329 6 1958 1 \n", "3 244000 2110 7 1968 2 \n", "4 189900 1629 5 1997 2 \n", "\n", " total_rooms full_baths lot_area neighborhood year_sold \n", "0 7 1 31770 NAmes 2010 \n", "1 5 1 11622 NAmes 2010 \n", "2 6 1 14267 NAmes 2010 \n", "3 8 2 11160 NAmes 2010 \n", "4 6 2 13830 Gilbert 2010 " ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "homes = pd.read_csv(path_data + 'ames_housing.csv')\n", "print(len(homes), 'homes')\n", "homes.head()" ] }, { "cell_type": "markdown", "id": "92789008", "metadata": {}, "source": [ "## More Features, Better Fit?\n", "\n", "```{index} training error\n", "```\n", "\n", "```{index} training RMSE\n", "```\n", "\n", "Imagine the analyst has records of only 100 sales to build the model. We take a random sample of 100 homes as the **training set**, the data used to fit a model, and set the other 2,412 - 100 = 2,312 homes aside. They play no role until the next section." ] }, { "cell_type": "code", "execution_count": 3, "id": "243a5012", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.450568Z", "iopub.status.busy": "2026-10-24T05:29:30.450510Z", "iopub.status.idle": "2026-10-24T05:29:30.456947Z", "shell.execute_reply": "2026-10-24T05:29:30.456611Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Training homes: 100\n", "Set aside for later: 2312\n" ] } ], "source": [ "train_rows, other_rows = train_test_split(homes.index, train_size=100, random_state=1)\n", "train = homes.loc[train_rows]\n", "print('Training homes:', len(train))\n", "print('Set aside for later:', len(other_rows))" ] }, { "cell_type": "markdown", "id": "086c8499", "metadata": {}, "source": [ "The first column is the most useful predictor: living area. We add the others one at a time, in the order of how an analyst might rank them (size, then quality, then age, then the amenities), refit the model each time, and record the **training RMSE**, the root mean squared error on the same homes used to fit the model. (RMSE was defined in {doc}`1301-linear-regression`.)" ] }, { "cell_type": "code", "execution_count": 4, "id": "526d7ab9", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.457941Z", "iopub.status.busy": "2026-10-24T05:29:30.457881Z", "iopub.status.idle": "2026-10-24T05:29:30.471139Z", "shell.execute_reply": "2026-10-24T05:29:30.470809Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
featureslast feature addedtraining RMSE ($1,000s)
01living_area44.5
12overall_quality36.8
23year_built32.8
34garage_cars32.0
45full_baths31.1
56lot_area24.5
67total_rooms24.4
\n", "
" ], "text/plain": [ " features last feature added training RMSE ($1,000s)\n", "0 1 living_area 44.5\n", "1 2 overall_quality 36.8\n", "2 3 year_built 32.8\n", "3 4 garage_cars 32.0\n", "4 5 full_baths 31.1\n", "5 6 lot_area 24.5\n", "6 7 total_rooms 24.4" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "numeric = ['living_area', 'overall_quality', 'year_built', 'garage_cars',\n", " 'full_baths', 'lot_area', 'total_rooms']\n", "\n", "results = []\n", "for k in range(1, len(numeric) + 1):\n", " chosen = numeric[:k]\n", " model = LinearRegression().fit(train[chosen], train['sale_price'])\n", " train_rmse = rmse(train['sale_price'], model.predict(train[chosen]))\n", " results.append([k, chosen[-1], round(train_rmse / 1000, 1)])\n", "\n", "pd.DataFrame(results, columns=['features', 'last feature added', 'training RMSE ($1,000s)'])" ] }, { "cell_type": "markdown", "id": "8c94ee64", "metadata": {}, "source": [ "Each added feature lowers the training RMSE, from about \\$44,500 with living area alone to about \\$24,400 with all seven. This is not a coincidence. Least squares picks the coefficients that make the squared errors on the training data as small as possible. A model with an extra feature can always set that feature's coefficient to zero and reproduce the smaller model, so the best it can do is the same or better. **Training error never rises when a feature is added.**\n", "\n", "That is true even when the feature is useless. To see this, we add columns of pure random noise: numbers drawn from a normal distribution that have nothing to do with the sales." ] }, { "cell_type": "code", "execution_count": 5, "id": "b3d8b2ad", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.472052Z", "iopub.status.busy": "2026-10-24T05:29:30.471997Z", "iopub.status.idle": "2026-10-24T05:29:30.479306Z", "shell.execute_reply": "2026-10-24T05:29:30.479022Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 noise columns: training RMSE = 24.39 ($1,000s)\n", "2 noise columns: training RMSE = 24.06 ($1,000s)\n", "5 noise columns: training RMSE = 23.56 ($1,000s)\n", "10 noise columns: training RMSE = 20.64 ($1,000s)\n" ] } ], "source": [ "rng = np.random.default_rng(0)\n", "noise = pd.DataFrame(rng.normal(size=(len(homes), 10)),\n", " columns=['noise' + str(i) for i in range(1, 11)],\n", " index=homes.index)\n", "\n", "for j in [0, 2, 5, 10]:\n", " chosen = numeric + list(noise.columns[:j])\n", " X_train = pd.concat([train[numeric], noise.loc[train_rows]], axis=1)[chosen]\n", " model = LinearRegression().fit(X_train, train['sale_price'])\n", " print(j, 'noise columns: training RMSE =',\n", " round(rmse(train['sale_price'], model.predict(X_train)) / 1000, 2), '($1,000s)')" ] }, { "cell_type": "markdown", "id": "98a498a5", "metadata": {}, "source": [ "Ten columns of random numbers lower the training RMSE from \\$24,390 to \\$20,640. The model \"found\" patterns in the noise because with only 100 homes, some combination of random columns always lines up a little with the prices. A manager who judged the model by its training error would conclude that the noise columns are valuable. They are not, and the next section shows how to see that.\n", "\n", "```{admonition} Key idea\n", "Training error rewards every feature, useful or not. It cannot tell you which features to keep.\n", "```" ] }, { "cell_type": "markdown", "id": "b4a610e1", "metadata": {}, "source": [ "## Overfitting\n", "\n", "```{index} overfitting\n", "```\n", "\n", "```{index} test set\n", "```\n", "\n", "What we want to know is how well the model prices homes it has **not** seen. The 2,312 homes we set aside are the **test set**: they were not used to fit the model, so the RMSE on them estimates how the model will do on new listings (see {doc}`1301-linear-regression` for the train/test idea).\n", "\n", "We also let the model become much more flexible. Beyond the seven numeric features, we add:\n", "\n", "- one column for each of the 28 neighborhoods (1 if the home is in that neighborhood, 0 otherwise),\n", "- the square of each numeric feature, which lets the model bend,\n", "- the product of each pair of numeric features, which lets the effect of one feature depend on another (an **interaction**).\n", "\n", "All numeric features are standardized first (as in {doc}`../12-regression/1201-correlation`) so squares and products are on a comparable scale." ] }, { "cell_type": "code", "execution_count": 6, "id": "2d1b09ca", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.480279Z", "iopub.status.busy": "2026-10-24T05:29:30.480222Z", "iopub.status.idle": "2026-10-24T05:29:30.485370Z", "shell.execute_reply": "2026-10-24T05:29:30.485034Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "numeric columns: 7\n", "neighborhood columns: 28\n", "squared columns: 7\n", "product columns: 21\n", "all columns: 63\n" ] } ], "source": [ "# Standardize the numeric columns so squares and products are on one scale\n", "Z = (homes[numeric] - homes[numeric].mean()) / homes[numeric].std()\n", "\n", "neighborhoods = pd.get_dummies(homes['neighborhood']).astype(float)\n", "squares = pd.DataFrame({a + '^2': Z[a] ** 2 for a in numeric})\n", "products = pd.DataFrame({a + ' x ' + b: Z[a] * Z[b] for a, b in combinations(numeric, 2)})\n", "\n", "X_all = pd.concat([Z, neighborhoods, squares, products], axis=1)\n", "n_numeric, n_nbhd, n_sq = len(numeric), neighborhoods.shape[1], squares.shape[1]\n", "\n", "print('numeric columns: ', n_numeric)\n", "print('neighborhood columns:', n_nbhd)\n", "print('squared columns: ', n_sq)\n", "print('product columns: ', products.shape[1])\n", "print('all columns: ', X_all.shape[1])" ] }, { "cell_type": "markdown", "id": "dd458901", "metadata": {}, "source": [ "The columns are ordered from the simplest to the most flexible, so using the first $k$ columns gives a sequence of larger and larger models. For each, we fit on the 100 training homes and compute the RMSE on the training homes and on the test set." ] }, { "cell_type": "code", "execution_count": 7, "id": "6db17c69", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.486295Z", "iopub.status.busy": "2026-10-24T05:29:30.486240Z", "iopub.status.idle": "2026-10-24T05:29:30.508263Z", "shell.execute_reply": "2026-10-24T05:29:30.507909Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
columnstraining RMSEtest RMSE
0144.548.1
1236.835.4
2332.834.0
3432.033.3
4531.133.0
5624.535.0
6724.434.8
73518.735.8
84215.6147.0
96312.6181.6
\n", "
" ], "text/plain": [ " columns training RMSE test RMSE\n", "0 1 44.5 48.1\n", "1 2 36.8 35.4\n", "2 3 32.8 34.0\n", "3 4 32.0 33.3\n", "4 5 31.1 33.0\n", "5 6 24.5 35.0\n", "6 7 24.4 34.8\n", "7 35 18.7 35.8\n", "8 42 15.6 147.0\n", "9 63 12.6 181.6" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sizes = [1, 2, 3, 4, 5, 6, 7,\n", " n_numeric + n_nbhd,\n", " n_numeric + n_nbhd + n_sq,\n", " X_all.shape[1]]\n", "\n", "test_rows = other_rows\n", "curve = []\n", "for k in sizes:\n", " cols = list(X_all.columns[:k])\n", " model = LinearRegression().fit(X_all.loc[train_rows, cols], homes.loc[train_rows, 'sale_price'])\n", " tr = rmse(homes.loc[train_rows, 'sale_price'], model.predict(X_all.loc[train_rows, cols]))\n", " te = rmse(homes.loc[test_rows, 'sale_price'], model.predict(X_all.loc[test_rows, cols]))\n", " curve.append([k, round(tr / 1000, 1), round(te / 1000, 1)])\n", "\n", "curve = pd.DataFrame(curve, columns=['columns', 'training RMSE', 'test RMSE'])\n", "curve" ] }, { "cell_type": "markdown", "id": "280a2fbb", "metadata": {}, "source": [ "Training RMSE falls at every step, from about \\$44,500 to about \\$12,600. Test RMSE tells a different story. It falls while the first few features are added, reaches its lowest value of about \\$33,000 with 5 columns, and then rises: \\$35,000 at 6 columns, and far worse for the most flexible models, about \\$147,000 at 42 columns and \\$182,000 at 63. The 63-column model fits 100 homes with 63 coefficients and a constant, so it chases the quirks of those particular homes. When it prices other homes, its typical error is as large as the average price of a home (about \\$175,000)." ] }, { "cell_type": "code", "execution_count": 8, "id": "98bdaf0b", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.509223Z", "iopub.status.busy": "2026-10-24T05:29:30.509165Z", "iopub.status.idle": "2026-10-24T05:29:30.562971Z", "shell.execute_reply": "2026-10-24T05:29:30.562638Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAGMCAYAAACI+5UEAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAfnlJREFUeJztnQd4U9Ubxr/uRWkLlBbKRkC27CEgS0ABGTJENu6NA3Ey/oI4UHHiZCiggoICKqKITNnIkL13gW5auvN/3tPekKRpm7SZ7ft7nts0J3ecO3Lz3u98w0On0+mEEEIIIYS4NJ7O7gAhhBBCCCkcijZCCCGEEDeAoo0QQgghxA2gaCOEEEIIcQMo2gghhBBC3ACKNkIIIYQQN4CijRBCCCHEDaBoI4QQQghxAyjaSJE4deqUTJkyRfbu3euU5YlIdna2OobLly8vFYdjzZo18u6776p93rVrl1XLrly5Ui2HY2ZtGyHE9dm8ebP67l6+fFlKMh6siFAy2L9/v/zwww8WzVujRg0ZM2ZMsbb3999/S5cuXeSbb76RESNGOHx5IpKZmSk+Pj5y3333yZdfflmiDwn2cc6cOfr3X3zxhdx///0WL//www/LZ599JhkZGeLt7W1Vm6vw008/yb///qt+mErStlyFAwcOyOLFi+Xee++VunXrOrs7xEpmzpwpEyZMkH379kmjRo1K7PFzrbsSKZZomzp1qkXz3nbbbcUWbRB+kydPliZNmjhleVJ6OHHihBJsuGaGDBkiAQEB0rx5c7tsq0+fPhIZGSmenq43CAEhNX/+fIeJNkdty5VEG+6hLVu2pGgjLgtFWwkBTxYQQYZ89913cvjwYXn22WelTJky+nb8+BUXrKM4N/TiLk9KD7iGwYwZM+See+6x67Yg2jARQogrQtFWgkSbqUkYwxv4wXvuueeU9QAcOnRIiTm81qpVS1asWKHm6dq1q7Rt21bNA38hTFevXpXatWvLnXfeKUFBQXl80ubNmycDBw7UW8u0deOHFcvB1+ro0aNStWpVNR8sJLZaXuPgwYPy+++/i06nU8Ott9xyixoqTE1Nlccff9yiY5ecnKzWgW35+flJx44dpUWLFkbzGPbN9LhdunRJP5QUHR2t+n3x4kV58sknJTQ0VPVl1apVan6sv02bNtKuXTuj9cOH6n//+5+yIN11113KP+Off/6R8PBwGTVqlFgDjsWnn36qfDswrFilShX9Mf/zzz9VO64HnFftujDFknkxHI/j8sorr8jp06eVP1hKSop06NAhz/6B48ePy7p16+TKlSvqnOI447Ug4MOGYwFwXLE9X19feemll1Rbenq6/PXXX+o6AK1atVLbLyrYhx07dsikSZP01jbD/cR5hRUqMTFR7WOnTp3MHn+cb/hrli9fXgYMGCBeXl4ya9Ys6dy5s5qsBecT1xgwfNiB5bFBgwb695ae44LOhaXbyg9rr50LFy6o7xO+O88884w6TjfffLMMGjRItm3bJps2bZLg4GD9cPj58+fV8cX8FStWlJ49e+a5jnDsly5dqr47lStXVtfOsWPHpFevXmattLiGMDQKFi1apK4BgPtJ//79ZePGjWqf8J0uV66c0bIfffSR+Pv76/tX0LbxignXL/b7559/lqSkJHUt4RyYgnnXr1+vP0+43rTvc0HgHojvxqOPPmrUjns6+tu9e3f99wT7rfUJ5w7fgevXr6tt4V5l6NKCafz48eq+pjF9+nTlTvDCCy+o46C5bqAdVsvevXvr501LS1O+qdqDGNbfvn37PP037BPON/qE8/3cc89JYGCgxMfHq2OHe2/9+vWNtlHigU8bKZn069dPh1N88eJFfduyZctU26effqpr2LCh+h/Tu+++q7t8+bKudevW+jZtqlSpku7ff/81WvfatWvVZ998802edX/55Ze6Zs2aGa3j5ptv1sXExNhsefDWW2/pPD09jeabMWOGrk2bNrp69epZdIxWr16tq1ixYp59HjlypC4jI8Oi4zZ69Gj1/++//64LCgrSf3bixAl13KpVq5Zn/b169dIlJibq149tof2+++7TjR07Vj9f165d8+274TIaaWlpuqFDh+q8vLx0n3/+uWrLzs7WTZw4UbUZ9iEwMFC3aNEio3VaMy+24+fnp/vpp5/Uq+H86BPWpTFt2rQ858rX11c3c+bMAs9PVFRUnmOHYwyWL1+ui4iIyPP5nXfeqbt+/brReh566CH1meE5tbRN209cK4bnF9MTTzxhtJ1r167punTpYjRPuXLldD/88IP6f/Lkyfp58b3E+++//15XGE2bNs2zn5iWLFli9Xkr7FwUtq38KMq1s2LFCl1AQIB+3gsXLqh2fP7444/r23FfArNmzVJ9Ne37O++8Y7R+3Ffw2dy5c3V16tTRz6t9J0x5+eWXze4zvtsA9xW8P3r0aJ5lca/BPceSbd99993qGvrtt9/UcTHc1tNPP2203ilTppg9TzgGhYH+4D5lysGDB/X3SQ2tTz/++GOeY4vvg/Y9xvccbZhPY//+/fp516xZo2/ftGlTnuONe6i5e23fvn11qampRv3U+oRt+vv76+e9fPmybseOHXnW06lTJ93//vc/9f++fft0JRla2kopr776qnh4eKinw6ioKGVli4mJke3bt6snLFjt8NSEJ0486Q0bNkz+++8/tUxh4OkIloUHH3xQPRXhiQhP1W+++aaabLE8+jRx4kT1P56Eb7rpJuXX9/LLL6unQFioCuPIkSNqWVjC8LTeuHFj9YSJp34ESOBpX7PoFHTc9uzZoz5DQEW9evWUz2DZsmWVVa1v375y9uxZ9TR4++23y7Vr19RTOCwFTzzxhLI2GoKn8ri4OBk8eLDyq6lTp45YCqw/2J8tW7bIsmXL1LbBO++8o44bjguseHj6h2UM/YBvI6wJ6J+182pP1GiHtRSWAjwNw4ry1VdfKWvlI488otpgUYGlFBbTatWqqadnnENYEwsClhdY2n788Ud1DeKYwNIGdu7cqc4XLFk4/7BY/Prrr2rCPpi6CxQH7CesrPhe3HrrrRIbGyvffvutfPjhhzJ27Fhp1qyZmu/555+XtWvXSoUKFdS5wD6vXr1aBTiYAisBfKjuvvtuZcUqCCwPCxiuNcP90ixflp43S85FYdvKj6JcO7BGNWzYUN1zYE3T3Dhg+cL9COcWy8CdApYuWHkArHfoDyysv/zyi3IBwXtYswzB+UCwDu4llSpVytcXEhZz3A+WLFmiv84A+lxUzG0bFn1YwBDs0LRpU2Vhw37iWnrvvffUtYT7EKxwsHLi/ofzBCvbuXPn9OfpqaeeEluCPo0bN05dx7i+cW3iO4egHHyPH3jgAWUhxn0Z5wZ9Avgf1xLud7Cg4Thq7UB7D/BbAksblsVICr6vOHe43yKIAPdu0z6NHj1a3VtgFcT14enpqb4vsOLi+GH9+C7inqNZ20s8zlaNxDmWNlh/rl69ajQ/LFmmFjUwadIktcy2bdssspTVr19fl5CQYLTe0NBQZT2z1fKDBw9W83733XdGfcVTKNotsbTdf//9al489RqSkpKithUZGWnRcdMsbbDOZWVl6du//vpr1d67d29denq6vv3MmTO68uXLK4tEdHS0kdXMw8PD6Im1IAwtbefPn9c1adJErXfz5s36ebDdsLAw9cR/5coVo+XxxIo+PPXUU1bPC2ANwfbHjRtnZFXbvn27ztvbW1lHwc6dO9V8n3zyidE6sb0tW7YUup/ffvutWh5WGUN2796ti42NNWqDhQ2Wolq1atnU0oa2Z5991midCxcuVO1vvvmm3soGqwCsf5cuXdLPByvCbbfdVixLm+F1Zoo1583Sc5HftvKjqNcOvoOG1w7QrLawpBoCCyrav/rqK6P2OXPmqPbu3bvnsXbVrVtXFx8fb9E+wJJo7jorqqXN3LZhQcJnsEgaMn/+fNWuWQy3bt1q1jJo6XfGWksb2h544AGjc4Ht4Lw1atRI39ayZUt1f9bo37+/rlu3brp7771X17ZtW307RghwrzRk165duri4OKO25ORktX4cK3PH6eGHHzbq0+LFi1U7Ps/MzDTaL81iS0sbKZHA1wH+NobAVyMkJET5LcAvIyEhQbKyspQfD4ClDT5DlqwbT16G68XTGpa3tG+FLQ9fF1hXhg4darQsrFeWRtFiP2EVgGUKkyF4qtu9e7fadzwlF3TcNPCkaBh1iD4C+EfhiVsDT82wQE2bNk09fcJqoIGnXMOnU0uA7wee2GEBhO8NLIQaOI+w3MGaAV8WU3C+sZ/WzmvO+qih+bHAQoprCOvDkzWehmG5wDWEdBo4JoY+M9aCdaG/WC8iTGF1gz8ZrBOwEuG9OT/IogKfHUN69OihXrXvB1INwGoLv5uIiAj9fLC4Yln4kBkCPy9bBONYc97sdS6Keu3A6mfOeo++aZZiw+8T/ElhETIE1qnXX39dtm7dmmc98EHDtp1BQdsu7FqCtQ3WRfh24X98p2xxngoC9ynDc4HtwHIJyzVGCGAFxb3prbfeUpY4+BPimkaaDVzv33//vbL243qHddw0aAhWPFjFYM08efKk/vsKf2nc2+EXZ3ifNHd9bM09x2iH1U8D9zxsb+7cuVLS4fBoKQU3blMwlAEBgdB3c+ALaQnmHGXxhceX0lbLwxHVNFgAQDRBFMEMXxhwysV6ChJ52GdD0WbuuOX3GX7EQM2aNfPMq7Vp81iy/vzAsDCGWLThWdN9BBBzmAo6r9bMq4EbKobYTKlevbp6xfHF/1gffliHDx+uhn4wVIRhkscee6zIwgpDShi2QSCJOTD8YivRhh8iDHkaog3ladcl9hUUdDzsgTXnDfthj3NRlGsHIsTcscrve4DvSn4pgvB9wsMLHjINf8yL8n2yFfltGyLF0JHf3LWE8wDhg/OE4VqIJNzvcJ7w4Kg5/NsKHDO4e5i7biGs8PBlKNow/Il7Dc4J2vAAgmMPEYf58PBi+vC5YMECeeihh1SwUn7fV8MgD1yrpn2Kz/2Omfs+2SIrgjtA0VZKMZc4FBYCCDZYbfBkh6dEfJnx9Pf555+7VJZ49O3MmTN52tFH+JDhKbAwwsLClH8UrF75YWpVKyjhquln2o0ZT5WmPnZo0/pg6frzA75T+NHC0yeWhz+ZZvHT1o8bqLnoNKBF9lkzrwZu6DgPpjdMPAAYHgMs98EHH6gJN3r8sON6g78T/rc2Lxpu8PArxA8FogyxfVjYICL/+OMP9YPn6OtVs6qYuy6142EPrD1vtj4XRekDwL0lPx9Zc98DXEuIbjQH2mEdNxRs+a2nKGhi1vQBAdc/fAJNRZgtto2HRfhMYoKFasOGDeo8wV8Vvm0F+Rejv/BfNAX3RnPge4T9MH1g1q5b7drGucU9E6IN82NEBL8VOO64B6Ed5wEYijaIPjxg4XjBXxeiS/u+ws8PIx2m39f8rgGtX6b3zvyujZIGRRvRg+EdOOBqKRY0IAZcDZjuYWbHhJuAxieffKJ+iCwRbQgYQMABUoXgf0NgqUOqAVPrirV9/Pjjj9UwKH4QtZsQbnZw9MaNDjc8W4ChWQwtIDgDT+tIjIr1Y5gJN1ykDIDzsumNDs7NmrXPmnkNee2111SKAe1HBEO+cDDGkAXWBwsIbsiaczfWi6EvDE8jpQeGNjHUbQ1IWYEndgzpGFpK8aOKlBzOAFYgWEBwLDBMr12DcKiGdaK4wPIAcA4Mz401583Sc5HftvKjqNeOtd8nDNXh2oaDugbeI10P0lgUB8N9NkWzCEJgwAFeA0PBGDq0NdgfoAUiwQLVr18/JYogtiFaCrIsob8Qdgjg0twlcF9AsEh+IOUQAg+07zGGIhEwheAbzRIIodW6dWsVdIBzivumJpS7deum2iHkcH0ZWsmwP7C+YRtwp9DAsUPAg6W0yR0aRjoRpGDStq2lZCoNULQRPfiiIZIHUY4QE/A5QKSSKz7BYCgHgg1+DBgm06JHcVPFDc6SKNcXX3xRrQM3G0RG4eaEp0j8cOGGB/8ZU183a4AFCGIKxxQ/avhR0aJH8cOAHx5LxKU10Wro/9NPP60i8zAcgfe4SeIJHftzxx13qOFjiBvkSsKP9dtvv638ZqyZVwM3TewPBD+ewhHVhWOK7cOnByDnF6IjkY8JxxjnB5F6EFcQsuasFIWB/mFZ/AjBsoVIRWwbYhFWA2eAYS/4V82ePVvtJyIfIeIQPaoNHxpelxjygnhH1GNh0aPaPgMMkcH3Efuv5U6z9LxZei4K2pY5inLtWAuua4g2RKLCJw9RpfixRg4vLdK4OGj7jP3A9QyBouVpw7WN84t7BnJYQhThWOI+YQ+fOURGI8IUxx7nCSIY+wo/UTycFbZN+KJ9/fXXqt+4R+K6w3WoDS+agnXCfw6+iYjUxLWJY4zvkvY91oAFDeILD5+GDyNoh5Uf9wTTEnMYKkY75scoA6yIsATi3MH6Zil9+/ZVVjr0Db9ReODWokchLPGbVeIpNAyFlMjoUbyasnfvXhV9aJj/pmrVqvocU++9955F0Z/m1o2+YN22Wh4gAso0j9Hrr7+ua9Cgga5FixYWHaN169apfTTNHYQowJdeesmivhUUaYcIR2vztFlKfst89NFHKgp14MCB+qjVqVOn5snBpEXDmkarWjqvlmsLx8R0fsM8bYhIRiSp6fqQpwp574oaPfr222/nOf8DBgzQvfDCC3mue1vkaTMFkaqY97HHHtO34ZwiZ5Rpnjbk7ML/6LPhtYG2QYMG6SwBUceI0Mwvd5ol583Sc1HYtvLD2mvHHFqeNnPgHmS6fh8fnzz5/rQIzj/++ENnDZ07dzabpw1gG4afIUIaUZ/5RY+a27aWf8yUpKQktYxhlC/Wa3ocsewXX3xR6H4gil2LWNamKlWq6H7++Wez0aNYL84vjqXhMg8++GCe6N6///5b//mePXv07YiY1toR5WnKG2+8oe5LhuvHtT9hwgT1v2HUcX7HSYtON83T1qFDB+ZpI+4PnrDwpGhYwgqmcuReMoww1MBTMIZPkCsMDsp4osGTDYahsIxWMSG/2qEFrRt9Mcx8XdzltQgs5IPSKiLgSQ/LYrgO1kJLQH4o7DPM+loeOvhmwPJm6BRbUN/wJJ7fUAWOPywNv/32m3qFRQLH0XRf4EeE9VtTUzO/ZWCFxJMsnpoRcYendQwjwgkY/l7wa8GTOxyJsf+mfkDWzKvtP6wAsHLhSde0IgKGk5BDCVZbRHXCHw1WGOTGs2ToTSvRZlrEG1YdWBTg/Izt4skbwzXIk4ahLsPr3lxNUUvbYDE1V4Aa1if0C8NFGvDnwfZxvmGtwfA6LG7IMQYMh4y0aGjTCOj8wDHDcYblFn6msGYaWr4sOW+WnovCtpUfll47+R1TzR3D3PcMIE8bljWsiIBrwLQiAu4rODea9cxScC/BdYx9h4uEYZ425ILDfQHnF5YpBG1h+BIBFoaBAQVtG9ZKc99x3BcM77GYRztP+B5r5wn7aollGtcv7mlaFRHcD2A1xVA9tmOuagiOKyI8Yc3EUKZpRQQNfLcR9YzzaWg5RQQpLN84Hhi9MAWuG7DAwjcP69e+r/huwKqJqbDjBLAc7qWwDsMiiGsF31vc62AZtOXohSviAeXq7E4QUhSQQkBLamoo5JDg84033tAn3yX2AUIaN07cgEkO+HGFmDD00YMTNn4A8eOLoSFN4COxLHzz4DtJiLPQRLA9fPOI7aFPG3Fb4OwNnxk8+eHJF0+lqOgAC4u19ToJsQUQaLBIwdKLV6RigdUNFgHkrjO0yMKKgRxXhBBiKRRtxG2Bmfz999/XF7cGiFyCA75hbjVCHAWuPwyva+W0NGARRlSp6TA2IYRYA0UbcVtmzZqlcv8gNB1DUvD5gM9EfhULiG0pyC+pNIs25DqE3xNe4RcFPzL4+JjzByTE2RTkP0ZcD/q0EUIIIYS4AdanviaEEEIIIQ7H5UQbHHYLq1GJrNtIpFnceQghhBBC3AWX8GlDhNWMGTNkzpw5Kr8SssUj9xNKABn6JyEzP4rnwvEc+bQaNmyoz4ZvzTwFgRIvyFGGfEuWZNUnhBBCCCkOyL6GlEGo7lJQ/V9vV8m3hY4ihxHKc0A0IdkjnMxRIkfbITg+o/A2RB0EFd4jYSAEGnbSknkKA9s2TdRICCGEEGJvkJS6SpUq7heI8OGHH6o6b1rCP2Q7Rl4jFKNu0aKFakNGb2SsRkZ0JK+0ZB5L8iwh4zQOHCLBCCGEEELsCSpJwGCE+rAF1ZZ1CUubOVB6w7DkC5KmotSHJsYAQumRSBUiDYLMknkKQxsShWCjaCOEEEKIoyjMLcslRdvmzZvliy++UD5tGlevXjWqBamB2n74zNJ5TEEeJUyGapcQQgghxNVwuejR/fv3qyLg9913n/Jp04A/mrmoUhTA1ZJWWjKPKQiAgClSm+jPRgghhBBXxKVEGzKII3M4Agc++eQTo88gphBcAAGmkZmZqSxomtOeJfOYAr85+LFpE3zZCCGEEEJcDZcZHj148KB07dpVWdk+++yzPOO6t912m0rH8ccff6jCy2DNmjVKoOEzS+cxBSlGMFlLVlZWofnkCLEEFLtniSNCCCFuIdqOHTumBBuiPJ9//nk5fvy4/rNatWqpYU8UYR4zZowqsowfOEyPPvqo3HPPPXLzzTereS2Zp7gg2BYJgBHhQYitQMRyZGQkcwMSQghx7ZQfSH776quv5pvDDYluAQIG4IO2fPlyJZ5gTXvllVfE399fP78l8xQEAhHg24ahUnPRoxcvXlSCrWLFihIYGMgfWVIscI2mpKSo6h0QbpUqVeIRJYQQF0GXlSUpO3ZK5pUr4h0eLoEtW4hHPj7yxaEw7eFSos2VKOjAYUj0yJEjSrAZVmogxBZVQSDc6taty6FSQghxARJXr5bo12dI5qVL+jbvyEiJeOlFKdujh1NEm0sFIrg6mg8bLGyE2BLtmqKfJCGEuIZgO//UeCPBBjKjo1U7PncGFG1FgDVJia3hNUUIIa4zJBr9+gz4r5j5MKcNn2M+R0PRRooMfLFQ2/XEiRN2XYYQQghxFMqHzcTCZoROpz7HfKUyerQ0kpWtk20nY+VyUqpUDPaX1jXLiZdnweUrisLff/8tH330UYHzTJ48WRo3bmz1upFK5ccff5QXXnhBRfnaa5miAnGogWAW+Ivdf//9Eh4ebna+mTNnSo0aNYw+e/LJJ+XChQvy9NNPy6233qoPdvnmm29k586dKg9g8+bNZeTIkapcmrltGzJ69Gjp27evTfeTEEKI7UDQgS3nsyUUbU5g1f6LMnXFAbmYkKpvqxTiL5P7NpBejWwbPYg0KEh5YijQEEQBMaIRERFRpHUHBQXJkiVLpHbt2nZdpqhAHCL9S+fOnZWT55w5c+S9996Tffv2Ge0z5oNPWf369eW1117Tt2M+LJOcnKwXYdevX5cOHToo8YmKHRCDiHBu0qSJbNmyRQWpmG7bkHr16tl9vwkhhBQxmn/LFon95muL5kc0qaOhaHOCYHtkwS4xHSm/lJCq2mePaG5T4Va9enU1acDqhuoQmgiJjY2VBx98UIkViKlDhw6p9xAb9957r74EWLVq1WTo0KHSsmVL/brgNP/dd98pS1NYWJh+XdOnT1cpV1DhIioqSsaPH6/qvxZ1GW1Y9YMPPpD//vtPWehQ5uyZZ56RadOmFZiDD/3V9vXuu+9WYm3BggXy7LPPGs03ePBgmT9/vkydOlXlBQRffvmlap83b55+vl9++UX+/fdflatPs9iNHTtWLWeaVsZw24QQQlyT7PR0SVz5i8TOny9phw8XvoCHh3hHRKj0H46GPm22UObpmRZNSakZMnn5f3kEm1pP7uuU5QfUfJaszxbZWiCGYBVC+TD8369fP2WdAxAssNINGDBAvL29lZBbtWpVnqFOCC/DdXXv3l3lskMVig0bNqjlkC6lqMuAPn36KPGEIUqIyPbt26vlUKLMUhBOXa5cOTl//nyez3rkhm+vzo0IwhAoxB3EoSGw2CFoQBN2GhCgAQEBFveFEEKIc8mMi5Ors2fLsa7d5OJLLynB5hEYKGEjRkjEpFeVOFOTIbnvkfbDHvnaCoOWtmJyPSNLGkz63SYnAxLsUmKqNJ5iWSjxgf/1lEBf25zCJ554QtVhNQRiTQPiDUOIb731lvTq1avAdaGqBdYHevbsKZUrV1bWqRYtWhRpmd9//10JOVTKgMUPwAqH4Udr2LFjh7KQmesHRBj8zb766iu1f0uXLlVWubZt2xrNB/GI4eXWrVvLiBEj1OcYLtUSQBvy8ccfy8qVK43aCrMMEkIIsS9pJ05I7PyvJeGnn0SXlqbaYDkrN3KEhA4eLF4hITltFSrkzdMWEWGXPG2WQtFGFLC0mXLu3DmZO3euEkvXrl1T7zEVRqdOnfT/I8M/LFCwbhUk2gpa5p9//lElzjTBBlCj1hLRpgknWMg2b96shl2HDRtmdt5x48YpvzZY7yDeTK1sAD5r+/fvV+uFVQ4iFhZBlE/78MMPjerYtmnTJo9Pm2kQBCGEEAeNim3dKrFz58m1dev07f4NG0q5sWOlbM8e4uHjY7QMhFkwRqEcUBHBUijaikmAj5eyeFkCokXHzN1e6HzzxrZS0aSWbNtWmGZgRkoO+J3dcccdqi4sSixBPMExvzAMhYtmxcrOzi7yMhhKxfYNMX2fH5pwio6OVuITr/nlRMOwMIZf4d+3ceNGWbRokdn5ILymTJmiJgyjwkcPUanwHXz55Zf189GnjRBCXMBf7ZdfJXbevBv+ah4eUqZbVyk/ZowEtGhRYJ5MCLSgNq3FVaBoKyY42ZYOUXasE66iRBF0YM4bDZdNZIi/ms8e6T+sYdmyZcrhH3VhNSB6nAEsbL/++qtR26lTpyxa1lA4YdgTqU0+//xzeeihh8zOD+va8OHDZeDAgcqqhpQehYlNDKtCzGL4lRBCiGv4q8V//73ELlwoWVdyfJ89AgIkdOBAKTdqpPgaBOi5EwxEcCAQYkjrAUwlmfYenztbsAFfX19l4YIlCcAXbPbs2U7pC3zrTp8+rYSkZuZGTjVrQZoRDI9OmjRJDfeaA2INUbQY9jQHLHBr1641ajt79qwaMm3UqJHVfSKEEGI70k6clItTpsixLl3lyqz3lWCDH1r4s89Inb/XSuSrr7itYAO0tDkYpPNAWg/TPG2RdsrTVlSQLPaTTz5RPl5ISotEsg0bNlTBAY4GYgspQRAM0a5dO+VzpqUxQSSpNUycOFFZ2iD6MLxpCtJ2FJSmA1GiTz31lLKuIaAAljjkZ4NPoGkgh7lABAw3m/OVI4QQYmN/tQYNcvzVevXM46/mrlC0OQEIs9sbRDqkIoIpyCdmWPAekZCwLCE3mqnP2J49e5QfG6xSCAQAEG/5JcrNb11InYFoy6Iuo0WXIs8a8shh2NbHx0cNmWrJbM2B9bZq1SpP2g+kLTFMFYL5TKNENSAKDT+HcP3zzz9VlQTkjMPnN910k1GQhLZOc2BeQgghxUeXni4Jv8Jfbb6kHTp0w1+tK/zVRktAy5Ylrq6zh84Wyb5KEIgyxA97QkJCHuf81NRUOXnypHJYN02kSuwLhBICCpAvDgEKDz/8sPz1119y9OjREvGl5LVFCCHW+KstlriFC/WlpJS/2oABOf5qJuUI3V17GEJLG3ELEM2KCE0M1x47dkwftVkSBBshhBDL/NViv54vCT/9LLrUHPci74oVVTLcsCGDxcvCrALuDEUbcQtQ6qp///7Kpw5+ZXD6ZwUCQggpDf5q21TKjmt//61v92tQX8qr/Go9xcPXV0oLFG3EbYD/mlZuihBCSMn2V0v87TeJgb/awYM3/NW6dJFyY0ZLYKtWpXKkhaKNEEIIIa7rr+bvL6EDB0jYyJHil1sbu7RC0UYIIYQQp5J2Ev5qX0vCsp+M/dWGD5ewoUNKhb+aJVC0EUIIIcQu6LKy8q3dqfzVtm3P8VczSFruVx/+amOkbK9epcpfzRIo2gghhBBicxJXr5bo12dI5qVLN0RHZKREPD9BdBkZEjN/vqQdyPVXE8n1Vxsjga1Lp7+aJVC0EUIIIcTmgu38U+NhTjNqh4A7/8yz+vfwVwsZ0F/KjRpV6v3VLIGijRBCCCE2HRKFhc1UsBnh6SkVnnhCwu4ZKt5hYTz6FkLRVsJBxYB1BrXYzHHnnXdK5cqVi7T+w4cPq9xpQ4cOzXceFHv/448/1P8weSN1R7NmzaRKlSpm5ytXrpwq3G5ISkqKLFq0SP0/duxYo5qjO3bsUAl3sRxKX6EEl7ltmzJkyJACM08TQgixHuXDZjAkapbsbAls3pyCzUoo2lzQOdOWXLlyRRU01/jll19U/U+UhNK49dZbiyzaNmzYINOmTStQtO3evVseeOABVYTe19dXCan169fLG2+8IU8//XSe+VBX9NSpU0Z9+v7771WCXTiujhgxQom2uLg4JTgxb6dOnVSNVNQDxXwvvfSS2W0b0q9fvyLtMyGEkPzRUnXYaj5yA4o2V3LOfOlFKWvj5LHt27dXkwbEGixcX375pdF8Bw8elF27dqkC7u3atVM10Aw5f/68bN26VdX+xPoqVKigxNfGjRuVWNLWB0tXkyZNzPblgw8+0FvBINhQBB4CLDw83Gg+CLB58+bphRfA+rt166ZqkGpMnTpVFX6HlQ1CVLPIoZB8QdsmhBBiP2CIsOV85AaeBv8TBzpnmpqOM6OjVTs+dyQovn7fffdJly5dZPny5fL2229L3bp1lRjTQI3Pm2++WebOnSvz589Xou23335Tlq7jx4+rOqCw5mGCuLMECLDMzEw5dOhQns/QH2wLVjVNUGIIdtiwYUbzYdnGjRvrBRsIDAyUQYMGFeOIEEIIKQ4YOYIhAhUMzOLhoT7HfMQ6aGkrJhAWuuvXLXfOnDbdvHMm2jxEoqe/LkHt2lk0VOoREFDssOhPPvlENm3aJEeOHNH7d7377rsyZswYZcECM2fOlEmTJsmECRPUe1jW9u3bJ7fccouMHj1azp49m8dyVxgQYsDcsGyfPn1k/Pjxagj1tttuU+uGEDP1P4NFDtY2CE0MddapUyff4/HNN98Y1Sr19/dXVj5CCCG2Bb9fGDlS0aN5Psy5R+Nze7gElXQo2ooJBNvh5jZ6WtDlWNyOtGpt0ez1du0Uj8DAYm1ywYIFSuwsXbo0R4DqdJKenq4saLCaRUVFKbGEodOEhAQ1bFqmTBk1hGotmnDCsOrHH38sjz32mNSuXTvPfPBpGzVqlMyZM0dZ9bDcDz/8IJcvXzaaDyISlkKsC0OtKCTfs2dPmTx5srIMGrJ9+3Yjn7bg4GCKNkIIsRPK1ef9WXJ+/NMq6EDDOyLCLq5ApQWKtlLOmTNnxNPT02g4VBuihCACmsCKjIyU5s2by1133SWPP/640bCkJUA4wScOwQIYxnzkkUfynXfcuHHSsmVLZWmDUIRVDcLNVNy98sorakIwAvbhrbfeUkLvwIEDqr8a9GkjhBDHEgR/6tzfkchp08S3alW7Bd2VFijaigmGKGHxsoSUHTvk7IMPFTpf1c8/k8CWLS3adnGBcz6sZu+8806+89SvX1/++usvZWlbu3atEkkQYKYiqjAMhRNEGYZBNQFnbptNmzaVJ598Ul5++eVC112jRg01dezYUb2uWbNGhg8fblX/CCGE2I70kyfVq1d4BQkbdDcPrQ1gIEIxgQ+VZ2CgRVPQrbda5JyJ+SxZny3KfPTu3VsNkSKowBAMj5r+D4tX//79VQoNpNLQhhmvW+jTZwj85hITE5W/XH7Aj+6ee+5R/nUF+cUZou0HolsJIYQ4j/RTp9SrX42aPA02gpY2ZzlnQnAZBiQ4yTnz1Vdflb///lslu0XSWli9kNoD4gfWKoBgAwQMtGnTRgk0DJciFxrAEGZMTIy88MILctNNNxWY8sMQWNxgQYP/2cMPP6wS7prSq1cvNeUHAhSQOLdHjx7Kunbx4kUVddq1a1c1FRSIUNykwoQQQiyztPnWqMFDZSNoaXMwcL6Men+WcsY0BO/Rbm/nTFjWkN5D35+yZWXz5s0yY8YMJdQgwBBVaZgPDVGcqB5w4cIFSU5Olq+//lpFbQIEEmDIFAEMEHvmUn5AUMFHzs/Pz6gdfnJIeouKBobzwe/NHDVr1jT6HEO6ixcvVsILlRmQcPfzzz9XQg7+bobr3LNnjz4tiTZhuJcQQoh9SNNEW01a2myFh05LhmUF+GFGAlP8mCMSMD4+XkXuVa9eXTmMu7MFA0N2GAbED7ppionU1FQ5efKkEg9IGeEOFRGIe2DLa4sQQlyBEwMGStrBg1Jl9icSbGAsINZpjyIPj8KSgkz2K1askKysLOXPFBERoTYAKwzSQmAYCpYQOJm/+OKLariM5AUCLagNjw0hhJCShy47+4ZPGy1tjh8exZAZogwvXbok77//vkqpAEWIguQ7d+5Ur1CKiAZ87733lH9R27ZtmQuLEEIIKWUg56hKPO/tLT5RUc7uTonBYktbbGysSvPQokX+iWQRzdigQQM1IY8XfJXg6E4IIYSQUhiEULWqeOT6GBMHijZzRbgLA5GFqFFJCCGEkNJDWu7QKIMQXCx6FHEMyOO1bds22/SIEEIIIW5N+slc0cZ0H64j2hYuXCjVqlVT+bmQw0ujb9++yuetpFKEgFtCeE0RQkrf8GhN5mhzCdGGAuNIutq5c2eV0NSQgQMHqpJFJQ0t91dKSoqzu0JKGNo1pV1jhBBSEkQbI0ddIE+b5q+GiNLx48fnrMjDQ2+BQiRpt27dVDHykpYrBVGxyEuHDP6BNiolRUov+M5AsF2+fFlViahUqZKzu0QIIcUiOy1NDt/STFX9qbNpo3iXL88j6ow8bYbs379fWdrMUaVKFZUapCQSidqhIupHlhBbAcGmXVuEEOLOpJ8+rQSbZ3CweJUr5+zulCiKLNp8fX3l2rVrqhICMLQ4wcIGxVgSwX7CGgJLW0ZGhrO7Q0oAGBJFCS5CCClRQQg1a3I0ylVEW6tWrZTf2ttvv51HtH300UfSvn17KcngR5Y/tIQQQogxNyohMAjBZUQbSlT16NFDjh07JoMGDVJtP//8s3z//feqiDfqkhJCCCGktEaOslC8y4i27t27y/z58+XJJ5+Un376SbX1799fDYt+/fXXJd7SRgghhJACRBtztLmOaAMjR45U6T02btyoAg8qVKggnTp1UoXkCSGEEFL6IuJZDcFFRRsICgqSnj172qY3hBBCCHFbsuLjJTshAY7u4lu9urO7U+IodhkrQgghhBDDoVGfSpXE09+fB8VZljYk07WWHTt2WL0MIYQQQtwT+rO5iGiDvxohhBBCSH4wctRFRNuqVavs2xNCCCGEuDUMQrAv9GkjhBBCiG2rITDdh+tFj169elUl1EWCXXMlnWbOnFmc1RNCCCHETdBlZkr6mTPqf1ZDcDHRtmfPHpXqIzo6Wl/GCvlZgFbeiaKNEEIIKR1kXLggkpEhHv7+4l2pkrO7UyIp8vDoK6+8Iu3atZO4uDj1Pi0tTY4eParaUdYqMTHRlv0khBBCiDsEIVSvLh6e9L6yB0U+qps3b5bXX39dQkND1Xtvb2+56aab5LXXXpNevXrJc889Z9X6MjMz5ccff5R+/fqpyRyw5H311Vfq87vuuks+++wzyc7OtnoeQgghhNiWNJavcl3RBgsbRJom2JKTk/WfDRkyRH744Qer1oc8cAsWLFAlsDZt2mR2nueff14mTpwoffv2lQEDBsirr74q48ePt3oeQgghhNgpCKFmDR5aV/Npg0XLx8dH/V+5cmXZvXu3dOzYUb0/d+6cpKamWrW+NWvWSPny5WXWrFlm04tcvHhR3nvvPVm4cKEMHTpUtZUtW1YJxGeffVaqV69u0TyEEEIIsT3pp3JEm1/Nmjy8dsImg849evSQhx56SJYuXSorV66UESNGSOvWra1aBwRbQfz1119qmLNPnz76tt69e4unp6cSfJbOQwghhBDbw8S6LmxpGz16tP7/qVOnSufOneXuu+9W76tUqaL8ymzJ6dOnlf8cCtRr+Pv7K7GHzyydxxQEUGDSYAAFIYQQYh1Z15Il8/Jl9T9ztLmgaJs3b57+fwyP7t+/X7Zs2aLet2jRwkg42YL09HQJCAjI0x4YGKg+s3QeU2bMmKFEJyGEEEKKNzTqVb68eJUty8NoJ2wWk+vr6yudOnVSk60FGwgLC5PY2Ng87TExMeozS+cx5cUXX5SEhAT9dPbsWZv3nRBCCCkNoo1BCC4q2v7880+VUsMcaLe1D1mzZs1UcMOBAwf0bcgLh+FMfGbpPKb4+fmpYAXDiRBCCCFF8Gdj+SrXFG3Ix4b0GuaYMGGCTJs2TWzJrbfeKnXq1JHp06fr2/B/jRo1lD+dpfMQQgghxD6ijZGjLurTtmvXLrnlllvMfoZ2pACxBgxT/vPPPypdCIYpNZE1d+5cqVmzpiqNtWTJEpU0t1q1aqpsFiJFly1bpk89Ysk8hBBCCLEtaadyLW1M9+Gaog0CCXnRYNkyBe1aHVJLuffee1UtU1MqVqyo/79p06Zy/PhxVfcU68d7JPY1xJJ5CCGEEGIb8FubfionQ4NvDeZosydFVjPIw/bOO+/Ip59+muezd999V1q1amXV+ho3bmyxWGzevHmx5yGEEEJI8UGqD11KCn58xbdKFA+pK4o2+LMhqS6iLWEli4qKkvPnz8uiRYvkt99+k99//922PSWEEEKI6wYhVKkiHr6+zu5OiabIoq179+4qgS7qev7666/6dkRfov3222+3VR8JIYQQ4vLpPjg0am+K5ew1duxYGTRokGzcuFEuX74s4eHhqv4oir4TQgghpOTD8lWOo9ge+hBod9xxh216QwghhBC3Io052lw/TxtSaRim9UD0yIcffihDhw41G5xACCGEkJJH+klWQ3B5Sxvyp6HW6BdffKHef/755/Lkk0+qclGLFy9WpaxGjhxpy74SQgghxIXITk+XjPPn1f9MrOvClrbZs2fLY489pn8/Z84cZWVD7c/PPvtMfU4IIYSQkkvGmTMYehPPoCDxqlDB2d0p8RRZtB0+fFjq1aun/k9KSpIdO3bIuHHj1HuIt0OHDtmul4QQQghxXX+2mjVVFSLioqINhdZhVQN///23Ollt27Y18nkjhBBCSGnwZ2O6D5cWbS1atJBJkyapGqSvv/66EmzI0QYOHjwojRo1smU/CSGEEOKyOdpqOLsrpYIiByJMmTJF1QqFLxuKsa9cuVL/2fz582X06NG26iMhhBBCXDhHm18NijaXFm3t2rWTI0eOyM6dO6VBgwZS08A0is8GDx5sqz4SQgghxAVhYl3H4qFDgjWiJzExUUJCQiQhIUE/3EsIIYQQYzLj4uRou/bq/3q7dopnYCAPkZ21R5F92gghhBBSetH82bwrVaJgcxAUbYQQQggpeuRojeo8eg6Coo0QQgghRQ9CYLoPh0HRRgghhJCip/uowRxtjoKijRBCCCFWk35Kq4bAdB+OgqKNEEIIIVahy8qS9NNn1P+shuAGedoqFFIYFmWuIiMjpXPnzvL0009LlSpViropQgghhLgQGRcuiC49XTx8fcWnUiVnd6fUUGRLW//+/aVZs2aSlpamSlrdcccd6jU1NVW1d+3aVeUa+fjjj6Vp06aqwDwhhBBCSpA/W/Xq4uHl5ezulBqKLNrGjRsnwcHBcubMGfn999/lm2++Ua94HxgYKI888oisXbtWjh8/rqolvPzyy7btOSGEEEKcWwmB5avcQ7SNHz9e3njjDQkLCzNqL1eunLz11ltqSBRERUXJRx99JOvWrSt+bwkhhBDidNI00cZ0H+4h2vbu3Svh4eFmP0M7Ptdo1KiRpKSkFHVThBBCCHHF4VGKNvcQbQgy+Pbbb81+tnDhQvW5xunTp6Vu3bpF3RQhhBBCXLAagh/TfbhH9Oijjz6qhkhhUevdu7dUrFhRLl++LCtWrJC5c+fK9OnTjUTcsGHDbNVnQgghhDiJ7JQUybx0Sf1PnzY3EW0TJkxQVelnzpwpn332mVGqj4kTJ6rPNXr16iXt2rUrfm8JIYQQ4hJDo15hYeIVGsqz4Q6izcPDQ6ZNm6YCDjZv3iwxMTFSvnx5ad++vXo1pFOnTrboKyGEEEKcDP3Z3FC0aUCgIYFubGysihxFGhBCCCGElPTIUZavcqsyVuvXr5e2bdtKSEiI1KhRQ73i/YYNG2zXQ0IIIYS4XBAC/dncyNKGIdHu3btLRESEPPjggypaNDo6WlauXKnaIejatGlj294SQgghxCUS6/ox3YfD8dDpdLqiLNijRw9VpmrRokXi6+urb09PT5d7771XkpKSVIUEdwPBFbAYJiQkqP0jhBBCSA6QDEdatpLs5GSp9ctK8atdm4fGgdqjyJa2LVu2yJ49e4wEG8B7VERA/VFCCCGElBwyr1xRgk08PcWnalVnd6fUUWSfNljU8gs6QDs+J4QQQkjJ82fzqVJFPE2MNsSFRVv9+vWN8rMZ8uWXX6rPCSGEEFIS030wctQZFHl49IknnpD7779fVUQYOHCgPhBh2bJlsnjxYpkzZ45te0oIIYQQ1whCqFGTZ8KdRNu4cePk/Pnz8vrrryuRpuHv7y+vvfaajB492lZ9JIQQQogLiTZa2twwue6rr74qjz32mGzatEkl19UqIiDJLiGEEEJKFmmnckUbLW3uWREBAq1v37626Q0hhBBCXBJderpknDuv/vdljjb3q4hACCGEkNJB+rlzIllZ4hkYKN4Vw53dnVKJxZa2li1bWr3yHTt2WL0MIYQQQlzYn61GDfHw8HB2d0olFou2ChUq2LcnhBBCCHGDIARGjrq8aFu1apV9e0IIIYQQlyVNn6ONos1Z0KeNEEIIIRZXQ8DwKHEOFG2EEEIIKRTmaHM+FG2EEEIIKZCshATJio1V//vR0uY0KNoIIYQQYlHNUe+ICPEMCuLRchIUbYQQQggpkDSDdB/EeVC0EUIIIcSyIISaFG0lTrR99NFHsnTpUsnOzrbH6gkhhBDihOFRP6b7KHmi7YknnpC7775bmjZtKsuXL7fHJgghhBDiIJhYtwSLtiVLlsjcuXNV6avHH3/cHpsghBBCiAPQZWdL+unT6n/6tLlJRQRrGDRokHodM2aMPVZPCCGEEAeRceGi6NLSxMPHR3yionjcnQgDEQghhBBSqD+bT/Vq4uHlxSPljqINQQa7d+/Wv9fpdPLhhx/K0KFD5dNPP7VV/wghhBDiCv5sTPfhvsOj8FnbsmWLfPHFF+r9559/Lk8++aSEhYXJ4sWLJSgoSEaOHGnLvhJCCCHESaKNkaNubGmbPXu2PPbYY/r3c+bMUVa22NhY+eyzz9TnhBBCCCkZw6O+NWo6uyulniKLtsOHD0u9evXU/0lJSbJjxw4ZN26ceg/xdujQoVJ/cAkhhBB3J+1U7vAoc7S57/Con5+fsqpFRUXJ33//LR4eHtK2bVv95/ZIrPvDDz9IamqqUVuTJk3UZEh8fLxs2LBB+dl17NhRDdkSQgghxDqyr1+XzAsX1f+shuDGoq1FixYyadIkNUT6+uuvK8FWtmxZ9dnBgwelUaNGYmsefvhhqV+/vlSvXl3fFhgYaCTaICAHDBigrIBeXl7Krw5543r06GHz/hBCCCElGS0/m1dIiHjTAOK+om3KlCnSs2dP5cvm4+MjK1eu1H82f/58GT16tNgDiMR77rnH7Gfp6ekyYsQIGTVqlLz//vuq7dlnn1XvT548KQEBAXbpEyGEEFKi/dk4NOrePm3t2rWTI0eOKLEG/zZDSxY+g3iyB//9958aJt2+fbtkZGQYfbZu3To5f/68jB8/Xt+G/6Ojo2XNmjV26Q8hhBBSUmH5qhJUESEyMlJ69+6dpx2WLXuxbNky2bdvnwp8wHAs0otoQ7Fox3BpTYMngqpVq0poaKj6rE+fPnnWl5aWpiaNxMREsSe6rCxJ2bFTMq9cEe/wcAls2YLJCgkhhLgkaczRVjJE25UrV5QQQjACHP4rVKigxFN4eLjYi2+//VZuv/129T8CEvr3768iVffv368CIRISEswGHZQrV059Zo4ZM2bI1KlTxREkrl4t0a/PkMxLl/Rt3pGREvHSi1KWPneEEEJcjPST2vBoDWd3hRRleHTbtm3SuXNniYiIkG7dusngwYNlyJAh0rVrV9XWpUsXZQWzB5pgA/7+/vLiiy/KgQMH5MSJE6oNPmvXrl3LsxxSkuTnz4Z1QNBp09mzZ+0m2M4/Nd5IsIHM6GjVjs8JIYQQVwEGGc2njYl13dDStmnTJiXUqlSponzFEMmJoUctzQaiRpcvX67SbPz111/Kt82eoOoCuHr1qtSuXVtq1aqlhBf6ovULIg7WQHyWX+oSTPYeEoWFTXQ6Mx/qRDw81OfB3bpxqJQQQohLkBUTI9lJSeo3yqdaNWd3h1gr2l566SWVQBc1RpFOwxxvv/22PPHEE8qChfQbtgLBBOXLlxdv7xtdRkAChFvDhg3VewhKWOCQ4uOBBx5Qbfjf09PTqSk/lA+biYXNCJ1OfY75gtq0dmTXCCGEkAKDEHyiosTTzsYNYgfRhmHPn376KV/BBvDZ9OnTVQCALYH/3IQJE6Rfv35SqVIllTwXguzjjz+WMmXKqHkg6v73v/8pKyCiSNGXt956S1599VU1dOssEHRgy/kIIYQQhwUhMN2He4o2WLlSUlIKrTCQnJxsZBGzBd27d5cff/xRFi1aJLt27ZIGDRqoAIQ6deoYzQdh17hxYzVMi/F4BC+Yixp1JIgStYTYr+eLh7eXlOnaVTx9fe3eL0IIISQ/0k/lJNZlEILr4KGDsrGQu+66S7KysuSbb75REZnmgP8YqhAg4S6scu4GUn6EhIQo3zitwoMtfNqOdeuugg7M+rWZ4BkSIiG9e0vIwIHi37CBiowlhBBCHMnZRx6Va2vXSsSkV6Xcvffy4LuA9rDKHIahxg4dOqhAhNtuu02VioLVDboPzv9IsosEtxiu3Lhxoy32o0Tg4eWl0nogShQOnUbCLVeQRbz8shoeTfj5Z+XfFrdokZr86tRR4i3krr7iXb6883aCEEJIqfRpY+Som1rawJkzZ5TfGJLcwqpmCHzKBg4cqHzIbO3T5s6WNo11X34v3p/MknIp8fq2mMBQyXp0vNx2/1C9VS75ny2SsHSpJP35p+jS03Nm9PaWMp06SejAAerVg8OnhBBC7IQuI0MONWsukpkpN/29VnwiI3msXUB7WC3aDIGzf0xMjF6wRUVFibtjL9G2av9FeWTBLvHQZUvDqyekXFqSxPoFy4EKtSTbw1Nmj2guvRpVMlomKzFREn/9VeKXLZPUPXv17V5hYcryFjJggPjffLPN+kgIIYRoQQgn7rhTPAICpN6unXTTKQmirSRiD9GWla2TDm/+JRcTUs1+jgHSyBB/2Tixq3h5mvdfSzt2TIm3hOXLJevKVX27X4P6EjpgoJTt01u8CwkQIYQQQiwh6a+1cu7RR8Wvfn2ptWwpD5o7+rQ5u4yVu7LtZGy+gg1ANeNzzNeutnm/Nb+bbpKICROk4tNPy7WNGyVh2U+S9NdfknbgoEQfmC7Rb70lwV26SMiA/lKmY0fxsHH0LiGEkNLoz8byVa6Ed1HKWD3//POyfv16JdYMQZQjAhSQYLdly5a27Kdbczkpf8FmyKfrjklcSrq0rB4mFcv6m50HYiy4c2c1ZcbFSeLKXyQBw6cHDkjS6tVq8gqvICF33SWhAwYosUcIIYRYg1a+yrdGTR44F8Kq4VHDMlZI/5FfGSv4ujmijJW7DI/+czxGhn2xxaplqpULlJY1wqRVjXLSqkaY1A4vU6BPQerhwyp4IWH5CsmKi9O3+zdpooIXyt5xh3iFhBRrPwghhJQOTo8YKSk7dkjlt9+SkL59nd2dEk+iPXzaYEVDyaiCylghjxvKWKGQuy3LWJUEn7ZLCalqKNQcoYE+0qdJJdl1Ol4OXkrMk84tLNBHWlTPEXAQc42iQsTPO+85QLTptfXrJX7ZT3Jt3ToV+QMQbRrcvbsKXghq3441TgkhhOTLkQ4dJevqVamxZIkENG7EI+WOog11Ps+dO1doRYS4uDiV8gPF2t0Ne0ePAsMDrtnODKNHE1MzZPeZeNlxKla2n4qVf8/GS2pGttH6fL095ZYqoXprXPNqYRIS6GM0T2ZMjCSsWCEJS5dJ2pEj+nbvyEgJ6ddPQvr3Y/4dQgghRmQlJcmRVjl1sOvu2C5euaUiiZuJNqwQFrTCUntA2CEoAUOm7oY987RBuE1dccAoKKFSiL9M7tsgT7oPQ9Izs+W/Cwmy83ScEnE7TsVJTHJu/rZcMHJaLyJYWlTPEXEQc1GhAWpIFac49b8DyvctceVKyUpI0C8X0Ly5Cl5Qw6f8YhJCSKnn+r59cmrwEFWCsc6G9aX+eLitaGMZK9sMlSJKFMEJFYP9pXXNcvmm+cgPnLKTV5OVeFMi7nScem8KBGHLXJ+4ltXLSb3IYPHIzJBrf62V+GVLJXnDRpHsHAueh7+/lO3ZQw2fBrZuLR6enjbYW0IIIe4GUktdeH6iBLZqJdW/+drZ3SkVJNpDtB06dEiVsULReEvKWNWtW1fcDXta2uzJlaQ0ZYlTQ6qn4+S/8wmSmW18aoP9vKUZLHHV4RdXThr5pUv6ql8kfukyST9xQj+fT+XKEtK/v4QMHCC+Vao4YW8IIYQ4i8vvvy8xsz+V0CFDpNL/pvJEuHNyXZaxcg9S0jOVL5xmjYOP3LW0nKAEDW9PDxXQ0Kp6qLRLj5aaO9ZK+h+/S3ZSkn4ePGmh9imscJ6BgU7YE0IIIY7k3NNPS9Jvq6TixIlSfuwYHnwHwDJWdj5w7gaGZQ9dStSLOEzRiWl55qsX6iP9k49Ji4ObJPi/3fri9hBswb16qfQhAS1asKQJIYSUUE70HyBphw5JlU9nq5ygxP6wjJWdD5y7AwPrubjrsuN0TmADpsPRNyxsIDwlTvpG/yu3n9khoXHR+nafatUkdEB/FYGKoVRCCCElA112thxu3kJ0qalSe9Vv4luDFRFKhGgrqWWsSotoM0d8SrrsOgNLXI5v3J6zCZKela2sbQ1jT8ntp7dLpwt7JCAzx0Kn8/AQ/9ZtpPyguyX49u7i6W++igMhhBD3IOPCBTnWtZuIj4/cvHsXSyK6e+1RlrEquYQG+krXmyPUBNIys2T/+YRcERch86rWldmJ/aTDhX1y+5nt0vTqcUnbukUubN0i6X4Bcr1DV6kybIhE3dqKw6eEEOLO5auqVqVgc0GsEm2GZazGjx+fbxmrjh07um0ZK3IDVFxAFQZMclttyc7WyfEr12T7qZay7VQ/+X7/UWmwd4MaPo24Hie+a36RpDW/yNqQCDnTqqsE9ukjTW+pI3UqlhFPK9OaEEIIcTxpuYXifWuy5qgrwjJWJpTm4dGiEJ2YKjtOxMjpvzZK2XWrpOGxneKflaE+yxIP2RVRTzbWbiu69h2l2U0VVeLfxlEh4u9jvgwaIYQQ53Fp2nSJW7BAyt03TiImTOCpcBAsY2XnA0fyOX4x8XLou2WS8csKKXfioL49ySdA1lZpJn9Uby1nylWVJlVDpQVKcClLXpiEBfnykBJCiJM5c9/9krxpk1Sa9pqEDhrk7O6UGhLt4dPm7e2tEusWVns0OTlZzUtKH2XLh0rrx8aKPDZW+UbELl0msct+kuArl+Wuk5vVdLJspPxRrZV8X6W5fOYfrJbDEKpWvQHWuCphOSW47FlZghBCSD4+bRwedUlYxsoEWtpsjy4rS5L/2SIJS5dK0p9/ii49p25qtqeX7K/aUH6KaC7bI2+WTM8bQj+irJ8qvdUyV8TdHBks3l6exarhSgghJH+yU1PlcLPmKmNAnU0bxbt8eR4uB8EyVnY+cKRoZCUmSuKvv0r8smWSumfvjfayIXLylo7ya1QL+TO9rGRkGWeiCfL1kubVwyQkwEdW7r2YZ72ajW32iOYUboQQUoSH6/hlP8mlV14Rj8BAqbttq3hyxMxhsIyVnQ8cKT5px44p8YbixFlXrurbfevXl4TbesqO2i1l89Us2XkqTpJMSnB56rKl4dUTUi4tSWL9guW/CrVE5+EpkSH+snFiVw6VEkKIhSSuXi3Rr8+QzEuX9G3ekZES8dKLUrZHDx5HB8AyVnY+cMR26DIz5drGjZKwdJkkrV0rkpETfYrkjsFdukhw//5yoW5TWbL7oszZdEraX9gnD+/9ScJTE/TruOIfIp826S+bKzeWx7vUluFtq0ulkACeJkIIKUSwnX9qvL5koZ5cn+Ko92dRuDkAlrGy84Ej9iEzLk4SV/4iCRg+PXBA3+4VXkGutO0qi0+mykP7l6s2w7CD7Nz301qPVsINVCsXqAIUMLWpWU69Lyi4gbjvsE7Kjp2SeeWKeIeHS2DLFuLh5eX0dRHi6uB6P9atu5GFzQgPD/GOiJCb1vzJ70FJFm0fffSRVK5cWfr37y+enjnO4+4CRZvrkHr4sApeSFi+QrLi4vTtuGDNSS8It6sBoTJ12HQ5FXddsk2ubAQ3tK5ZXi/ibgpn0l93x5bDOhwiIsVB/ZRmZyshJJmZOa9ZWepVl4n/c9owsqDmM21T82bnbctEu8kyWpvhMoZt+v9z+2P4OV6zc5bJuHxFru/YUei+VZs/X4LatOYFUlJFm2bNQC3S6dOny1133SXuAkWb64Fo02vr18vVr+ZI6u7dhc4f8dpr4tm+g+xJ1Mm2MwkqNYi+jqoBYYE+KjI1R8SVl/qVbkSoktI1rMMhoqKTIwDyipN8BUsB4kSytWVMxUdOe542bRlzgkWJE5Nl1KvJMrltRiInV3QZLaPtj9HnN9rRVlKpPHOmhPTp7exulGicKtp++OEHuXbtmqxbt07WrFkjZ86cEXeBos11SVj5i1x47jnLF/DwEK+QEPEqX148w8Ikyb+MXPIMlJOZvnLoupdc8Q6UeL8yEu8XLAl+QSJlgqV5jRuWuMZVQlQpL3vBoTg7DutgSD00VCq9+YZ4+geIp5+vePj5iYefv8H/OZN4esrx7rdbNUSkbpsmP9o3BEnuj7qyjNxoyyNO8rWMmBEslixjKD40wWLallXYMtlmBYt5QZMrVGz/E1Ly8PYWD4w64dXLK+c6yv1fvDzFwyv3f298prV75W3D/55m2vA/7lWeXnnbvEyW0dq8vCX97BmJ+/qbQrtPS5v9oU+bnQ8ccTzJW7fJmdGjC53PMzhYsq9ds/rHJMvDUxJ8g5SQS/Aro0Seb4XyUi4qQqJqVJYadapIYMVw8S5XLkcIBgUV2UeOQ3FFB4Ip/ocf5dKrr4pDgcjT6fRCjRSCp6eRODEWIibixSJBk/u5Eh3aMoZCxGQZA3GiFzSmbWoZT5M2zJe7LsNl9H0yWMawTW3D80abfh2eLutLq3/4iY42f7+kT5vDoGiz84Ejrn2DAVnx8ZIZEyNZsXGSFRsjmTGxkhUXm/OqvY+NlczYWMlOTLS+Q76+SsBpIk694n15vJbPfb3xuae/v0sOxbmDxS/z6lVJ3rxZlde5tmmzZF29kSKmIHyiosTD3190aWlqyk5P1/9vNwuRj48SFEY/3GYsKEZCI9eCYn4ZCATv/JcxbctvGU1UGAoWvYjRRFVOm7EQyVlH3s8LEGDKukNXA3dAfz8Cht8JRo86FIo2Ox84UrJuMPCby4yL14u5zJircvlstFw4dUHizkdL6pWr4p+cKCFp1yQ07ZoEZOVUdbAGJKz0CguTzMuXb6Q1MQMEXrX588Q7NFQ8y5YVT1/71WV1tMXPUoEIcXV9164ckbZxk6QdvFHHVoFhTQivIg7rqOHNjAy1neTN/8j5J5+0yK8nsEVzEwuLV17BQoibQcu/83G6aKtQoYJctfBp2JWgaHN9nHWDOReXooIaMO0+ekniLkQrAQchF5KWLKHp16S6Z5rU9E6TiOxUKZuaJJ4J8ZIVEyO6AkRaYcDvyrNssHgFlxWv4GA1/OtVFq9lb7wGlzF4j8/LimeZnPlgaTI3PONoi19B5y349tsl/eRJSd4IS9pGSdm2XXTXrxst79+ggQTdeqsEdeggAU0ay/E77rTJsA6HiAhxD4t7Scbpog0/EnZatV2haHMPXOEGg2L120/GybaTMbL1ZKwcupSUZ56o0ABpUyNM2lbylxZlRYI3rZGYjz4udN1qSC/1Rm3VYuHjo8SeEnxlcwSeR5lgSd6wIY8wMgRDu1VmzxavoMAc4efrq4Z41f8YArTCTydfgZiLZ1ioZMfFG28/vIKUaZ8j0oLat8tTB9GWVlcOERFCSpxo+/fffy3uQLNmzSjaSKkiPiVddpyKk22nYpWI238+QbJMksV1uHZaXv7zw0LXVWXePCnTqqVkJydLVmKSZCclSlYSXpNuvFevSbntOe+zkhIl26AduZ3sgoeHEm+ecM6HiPPzFU8/f6M2T/+caE3x85Wk31aJLiWl4HX6+EhQq5YSdGsHCepwq/jVrVuoMGSeNkJIScAuos3aCBha2khpJjktU3adgSUuR8T9ezZeMjMyZd7v06V8aoJ4FpAg+K9XP5WmNcpLWKCvhAb66F/L+vuIp6eHxd+/7OSUGwLvGgRfohJ01zb/I4k//VToOjxDQtT3PhvO+7D82dF6XvWLz6VMx45WL8eKCISQ0iLavK1ZaVBQkDzxxBPqtTBedXQ4PiEuRpCft3SsE64mkJqRJZ+sPSafXuwvr2ybrwSap5lSXJ817iebt50TwWQC9FpIQI6ICzEQc+o1wEdCg3xV0mD1OebD+/LhEhAZafTQ5R0RaZFoq/LBB3pHfohA+OZBvKlIzFwhl52KaEy8oh3RmcZt1//dI0mrVhW6rayEIkTw4ph5edksW7st10UIIbbGKtHWu3dvueWWW2To0KGFzkvRRogx/j5e0q52BfmgcmNVI9W06D0sbEqwVW4sbWuVE29PT4m/ni5xyRlq6DU5PUuV5opLyVCTNfh6eypRpxd5fl4yPCBEyl3P3+IXFxgqdZo317dB9MGvDalOrM2vZ4log28iIYQQG4m2IUOGyKJFiywSbYSQvKDaQqUQf/lHGsuWSg2l4dUTUi4tSWL9guW/CrVE5+GpPl94f1vxMhkGTcvMkoTrEHAZEpecroRbAkSdEnHpEg9xl/seIk97zcjSSXpmtlxOSlOTRmLjgi1+nzTqJ9M/2Sz1IoIlMiRA9SsyxF//WiHIz6KhWgSJIEo049IlszVjMeDqExmp5iOEEGIj0XbnnXdKNELsLeCff/6xZtWElAogxCb3bSCPLNilBNq+8Jv0n2mCBp+bCjaAkloVgzHlJOm1BAxppqRn5Yg6TdylZMj6I1dkiUihFj+5mCQHL+aNigXenh4SUfaGiMt5DZDIsjfeVwz2E28vL7kw4mEJnzklX4GIz+swvQAhhBSI3VJ+uCtM+UEcwar9F2XqigNyMeFGWg+IHAi2Xo0q2X37/xyPkWFfbFH/e+qy81j8sj1ypNXjXWpLsL+P6uelhFS5mIjX63IlKU0N1RYGtGeFMr7K6tfq7N48AvFyQKh83rifnKjfSjZO7GpWrBJCSEkn0dl52twVijbiKJAOBJGlyPcG6xmGTh0lWrDtDm/+pYSYuRsAegFrWX5CKjMrW65cS7sh5tTrdaP30Ympkmmi7AoSiN8+0Fba1TbOxUYIIaWBRHtEjxZGRkaGeHl5iSdLuRBSKBBDzhIphsO0kGSG0qqwYVrg7QXfO/i5BeS7jexsnVxNTpPvt5+Vd1YfyWkzGRI2BOKVEEJI/lhdKC89PV1WrFhh1JaamiojR46UMmXKKIX49NNPS1ZWlrWrJoQ4EAzDzh7RXFnUDMF7tBd3mBZBCrAgtqxezqL5rfHVI4SQ0ojVlrYvvvhC/vvvP+nbt6++bebMmbJgwQIZPHiw1K5dW2bPni1169aVRx55xNb9JYTYEAiz2xtE2nWYVouYzW8oFuBzzEcIIcSGom3OnDnyySef5Gm7/fbbZfHixfoSVu+++y5FGyFugL2HaQsaitXoVCecQQiEEGIr0fbRRx+p171798qGDRtk+/bt+qHRkydPStOmTfXzpKSkyL59+/TvweOPP27ppgghJXQo1jRiNtjPW5LSMmXxzrNyW71wubOx/SNnCSGkxIs2b29v/as2AQyVgpYtW+rb/P39JTMzMyeDeu5ECCndmBuKbVUjTCYv/08Wbj0j47//V8oH+UqbWowgJYQQm6T8aNOmjTz55JMyfPhwfZWEdevWySVkO88VZ4cPH5aBAwfqBZ07wZQfhDgWpB95ZMFOWX0gWoL9veWHh9tLvchgngZCSKkh0cKUH1ZHjz766KMyduxYJco6deokS5Yskccee8zImrZx40b1GSGEWOLz9sGwZtKyepgkpWbK6Dnb5EL8dR44QggprmgbPXq0CjKAHxusaxMmTJAXXnjBaJ5vv/1WiTtCCLEEfx8v+XJ0S7mpYhm5lJiqhBvqphJCCLkBKyKYwOFRQpzH+fjrcvcnm5Vwg7/bN/e1UYKOEEJKMnYbHiWEEHsRFRog88e1Vr5t20/FyVPf7VY+b4QQQijaCCEuBoIQvhjVUny9POX3/6Jl8vL9whLJhBBihWhD4MH+/futOmaY/+677+ZxJoRYRdta5WXWPbcI4psWbDkjH689xiNICCn1WCzaEB2KBLrdu3eX+fPny9mzZ83Od+bMGZk3b55069ZNzU8IIUUBiXan9G2o/p+5+ogs3m7+nkMIIaUFi5Pr/vjjj7J+/Xp5/fXXVcoPDFeEh4dLpUqVJDg4WDnRXbx4Ua5evSqenp7Ss2dPlb+tQ4cO9t0DQkiJZXT7GiooYfbfx+XFZfukQrCvdL05wtndIoQQ94keRbqPlStXqnxssKzFx8dLaGioVK9eXTp27Ci9e/eWGjVqiDvC6FFCXAvcop5dskeW7jovAT5esuiBNtKsWpizu0UIIQ7XHkz5UcQDRwhxHBlZ2XL//B2y7sgVCQv0kR8faS+1wsvwFBBCSgRM+UEIKTH4eHnKJ8ObS5MqIRKXkiGj5mxT9UsJIaQ0wTxthBC3IMjPW+aMaSXVywfKubjrMnbudklKzXB2twghxGFQtBFC3IYKZfzk63GtpUIZX/nvQqI8vGCnXE/Pkn+Ox8jP/55Xr8VJxotlbbUuQgixNfRpM4E+bYS4PnvPxcs9n2+RlPQs8ffxlNSMbP1nlUL8ZXLfBtKrUSWr1rlq/0WZuuKAXExILfa6CCHEGujTRggpsTSpEir3daip/jcUbOBSQqo8smCXEmGWgnmxjKFgK+q6CCHEXnB4lBDidmDY8oed58x+psudJv38n5yPuy4x19IkMTVDUjOyzA53og0WNnMDoVobPudQKSHEbZLramWpGjVqZNG8P/zwgwwaNKio/SKEkHzZdjI2j1XMlMtJaXLrm3/laUdpLESjorapj5eHIFNl/PX8Axog3LCtLSdi5NabKvCsEELcQ7Q1btw4T+HmChUqqCoIpgwePJhFngkhdsHSdB+eHiKmxjXcwtIzs9VkDaO+2irVKwRJjfK5U4VAqa7+D5So0ADx9uLABSHEhUSbOWJiYmzTE0IIsZCKwf4Wzbfw/rbSpmY5ycjOlowsnWRkZhv/n5UtO0/HyQtL9xW6riydyIkryWoyxdvTQ6qEBehFHF5rVghS6UmqhAWKrzcFHSHEBUQbIYQ4mtY1y6nITgQKmPNF8xCRyBB/NZ+np4f4eXqJH+52fnnnRWWF99ccLXRd3z7QVuWHOxWTLKdjkuVUTIp6PR2TImmZ2eo9pnVmrH1RYQHKOgcRl/MKUZcj6Px9vGx1WAghJRyrUn54eMD/Q1doW0Htrg5TfhDiHmgRn0BnIrLA7BHNLU7VUZx1ZWfrJDopVU5ezRFwStRdzX2NSZHrGVn5bhf+dZVDYKG7MdRaXRt6LRckAb4UdISUBhLtUXuUoo0Q4krYMreaPfK04fZ6JSkt1wqXa6EzEHTX0jILXD6yrP8N61wFzUqXI+zKKNMhIaQkYDfRlpFhHGXl4+OTp01rp6WNEGJvkIoD0aQIToCvG4ZEvTAm6eR1FQbujzHJ6XohZzjkCqtdYmpmodUhbvjPaZa6HHFX1t/HLn0mhLiZaLMGijZCCCkaccnpeouc6WtscnqBy5YL8jXwn8t5raEiXwMlNNCXp4QQNxVtVtnXhw4daou+EUIIKYSwIF81NasWluezhOsZciYmRU4q/7kbFjq8Xr2WpkQdpt1n4vMsGxLgc8N3ztCHrnyQlA/ytfrhnBDiOFh71AQGIhBC3Bn4yZ02sMydMhB10YlpBS4b7OethleNBF3u/+HBfhR0hLjT8GhpgKKNEFJSSUnPlDOxKXl86CDwLiRcV4mH8yPQ10uqlQvMzT9nbKWLCPZXqVUIIS40PFoQmzZtUlNoaKgMGTJEvRJCCHEdAn295ebIsmoyBbVZz0LQ6Ydab1jrUMM1JT1LDl1KUpMpft6eedKWaMmFK4UE2C2Yg5DShlWWtjNnzsi9994rvXr1kldeeUXf/sYbb8iLL76of1+zZk3ZunWrhIeHi7tBSxshhBiDkl/n4iDmUnLz0d2w0p2Nu66ibvMDNV6rltOSC7P8FyEOGx6FUPvpp59k27ZtEhgYqNrOnTsntWvXlo4dO8rTTz+t3kPAjRs3TmbOnCnuBkUbIYRYDkqBXYi/fsNCl5uHDhMsdygZlh8o/1W1XGCeSFeW/yKljUR7DI/+8ssvSohpgg2sXLlSvV+6dKl+Q2XKlJHp06e7pWgjhBBiOT5eGBrNsaKJGI+uwAIHQXcjXUne8l+w3GESuZJv+S/j1CUs/0VKL1ZZ2iDKoqOjJSAgQN82evRoSUtLk++++07fduXKFalRo4YkJ+ctrOzq0NJGCCH2B+W/LiWmGuegK0L5Ly3/nBbpimAJlv8i7oZdLG2+vr5qhYaibceOHTJy5Eij+YKC8MRlHy5duiSHDh2SihUrSv369c2GoKempsr27dtVct/WrVuLv7+/3fpDCCHEehBtWjk0QE3taxdc/gtpSwyTCyOtyfn462rafDym0PJfNXMtgWgLYvkv4sZYJdpuuukmWbhwoTz77LPq/d69e+XAgQPSrl07o/ng84ZgBFty8uRJeeKJJ2T37t1St25dOXz4sERGRsqSJUuUT53htvv166cUq5eXl7L6LVu2TG699Vab9ocQQoh9wMN4xbL+akIpMUvKf0HQYZg1KTVTWfAwbT0Zm2fdyDeXJ7kwy3+Rkjg8+tVXX8mDDz4ow4YNk6ioKCXgvL295fjx40ogadxxxx1StWpV+fzzz23WUYixq1evyp133qneY0gW28nKypJ169aptszMTCXounTpovoKHn74YeWLhz7CUlgYHB4lhBD3BD9n8SkZRS7/hYoQN4IhjCNdWf6LuF30KAQSRNvcuXPVl6NSpUryww8/SPv27fXzvP/++zJ+/HhZv369iii1Jx9//LFMmDBBUlJS1Pu1a9dK165dlRUO4k2z0NWqVUt+++03laqkMCjaCCGkZILyX/pgCDPlvwpCK/8FHzpjK12gqvXK8l/EZSsiIBghLi5ODZfC0mbI0aNHlbirV6+e3S9iWPwg0Hbt2qXez5o1S6UbuX79utF8OBAvvfSSTJw4Mc86YLHDZHjgYCUs7MARQggpmeW/THPRWVv+KyfKNceHLrwMy38RJ1dEiIiIUJM56tSpI45g8eLF8v3338uKFSv0bfHx8VKunLH/Ayhfvrz6zBwzZsyQqVOn2rWvhBBCXJsyft7SsHKImiwt/4UAiQsJqZKUlin7zyeqyVz5L1PLnFYxomKwH8t/EauwWRkrR7J69WoZNWqUvPPOO9K7d299O3zWtKFSQ9CWnz8bLHPPPPNMHksbIYQQYovyXwcvJqrJFH8fT6leLjcHXa5lTstJx/JfpNiibdCgQdbMrvzdbM0ff/wh/fv3l9dee01VYDAEEaswLSYlJUlwcLBqw1BpTExMvtGsfn5+aiKEEEKsxd/HS+pEBKspv/JfOWlLUvKU/0rNyJbD0Ulqyq/8V04NV+NI18qh/uLt5WnTk4VEyNtOxsrlpFSpGJwTtcuasW4u2n788Ufx8fGxKArTHqxZs0al85gyZYoKQDClW7duyscOpba03HH4H2573bt3d0KPCSGElFZ8vT2lVngZNRVW/ivHjy5H4MFyl56VLcevJKupsPJfStCpJMNBUiUsQFWpsIZV+y/K1BUH5GJCqr6tUoi/TO7bQHo1qlTEvSdOF21ImpuRkaGiMMeOHateDVN92BMk8b3rrrukU6dO0qBBA1U+SwP9gFhDwt0XXnhB5XODdQ19mzx5shr+rFKlikP6SQghhNij/JdKMgxBV0D5L1jHokIDjIZaNWsdLHd+3l55BNsjC3aJaUTipYRU1T57RHMKNxfCqujRa9euqWS2SPmxYcMGldwWFq0xY8YoIWVPkLIDKT7MgYAEwyoMKKm1fPlyZWGDz9vw4cMtjmRlyg9CCCElufyXln+uWrkA+WzdCYlLyTA/P6pLhPjLxoldOVTqzik/wLFjx2TevHny9ddfy9mzZ1W5KFjf7rnnHgkNDRV3haKNEEKIO6KV/zIcajUt/1UUvn2grbSrXd7m/SUOFG0a2dnZ8ueffyrrG/zHgGmeNHeCoo0QQkhJw1z5r43HrsquM+bTYRny/j23SL9bohzSz9JKoj3ztJlSTN1HCCGEEDsCF6EKZfzU1KJ6Tj7TdrUryLAvthS6LKJJiWtQ5JhhVCJAjrNq1aqpQIATJ07Iu+++KxcvXrRtDwkhhBBic5DWA1GiBXl8l/X3VvMRNxRtqCrw2WefSbt27eTmm29W/mxw8j9w4IBs3bpVHnnkEbf2ZyOEEEJKC4g0RVoPkJ9wS0zNlP+t+E8ys7Id2jdiA5+2gIAANRSK1BuIGO3Zs6fDUn44Cvq0EUIIKU3kl6etba3ysmz3efW+c71w+XBYMwn293FiT0sudglEwJg4RJppkfj8SE29cQG4CxRthBBCShv5VUSAoBv//b+qekO9iGD5akxLqRIW6OzuljjsEogwdOhQW/SNEEIIIS4EBJq5tB6oiLA4NEDun79Dldvq//Em+XxUS2leLcwp/SztFDvlR0mDljZCCCHEmIsJ1+W+eTvkwMVEVZ7rncFNpW/TyjxMDtYetq04SwghhJASR6WQAFnycDvpXr+iKqP1xLe75cM1R5nyy8HYRbRt27aNBdoJIYSQEkSQn7d8NrKl3N+hpnr/zh9H5NnFeyQtM//SWcTJog2jqdu3b5dly5bJzp078+RuGzhwoLRp00Y2bdpky34SQgghxAV8317p00CmD2ik/l+6+7yM+HKrxCanO7trpQKrRFtcXJy0b99e1RmFOGvZsqX06NFDla368MMPpUmTJrJy5Up56KGHVG1SQgghhJQ8hrepLvPGtpJgf2/ZfipOBSgcu3zN2d0q8VgViPDUU0/Jp59+KqNHj5Y6deooyxoS7CJfG8QahNzbb78ttWrVEneFgQiEEEKIZRy7nCRj522Xs7HXlYCbPbyFdKhTgYfPFfK0QYy99tprqgqCxjfffCOjRo2Sl19+WaZNmybuDkUbIYQQYjkx19LkoW92yo7TcWrIdFr/RjKsdTUeQmeLNl9fX7l69arRCrUNXblyRSpUcH91TdFGCCGEWEdqRpa88ONe+enfC+r9Ax1rygt31Ff/m0vaSxxUEcHc7Pm1uyMUbYQQQoj1QAd8+NcxefePI+p9k6gQiU5KlejENKPyWKh3iqS9xAF52jIzM/NMBbUTQgghpOQDA86T3erIB8Oaibenh+w9n2Ak2MClhFR5ZMEuVR6LWI/Vos3HxyfPVFA7IYQQQkoPvRtXkpAA8xpAG5NDgXrUOyXWwdqjhBBCCLEZ8GGLKSBvG6TaxYRUNZ+5eqfERqLtu+++s2Z2QgghhJQyEHRgCefjU0SEos0aWHuUEEIIITYDUaKW8L8VB+T9P4+ymoIVULQRQgghxGYgrQeiRAtK7IGsH4mpmfLen0ek3Yw18tKyfXL8CisqFAZFGyGEEEJsBvKwIa0HMBVuHrkTIkwxNakSImmZ2bJo6xnp9s46uX/+dtlyIqbEpBGzNVblaSsNME8bIYQQUnyQ1gNRogg6yC9PGyQIAhK+2HBC/jx4WT8fxNx9HWrKnY0riY9XybcvJdojuW5pgKKNEEIIsQ1I62FpRQQMj3618aT8uPOcsr6BqNAAGdO+htzTuqoE+5fcVGIUbXY+cIQQQgixTy3TBVvOyNf/nNKnDinj5y3DWleVMbfWVEKupEHRZucDRwghhBD71jP9afd5+XLjSTl2OSdIAVY6JO99oGMtaVwlpMQcfoo2Ox84QgghhNif7GydrDtyRfm9bT4eo29vU7OcEm9db64onm5ehJ6izc4HjhBCCCGOZf/5BOX3tmLPBcnMLYNVKzxIBS3c3byK+Pt4ueUpoWiz84EjhBBCiHO4mHBd5m0+pVKFJKVmqrZyQb4yom11GdWuulQo4+dWp4aizc4HjhBCCCHO5VpapizeflbmbDop5+KuqzZfb08Z2CxK7u9YU26qGOwWp4iizc4HjhBCCCGuQWZWtvz+X7R8vuGE7Dkbr2/vUi9c+b2hML2Hh+v6vVG02fnAEUIIIcS10Ol0suN0nHyx/oT8cTBatEy0DSuXVZa3Pk0qu2SyXoo2Ox84QgghhLguJ68my5yNJ2XJzrOSmpGTrDeyrL+MvRXJeqtJSIDrJOulaLPzgSOEEEKI6xOXnC4Lt56WeZtPy9VraaotyNdLhraqpgRc1XKBzu4iRZu9DxwhhBBC3Ie0zCz5+d8L8tWGk3I4Okm1Ib3bHbnJem+pGlqsMlzFgZY2Ox84QgghhLin39v6o1flyw0nZMPRq/r2VjXC5P6OtaR7/QglzCwpeG8rKNrsfOAIIYQQ4t4cvJgoX244Kcv3nJeMrJyohRrlA6VtrfLy/fazkhvHoEezsc0e0dymwo2izc4HjhBCCCElg+jEVJm/+ZQs3HpGEq5nFDgvhFtkiL9snNjVZkOllmoP14t7JYQQQghxIBFl/eX5XjfL5he6ypj2NQqcF9Y3DJnC183RULQRQgghhIhIkJ+3NKuWNyDBHAhOcDQUbYQQQgghuSBK1Jbz2RKKNkIIIYSQXJDWA1Gi+XmroR2fYz5HQ9FGCCGEEJILgguQ1gOYCjftPT63R762wqBoI4QQQggxAOk8kNYDUaKG4L2t031Yg7dTtkoIIYQQ4sL0alRJbm8Q6ZCKCJZC0UYIIYQQYgYItHa1y4urwOFRQgghhBA3gKKNEEIIIcQNoGgjhBBCCHEDKNoIIYQQQtwAijZCCCGEEDeAoo0QQgghxA2gaCOEEEIIcQMo2gghhBBC3ACKNkIIIYQQN4CijRBCCCHEDaBoI4QQQghxAyjaCCGEEELcALcSbVeuXJE333xTGjZsKE2aNDE7T2pqqrz44otqngYNGsiECRMkJSXF4X0lhBBCCCm1oq1Hjx5y9epVuf322+XChQtm57n//vtl8eLF8vHHH8tnn30mP//8s4waNcrhfSWEEEIIsSXe4kbs2LFDvLy8ZNasWWY/P378uCxcuFB+++036dy5s2r75JNPlMg7cOCAsrwRQgghhLgjbmVpg2AriPXr14unp6d07dpV39alSxfx8/NTnxFCCCGEuCtuZWkrjHPnzkm5cuXE19fXSOiVL19ezp8/b3aZtLQ0NWkkJCSo18TERAf0mBBCCCGlncRczaHT6UqPaMvOzhZv77y75OPjI1lZWWaXmTFjhkydOjVPe9WqVe3SR0IIIYQQcyQlJUlISIiUCtEWHh4uMTExSql6eHjo2xG8gM/MgUjTZ555xkj4xcbGKuuc4TpsraghCs+ePStly5a1yzZKAyXhOLrKPjiyH/bYli3X6SrnhBBHY69r392/U4kO6D90CwRb5cqVC5yvRIm21q1bS0ZGhmzfvl39D3bv3i3Jycn696bA3w2TIaGhoQ7pL06+O17ArkZJOI6usg+O7Ic9tmXLdbrKOSHE0djr2nf371RZO/e/IAubWwYiFEbLli2lTZs28tJLLynFCrEGS9ott9wi7du3d3b3CCGEEEKKjFuJtjFjxkiVKlVkypQpEhcXp/7HdPjwYf08S5YskczMTDW8iaCEa9euybJly+w21EkIIYQQ4gjcanj0nXfekevXr+dpj4yM1P+Pcee///5bjUFjjNgSc6OjwXDs5MmT8wzLktJ3HF1lHxzZD3tsy5brdJVzQoijsde17+7fKT8X6r+HrrD4UkIIIYQQ4nTcaniUEEIIIaS0QtFGCCGEEOIGULQRQgghhLgBFG0O5MyZM/LKK69InTp1pHfv3o7cdIni1KlT8tRTT0nz5s2lWbNm8sgjj+RbpswVQQLnRYsWSc+ePaVu3bqqPu7cuXOd2idEWSNdDqKxUQ7O1nzyySf6aG9tqlevnk3WfejQIRkxYoRaX+fOnWXVqlVWr2PgwIF5+oepf//+NukjIa4CsivMmTNHunfvru4/eP3uu++M5vn111/Nfh+uXLli0TZQaQjz49WQgwcPyoMPPihNmjRRKbqeffZZlfzeVUhLS5M33nhDWrVqJY0aNVLpw1JTU/WfIw/s559/Lt26dVPHrkePHvLjjz86tI9uFT3qzuBiuO2222Ts2LHqgjhy5Iizu+S2DBgwQMaNG6dSwKA8Gb5YnTp1UomU3SFx4/z582XLli0yceJEFe28efNmefTRR1XEM8SoM3jsscdUzV6IX9zUbQ32DWl4fvnlF32bp2fxnxn//fdfde4h2nDzxA0WZemQTBspfyzls88+M6pBHB8frx4KIGQJKUngAer48ePKgIDs+2vXrpXRo0er6x+vICUlReU6/e+//4yWxXe4MDZu3ChfffWVSrOl1fLW1jlkyBB1j3v88cfVgyKqEUH44H5oWDPcWQ/T/fr1k9OnT8v7778vtWrVkp9//lk+/fRTGT9+vJpn1qxZcvHiRRVJiqwVq1evlnvuuUcWLFggQ4cOdUxHET1KHENWVpZ6feqpp3QtWrTgYS8imZmZRu8vX76MCGjdsmXL3Oo6MGTw4MG6O+64wyn9+eabb9T1uHr1anUcT548afNtzJgxwy7XfLt27cwet+zs7GKt9+OPP9Z5e3vrLly4UKz1EOIO95/bb79dN2zYMP37JUuW6EJCQqxed0xMjK5GjRq6DRs26OrVq6ebOHGi0XfSdNt79+5V95x//vlH52zmz5+vvvPHjh0zajfss7lj17FjR93YsWN1joLDow7EFpYFIuLl5WV0GGCyBj4+Pm55HeCp959//lHDpI7m6NGj8txzz8nChQvtfvxgXa5fv76qUPLAAw8Ue0gbw7g4brBem1LcZNqwFPTp00cqVapUrPUQ4ur3nwMHDqhRCtP7DyoKYRgTw4TDhg1T8xXGfffdJ8OHD5cOHTqY/U6abtuV7t2LFy9WVvvatWsbtRv22bT/e/bskX379jn03k0VQdwelCqDmR/+TO4EfLAqVqyofCNgYp8wYYJDt5+enq62i+FEW/mX5UdQUJC8/PLL6sb40UcfqSEIDD9a6iNjDs3FADd8+AfCVxSvRfFpMx1y3bVrlxKWhJRE4FYCn7Pw8HAlzDBkaXi9e3t7q6FL+NpiwncM31fT4VJD8L3GgxgqFlkCUsTintCgQQP1IOdsjhw5ovoCtxXcD+FiMWnSJCVeDcEwL45dhQoVpEWLFso9Z+TIkY7rqMNsekQPh0dtxxtvvKHz9/fXrV+/3u2usPPnz+uOHDmi++qrr3RlypTRffDBBw7d/tNPP62766679O/Xrl1rt+FR0+HK5ORkXUREhG7y5MlFXueqVatUf6tVq6b7+eefdYcOHVLDsB4eHro//vijyOt9/PHHdVWrVjU7FEJISeHs2bPqOwNXgICAADU8mN/3Fe+bNWumGz58uNl17dmzR1e+fHl1P9MwHR41Zfz48WoIFkOkrgCGdX19fXUTJkzQHThwQLmL1KpVSzdo0KA8xwLH7uDBg7r3339f/f589913DusnRZsToGizDe+++676wuDH292BeIGIcSS1a9dWN9qoqCg1VahQQYmgyMhI3SuvvGL37Xfv3l03ZMiQIi+/a9cu1V9cB6Y+JiNGjCjSOlNTU3VhYWHFEpOEuONvEkRWQTz88MO65s2bm/3s7bffVvdi7V6CCf5hwcHBuurVq+eZH2KubNmyuq1bt+pchdatW+vq1q1r1LZw4UL1EJiUlJTvcg899JCuadOmOkfB6FHiliC6B2bpZcuWqSExdycwMFCZ3R0JorwMI0XhH4boLhzTm2++2a7bxgPjyZMn5aabbiryOjCUERwcLGXKlDFqR5u5GsWWsHTpUhXxhuhkQkoLltx/Tpw4kW9E9sMPP6xcLQyBuwruzbhPG4L3s2fPVpGXGIJ0Fdq2bavuiab3EtyrEJVuep9x2r3bYfKQ6KGlrXh8+OGH6qnut99+c8uratq0abrt27frh9/wf+XKlXXjxo1zar/sOTyKa/7o0aPq/+vXr+uee+45nZeXl27Lli3FWi+GMvCUe/HiRf0+YIjj66+/LtL6unbt6rQoXkIcwauvvqqGM7UhUER6wuL+zDPP6OeZNGmSbvfu3WoeROvjnot7gzXDgOaGR2HBh4WtuN97e4DhTj8/P92PP/6o3sfGxupuu+02Xfv27fXzvPjii7r9+/frj91ff/2lhnhfeuklnaOgaHMg3bp1U2Zj+C/5+Pjozcjw7yGWgRsIzNWmpnhMn376qVscRoS3d+7cWQ0dYMKw5AsvvKBLSUkpsaINN8JGjRqpGxxujPCPKY7fmUZ6eroanoBPDoY18ePz1ltvFWldJ06cUNfW0qVLi90vQlwVfM+RKgfiKSgoSLlDTJkyRX2XNNatW6efB/faOnXqqKFCazAVbfDhxf0F2zS9d7vKd27FihXKRxb3Kex3nz59lP+aBu5ZGEbFcQkMDFQP23gIz8jIcFgfPfDHcXa90k10dLQ+xNmQqKioYqcoKE3kl7E/NDQ0XxO2KwKTO64HmOBdASTXRDQn0lyYplWxFUjY6e/vb/MQf/QdUV7WJNQ1BUOqMTExKhKZ6XlISQfXOxLKIrI7PzDsh9+mgICAIv3eYTkt4TkiVpGY1hz43mKY0RXQ6XTqPoDfE0TR5ndcMF9Bx85eULQRQgghhLgBzNNGCCGEEOIGULQRQgghhLgBFG2EEEIIIW4ARRshhBBCiBtA0UYIIYQQ4gZQtBFCCCGEuAEUbYQQQgghbgBFGyEmVKhQQe6//363OS4zZsyQmjVrqoS49q4ZqoGaoaa1Bksijt5Pd7v2igNqU6Leo73PCxLINm/eXCZOnFjkdbg7Rd1XJMStX7++TJkyxS79ItZD0UYcSsuWLVWGbXNF3q9du6Y+e/zxx3lWLGTlypWqAPOsWbNUdYVDhw7x2Lk4VapUkREjRji7G6WGefPmybFjx+T55593dlfcDjwIQrC9/fbb+VaiIY6Foo04hdWrV6uJFI9169ap8i/9+vVj6SU7gB/77777zh6rJg4ApYbeeOMNGTZsmJQvX57HvAgMGjRIldp7//33efxcAIo24nDq1q0rtWvXVk++GLogRefq1atFqgtISGlgzZo1cvToURk5cqSzu+LW1jYMrc6dO1fV+CXOhaKNOBwUC4cf1p49e+Sbb74ptKg6hkzN+VR06NBBTeZ8gv755x9p06aNskI1a9ZMvQebNm1SfjQQOhCPGF7Mj82bN0vr1q3VvLVq1ZL33nvP7HywxLRr104VD8bUtWtXtR1z/dq6davqM9b5yiuv5LttiFkMScCfxM/PTy0/ePBgOXLkiF6s4bhg6AfFjfE/JgyTFsT333+vto8nZxRp7tOnj+zatcvi7brSedLWsXv3bv0xrV69unzwwQd55l2xYoXceuutEhYWpqZOnTrJ8uXLpSi+QNZs15DMzEx1jM6fPy8LFy7UnzO4DJhi6botufbMYYvjb+m1kp6eLi+88IJUrlxZ9bF79+4FXk9F3Sdz4LxjHyz1nbN0nyBeJkyYIJUqVTLaJ0v99Cy9Hgv7vuIz7TqCuEJ/hg8fLmfOnLFofy091mjHfQb3ROJkdIQ4kBYtWugaNmyo/m/Tpo2uatWquuvXr6v3SUlJOlySjz32mH5+fIa2yZMn51nXrbfeqiZDypcvr7v99tt1gwcP1h09elR3+fJl3ZAhQ3ShoaG6zZs36/r166c7fPiw7sqVK7p7771X5+/vr7t06VKedXTv3l3Xp08f3YEDB3QxMTG69957T+fl5aX73//+ZzTvlClTdN7e3ro333xTd+HCBbWuZ599Vufr66vbtGmT2XX+999/unPnzumWLVuW73EaO3aszs/PT/fFF1/o4uLidHv37tW1bdtWFxYWpjt27Jh+vtGjR6t1W8Irr7yi9gHH8uTJk2q9v/76q27UqFFWb7d27dq6oUOHOv089ezZUzdo0CB1ntBf7CP68ccff+jn27lzp9rvSZMmqe3gOsO5wTbi4+MLPGam+2nNdvMjKipKN3z4cLOfWbNuS6+9/LZT3ONv6bWC5cuUKaNbtGiROt7btm3T9ejRQ90LcA8oyj6ZOy/maN68ua5du3ZmPzO3Dkv3CcsFBwfrvv/+e11CQoJu+/btujvvvNPsPpli6fVoyffVEHwHcWyx/caNG+vS0tIK3Fdrrp+LFy+qaxDLEOdC0UacJtrWr1+vbgRvvPGGTUVbSEiI0c3vzJkzah2VKlXSxcbG6tvPnz+v2t9555086wgKCtJdvXrVqH3MmDHqx0tbx/Hjx9VNFTc6U/BDcdtttxmtMyAgQP0IFsaePXtUv1566aU8N05s3/AH31LRduTIEZ2np6fuySeftMl2bSHabHGeypYtq37MNLKysnRVqlRRYkRj1qxZann8uFpLfqLNku0WVbRZsm5rrr38tlOc42/ptaLNN2PGDKP5du/erdoNBY41+2SpaKtQoYJu4MCBZj8zXYel+6T1HWLHkH379uk8PDwKFW2WXI+WfF/zY8eOHWr9f//9d777au31k52drfoDUUucC4dHidPo2LGjcqDHUClM77aiffv2EhISon9ftWpVNfzTqFEjNRShgeGasmXLyokTJ8yuw9RxuX///moYUBsi+O2331RIPIZPTOnWrZsaZsCQmOE6MdxSGH/99Zd6HThwoFF7ZGSkGlKBn461rFq1Sg393HvvvQ7drr3PE/oVGhqqf+/p6amGtgznbdq0qXpFxCb2EeewuFiyXXuu29prz9bH39JrRXu96667jOa75ZZbpEaNGkZtttgnU+Lj49XQoiVYuk9r165Vr3379jWaD8cNqXcKw5Lr0ZLvKzh48KAavsc58vb2NhpuRxBNflh7rLHeMmXKqONJnAtFG3Eqb775piQnJ8trr71m9bKwFJsDfh2m4MadX7u5G1FERES+bfAnA5cuXVKvuKHjhgmfEvzAYpo2bZq66SUmJuqXj4qKsmi/NAGLHwtT0KZt3xouX75caB/ssV17nydz80JgGM4LP6Ovv/5a+fngBwlCBT46P//8cxH2xvLt2nPd1l57tj7+ll4r2nwFfZ9suU+mQPxauoy1+1SxYsU885lrM8WS69GS7yvmwbG6ePGi8pFLSEhQ3zUt7Q9SAOWHtcca60VKJkMxT5wDRRtxKvXq1ZP77rtPZs+ebdZKAWdgf39/SUpKyvMZHLrNgadCa9rNER0dnW+bZoHTrGYIqMBNDk+ueDrGlOt6oJyHDQMwLEFbJr8+WGKtMyU8PLzAY1bc7TrrPFk6L6IH//33X7ly5YpyvsZysJzColEUrOmjPdZt7bVnzXYs2b6l14r2XSno+2TLfTIFQRyaQLH1PmnCyhBzbUW5Hi35viI4JC4uTj788ENp0aKFCiQAJ0+eLHT71h5r7D8+w/EkzoWijTidqVOniq+vr7z44ot5PsPNDEMO+/fvN2rft2+fxRFSRQFDoLGxsUZteBKGMMGwErjzzjvVkykivGwJnr7BsmXLjNpx48Swhfa5Ndxxxx2qr4sWLbLLdp11nqwFP1YDBgzQn7P169c7vA/4cS1u6gR7XXuWYum1AgsSgCXIkL1798qpU6fsvk+IroQ4smRY1dJ96tKli3r95ZdfjOY7cOCARYLJkuvRku+r4QOTIbDiFYa1x3rbtm16lxbiXCjaiNPBMMlzzz0nv/76q9nPH3jgAfnzzz9VehCY6Lds2SKTJk1SfjH2AmH7Y8eOlcOHD6unWaRcwM0QpXC0J9A6derI5MmTVfJODCngR+j69evKzwSJKItajqhJkyYyatQomTlzpsqNhGEpiKG7775biVts01rQV4jijz/+WInk06dPq/XCt2X06NE22a4zzpMlvP7666pqBNJoYCgew1ywTmhDVY4Gvk87duwoVoZ5e117lmLptQL/LfhcTZ8+XQkEDLlh3/E9gnXI3vsEvzOsA9eirfYJ1/PQoUOVS8eSJUvUPiENB/YJ5bJscT1a8n2FiIQPIraLIVJYFJFyxxKBau2xhu8dLIzaAytxIk4OhCClOHrUkGvXrqmoNdPoUZCRkaGinMLDw3WBgYG6O+64Q0W65ReVeN999+VZf0REhIq0tCSST1sHolvRX6QAqFGjhm7mzJlm9wmpO7p06aKi8RAhiv1Df0+fPl1ov/IjMzNTRdXWq1dP5+Pjo1IOIAru4MGDRvNZk/IDLFiwQEW34ThiOaQg2bVrl9XbNRe956zzZMrdd9+t+qeBSEhE+iH9A6KCEVHYtWtX3YoVK4ocPWrJdguKDOzYsaM6RrjecY0Vdd2WXHvmsMXxt/RaSU1N1U2YMEEXGRmp+ti5c2c1DyIUzUVaWrJPlkaPIuqxbt26ugcffDDPZ+bWYc0+Pffcc+p4oY+4npAiBelBOnXqVGCfrLkeC/u+rlmzRteqVSvVB0QYv/rqqyo1Ca6r2bNnF3q8LDnWOCbYT+wvcT4e+ONM0UgIIYTYCySgfvLJJ9XQpb1LWSECF9YoZw1b2wPsy7hx49SoA+rmEufC4VFCCCElFgx5orIFItXtCfy+MOR92223SUkBwQcYckX1Bwo214CWNkIIIcQKvv32W+VDhohPBBOgPN3DDz+sPkPgA3KaEWIPaGkjhBBCrKB3795y4cIF6dWrlxJtSOGBSFVEf1KwEXtCSxshhBBCiBtASxshhBBCiBtA0UYIIYQQ4gZQtBFCCCGEuAEUbYQQQgghbgBFGyGEEEKIG0DRRgghhBDiBlC0EUIIIYS4ARRthBBCCCFuAEUbIYQQQoi4Pv8HU/g/tBM+LsQAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(7, 4))\n", "ax.plot(curve['columns'], curve['training RMSE'], marker='o', color='tab:blue', label='Training RMSE')\n", "ax.plot(curve['columns'], curve['test RMSE'], marker='o', color='tab:red', label='Test RMSE')\n", "ax.set_xscale('log')\n", "ax.set_yscale('log')\n", "ax.set_xticks(curve['columns'])\n", "ax.set_xticklabels(curve['columns'])\n", "ax.set_yticks([10, 20, 50, 100, 200])\n", "ax.set_yticklabels([10, 20, 50, 100, 200])\n", "ax.minorticks_off()\n", "ax.set_xlabel('Number of columns in the model (log scale)', fontsize=12)\n", "ax.set_ylabel('RMSE ($1,000s, log scale)', fontsize=12)\n", "ax.set_title('Training error keeps falling; test error turns upward', fontsize=13)\n", "ax.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "1c536644", "metadata": {}, "source": [ "This pattern is **overfitting**: a model fits its training data closely, including its random quirks, and so predicts new data worse. The gap between the two curves is the warning sign. A model whose training error is far below its test error has learned the sample, not the market." ] }, { "cell_type": "code", "execution_count": 9, "id": "3956e2a3", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.564111Z", "iopub.status.busy": "2026-10-24T05:29:30.564047Z", "iopub.status.idle": "2026-10-24T05:29:30.566279Z", "shell.execute_reply": "2026-10-24T05:29:30.565938Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Lowest test RMSE: 33.0 with 5 columns\n", "At 63 columns: training RMSE 12.6 and test RMSE 181.6\n" ] } ], "source": [ "best = curve.loc[curve['test RMSE'].idxmin()]\n", "print('Lowest test RMSE:', best['test RMSE'], 'with', int(best['columns']), 'columns')\n", "print('At', int(curve['columns'].iloc[-1]), 'columns: training RMSE', curve['training RMSE'].iloc[-1],\n", " 'and test RMSE', curve['test RMSE'].iloc[-1])" ] }, { "cell_type": "markdown", "id": "15bb60c9", "metadata": {}, "source": [ "The lesson is not that big models are always bad. The more training data you have compared with the number of columns, the more complexity a model can support. With 100 homes, five or so columns is the limit; with thousands of homes, as in {doc}`1302-multiple-regression`, seven numeric features are no problem. The right amount of complexity depends on how much data you have." ] }, { "cell_type": "code", "execution_count": 10, "id": "bffbddba", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.567218Z", "iopub.status.busy": "2026-10-24T05:29:30.567162Z", "iopub.status.idle": "2026-10-24T05:29:30.568908Z", "shell.execute_reply": "2026-10-24T05:29:30.568656Z" }, "tags": [ "thebe-interactive" ] }, "outputs": [], "source": [ "### Exercise: The Training Gap\n", "# 1. Load ames_housing.csv. Split homes.index with train_test_split(train_size=100, random_state=2); the other rows are the test set.\n", "# 2. Define two feature lists: simple = living_area and overall_quality; full = all seven numeric features, in this order: living_area, overall_quality, year_built, garage_cars, full_baths, lot_area, total_rooms.\n", "# 3. For each list, fit LinearRegression on the training rows and compute the training RMSE and the test RMSE in $1,000s, rounded to 1 decimal.\n", "# 4. Print both RMSEs for each model and the gap (test minus training). In a comment, say which model looks more overfit.\n", "# Print exactly:\n", "# simple train 36.6 test 35.1 gap -1.4\n", "# full train 28.1 test 33.4 gap 5.3\n", "path_data = '../../data/'\n", "import numpy as np\n", "import pandas as pd\n", "from sklearn.linear_model import LinearRegression\n", "from sklearn.model_selection import train_test_split\n", "\n", "\n", "### Your code starts here.\n", "\n", "\n", "\n", "\n", "### Your code ends here." ] }, { "cell_type": "code", "execution_count": 11, "id": "55d00264", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.569746Z", "iopub.status.busy": "2026-10-24T05:29:30.569691Z", "iopub.status.idle": "2026-10-24T05:29:30.577103Z", "shell.execute_reply": "2026-10-24T05:29:30.576809Z" }, "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "simple train 36.6 test 35.1 gap -1.4\n", "full train 28.1 test 33.4 gap 5.3\n" ] } ], "source": [ "### Solution\n", "path_data = '../../data/'\n", "import numpy as np\n", "import pandas as pd\n", "from sklearn.linear_model import LinearRegression\n", "from sklearn.model_selection import train_test_split\n", "\n", "homes = pd.read_csv(path_data + 'ames_housing.csv')\n", "train_rows, test_rows = train_test_split(homes.index, train_size=100, random_state=2)\n", "train, test = homes.loc[train_rows], homes.loc[test_rows]\n", "\n", "simple = ['living_area', 'overall_quality']\n", "full = ['living_area', 'overall_quality', 'year_built', 'garage_cars',\n", " 'full_baths', 'lot_area', 'total_rooms']\n", "\n", "def rmse_k(actual, predicted):\n", " return np.sqrt(np.mean((actual - predicted) ** 2)) / 1000\n", "\n", "for name, features in [('simple', simple), ('full', full)]:\n", " model = LinearRegression().fit(train[features], train['sale_price'])\n", " tr = rmse_k(train['sale_price'], model.predict(train[features]))\n", " te = rmse_k(test['sale_price'], model.predict(test[features]))\n", " print(name, 'train', round(tr, 1), 'test', round(te, 1), 'gap', round(te - tr, 1))\n", "# The full model has the smaller training RMSE but the larger gap, so it is more overfit." ] }, { "cell_type": "markdown", "id": "395074a5", "metadata": {}, "source": [ "## Validation\n", "\n", "```{index} validation set\n", "```\n", "\n", "```{index} model selection\n", "```\n", "\n", "We used the test set above to *show* overfitting. But suppose you used it to *choose* a model: try ten models, pick the one with the lowest test RMSE, and report that number. The number is now too optimistic, for the same reason as before. You picked the winner because it did well on those particular homes, so those homes are no longer an honest test. The test set is meant to be used once, at the end.\n", "\n", "The fix is a third group of homes. We split the data three ways:\n", "\n", "- the **training set** fits each candidate model,\n", "- the **validation set** compares the candidates and chooses one,\n", "- the **test set** is touched once, to report the chosen model's error.\n", "\n", "Here is a split with 150 homes for training, 150 for validation, and the remaining homes for testing." ] }, { "cell_type": "code", "execution_count": 12, "id": "e06e8c72", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.578095Z", "iopub.status.busy": "2026-10-24T05:29:30.578027Z", "iopub.status.idle": "2026-10-24T05:29:30.580333Z", "shell.execute_reply": "2026-10-24T05:29:30.580036Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Training: 150\n", "Validation: 150\n", "Test: 2112\n" ] } ], "source": [ "train_rows3, rest = train_test_split(homes.index, train_size=150, random_state=7)\n", "valid_rows, test_rows3 = train_test_split(rest, train_size=150, random_state=7)\n", "\n", "print('Training: ', len(train_rows3))\n", "print('Validation:', len(valid_rows))\n", "print('Test: ', len(test_rows3))" ] }, { "cell_type": "markdown", "id": "b2d4865d", "metadata": {}, "source": [ "We fit six candidate models on the training set and compare their RMSE on the validation set." ] }, { "cell_type": "code", "execution_count": 13, "id": "7547238b", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.581241Z", "iopub.status.busy": "2026-10-24T05:29:30.581174Z", "iopub.status.idle": "2026-10-24T05:29:30.595929Z", "shell.execute_reply": "2026-10-24T05:29:30.595614Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
columnstraining RMSEvalidation RMSE
0143.650.0
1232.735.7
2430.135.7
3726.933.5
43522.332.8
56313.836.1
\n", "
" ], "text/plain": [ " columns training RMSE validation RMSE\n", "0 1 43.6 50.0\n", "1 2 32.7 35.7\n", "2 4 30.1 35.7\n", "3 7 26.9 33.5\n", "4 35 22.3 32.8\n", "5 63 13.8 36.1" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "candidates = [1, 2, 4, 7, n_numeric + n_nbhd, X_all.shape[1]]\n", "y = homes['sale_price']\n", "\n", "rows = []\n", "for k in candidates:\n", " cols = list(X_all.columns[:k])\n", " model = LinearRegression().fit(X_all.loc[train_rows3, cols], y[train_rows3])\n", " rows.append([k,\n", " round(rmse(y[train_rows3], model.predict(X_all.loc[train_rows3, cols])) / 1000, 1),\n", " round(rmse(y[valid_rows], model.predict(X_all.loc[valid_rows, cols])) / 1000, 1)])\n", "\n", "val_table = pd.DataFrame(rows, columns=['columns', 'training RMSE', 'validation RMSE'])\n", "val_table" ] }, { "cell_type": "markdown", "id": "28419c3c", "metadata": {}, "source": [ "The training RMSE falls steadily from \\$43,600 to \\$13,800 as columns are added, while the validation RMSE is lowest, about \\$32,800, for the 35-column model (the seven numeric features plus neighborhood). The 63-column model has the lowest training error but a validation RMSE of \\$36,100, a sign of overfitting. We choose the 35-column model and only now use the test set." ] }, { "cell_type": "code", "execution_count": 14, "id": "f896fe84", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.596999Z", "iopub.status.busy": "2026-10-24T05:29:30.596939Z", "iopub.status.idle": "2026-10-24T05:29:30.601028Z", "shell.execute_reply": "2026-10-24T05:29:30.600729Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Chosen by validation RMSE: 35 columns\n", "Test RMSE, reported once: 31.5 ($1,000s)\n" ] } ], "source": [ "best_k = int(val_table.loc[val_table['validation RMSE'].idxmin(), 'columns'])\n", "cols = list(X_all.columns[:best_k])\n", "final = LinearRegression().fit(X_all.loc[train_rows3, cols], y[train_rows3])\n", "test_rmse = rmse(y[test_rows3], final.predict(X_all.loc[test_rows3, cols]))\n", "\n", "print('Chosen by validation RMSE:', best_k, 'columns')\n", "print('Test RMSE, reported once: ', round(test_rmse / 1000, 1), '($1,000s)')" ] }, { "cell_type": "markdown", "id": "87ebe08b", "metadata": {}, "source": [ "The analyst can tell the broker that the chosen model should price new homes with a typical error of about \\$31,500. Because the test set played no part in choosing, that is an honest estimate.\n", "\n", "```{admonition} Rule\n", "Fit on the training set. Choose on the validation set. Report once on the test set.\n", "```\n", "\n", "```{index} cross-validation\n", "```\n", "\n", "### Cross-Validation\n", "\n", "With 300 homes, a validation set of 150 is a noisy judge: a different random split might choose a different model. **Cross-validation** uses the data more efficiently. In **5-fold cross-validation**, the labeled homes are divided into 5 equal groups (folds). Each fold takes a turn as the validation set while the model is fitted on the other four, and the five RMSEs are averaged. scikit-learn does the whole loop in one call, `cross_val_score`. We apply it to the 300 homes used for training and validation, still keeping the test set out of sight." ] }, { "cell_type": "code", "execution_count": 15, "id": "bbdbe371", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.601999Z", "iopub.status.busy": "2026-10-24T05:29:30.601946Z", "iopub.status.idle": "2026-10-24T05:29:30.639947Z", "shell.execute_reply": "2026-10-24T05:29:30.639544Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 columns: 5-fold RMSE = 46.9 ($1,000s)\n", "2 columns: 5-fold RMSE = 34.4 ($1,000s)\n", "4 columns: 5-fold RMSE = 33.1 ($1,000s)\n", "7 columns: 5-fold RMSE = 31.2 ($1,000s)\n", "35 columns: 5-fold RMSE = 29.8 ($1,000s)\n", "63 columns: 5-fold RMSE = 32.2 ($1,000s)\n" ] } ], "source": [ "pool = list(train_rows3) + list(valid_rows) # the 300 homes used for choosing\n", "\n", "for k in candidates:\n", " cols = list(X_all.columns[:k])\n", " scores = cross_val_score(LinearRegression(), X_all.loc[pool, cols], y[pool],\n", " cv=5, scoring='neg_root_mean_squared_error')\n", " print(k, 'columns: 5-fold RMSE =', round(-scores.mean() / 1000, 1), '($1,000s)')" ] }, { "cell_type": "markdown", "id": "998b2e83", "metadata": {}, "source": [ "Cross-validation agrees with the validation set about the main point: the 35-column model is best (about \\$29,800), and the 63-column model is worse (about \\$32,200). Because every home is used for validation once, the comparison is less affected by the luck of one split. The cost is that the model is fitted five times for each candidate." ] }, { "cell_type": "code", "execution_count": 16, "id": "3811b9df", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.640895Z", "iopub.status.busy": "2026-10-24T05:29:30.640831Z", "iopub.status.idle": "2026-10-24T05:29:30.642707Z", "shell.execute_reply": "2026-10-24T05:29:30.642381Z" }, "tags": [ "thebe-interactive" ] }, "outputs": [], "source": [ "### Exercise: Choose by Validation\n", "# 1. Load ames_housing.csv. Split homes.index into 150 training rows and 150 validation rows: first train_test_split(train_size=150, random_state=3) gives training rows and the rest; then split the rest with train_test_split(train_size=150, random_state=3) into validation rows and test rows.\n", "# 2. Candidate models use the first 1, 2, 4, and 7 numeric features, in this order: living_area, overall_quality, year_built, garage_cars, full_baths, lot_area, total_rooms.\n", "# 3. Fit each candidate on the training rows. Compute the validation RMSE in $1,000s, rounded to 1 decimal, and print it for each.\n", "# 4. Print the number of features of the candidate with the lowest validation RMSE, then its test RMSE (rounded to 1 decimal) once.\n", "# Print exactly:\n", "# 1 features: validation RMSE 56.9\n", "# 2 features: validation RMSE 41.5\n", "# 4 features: validation RMSE 37.7\n", "# 7 features: validation RMSE 35.4\n", "# Chosen: 7 features\n", "# Test RMSE: 32.2\n", "path_data = '../../data/'\n", "import numpy as np\n", "import pandas as pd\n", "from sklearn.linear_model import LinearRegression\n", "from sklearn.model_selection import train_test_split\n", "\n", "\n", "### Your code starts here.\n", "\n", "\n", "\n", "\n", "### Your code ends here." ] }, { "cell_type": "code", "execution_count": 17, "id": "82591b67", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.643546Z", "iopub.status.busy": "2026-10-24T05:29:30.643493Z", "iopub.status.idle": "2026-10-24T05:29:30.653764Z", "shell.execute_reply": "2026-10-24T05:29:30.653439Z" }, "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1" ] }, { "name": "stdout", "output_type": "stream", "text": [ " features: validation RMSE 56.9\n", "2 features: validation RMSE 41.5\n", "4 features: validation RMSE 37.7\n", "7 features: validation RMSE 35.4\n", "Chosen: 7 features\n", "Test RMSE: 32.2\n" ] } ], "source": [ "### Solution\n", "path_data = '../../data/'\n", "import numpy as np\n", "import pandas as pd\n", "from sklearn.linear_model import LinearRegression\n", "from sklearn.model_selection import train_test_split\n", "\n", "homes = pd.read_csv(path_data + 'ames_housing.csv')\n", "train_rows, rest = train_test_split(homes.index, train_size=150, random_state=3)\n", "valid_rows, test_rows = train_test_split(rest, train_size=150, random_state=3)\n", "\n", "numeric = ['living_area', 'overall_quality', 'year_built', 'garage_cars',\n", " 'full_baths', 'lot_area', 'total_rooms']\n", "y = homes['sale_price']\n", "\n", "def rmse_k(actual, predicted):\n", " return np.sqrt(np.mean((actual - predicted) ** 2)) / 1000\n", "\n", "valid_scores = {}\n", "models = {}\n", "for k in [1, 2, 4, 7]:\n", " features = numeric[:k]\n", " models[k] = LinearRegression().fit(homes.loc[train_rows, features], y[train_rows])\n", " valid_scores[k] = rmse_k(y[valid_rows], models[k].predict(homes.loc[valid_rows, features]))\n", " print(k, 'features: validation RMSE', round(valid_scores[k], 1))\n", "\n", "best_k = min(valid_scores, key=valid_scores.get)\n", "test_pred = models[best_k].predict(homes.loc[test_rows, numeric[:best_k]])\n", "print('Chosen:', best_k, 'features')\n", "print('Test RMSE:', round(rmse_k(y[test_rows], test_pred), 1))\n", "# The test set was used only once, after the choice was made." ] }, { "cell_type": "markdown", "id": "523707e7", "metadata": {}, "source": [ "## Interpretability\n", "\n", "```{index} interpretability\n", "```\n", "\n", "Accuracy is not the only thing a business needs from a model. A model that a manager can explain is one the manager can check, trust, and defend. A model's **interpretability** is how easily people can understand why it gives the predictions it does.\n", "\n", "Compare two models for the broker, both fitted on the 300 homes used for training and validation and judged on the test set. One uses three features: living area, overall quality, and year built. The other is the 35-column model, with the seven numeric features and 28 neighborhood indicators." ] }, { "cell_type": "code", "execution_count": 18, "id": "58ecb964", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.654671Z", "iopub.status.busy": "2026-10-24T05:29:30.654616Z", "iopub.status.idle": "2026-10-24T05:29:30.660011Z", "shell.execute_reply": "2026-10-24T05:29:30.659683Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3-feature model: test RMSE 33.0 ($1,000s)\n", "35-column model: test RMSE 30.5 ($1,000s)\n" ] } ], "source": [ "small = ['living_area', 'overall_quality', 'year_built']\n", "big = list(X_all.columns[:n_numeric + n_nbhd])\n", "\n", "simple_model = LinearRegression().fit(homes.loc[pool, small], y[pool])\n", "big_model = LinearRegression().fit(X_all.loc[pool, big], y[pool])\n", "\n", "simple_test = rmse(y[test_rows3], simple_model.predict(homes.loc[test_rows3, small]))\n", "big_test = rmse(y[test_rows3], big_model.predict(X_all.loc[test_rows3, big]))\n", "\n", "print('3-feature model: test RMSE', round(simple_test / 1000, 1), '($1,000s)')\n", "print('35-column model: test RMSE', round(big_test / 1000, 1), '($1,000s)')" ] }, { "cell_type": "markdown", "id": "cfdd8018", "metadata": {}, "source": [ "The larger model is more accurate on the test set, with an RMSE of about \\$30,500 against \\$33,000: a gain of about \\$2,500 on homes that sell for about \\$175,000 on average, or roughly 1.4 percent. Is that worth having 35 coefficients instead of three?\n", "\n", "It depends on the use. If an automated tool prices thousands of listings and only accuracy matters, the larger model may win. If agents need to explain a price to a seller, the three-feature model has a real advantage, because its whole logic fits in one sentence. Here are its coefficients." ] }, { "cell_type": "code", "execution_count": 19, "id": "a84ef381", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.660906Z", "iopub.status.busy": "2026-10-24T05:29:30.660847Z", "iopub.status.idle": "2026-10-24T05:29:30.662327Z", "shell.execute_reply": "2026-10-24T05:29:30.661987Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "living_area 64.9\n", "overall_quality 24275.8\n", "year_built 397.6\n", "intercept -851224\n" ] } ], "source": [ "for name, value in zip(small, simple_model.coef_):\n", " print(f'{name:16s} {value:10.1f}')\n", "print('intercept ', round(simple_model.intercept_))" ] }, { "cell_type": "markdown", "id": "82189cd4", "metadata": {}, "source": [ "An agent can say: \"Each extra square foot adds about \\$65, each step of overall quality adds about \\$24,300, and each year of newer construction adds about \\$400, other things held fixed.\" Each statement can be checked against what the agent sees in the home. Reading 35 coefficients in the same way would be much harder.\n", "\n", "### Report the Uncertainty\n", "\n", "```{index} coefficient; uncertainty\n", "```\n", "\n", "A coefficient from 300 homes is an estimate, so it comes with sampling variability, just like the slope in {doc}`../12-regression/1205-slope-inference`. We use the same remedy: bootstrap the homes, refit the model each time, and take the middle 95 percent of the coefficients." ] }, { "cell_type": "code", "execution_count": 20, "id": "29e3e84b", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:30.663129Z", "iopub.status.busy": "2026-10-24T05:29:30.663077Z", "iopub.status.idle": "2026-10-24T05:29:31.238863Z", "shell.execute_reply": "2026-10-24T05:29:31.238462Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
estimate95% low95% high
living_area64.954.575.2
overall_quality24275.819412.928990.8
year_built397.6229.6568.9
\n", "
" ], "text/plain": [ " estimate 95% low 95% high\n", "living_area 64.9 54.5 75.2\n", "overall_quality 24275.8 19412.9 28990.8\n", "year_built 397.6 229.6 568.9" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "boot_rng = np.random.default_rng(13)\n", "pool_homes = homes.loc[pool]\n", "boot_coefs = []\n", "for i in range(1000):\n", " resample = pool_homes.sample(len(pool_homes), replace=True, random_state=boot_rng.integers(1_000_000))\n", " boot_model = LinearRegression().fit(resample[small], resample['sale_price'])\n", " boot_coefs.append(boot_model.coef_)\n", "boot_coefs = np.array(boot_coefs)\n", "\n", "low, high = np.percentile(boot_coefs, [2.5, 97.5], axis=0)\n", "interval_table = pd.DataFrame({'estimate': simple_model.coef_.round(1),\n", " '95% low': low.round(1),\n", " '95% high': high.round(1)}, index=small)\n", "interval_table" ] }, { "cell_type": "markdown", "id": "49dc2503", "metadata": {}, "source": [ "The interval for living area is about \\$54 to \\$75 per square foot, and for year built about \\$230 to \\$570 per year. A report that says only \"\\$65 per square foot\" suggests more precision than the data support. Giving the interval lets the manager see which coefficients are well pinned down and which are not.\n", "\n", "These intervals assume the model is a reasonable description of the data and the homes are like a random sample of the market, the same cautions as in {doc}`../12-regression/1205-slope-inference`. They also describe association: in these records, bigger and better homes sell for more, not that adding a square foot to a particular home would add \\$65." ] }, { "cell_type": "markdown", "id": "81a9e8cc", "metadata": {}, "source": [ "## Section Practice\n", "\n", "```{index} model choice\n", "```\n", "\n", "Four models for the brokerage were fitted on a training set and then scored on a validation set and a test set. The RMSEs are in thousands of dollars. Model A uses 3 features, model B uses 5, model C uses 7 features plus neighborhood, and model D adds all squares and interactions (63 columns).\n", "\n", "| Model | Training RMSE | Validation RMSE | Test RMSE |\n", "|---|---|---|---|\n", "| A (3 features) | 34.0 | 35.2 | 35.6 |\n", "| B (5 features) | 30.5 | 31.0 | 31.9 |\n", "| C (35 columns) | 22.0 | 31.8 | 33.5 |\n", "| D (63 columns) | 9.0 | 48.0 | 52.0 |\n", "\n", "Which model would you recommend? Model D has by far the lowest training RMSE, but a validation RMSE of 48.0 and a test RMSE of 52.0 show it has overfit. Model C also has a large gap between training and validation. The choice should rest on the validation column: B has the lowest validation RMSE (31.0), so B is chosen, and its test RMSE (31.9) is the number to report. B is also simpler than C, which makes it easier to explain." ] }, { "cell_type": "code", "execution_count": 21, "id": "8a4c6954", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:31.240035Z", "iopub.status.busy": "2026-10-24T05:29:31.239962Z", "iopub.status.idle": "2026-10-24T05:29:31.241870Z", "shell.execute_reply": "2026-10-24T05:29:31.241557Z" }, "tags": [ "thebe-interactive" ] }, "outputs": [], "source": [ "### Exercise: Choose a Model\n", "# 1. Build a table (DataFrame) with columns train, valid, and test, with one row per model A, B, C, D, using the RMSEs in the table above (in $1,000s).\n", "# 2. Find the model with the lowest validation RMSE and print its name.\n", "# 3. Print that model's test RMSE.\n", "# 4. Compute each model's gap (validation minus training) and print the model with the largest gap.\n", "# Print exactly:\n", "# Chosen model: B\n", "# Test RMSE: 31.9\n", "# Largest gap: D\n", "path_data = '../../data/'\n", "import numpy as np\n", "import pandas as pd\n", "from sklearn.linear_model import LinearRegression\n", "from sklearn.model_selection import train_test_split\n", "\n", "\n", "### Your code starts here.\n", "\n", "\n", "\n", "\n", "### Your code ends here." ] }, { "cell_type": "code", "execution_count": 22, "id": "fe5844c6", "metadata": { "execution": { "iopub.execute_input": "2026-10-24T05:29:31.242811Z", "iopub.status.busy": "2026-10-24T05:29:31.242754Z", "iopub.status.idle": "2026-10-24T05:29:31.245221Z", "shell.execute_reply": "2026-10-24T05:29:31.244913Z" }, "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Chosen model: B\n", "Test RMSE: 31.9\n", "Largest gap: D\n" ] } ], "source": [ "### Solution\n", "path_data = '../../data/'\n", "import numpy as np\n", "import pandas as pd\n", "from sklearn.linear_model import LinearRegression\n", "from sklearn.model_selection import train_test_split\n", "\n", "models = pd.DataFrame({\n", " 'train': [34.0, 30.5, 22.0, 9.0],\n", " 'valid': [35.2, 31.0, 31.8, 48.0],\n", " 'test': [35.6, 31.9, 33.5, 52.0],\n", "}, index=['A', 'B', 'C', 'D'])\n", "\n", "chosen = models['valid'].idxmin()\n", "print('Chosen model:', chosen)\n", "print('Test RMSE:', models.loc[chosen, 'test'])\n", "\n", "gap = models['valid'] - models['train']\n", "print('Largest gap:', gap.idxmax())\n", "# D fits its training data best but has the largest gap, so it is overfit." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.13.16" } }, "nbformat": 4, "nbformat_minor": 5 }