Files
FYS-STK4155/doc/LectureNotes/_build/html/_sources/exercisesweek36.ipynb
T
Morten Hjorth-Jensen fa98600d75 updating
2025-08-26 08:05:50 +02:00

384 lines
42 KiB
Plaintext

{
"cells": [
{
"cell_type": "markdown",
"id": "9dd8f8d0",
"metadata": {},
"source": [
"# Exercises week 36\n",
"\n",
"## Deriving and Implementing Ridge Regression"
]
},
{
"cell_type": "markdown",
"id": "e5cc5739",
"metadata": {},
"source": [
"## Learning goals\n",
"\n",
"After completing these exercises, you will know how to\n",
"- Take more derivatives of simple products between vectors and matrices\n",
"- Implement Ridge regression using the analytical expressions\n",
"- Scale data appropriately for linear regression\n",
"- Evaluate a model across two different hyperparameters"
]
},
{
"cell_type": "markdown",
"id": "427d9b9d",
"metadata": {},
"source": [
"## Exercise 1 - Choice of model and degrees of freedom"
]
},
{
"cell_type": "markdown",
"id": "ff6ebdb5",
"metadata": {},
"source": [
"**a)** How many degrees of freedom does an OLS model fit to the features $x, x^2, x^3$ and the intercept have?\n",
"\n",
"**b)** Why is it bad for a model to have too many degrees of freedom?\n",
"\n",
"**c)** Why is it bad for a model to have too few degrees of freedom?\n",
"\n",
"**d)** Read [chapter 3.4.1 of Hastie et al.'s book](https://link.springer.com/book/10.1007/978-0-387-84858-7). What is the expression for the effective degrees of freedom of the ridge regression fit?\n",
"\n",
"**e)** Why might we want to use Ridge regression instead of OLS?\n",
"\n",
"**f)** Why migth we want to use OLS instead of Ridge regression?"
]
},
{
"cell_type": "markdown",
"id": "a78b2ce2",
"metadata": {},
"source": [
"## Exercise 2 - Deriving the expression for Ridge Regression\n",
"\n",
"The aim here is to derive the expression for the optimal parameters using Ridge regression.\n",
"\n",
"The expression for the standard Mean Squared Error (MSE) which we used to define our cost function and the equations for the ordinary least squares (OLS) method, was given by the optimization problem"
]
},
{
"cell_type": "markdown",
"id": "93a7ad16",
"metadata": {},
"source": [
"$$\n",
"{\\displaystyle \\min_{\\boldsymbol{\\beta}\\in {\\mathbb{R}}^{p}}}\\frac{1}{n}\\left\\{\\left(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\right)^T\\left(\\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\right)\\right\\}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "6b182c01",
"metadata": {},
"source": [
"By minimizing the above equation with respect to the parameters $\\boldsymbol{\\beta}$ we could then obtain an analytical expression for the parameters $\\boldsymbol{\\hat\\beta_{OLS}}$.\n",
"\n",
"We can add a regularization parameter $\\lambda$ by\n",
"defining a new cost function to be optimized, that is"
]
},
{
"cell_type": "markdown",
"id": "fb065daf",
"metadata": {},
"source": [
"$$\n",
"{\\displaystyle \\min_{\\boldsymbol{\\beta}\\in\n",
"{\\mathbb{R}}^{p}}}\\frac{1}{n}\\vert\\vert \\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\vert\\vert_2^2+\\lambda\\vert\\vert \\boldsymbol{\\beta}\\vert\\vert_2^2\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "89ae11e8",
"metadata": {},
"source": [
"which leads to the Ridge regression minimization problem. (One can require as part of the optimization problem that $\\vert\\vert \\boldsymbol{\\beta}\\vert\\vert_2^2\\le t$, where $t$ is a finite number larger than zero. We will not implement that in this course.)"
]
},
{
"cell_type": "markdown",
"id": "456d58d4",
"metadata": {},
"source": [
"### a) Expression for Ridge regression\n",
"\n",
"Show that the optimal parameters"
]
},
{
"cell_type": "markdown",
"id": "a5677efd",
"metadata": {},
"source": [
"$$\n",
"\\hat{\\boldsymbol{\\beta}}_{\\mathrm{Ridge}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}+\\lambda\\boldsymbol{I}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y},\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "e7b08048",
"metadata": {},
"source": [
"with $\\boldsymbol{I}$ being a $p\\times p$ identity matrix.\n",
"\n",
"The ordinary least squares result is"
]
},
{
"cell_type": "markdown",
"id": "8bb0e6af",
"metadata": {},
"source": [
"$$\n",
"\\hat{\\boldsymbol{\\beta}}_{\\mathrm{OLS}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y},\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "2faaae78",
"metadata": {},
"source": [
"## Exercise 3 - Scaling data"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "8ed79466",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn.preprocessing import StandardScaler"
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "59d6f347",
"metadata": {},
"outputs": [],
"source": [
"n = 100\n",
"x = np.linspace(-3, 3, n)\n",
"y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2) + np.random.normal(0, 0.1)"
]
},
{
"cell_type": "markdown",
"id": "c9fc9fe0",
"metadata": {},
"source": [
"**a)** Adapt your function from last week to only include the intercept column if the boolean argument `intercept` is set to true."
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "5aa9c86d",
"metadata": {},
"outputs": [],
"source": [
"def polynomial_features(x, p, intercept=False):\n",
" n = len(x)\n",
" X = np.zeros((n, p + 1))\n",
" #X[:, 0] = ...\n",
" #X[:, 1] = ...\n",
" #X[:, 2] = ...\n",
" # could this be a loop?"
]
},
{
"cell_type": "code",
"execution_count": 41,
"id": "2ebf249a",
"metadata": {},
"outputs": [],
"source": [
"def polynomial_features(x, p, intercept=False):\n",
" n = len(x)\n",
" X = np.zeros((n, p))\n",
" X[:, 0] = x[:]\n",
" X[:, 1] = x**2\n",
" X[:, 2] = x**3\n",
" return X"
]
},
{
"cell_type": "markdown",
"id": "c5924634",
"metadata": {},
"source": [
"**b)** Split your data into training and test data(80/20 split)"
]
},
{
"cell_type": "code",
"execution_count": 42,
"id": "676541e7",
"metadata": {},
"outputs": [],
"source": [
"X = polynomial_features(x, 3)"
]
},
{
"cell_type": "code",
"execution_count": 78,
"id": "f895606f",
"metadata": {},
"outputs": [],
"source": [
"X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n",
"x_train = X_train[:, 0] # These are used for plotting\n",
"x_test = X_test[:, 0] # These are used for plotting"
]
},
{
"cell_type": "markdown",
"id": "a58f05ad",
"metadata": {},
"source": [
"**c)** Scale your design matrix with the sklearn standard scaler, though based on the mean and standard deviation of the training data only."
]
},
{
"cell_type": "code",
"execution_count": 79,
"id": "0c07995a",
"metadata": {},
"outputs": [],
"source": [
"scaler = StandardScaler()\n",
"scaler.fit(X_train)\n",
"X_train_s = scaler.transform(X_train)\n",
"X_test_s = scaler.transform(X_test)\n",
"y_offset = np.mean(y_train)"
]
},
{
"cell_type": "markdown",
"id": "63b86003",
"metadata": {},
"source": [
"## Exercise 4 - Implementing Ridge Regression"
]
},
{
"cell_type": "markdown",
"id": "a2f52ef2",
"metadata": {},
"source": [
"**a)** Implement a function for computing the optimal Ridge parameters using the expression from **2a)**."
]
},
{
"cell_type": "code",
"execution_count": 80,
"id": "69f19bde",
"metadata": {},
"outputs": [],
"source": [
"def Ridge_parameters(X, y):\n",
" # Assumes X is scaled and has no intercept column\n",
" return np.linalg.inv(X.T @ X) @ X.T @ y\n",
"\n",
"beta = Ridge_parameters(X_train_s, y_train)"
]
},
{
"cell_type": "markdown",
"id": "626c81ba",
"metadata": {},
"source": [
"**b)** Fit a model to the data, and plot the prediction using both the training and test x-values extracted before scaling, and the y_offset."
]
},
{
"cell_type": "code",
"execution_count": 82,
"id": "f9ba0655",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.PathCollection at 0x113e21950>"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXxtJREFUeJzt3Qd409X6B/Bv0r2hm0JpS8tGyixbQVmKuK8IKIjryhUX/3sv4ICLC/V6Fa8gXvdEwAWoyF4yy1Q2FAqFlm66IW2T/J9z0pSOtE1Kdr6f58nTX9JfmtNQkjfnvO97FFqtVgsiIiIiG1Ha6oGJiIiIBAYjREREZFMMRoiIiMimGIwQERGRTTEYISIiIptiMEJEREQ2xWCEiIiIbIrBCBEREdmUOxyARqNBRkYGAgICoFAobD0cIiIiMoLoq1pcXIyoqCgolUrHDkZEIBIdHW3rYRAREVEznD9/Hm3atHHsYETMiOh/mcDAQFsPh4iIiIxQVFQkJxP07+MOHYzol2ZEIMJghIiIyLE0lWLBBFYiIiKyKQYjREREZFMMRoiIiMimGIwQERGRTTEYISIiIptiMEJEREQ2xWCEiIiIbIrBCBEREdkUgxEiIiKyKQYjREREZFMMRoiIiMimGIwQERGRTTnERnlERESNUVWqcSCtAMmp+fByVyIxugWuax0EPy++zTkC/isREZFDyi6+gp/2p2P76TzsSc3H5Qp1re8rFUCHiAD0immJx6+PR9sQX5uNlcy8TLN161aMHTsWUVFRckvg5cuXN3kflUqF559/HjExMfDy8kJsbCw+/fRTUx+aiIhIWnMkEyPf2Yp5vx3H1pM5MhAJ9ffErd1b4eZukWgV5A2NFjieWYzFu9Nw87tbsXRPGrRara2HTuaYGSktLUViYiIeeugh3HXXXUbd595770VWVhY++eQTJCQk4OLFi9BoNKY+NBERubjL5Wq8/OtRGWAInSIDcE/vNhjcPhQdIwLkh2S9rKIr+ON8AT7eliqXb2b8cAjrjmZh3l3dERbgZcPfgupSaK8hTBT/6D/99BPuuOOOBs9ZvXo17rvvPpw5cwbBwcHNepyioiIEBQWhsLAQgYGBzR0uERE5sCMZhXjq2wM4nVMqr//1+nb4v5Ed4ene+CS/WqPFJ9vO4K01J1Gu1iDEzxNv3tMdN3WOsNLIXVeRke/fFq+mWblyJfr06YM333wTrVu3RocOHfD3v/8dly9fbnRZR/wCNS9EROS6tp3KxZ0Ld8hAJDzAC18/3A+zbuncZCAiuCkVeOz6eKyYNkjOpOSVluOxr/Zhy8kcq4ydmmbxYETMiGzbtg2HDx+Wsyjz58/H999/j7/97W8N3mfevHkyktJfoqOjLT1MIiKyU8cuFuHxr/fJWY2hHcOw+pnr5bKMqTq3CpQBye09ouRsyRPf7MeJzGKLjJnsLBgRuSFiOeebb75BUlISbrnlFrz99tv44osvGpwdmTVrlpzS0V/Onz9v6WESEZEdulh4GVM+24MSVSX6xQXjfw/0RrCfZ7N/npe7m1yiSYoLlj/zoc/3yKoccvJgpFWrVnJ5Rsxw6HXu3FlmNF+4cMHgfUTFjVhbqnkhIiLXUnSlAg9+ugeZRVfQPtwfHz7QRwYT10r8jP/d3xtxoX5IL7iMR7/chyt1yoLJyYKRQYMGISMjAyUlJdW3nTx5EkqlEm3atLH0wxMRkQMqr9Rg6tf7cCKrWFa+fDalL4J8Pcz281v6eeLTB/uiha+HrLiZvuwgNKIWmBwjGBFBxcGDB+VFSE1NlcdpaWnVSyyTJk2qPn/ChAkICQnBlClTcPToUdmn5B//+IcsDfbx8THn70JERE5i9orD2J6SBz9PN3z2YF+0aWn+hmViZkTMkHi4KbDqUCYWbTlt9scgCwUje/fuRc+ePeVFmD59ujyePXu2vC56iOgDE8Hf3x/r1q1DQUGBrKqZOHGibJr23//+19SHJiIiF7DxeBaW7DkvO6gunNgL3VpfXeY3t37tQvDKHd3k8X83nEJaXpnFHoss1GfEWthnhIjIdfJERr69VeaJPDokDs+P6WLxxxRvg/d/slvOxIhqHTETU7N5GjlBnxEiIiJjvfbrMRmIxIb4YvqIjlZ5TBF4vHR7N3i6KbH5RA5WH860yuPSVQxGiIjILvx+Kkcuzwhv3pMIH89rr5wxVnyYP/56Qzt5PPfno7Lsl6yHwQgREdmcePOf+cMheTx5QIzsA2JtTwxLQNtgXzkz8+76k1Z/fFfGYISIiGzujd+Oy54fbVr64J+jO9lkDN4ebph7e1d5/On2s7LzK1kHgxEiIrKpPWfz8dWuc/L4jbu7w8/L5A3lzWZYx3CM7hop28W/sPywTG4ly2MwQkRENiPe7F/99Zg8HtcnGoMSTN9zxtzm3NYFvp5u2HfuEjYcy7b1cFwCgxEiIrIZUbly8HyBfPP/v1EdbD0cqVWQDyYPjJXH721K4eyIFTAYISIim6hQa/DmmhPy+JEh7RAe4A178fDgOHh7KGWr+G0pubYejtNjMEJERDYhynhTc0sR4ueJx67XldXai1B/L0xIipHH721MsfVwnB6DESIisrpSVSXeXX9KHj91U3v42zBptSEiQBKN0JJT87H7TJ6th+PUGIwQEZHVffx7KnJLVIgJ8cX4pLawR5FB3vhLH93u8gs2cXbEkhiMEBGRVYkg5MOtuh1y/zGqIzzd7fet6PEb4uGmVOD3U7kyf4Qsw37/AoiIyCmJ3XFLy9VIbBOEMde1gj2LDvbFnT1by2POjlgOgxEiIrKai4WX8W1ymjyeeXNnh9gdd+rQeIhhrjuaxa6sFsJghIiIrOajramoUGvRv10wBsSHwBGITfT0MzgfbNEtL5F5MRghIiKryC8tr54V+dvQBDgSkTsi/HYoE3klKlsPx+kwGCEiIqv4fMdZXK5Qo1vrQAxpb/u276bo1joIidEtUK7W4Pt9F2w9HKfDYISIiCyuRFWJL3acrZ4VcYRckbomVpUgL05Og0bDFvHmxGCEiIgsbvHucyi8XIF2YX4Y1TUSjujWxFYI8HbHubwybD/NFvHmxGCEiIgsSlWplk3OavbtcES+nu64u5euCdo3u3S5L2QeDEaIiMiiftiXjuxiFVoFeeOOHrqeHY5qQj/dUs26Y1nIKrpi6+E4DQYjRERkMZVqDf5X1W310SHt7LrbqjE6RAQgKTYYao0WS/ect/VwnIZj/1UQEZFdW3U4U+ZYBPt54r6kaDiDif11syOiTFkEW3TtGIwQEZHFfLJNlysyeUCszLlwBqO7Rcrg6mLhFWw+kWPr4TgFBiNERGQR+9Muyc3lxNLM/VWzCc7Ay90N9/SuSmTdfc7Ww3EKDEaIiMgiPt+u6ytyW2IUQvy94EzGV/Uc2XwyB+fzy2w9HIfHYISIiMwus/AKVh26KI+nDIqFs4kL9cOAdiHQaoEVB9NtPRyHx2CEiIjM7utd51Cp0SIpLhhdo4LgjO7sqStTXn4wA1oRlVCzMRghIiKzulKhli3ThYeccFZEb/R1kTIfJiW7BEcyimw9HIfGYISIiMxq5cEMuUNv6xY+GN45As4q0NsDN3UKl8dcqrk2DEaIiMhsxHLFZ1Ub4k0aEAN3N+d+m7m9qqPsyj8yZCM0ah7n/ishIiKr2p2aj2MXi+Dj4Yb7+jpPOW9DhnUKQ6C3O7KKVNh9Js/Ww3FYDEaIiMhsPtuua3J2V6/WCPL1gLMTPUduua6VPF7OpRrrBSNbt27F2LFjERUVBYVCgeXLlxt93+3bt8Pd3R09evQw9WGJiMjOpRdcxrqjWfL4wYHOm7ja0FLNb4cyZfIuWSEYKS0tRWJiIhYuXGjS/QoKCjBp0iTcdNNNpj4kERE5gG93p0GkTQyMD0H7iAC4in5xwXJH4mJVJTafyLb1cFwjGLn55pvxyiuv4M477zTpfo8//jgmTJiAAQMGmPqQRERk58orNVhStYvt/f1j4EqUSgVu6xElj5cfyLD1cBySVXJGPvvsM5w5cwZz5swx6nyVSoWioqJaFyIisiMaNZD6O3Doe/l17eEM5Jao0MLXQy5V7DydJ6tLxEUci9JX8VUELTWvGzrHEatS7qhaqtl4PBuFlytsPRyHY/EtFE+dOoWZM2fi999/l/kixpg3bx7mzp1r6aEREZEpwce5HUBJFpB3Gtj3GVCsa/cufFn5EoAEFJRVYPqyP+RtIjARxG16SgXkUo6eoXNaB3rgnf5lSAqrBPwjgJiBgNIN9qxzq0B0jAjAiaxirD58EeNcoJLIYYIRtVotl2ZEYNGhQwej7zdr1ixMnz69+rqYGYmOjrbQKImIqFFHVwKrZwBFhpcgUjRRSK5MgBIaDFb8ia3aHvUCDL26kx51zxmlTMYc1ZeI2pp/9cbAKGDUPMA3RBcM2WmAcnvPKLy5+gRWHMxgMGJPwUhxcTH27t2LAwcOYNq0afI2jUYjm+KIWZK1a9fixhtvrHc/Ly8veSEiIusTyyTJqfnILr6CTpc2o8OWJ6BAw0snX6uHy683KvfjJY8vMFj1LjTNyAIQgcgij/n1vyGCoO8m175NBCij3wC63AZ7Mba7LhjZdSZPdqAN9vO09ZAchkWDkcDAQBw6dKjWbe+//z42btyI77//HnFxcZZ8eCIiMjH4OJtbhm+T05BZdEXOdGzzegFahRaKBu5bpvXCD+oh8vh+t/WIUuQhSXkcuzRdTBqDeKw5Hl/qjht6sBq0RReBZZOwJ2k+1B3Hyg353Iy5owVFB/uia1Sg3Kdm/dEs3NuXM/oWC0ZKSkqQkpJSfT01NRUHDx5EcHAw2rZtK5dY0tPT8eWXX0KpVKJbt2617h8eHg5vb+96txMRkfWJ/Ia5Px/FxcIr9b4ngoooRY3lEgN+Vg9AMfzQVpGF65W6D5/hKDB5HMY8Vk1ipkYs+bTZ/RIGbw1BRJAv5oztgtHddA3IbOXmbpEyGPnt8EUGI5YMRsSyy7Bhw6qv63M7Jk+ejM8//xwXL15EWpput0YiIrLnWZBSvLP+VIPnGhNU6JdoJrhtgFKhW8rJRguTx9WcAEZMhERBNxOzu7ALpn69H88M74DYUF+EB3jbZLZkdLdIvLX2JLan5KHoSoXcTI8sEIwMHTpU5nw0RAQkjfnXv/4lL0REZD+zIIY0FVT8oWmHQ9p28EQF/uK2Rc5UZCIEyZpOJo+tOQFMzUBG/670zvqT1beLRmTWni1JCA9AfJgfTueUYtPx7OrurNQ47k1DROTE9D08Xv75CB7/er/RgYgggooMbXC9Chi9b9S6jtq3KHejJYrl8dyKB5qVvNrUYzUnkMksvCJnS0QQZu3ZEWHNkUyrPq4jYzBCROSkxJvw4Dc2YvxHu/DJ9rMm318EFXMrJumO6wQJRVofmS8iTHDfAJVvJJ7z+CfWaJKqz2np61HdR0Sv7qqJ/pxaj2Xs+LRAhrbhmRht1WXmD4ewPSXXas3Ubq6aidl0PAeXy7lXjV00PSMiIvvLBzGWCC6mVjwjK12icDXBdIXHGFxWeSMhUIO+970FRewgvAolbq96fH3OhpBc47beMS2x79ylBs7pgVOXOqDDgVca7Gmip48rjJmJKbhcgYkf77baso2oqGndwkduHLj1VA5GddXNlFDDFNrGEkDshGh6FhQUhMLCQlkuTERE15YPYgpFVenttyPVsiuq1i8cN68EjmcWY/atXfDQ4DjLdXsVDc5K84C1s2oFKGJGRAQiNWdijPk9hEX397J4QPLyL0fxybZU3NWzNd4e57o71RcZ+f7NmREiIicIRERuhKU+WUZWzSgkVb2BH0y7hOOZO+DprsRdvSyQoCk6q8bpepdU6zJWBih/HDuOedsKsEfTCWoTMw3E8+MGDVauWIaRmigoAyIt1slVlPiKYGTdsSy5H494rqhhDEaIiBx4WWbX6TyZE2HOQCQy0Avjk9oiNtTPYIns4t269g23XtcKLXyt1GW0KkBJjBuCB9texLlmzALJVvNiuakiH/jRsp1ce7VtibAAL+QUq7DzTB5u6BBm1p/vbBiMEBG58LKMPsQwtj+H6J3x85+65ZLx/Wyz/4pYYhnRJbJefowYsdbkVvO6Tq6490uzBiRKpQIju0Tgm91pWH04k8FIExiMEBG58LKMfgnG2ByKFQfScaVCg/bh/ugT0xK2IoKlAfEh1dc7RgY0GJw13mpeC9HsvvyXf2K1qgfCA/3M1ixNVNWIYGTd0Uy8ckc3m7ert2cMRoiIHEFVUqemOBMrf86AAu2gbUZ3Bv3swUODYuXsgilvvKLeQby5CmIZR6GwnzdX/WyJWLZ6YvF+WUFjbKt50Vreq+wivv1uqdxTx1xVN/3aBSPIxwO5JeXYezYf/dpdDZ6oNmbUEBHZu6MrgfndgC9uhfLHR/B+xWxs83pKLj00Zybkg/t7YfbYrnJmwZRP6wfPF8gKGi93Je7u1Qb2Rvwug9qH4vW7r5NBl8LEVvP688zVLM3DTYmbOoXL4w3Hs6/pZzk7BiNERPYeiIichjp9NyKRL3MgGgtI9G/Izw7vgHfv64FvH+2PbTNubPYnfn3i6pjurRBUp5mZPRG/nyjfFYGXKa3m9efpl7/Ess+1Nkq7qXOE/LrhWNY1/Rxnx2UaIiJ7bV5WVIrR6/4BT5nVUJuY0BDvk3M8vsI6VR+Djb9MzQdpTM3E1QlJtklcbXaSa9F1UK37GJ5lWXJJpq66e+qIHBOxtBNeXIDjO6+g64DRzS7/vb5DKNyVCrlXjUi0FRVKVB+DESIiO62S6a88its9M43atVbkOui18PHAwom90L+dacswjVleI3FVdFF1BFeTXEMAr3/rZpjq1NzU7eRaXf6rzzFZtwDY3fzy3wBvD5k7InbxFUs1D5u7QZyT4DINEZGdVcnoK0JMzXXQL8uInIlBCaFmC0Rk4uou3RLNxH72lbhqNBFIiPLdwNqzRGJGRLS7F51c9eW/YgnMYPmvWDJrhhs76ZZqNh7nUk1DODNCRGQnSzNiRqTmIoKpuQ7mXJapaX/aJZzIKoa3hxJ32mHiqkkBSacx1VVJ037OwJridrKTa1PlvzLMWz1Td38Tl2xEEqtoD7/7TL5c7gr0tt98G1thMEJEZAeByOfbU+v1yBA5DBnaYPlJ3dAkh+yP4RuJ8bePw9Nm7I9Rl76c99buUbJU1aFVdXIVywK3KS/it6/3y9mkpsp/ZUBSlI5tG1bCrd31Jj3XIk8kPsxP5o38fjJXJgBTbVymISKy8dLM4Dc24uVfj9X7nshhmFsxqVZuw1UK+SbqdeubuL1nW5PLdI1VWFaBX//UlbhOsFHHVWtU3Ri7JLZs016M/2iX/DczpfS3uqqGSzUGMRghIrKTHBFDRC6DyGnIRHDtb4g9VczcwtyQH/ZfgKpSg06RAegZbdyykaMFJKLc+a9jBpq0JGZqL5Ibq/qNbD6Rc83lws6IyzRERDYo280svCxnQ4x5WxIByXpVH4wKOIMFYy2722zdxNXFyQ6euGoEMaMky3dF1YxIVjWi/Lcqi0Tm+YgS4qZmpUTr/EBvd+SXluPg+UvoHVMnuHRxDEaIiOx8cztF1ZLNbbffC6WZk1Mbs+fsJaRkl8DHww2392wNpyYCO1G+a0T5r564WfxbiuCy5j45hri7KTG0YzhW/pGBDceyGYzUwWUaIiI7WZJpiMhpELkN5q6Sacri3efk19sSo1yjAsSI8l9DtqcYt/RyU+eq1vDH2Bq+Ls6MEBHZoGzXWC+O6YwHB8VZfcfXS6XlWHU40ykTV40t/z15OgWzN+bKpRlDHW71Fmw6jR/2pzdZVn1DhzD57yjKpM/nlyE62NdCv4Tj4cwIEZENynabIkIPsXusLQIRfeJqeaUGXaMC0b1NEFxKVflv/I0P4lxAL6N2RzYmobWFr2d199pNJzg7UhODESIiG5TtNkYfeohP2rYIRDQabXVvkQlOnLjaFPHci38DoalnwNjN9fS7+K7nUk0tDEaIiCzAEXNE9LafzkVqbikCvNxxRw8nT1w1cQfgxtRMaG0qb2TX6TyUlVeadayOjDkjRER2kCMS7OeBF2/tishAb4t1UjXWVzt1iat3924DPy++Teh3AH5n3Uks2JTS4Hn63X49jqUDys4Gy6/jw/zRuoUP0gsuY9eZvOp9a1wd/8qIiMzcQ0RUVxg7I6IPOV678zqbzYTUJN4k1x/TdQm9v78LJa42QQSHYvPBhoKRWrv97oXuElh/t1+x5DW0Y5hcBttyIofBSBUGI0RENuwhYqnN7Zrr291psq/GwPgQJIQH2Ho4dkXMWImkYpGsWnPWS7/bbz363X7rdMoVVTUyGDmZY52BOwDmjBAR2Sg/RJTtilbk9hKIqCrVWLJHl7j6QP8YWw/HIRJam97tV/yBzAQ06upbByaEwsNNgbN5ZTibW2qVsds7BiNERFbOD7F12W5DVh/ORG5JOSICvTCiC5cPjElo1e/22/A/o263X9G3RM/fyx19qjqwcnZEh8s0RERW7CFi67JdYxJXJyTFyPbl1HhCq8gPksmqIj+kKSW1d+u9oWMYdp7Jw+YT2Zg8MBaujn9tRERW7CFi67LdhhzNKMLec5fgrlRgfFK0rYdj90QgKfaj6dO1s1HnHynyqdV/RCSxCiIguVJxdQnHVZkcjGzduhVjx45FVFSUzApevnx5o+f/+OOPGDFiBMLCwhAYGIgBAwZgzZo11zJmIiKHyxGZNiwB3z7a365yRGr6apduVmRUt0iEBzbdU4OqiPJdUTXTQFs0EX9kaEMw9meNDF71HVo7RgTI5bArFZpG+5K4CpODkdLSUiQmJmLhwoVGBy8iGFm1ahX27duHYcOGyWDmwIEDzRkvEZFD5Yjo80OeHdFBfpK2t6UZoehKBZYfSJfHTFxt5m6/Uu1/27q7/dZsGS8+zIuqGmEL80ZMzxm5+eab5cVY8+fXLnd67bXXsGLFCvz888/o2bOnqQ9PRGRT4lOsqT1E7DE/pKZle87jcoUaHSL80S+OW9s3e7ff1TOAooxau/2KQES/26+ITcRfgQhmRc7J0I7hWLb3ggxGXoRrs3oCq0ajQXFxMYKDG/6DV6lU8qJXVFRkpdERETVAlGae2wGPY8fQX1nY5E6u9thDxJBKtQafbT8rj6cMinPZfWjMtdvvkZ2r8b9fdyAbLQz+jdRsGS+aqIkgNSW7BBculaFNS9fdxdfqwchbb72FkpIS3HvvvQ2eM2/ePMydO9eq4yIiatDRldWfevsAWOIp8gCCMbdiUvWnXkM9ROytdNeQNUeyZNfVYD9P3NnTtfehuWZKN6T49cBKTdOnZhdfwQCfEPRq2wJ7zl6SsyMT+7nuEplVq2kWL14sg4xly5YhPFy3WZAhs2bNQmFhYfXl/Pnz1hwmEVHtQER00awx/S5EIl923RTdNx2hh0hDPtl2Rn69v19beHvU3keFTBceYFzy76msYuw8nYch7XV5I5tPuHbeiNVmRpYsWYJHHnkE3333HYYPH97ouV5eXvJCRGTzpRkxI2IgXVXEGSJBcY7HV1in6iOn4x0lR0Rvf9ol7E8rgKebEvcPcN1P5dZoGV/Xgk2n5SXU31Ne35GSi/JKDTzdXbPjhlV+62+//RZTpkyRX8eMGWONhyQiunaia2adGZGaRLwRpciTXTjtuYdIQz7Zliq/3tYjyuhP9GR6y/jGiI63Qmm5GnvPuW6Jr8kzIyLfIyXl6q6FqampOHjwoExIbdu2rVxiSU9Px5dfflm9NDN58mS8++676NevHzIzM+XtPj4+CAoKMufvQkRk1hLe06dT0MGIc/8+MAgVnfvLT8WOMCMiiIRJ0f5deGhQnK2H45Qt403dOHHziRwMjA+FKzJ5ZmTv3r2yJFdfljt9+nR5PHv2bHn94sWLSEvTbbQkfPjhh6isrMQTTzyBVq1aVV+efvppc/4eRERm7646e2OuUeeLLpz22kOkIV/sOCsDLrE7b5eoQFsPxykDEtHgTjS6Ew3vjLH2iC44dEUmz4wMHToUWm3DK2Gff/55reubN29u3siIiGzYXVW8ymWhk6yaEcmqhuMMha77pujC6UBKVJVYkqwrDHhkCGdFLN0yXlTOGONsXhlyS1QI9Xe9nEnXzJQhIjKiu6pIShXlu/K43mewquhk9Ou6LpwO5Lu951GsqkS7MD8M7dBwZSOZhyn5OL+7aDdWBiNE5PJEECLKLN9Zd7LeGr/oIzK14hlkok6jRjEjIrpuimZXDqRCralOXBVNzpQOtLTk6BU2xjzTz/10qHr/Gldi9aZnRET2RLzwN5VoKAISUb4rqmbCUYB7h/XB4Jtuc7gZEWHFwQxcuHRZlpTe06uNrYfjUhU2YvlPBCSNlfxertDg8a/34wMHqsoyB86MEJHLMmX3XbFks0vTBSs1A+HW7nqHDETEDND7m3TVkI8MaQcfT8f7HRy9wkaUfxtj7s9H5b+Xq2AwQkQuydTdd2t2VxXT7o7o10MXcSa3FC18PXA/d+e1WYWN2CqgKRer9q9xFQxGiMglmbL7ruBo3VXr0mi0WLgxpbqviL8XV+ltQfzthAYYVy2TbWQVjjPgXyMR2QXxwns6uxRl5ZUoK1fjcrkaqko12ob4IbFNEFr46tpmm2NGRAQiv5mYJOgIO/A2Zt2xLJzIKkaAlzsmD4y19XBcmrHVNeEu1BWXwQgR2UxaXhnWHMnE6iOZcp+URloYIS5UF5T0iQ2W7csDvT0skqxa17Rh8RiUEOZQ3VXrEr2h3tt4Sh6LQCTIx/Tnjqy7f02wn4fDLgc2B4MRIrL6G+Pao1l4d/0pHL1YVC/gCPB2h4+HG3w93eDupkRKdglSc0urL8sPZmDeqmP4S59oPDgwFrGhfiY3MzOGomo25NkRHR02CNHbfDIHh9OL5HP60GA2OXOE6pq4UD/88meGnB1x5EDYWAxGiMiqMyFzVh7Gpqrt0sULbL+4YIzuFomRXSIbrDQoKCvHHxcKcTCtQL5An8ouwec7zuKLnWdxU6dwTLuxPXpEtzBbsqqj54fUmxXZoJsVEUmrwX7mWe4iy+xf4+vpJpcp950rwL5zB+VtrRx8idAYCm1jvd3tRFFRkdxUr7CwEIGB3EOByNFcqVDjgy2n8f7m03KbdA83BR67vh0eGdwOLU18cxQvWdtScvHpttTqoEahACb2a4t/ju5kcPlGNDQb/9Euox/DmV78t53Kxf2f7IaXuxK/zxjmUnkIjkCfw5RdfAVnc8vwzvqT9c7Rh8OOtCO0qe/fnBkhIos6n1+GKZ/vkcstwuCEUMy9vSviw/yb9fMUCgWGtA+Tl9M5JbJC5McD6fh6VxrWHMnCv8Z2xS3XRcrz9C/221OM2/Bu0oAY3NytldNMi4sKmnm/HZPHE/q1ZSBix/vXqDVauTmjIWLGQPw1ilmUEV0ineJvsy4GI0RkMUcyCvHgZ3uQU6xCeIAXZo/tgjHXtaoOFK6VCGjeHtcD9/Rpgxd+Oix7aDyxeD+GdQzDG3d3l0mxpiSsikBEvDE4ixV/pONIRpGsoHnyxva2Hg5dQ6m5tkbvEWf6G9VjMEJEFiFmI/761T65Q2ynyAB88VASIgIt88l8YHwoVj09BIs2n5YXsXwz/O0tKLpSaVKyqjNVL4ilsbfW6Kb8pw6LZ66Incs2sqeIs/YeYdMzIjK7FQfT8eBnyTIQ6d8uGMseH2CxQETP28MNz47ogFVPD0aHCH+TAhFnSVatSST4phdcRlSQt2xyRvYt3MV7jzAYISKzWpKchqeXHESFWosx3VvJGZHm9ARproTwAMy6uel223piRsQREwMbk19aXt1t9f9GdpSBGjn2zr4KACF+nsgsvCwTsp1t3xou0xCR2fx+KgfPLz8sj0UPkNm3drHJFvVFVyqMbmjmDH1E6vrvhlMoVlWiS6tA3Nmzta2HQ2boPaIFkFdajmeX/eF0FV8CZ0aIyCxSsovxt2/2y09sd/VqLV8orR2IiMcWnxpPZekqd5oiOqs6WyByNrcUX+86J4+fH9PZJsEgWX5n38zCKzJwEc38nAFnRojomuWVqGT5bvGVSiTFBmPeXdeZrWLGUq3eHXn33ca8sfo4KjVaDO0YhkEJobYeDjUjIBnRJRK7z+Thr1/vk/+nXKHclzMjRHRNxGZ2omrmfP5lxIT44oMHesPL3bo5CvpW76bsOSM6vzr463c9649m4bfDmfL3mnlzJ1sPh5rJTanAwIRQJLZpuKtw3XJfR8dghIiaTXRDnfnDIew9d0nuKfPJ5L5WLyE1tdW7n6cuUBJ73Lyw/LDTJAIWXq7A88sPyeNHh7RDp0h2q3Z0MSG+LlPuy2CEiJpt2d7z+OlAuvwkt2hibySEN6+rqiWbRdVMVv320f7481+jZD6LWEX6Zncanvhmv+zJ4ejE5oFZRSq5wZoocSbHNzA+1GXKfRmMEFGzN7176eej8vgfozpicHvr5yfoWr3r9qdpSvuIANm5UgROUwbFYcH4XvB0U2L1kUzZE8XYChx7bTC3ZM95eSw6z7KU1zmM7hYJ90bWEhVOlPvEYISImhUE/N93B1FarpYJq2JZwNpEnojYy2PBptPN+vQoeqB8PqUv/L3csetMPsb9bxeyixxvurtUVYmZP/5ZvbeOM7wxkY4InK/vEOoSzfoYjBCRyT7+/Qz2nL0k8y/+c2+i1V8MTUlYbezTo0gSXPJYf4T6e+HYxSLctWgHzuQYVxZsL/695oRMHm7dwkfuWkzO5f7+MfKrW53qtIhALzwzvD1UlRqnaILG0l4iMol40/7PWt2eJ2Lju+hg45LsbJGwasynx26tg/Dj1IF44NPdOJdXhrsX7cAnD/ZFr7YtYe9EvswXO8/KY1FOLWZ5yLn0iwuBh5tCdjSeP66HzHU6m1uGb5PT8M76U9XnOXoTNM6MEJFJZbzPLj2IcrUGwzuH494+0XabsGpKq/e2Ib74YepAdG8ThEtlFZjw0S5ZJmvPMgouyyZzWi3wl95tcH2HMFsPiSzAz8sdPasCY7HXk5e7EvPXn0RmnSVFR2+CxmCEiIz23oYUHM8sluW78+7qbvXGZqaUMU4bloBtM240+pOiWKoR1TaiWdiVCg0e+2ovFu9Ogz26XK6W48stUckdkf91W1dbD4ksaEhV8zqx3UJDs4L628T3HXHJhsEIERklNbcU/9uqSxZ99Y5uCAvwsu4ANGog9Xd0zl2L/sqjUELT6Omi+6ipuSziU+hHk/rg3j5tIF7Pn/vpkOxoqrGjF3fR2+WfP/yJw+lFMigU4xXjJuc1uKpSbdup3EZnBR25CRr/gonIKC/9fESuW9/QIUyWHFrV0ZXA6hlAUQZEB40lnkCGNhhzKyZhjSap1qmKquWZ5laVeLgpZXlsZKA3/rsxBYs2n8aprGK8M64HAqy4+3BDFm05jZ//yJAln4sm9rJ6zg5Z33Wtg2RTwYZawztDEzTOjBBRkzYcy8KmEzkykU7XMExh3UBk2SQZiNQUiXws8piPUcpks5c7it9v+siOeGdcIjzdlVh/LBt3vb8D5/JKYUsij0VUzwhiaaZfuxCbjoesw91NiYHxxv9bO2ITNAYjRNQo0Z30pV90zc0eGhyHdmH+1l2aETMiBlbJZayhAOZ4fFW9ZGNswqqx7uzZBt/9dYAsozyVXYLbFmyXDcZsYcvJHDy15IBMWJ3Yr211ySe5hsHtdQnKnm6K6qDbmZqgMRghokZ9si1VlryGB3jhyRvbW/fBz+2oNyNS9wUsSpGHL2+qlMmnpiSsGisxugVWThssv4r9Xx74ZDf+veY4yisbz1kxp58OXMDDn+9BWblaVs3MGcuEVVdNYlVrr+7YW5O4LoLyd/sVw+3IDzK/SgbzzhqMbN26FWPHjkVUVJScyly+fHmT99m8eTN69eoFLy8vJCQk4PPPP2/ueInIii4WXsaCjSny+LlbOlu1j4WoCDh5WvfYTRkcqa5u9W4JEYHeWPpYf1lCK3JZF246jTvf346TWcWwtA+3nsazS/9ApUaL2xKj8PGkPnLpiFxv07w2LX3k/4snhsXLWcCaxvkfxOGWf0fS1snADw8DX9wKzO+mW+Z0ACb/RZeWliIxMRELFy406vzU1FSMGTMGw4YNw8GDB/HMM8/gkUcewZo1a5ozXiKyotdWHcflCjX6xLTE7T2irN7qffZGI5dE/CMsPSS538u//5Iok0Zb+nrgSEYRbn1vm+xGa4lqG/EzX/nlqPw3EB4eHCebXjEQcU0KhQKDq2ZHROm5mAUUs4Hv3tcDa0YVYF7lv+FzObP2nYou6vKtHCAgUWhFnVhz76xQ4KeffsIdd9zR4DkzZszAr7/+isOHD1ffdt9996GgoACrV6826nGKiooQFBSEwsJCBAZyW2wia9h3Lh93L9opOz7+PG2w7FRqzVbv2qpp521eT8lkVcOTHgogMAp45hCgtN7mcGIPmxk//CmTevXVDtNHdsDQDmFmSe49klEoNyHcXVWi+dwtnfDY9fHX/HPJsf3yZwamLT6AjhEBWPPs9bobxVKMmAFpcDnTNv9HTH3/tniIvXPnTgwfPrzWbaNGjZK3N0SlUslfoOaFiKxHfEZ5c7WuamNcn2irBSJ1W71roJTlu/K43semqjf90a9b/UU2PNAbnz7YF6/deZ3cn+dQeiGmfLYH93ywEztONz/BNa9EhVk/HsLY97bJQER023z73kQGIiQNjA+VHw5OZBVf3dSxibwqmWFSlK47z45ZPBjJzMxERETtKVRxXQQYly9fNnifefPmyUhKf4mOtn7LaSJXti0lV74Zerop8dRN7W3a6l30EZla8QwyUadCQHzau/dLoMttsAUxAzKhX1ts/ecwPHZ9Oxk47Dt3CRM+2o17P9iJz7enIi2vzKjA70RmMRZuSsHQtzbLPUdE4HVr91bY+PehuKtXG6v8PmT/gv080TVKN7uwXR/0lhi5bYGx59mIXTY9mzVrFqZPn159XQQuDEiIrEO8Ob5V1ctClI9GtfCx2mM31KxJBCTrVH2QpDyOcBTg3mF9MPim22wy7VxXiL+XTO59ZHCcDCgWJ6ch+Wy+vPzr56OID/PDsI7hiA31k/keXlWXosuVchZl++k85BSrqn+eeLMR1TKOWJ5Jljc4IUx23/39VK4sPTc6X8oKeVV2HYxERkYiK6t2RCaui7UjHx/DL3Ki6kZciMj61h7Nwh8XCuHr6Ya/DYu32vKMmBU5lVXS4DliyWaXpos8Ht+uv10EInWXbube3g2P3RCPX/7IwMbj2dh77hJO55TidE5qo/cVwYkIPsYmRuHuXm0sVhVEjm9I+1B8sOW0bA0vPjgoYgbqZglFsqrBXWuqckbEea4cjAwYMACrVq2qddu6devk7URkX0RQ8J+1ulmRKYNi5eZx1khYFXkixuzEe62t3q2hdQsf/PWGeHkRfUnEm8a2lBzkl5ZDVamR/UnERQQcfWJbyj10erVtKat1iJrSO6alDF6zi1WyEV+HiABg9Bu6qhn5P0RrF3lVFg9GSkpKkJKSUqt0V5TsBgcHo23btnKJJT09HV9++aX8/uOPP44FCxbgn//8Jx566CFs3LgRy5YtkxU2RGRfxJ4nJ7NKEOjtjseGxFu1cqYp5mr1bk1BPh4Y072VvBCZg7eHmwzGxTKNCHRlMCLypkT+VNX+TdXEjIgIRGyUV2XRYGTv3r2yZ4iePrdj8uTJspnZxYsXkZZ2ddvtuLg4GXg8++yzePfdd9GmTRt8/PHHsqKGiOxHhVqDt9edlMfiU32Qr4dVK2eaImZERCBi7g6rRI5mcEKoLhhJyZVbNEgi4Og0Rlc1I5JV/SOgjh6A5HOFyD6YLverEUGMvQbyJgcjQ4cOletUDTHUXVXc58CBA6aPjois5ru9F5CWX4ZQf088ODDWJpUzhkwbFo9BCWF2/UJKZE2D24cCvwG7zuTJJb/qRnhiKSZuyNXlz39vqfV/rJUdB/Rs5UdE8gVtwcZT8vhvQxPgZ4W278Zuc94+IsCird6JHE3nyECE+HnKvYoOni9ocPmzbrCfWXhF3i6+b28YjBARVhxMR0bhFZmwKnpnWJpYosmtUc7qbNuhE1mSUqnAwKrW8NtO6boAG7P8qb9NfF+cZ08YjBC5OPGitGjLaXn8yJA4i1d16PedefnXY42e58jboRNZaxff31NyTVr+FCGI+L44z57YZdMzIrKetUcycSanVFbQTLTwrIix1TOOWDlDZE2DRN4IgD/OF8gSclG5Zcryp7HnWQtnRohcmEhGf3+zblZk8sBYBHh72EX1jKicWXR/L7tMtCOyl3427cL85NYBO0/nmbysaW/LnwxGiFyYKA0Um7z5eLhhyqCqEkEbV8+8OKaz3B6dgQiRcUs1oqmenljWFMubCgdb/mQwQuTCxF4qwn1J0XITLksydlo4NMCLSzNERhjSPkx+FT1H9MT/HbG8KSgcaPmTwQiRixI7zO46kw8PNwUeHdLOosszYhq5sX1n7Hn6mMhe9Y8PgbtSgXN5ZbV2iBazimKZUyx3OsryJxNYiVzUos26WZG7erax2M68zrbvDJE98fdyl/saiR2if0/JwcSQmOrviYBjRJdIuTwqZiX1HVgF8eGg5m32MEvCYITIBZ3ILMb6Y9lQKETrd8vMijj7vjNE9tKNNVkEIydzMbHf1WBEEP+XRMPAxj4c2EtXVi7TELmgj38/I7/e3C0S7cL87WLfGXudPiayZ0OqSnx3nM5FpVrT4Hn23pWVMyNELianWIUVB3U7ez5ioVwR7jtDZB3d27SQPYKKrlTiz/RCuWxTV1NdWcX/PPF9saxjq/+HnBkhcjFf7TqHcrUGPdu2MPjCZQ7cd4bIOtyUCgyqbg1fuxurI3VlZTBC5EKuVKjxza5z8vhh/dbjFqmcKTbqfFbOEF27qyW+tfep0XOErqxcpiFysQ3x8krLZffG0V0jzfqzWTlDZNu8kQNpBSi+UlGvk7IjdGXlzAiRC7V+/2Rbqjx+cGAs3N3M99+/oeQ4Q1g5Q2Re0cG+iA3xRaVGK3sH1eUIXVkZjBC5CNGl8WRWCfw83TAuKdpsP5eVM0T2UeIrbDOwVOMIXVkZjBC5CP2syF/6RCPQjBviGV85k4BvH+3PfWeIrNQaviZ778rKnBEiFyASSreczJFNzqYMirVR5Yx/rQZMRGQ+A6qq0s7kluLCpTK0aelb7xxDXVl7x7SUW0OIfDJbdmRlMELkAj7drpsVGdE5AjEhfmZdosktVhl1LitniCwn0NsDPaJbyMBClPjel9TW4Hk1u7KKXK8b/r3JLjqycpmGyMkVlJXjx/3pZi/nFS9kg9/YiJd/PdboefaQHEfkSlU1vzewVGPPHVkZjBA5uWV7z0NVqUGXVoFmCwiMrZ6xl+Q4IlfKG9mWkitnLRvSVEdWQXy/sZ9hbgxGiJyYeDERHVeFSQNioBBJI1asnrGX5DgiV5DYJghBPh4ovFyBPy4UNHiePXZkZTBC5MS2nMzG+fzLcu+K23u0tmr1zItjOrNyhsiK3N2UGFzVGn7LCcPdWO21IyuDESIn9sUO3azIuL7R8PF0M8vPNPYFKjTAi0szRFZ2fQddMLK1gdbw9tqRlcEIkZNKzS2tLue9v3/MNf887jtDZP+u76DLG/njfIFMXjfEHjuyMhghclJfV+WKDOsYfs3lvPrKmfEf7cKCTacbPZfVM0S20yrIBx0i/CFyT0UiqyH22JGVwQiREyorr5RVNPrE1WvBfWeIHMv1VVU1jeWN2FtHVjY9I3JCyw9koPhKpdw8S//CZK19Z2zRMImIrrqhYxg+3pYq80bEBpkNVdEZ6sjKDqxEZBbixefLnWflscgVUV7DC4sp+84MSgi12QsZEV3VNzYY3h5KZBWp5OaYHSMD0JCaHVltics0RE5mz9lLOJ5ZDB8PN/yld7RV951hIEJke94ebujfLqS6vN8RMBghcjL6Jmd39IxCkO817M6rUSOh9CBuU+5Af+VRKKFp8FRWzhDZl+urlme3nmy6NbzDBiMLFy5EbGwsvL290a9fPyQnJzd6/vz589GxY0f4+PggOjoazz77LK5csV4zFSJXkVuiqt5TYmK/a0hcPboSmN8NXddNwH89F2CJ5yvY5vUURilr/19n5QyRfZf4Jqfmy4R2pwtGli5diunTp2POnDnYv38/EhMTMWrUKGRnG54KWrx4MWbOnCnPP3bsGD755BP5M5577jlzjJ+IahAVNBVqrdy9s1vroOYHIssmAUUZtW6ORD4WecyvDkhYOUNkv+LD/NC6hQ/K1RrsPmO9tu5WC0befvttPProo5gyZQq6dOmCDz74AL6+vvj0008Nnr9jxw4MGjQIEyZMkLMpI0eOxPjx45ucTSEi02g0WizenSaPJ/Zr28wfogZWz6ixXdZV+nhjjsdXcsmG+84Q2S+FQlE9OyKaHzpVMFJeXo59+/Zh+PDhV3+AUimv79y50+B9Bg4cKO+jDz7OnDmDVatW4ZZbbmnwcVQqFYqKimpdiKhxW07l4MIl3T40YxOjmvdDzu2oNyNSNyCJUuTh57FK7jtDZOduqApGtjpAMGJSaW9ubi7UajUiIiJq3S6uHz9+3OB9xIyIuN/gwYNlyWFlZSUef/zxRpdp5s2bh7lz55oyNCKX980u3azIPb2jZTZ9s5RkGXVa18DLV6dKiMguDUzQVbidyS3F+fwyRAf7wmWraTZv3ozXXnsN77//vswx+fHHH/Hrr7/i5ZdfbvA+s2bNQmFhYfXl/HldJ0kiMiy94DI2HtcFEhP7t232vjPbMo0MYvxrfyAhIvsT6O2B3m1byuPNdj47YtLMSGhoKNzc3JCVVfvTk7geGRlp8D4vvvgiHnjgATzyyCPy+nXXXYfS0lI89thjeP755+UyT11eXl7yQkTGWZqcJveiGNAuBPFh/ibdV1TfiC6rormZEu7Y5hWMSEV+A59UFEBgFBAz0EwjJyJLGtopDMln87HpeDYeMMOGmXYxM+Lp6YnevXtjw4YN1bdpNBp5fcCAAQbvU1ZWVi/gEAGNIJZtiOjaVKg1WLLnfLNmReruO6OBEnMrJsn8VRHc1Fa1LDP6dUDZzGUgIrKqYR3D5dcdp3NxpUINp1mmEWW9H330Eb744gtZqjt16lQ50yGqa4RJkybJZRa9sWPHYtGiRViyZAlSU1Oxbt06OVsibtcHJUTUfOuPZiG7WIVQfy+M7GJ4htKUfWfWaJIwteIZZKJO7xAxI3Lvl0CX28wzcCKyuE6RAbIX0JUKDXadyYO9MnlvmnHjxiEnJwezZ89GZmYmevTogdWrV1cntaalpdWaCXnhhRdkiZH4mp6ejrCwMBmIvPrqq+b9TYhc1DdV5bzj+raBp7vSLPvOiIBknaoPkpTH8dKNoegQn6BbmuGMCJFDUSgUGNoxHN8mp2HziRx5bI8UWgdYKxGlvUFBQTKZNTAw0NbDIbIbZ3NLMfStzRCbcm79xzCTsuVXHEzH00sONnneu/f1wO09Wl/jSInIVtYeycRjX+1D22BfbPnH0AZ38bXl+zf3piFyYOLTjjC0Q5hJgYhYosktVhl1LvedIXJsgxJC4eGmQFp+GVJzS+EUyzREZB9UlWp8t++CPB6f1LZZ1TONEZ+dRJdV7jtD5Nj8vNzRLy4E21JyselEDtqZWHFnDZwZIXJQa45kIb+0HJGB3rixU3izqmcawn1niJzL0I66bqybTxjeR87WGIwQOahvqxJX7+0bDXc3ZbOrZwzhvjNEzmVY1QcWsWleqcr+dvHlMg2RAzqTU4KdZ/JkR/ZxfaOvuXqmphfHdMaDg+I4I0LkRNqF+skEVpE3sj0lFyO7Gt8GwBo4M0LkyImrHcPlNuHGyC5uOhARQgO8GIgQORmFQoFhVUs1Im/E3jAYIXIwoovi91WJqxOMSFzV7ztzKqvYqJ/P6hki516q2Xwi2+46oHOZhsjBrDmSiUtlFbKroj4p7VorZwRWzxA5t/7tQuDtoZSvByeyitEp0n76dnFmhMjBLK7uuNp44qqxlTMCq2eInJ+3hxsGxofK403H7WuphsEIkQNJyS7B7tT8JhNXTamcEVg9Q+QahunzRo7bV4kvl2mIHDBxVfQVaRXkc82VM9OGJcjujGJphjMiRC6SN7LiCPaey8el0nK09POEPeDMCJEDJa7+sL8qcbVfW7NUzrSP8MeA+BAGIkQuok1LX3RuFQiNVlTV2M/sCIMRIgex+nAmCsoqZCnvDR0a7rjKfWeIqDEjOuteP9Yfy4K9YDBC5ICJqw3NZIik1cFvbMTLvx5r9GeJe4tqHFbOELmemzpHyK9bTuTIPa7sAYMRIgcgeoQkn82XQci9fQwnrnLfGSIyxnWtgxAe4IXScjV2ncmHPWAwQuQAFtdIXBWVL3Vx3xkiMpZSqaieHVl/1D6WahiMEDlA4uqP+9MbTVw1Zd+ZbTNuZCBC5OJGdNHljWw4lmUX3VgZjBDZuVWHLqLwsi5x9fr2hjuuct8ZIjKFaH7m4+GGjMIrOHqxCLbGYITIQRJXxyfVT1zlvjNE1NxurEPa67qxrj9q+xJfNj0jsmMns4qx99wlGYT8pU7iKvedIaJrMbxLBNYezZIlvk8Pbw9b4swIkQPMigzvHI6IwKuzGtx3hoiulUiIVyiAQ+mFuFh4GbbEYITITl0uF4mr+o6rMdW3c98ZIjKHUH8v9IxuIY83HLPtUg2XaYjs1C9/ZqDoSiXatPTBkATd2q7AfWeIyJxLNfvTCuRSzf39r37osTbOjBDZqW+qlmhEOa/oC6DHfWeIyFxGVPUb2ZGSh1JVJWyFwQiRHTqcXoiD5wvg4Va/46qxFTGsnCGipiSE+yMmxBflag1+P5ULW2EwQmTHHVdFnodY15U0aiD1d/Qr3YhbAlLgBo3B+3LfGSIylkKhwKiukegb2xJeHrYLCZgzQmRnSlSVWHFA13F1or7j6tGVwOoZQFGG/ATxPoAMr2C8VDEJqzVJ1fdl5QwRmWrWzZ1kUGJLnBkhsjPLD6TLDaziw/zQT8xuiEBk2SQZiNTUSnEJizznY5Qyufo2Vs4QkalsHYgInBkhsiNijwh94urEfjFQaDW6GREDhbwKaKGFAv9tsRSrRzyM8EA/Vs4QkUNiMEJkRw6cL8Cxi0Xwclfi7l5tgHM76s2I1A1IvMou4vYW54C4IVYdKxGRuXCZhsiOfLNLNysyNjEKQb4eQImR23sbex4RkR1iMEJkJwrKymWjM33iqui0eqTIx7g7++t6BRAROSIu0xDZie/3XYCqUoMurQKRWXgZf/tmP7IKNdjmFYxI5MNwKogCCIwCYgZaf8BERLacGVm4cCFiY2Ph7e2Nfv36ITn5aja/IQUFBXjiiSfQqlUreHl5oUOHDli1alVzx0zkdDSaq4mrPaKD8LdvDsiW7xooMbdiku6cejmsVdHJ6NcBpZuVR0xEZMNgZOnSpZg+fTrmzJmD/fv3IzExEaNGjUJ2tuFNdsrLyzFixAicPXsW33//PU6cOIGPPvoIrVu3Nsf4iZzCtpRcpOaWwt/LHRuOZ9eqnVmjScLUimeQiTpNzMSMyL1fAl1us/ZwiYjMSqEVtYQmEDMhffv2xYIFC+R1jUaD6OhoPPnkk5g5c2a98z/44AP8+9//xvHjx+Hh4dGsQRYVFSEoKAiFhYUIDAxs1s8gsmePfLFXblQlOiGuOZJp8BwlNEhSHkc4CvDXMQPRdcBozogQkV0z9v3bpJkRMcuxb98+DB8+/OoPUCrl9Z07dxq8z8qVKzFgwAC5TBMREYFu3brhtddeg1qtbvBxVCqV/AVqXoic1YVLZdh4XFcN0ztGt523IWLJZpemC1ZqBiLFrwcDESJyGiYFI7m5uTKIEEFFTeJ6ZqbhT3NnzpyRyzPifiJP5MUXX8R//vMfvPLKKw0+zrx582Qkpb+ImRciZ7V4d5rMB+kWFYiiy8btmslN8IjImVi8tFcs44SHh+PDDz9E7969MW7cODz//PNy+aYhs2bNklM6+sv58+ctPUwim1BVqvHVznPy+HBGERZsSmn0fG6CR0Rw9dLe0NBQuLm5ISurdoMlcT0yMtLgfUQFjcgVEffT69y5s5xJEcs+np6e9e4jKm7EhcjZvb7qGIpVxs2GcBM8InJWJs2MiMBBzG5s2LCh1syHuC7yQgwZNGgQUlJS5Hl6J0+elEGKoUCEyFWIpmZfV5XzGoOb4BGRszK56Zko6508eTL69OmDpKQkzJ8/H6WlpZgyZYr8/qRJk2TZrsj7EKZOnSorb55++mlZcXPq1CmZwPrUU0+Z/7chciBL95xHhbrpYrZpw+IxKCGMm+ARkdMyORgROR85OTmYPXu2XGrp0aMHVq9eXZ3UmpaWJits9ETy6Zo1a/Dss8+ie/fuMlARgcmMGWInUiLXteJgulHntY8IwID4EIuPh4jIYfqM2AL7jJCzyStRod9rG1BZv61qPd8+2p/BCBE5JIv0GSGia7f68EUMfWtzk4EIK2eIyFVwozwiKwcij3+9v8nzWDlDRK6EwQiRFatn5v581OjKGRGIsHKGiFwBgxEiK0lOzZc78TblxTGd8eCgOM6IEJHLYM4IkZVkFzcdiAihAV4MRIjIpXBmhMgKyzNiVuRUVrFR53PfGSJyNQxGiCycsCryRIxZnlFU5YqweoaIXA2DESILBiJTv94PYxr5sHqGiFwZgxEiC1bOGNtRkNUzROTKGIwQ2bByZtqwBAxKCOW+M0Tk0hiMEFlgVmR7So5R57aP8GerdyJyeQxGiGyUsCqwcoaIiMEIkc0SVlk5Q0Skw6ZnRFZOWGXlDBFRbZwZIbJiwqrAyhkiotoYjBBZsdW7qJ55dkQHzogQEdXAYITIiq3eRRkvAxEiotoYjBA1E1u9ExGZB4MRIgtXzugxYZWIyDAGI0QWbvXeigmrRESNYjBCZKHKGV9PN3xwf2/miRARNYHBCJGFKmdu6hyO6zuEWXw8RESOjk3PiExcosktVhl17p092lh8PEREzoAzI0QWqJ7x93LHDR05K0JEZAwGI0QWqJ6ZexsrZ4iIjMVlGiIzV8+M7hqBu3tHW3hURETOg8EIkRn3nfH2UOLVO6+z+JiIiJwJgxGihmjUQOrv8Dj2A/orj0IJTZN3mdgvBiH+XlYZHhGRs2DOCJEhR1cCq2cARRnoA2CJJ5ChDcbciklYo0kyeBd3pQKPXd/O6kMlInJ0nBkhMhSILJskA5GaIpGPRR7zMUqZbPBu4/pGIyLQ20qDJCJyHgxGiOouzYgZEQPpqvrimDkeX9VbshGVM4/fEG+tURIRORUGI0Q1ndtRb0akbkASpchDkvK4vO7hpotQxidFIzrY12rDJCKCqwcjCxcuRGxsLLy9vdGvXz8kJxuetq5ryZIlUCgUuOOOO5rzsESWV5Jl1Gl/HxiEGaM7okKthZe7Ek/e2N7iQyMiclYmByNLly7F9OnTMWfOHOzfvx+JiYkYNWoUsrOzG73f2bNn8fe//x1Dhgy5lvESWayXyM7TediW6WbU+b27dMKvhy7K48kDY5krQkRkzWDk7bffxqOPPoopU6agS5cu+OCDD+Dr64tPP/20wfuo1WpMnDgRc+fORbt2rDYg++uuOviNjRj/0S5M2uAuq2YaLuJVAIGtsbokHofTi2Tbd+aKEBFZMRgpLy/Hvn37MHz48Ks/QKmU13fu3Nng/V566SWEh4fj4YcfNupxVCoVioqKal2ILNnmXd/UTAOlLN8V+auaejmsuvwQ9cjX8Z/1KfL4ocFxCPbztPawiYhcNxjJzc2VsxwRERG1bhfXMzMzDd5n27Zt+OSTT/DRRx8Z/Tjz5s1DUFBQ9SU6mq21yXpt3kUfkakVzyATwbW/ERgF3Psllqt6ISW7BEE+HnhkSJw1h0xE5JQs2vSsuLgYDzzwgAxEQkNDjb7frFmzZF6KnpgZYUBC1mzzLgKSdao+smrmpRtD0SE+AYgZiHKNAvPf3izPEcszgd4eVh41EZGLByMioHBzc0NWVu2KA3E9MjKy3vmnT5+Wiatjx46tvk2j0a3Gu7u748SJE4iPr7/e7uXlJS9ElpwV2Z6S0+g5Yslml6YLjoX2QIe41vK2pbvP4nz+ZYT6e2HywBgrjZaIyLmZFIx4enqid+/e2LBhQ3V5rgguxPVp06bVO79Tp044dOhQrdteeOEFOWPy7rvvcraDbJYnIpZnjN38LjxAVylTeLkC76w/JY+fuikBvp7cTYGIyBxMfjUVyyeTJ09Gnz59kJSUhPnz56O0tFRW1wiTJk1C69atZd6H6EPSrVu3Wvdv0aKF/Fr3diJrJqzW769an0hXjQzyRlKcLndk4aYU5JeWIyHcHxOS2lp8rERErsLkYGTcuHHIycnB7NmzZdJqjx49sHr16uqk1rS0NFlhQ+QoCauGVHV+x5yxXWSr93N5pfhse6q87fkxneHuxr9xIiJzUWi1WmNem21KJLCKqprCwkIEBgbaejjkoERTM9FLxBitgrxlIDK6Wyt5ferX+/Db4UwMaR+KLx9Kkp2EiYjIPO/fXPQmpydmRETlzG+HdR1TmzJtWAKeHdFBzogIu8/kyUBEXH1hTBcGIkREZsZghJyaqcmqwqCE0OpARKPR4pVfj8nj8Ult0TEywGJjJSJyVQxGyGmZkqxqKGFV+OlAOg6lFyLAy13OlhARkfkxC4/g6smqhhJWhRJVJd5cc1weP3FjguwtQkRE5seZEXK57qqGRNZJWBXeXnsSWUUqxIT44sGBsRYaKRERMRghJ+2ummvUuZMGxODmbq3k0ox+RkQ4nF6Iz3foSnlfvr0bvD3cLDZeIiJXx2CEXDphVQQiA+JD6gUzz/90SO7aOzYxCtd3CLPQaImISGAwQk7jWrqr1vTN7nP444IuafXFMZ0tMlYiIrqKCawEV++uWlNW0RX8e/UJefzP0R0RHqjbl4aIiCyHMyPkcgmrhpJV9V7+5SiKVZVIbBOECf24Ky8RkTUwGCEX664aj2dHdKw3IyJsPpGNX/68KDutvnrndQbPISIi82MwQi7WXTXMYJBRWFaBmT8cksdTBsWhW+sgs46ViIgaxmCEXLa7ak3/+vkIMouuIC7UD/83kp1WiYisiQms5JLdVWv67dBF2fZdfOs/9ybC15MxOhGRNfFVl1yyu6pedvEVPPeTbnlm6tB49Grb0qxjJSKipjEYIYdiju6qelqtFs/9eAiXyirQuVUgnr6JyzNERLbAYIRcqrtqTd/tu4D1x7Lh6abEO+MS4enOVUsiIltgMEIu1V1VLzW3FC/9fFQeTx/ZAZ0iA804WiIiMgU/CpLLdFfVu1yuxtSv96FEVYmk2GA8OqSdWcdLRESm4cwI2X1DM5EjYo7uqnqzVxzG8cxihPp74r0JPdncjIjIxhiMkNM0NGusu6resj3nZa6IOOW/43signvPEBHZHIMRcviGZk11V9U7klGIF1cclsf/N7IjBsaHXuNIiYjIHBiMkEM3NDM2YbXwcgX+9s1+qCo1GNYxDFNviDfLeImI6NoxgZXsh0aN4ztXoW/xRvRXHoUSGrMkrFaqNXh26UGcyytD6xY+eGdcDyiZJ0JEZDc4M0L24ehKYPUMdC3KwH89dTdlaIMxt2IS1miSmp2wKhqbiZmWjcez4eWuxKL7e6GFb9UDEBGRXWAwQvYRiCybJEKHWjdHIh+LPOZjasUz9QISkawqckQa6q6q98m2VHy16xwUCmD+uB7o3qaFxX4NIiJqHi7TkG1p1HJGpG4gIuhjjDkeX1Uv2YibWgV5y6oZ0V21sUBEJMK+uuqYPH7+ls64+bqGy32JiMh2GIyQbZ3bARRlNPhtEWtEKfKQpDxuVH6I3oG0S3h6yUFotcAD/WPw8OA4Mw+ciIjMhcs0ZNOGZh7HjqGPEeeHo8CohmbC6ZwSPPrlXlk5c2OncHkfhVinISIiu8RghGza0Ky/shBLjMgn/euYgXhnwI1NzoiIPWfGf7gLuSXl6BoViPfG94S7GycAiYjsGYMRsmlDs2RNJ1k1I5JVDccZCiAwCl0HjL6aRNKAc3m6QCS7WIWOEQH48qEk+HnxT5yIyN7xIyNZZUlm5+k8/HQgHc/9dKhWqqoGSlm+K4/r5bBWBR+jXweUbo0+xvn8MhmIZBZdQftwf3zzaD+E+HuZ+TchIiJL4MdGsvkeM6JsV5TvzvH4ElHIv/qNwChdINLltkYf48KlMtz34S5kFF5BuzA/GYiEMhAhInLumZGFCxciNjYW3t7e6NevH5KTkxs896OPPsKQIUPQsmVLeRk+fHij55PzLckYs9mdCEgGq/6L+8pfwN4+/wYm/wI8c6jJQOREZjH+8sFOpBdcRlyoH759tD/CA7j5HRGRUwcjS5cuxfTp0zFnzhzs378fiYmJGDVqFLKzsw2ev3nzZowfPx6bNm3Czp07ER0djZEjRyI9Pd0c4ycn2mNGLNns0nRBRee7gbghTS7N7Didi3s+2CGDnfgwXSDCXXiJiByPQiv6ZZtAzIT07dsXCxYskNc1Go0MMJ588knMnDmzyfur1Wo5QyLuP2mSLlegKUVFRQgKCkJhYSECAwNNGS7ZKBD5fHsqXv5V13DM1A3vts1oumpmxcF0/P27P1Ch1qJvbEt8NKkP27wTEdkZY9+/TcoZKS8vx759+zBr1qzq25RKpVx6EbMexigrK0NFRQWCgxveYVWlUslLzV+GnCdHxBBjG5qJ2Pl/W8/g9d+Oy+u3XBeJt+/tAW+PxmdRiIjISZZpcnNz5cxGRERErdvF9czMTKN+xowZMxAVFSUDmIbMmzdPRlL6i5h5IefKEalLzIiITewaa2hWqqrEM0sPVgcioqvqgvG9GIgQETk4q1bTvP7661iyZInMIxHJrw0RMy8iL6XmzAgDEufKERFzH8F+nnhhTGdEBvk0ueHdqaxiTP1mP1KyS+R54n5TBrHFOxGRywUjoaGhcHNzQ1ZWVq3bxfXIyMhG7/vWW2/JYGT9+vXo3r17o+d6eXnJCzlOW/ftKTlGz4joQ45X7+zWZGt34acDF/Dcj4dxuUKNiEAvLJjQC31jG17mIyIiJw5GPD090bt3b2zYsAF33HFHdQKruD5t2rQG7/fmm2/i1VdfxZo1a9CnjzE7kZAz54cYu8dM0ZUKvPLLUSzbe0FeH5wQivn39WAPESIiV1+mEcsnkydPlkFFUlIS5s+fj9LSUkyZMkV+X1TItG7dWuZ9CG+88QZmz56NxYsXy94k+twSf39/eSHnaOturBfHdMaDg+KarJZZfzQLzy8/hKwiFcQed0/d2B5P3dS+yfsREZELBCPjxo1DTk6ODDBEYNGjRw+sXr26Oqk1LS1NVtjoLVq0SFbh3HPPPbV+juhT8q9//cscvwNZeVlm15k8zPyhdlt3Y8t2mwpE8kpUcrZl5R8Z8npsiC9ev7s7+rcLMcPoiYjIKfqM2AL7jDhH2W5j1TKVag2W7j2P/6w9ifzScrkn3qND2uHZER1YLUNE5KAs0meEXFdzl2WayhERsfDG49mY99txWSkjdIoMwBt3d0didAszjJyIiOwdgxFqmEYNnNsBTXEmVv6cAQXaQWtCa5ppwxIwKCG0wbLdg+cL8Ppvx7DrjG5zvJa+HjIvZGK/GHi6c0NpIiJXwWCEDDu6Elg9AyjKkOHH+wAyvIIxt2KS3NTOmPwQscRSNwgRMyG/n8rFB1tOY8fpPHmbCDweGhSHqUPjEeTjYdFfi4iI7A+DEarXM8TtxM/om/yMCB2q8z2ESORjkcd8TK14psGApKG27hVqDVYduoj/bTmDoxd17f3dlQrc3qM1nh3RHm1a+lr0dyMiIvvFYIRqJadmFZZhm9dcaMWCTJ2VFXFdowXmeHyFdao+cpfdpvJDRB7Isr3n8eP+C8gtKZe3+Xq64b6+bfHwkDi0buFjnV+QiIjsFoMRqpWc2l95HFEKXQ6HISIgiUIekpTHsUvTpfr2Fr4eWDi+F/rHh+BSWTm+TU7D9/suYN+5S9XniGZlkwbEyAt32CUiIrh6MKJfksguvoLwAO8m90ZxleTUcBQYdXf9efpn7J+jOuJEVjHe3XgKe8/myxkUQTynwzqGY1zfaAztGAYPNyamEhFRbS4ZjBjql9HKyBblzp6cmg3jymkztLq9YXw83WTS6XM/Ha71/etaB+GW61rh7l6tER7Y8KaIRERELtf0rKF+GeITvhIafDtSjaSwSsA/AogZCCjdXCI5VT+T8beKpzDb42uZrKqfKNJoFTirjcBBbQIOauKxH51xRNMWNf9yRMt2sXndqK6RGNU1ggmpREQENj0zcZv7kcpkzPH4ElFba+RLBEYBo98AutwGV0lOfdr9Jzxb/jeMcduFU5o2OKltg2PaGBTBr97PjQv1w8D4ENlLZEC7ELT0Yx4IERGZzqWCETEzYKiV+ShlsixZrafoIrTLJuHkDQtxvOVQx8ktqcoHQUkWknPc8cRaN6ihRD/FcfigHMe1bZGlCcZFbTDStaFy1iNNG4Fz2ojqoGN35dXkVMEL5ejmkYHEhLbo0aM3+sS0RBQrYYiIyAxcKhgRyap1iaUZMSMij+vFGGIOAQjY/CKeVb0rS1ltnVsiVtVUlRqUlatRqqpE0eUKJJ/NR2bhFZkcGlp6AsVH1qLwigaXtP7IQxDaK4JlwLFH2wk9VR82+Rg+UOEyvKCAFu0UGZjW0wO39oqHR7vbnGbZioiI7IdLBSNiZqOuOMVF7NB0hVYr3nqvRiPiWFv9VYEYRSZStVFyZuXxr/ejf7tgdIwIQHSwLjdC5E+otVpoxEUjvuqWhcR18VV+T6NFpbiodV/VGo08LldrZFOwCrVWfhXBRnml/qtafr1SocblcjUuV6ir8zsME8HCzY0+DyEoRKQiH60U+fJrjCILbRXZ1V8fqviHLNuNDPLBP8aOdY2kXiIishmXSmAVQcHgNzbKWQS7/6WNIGZ1NDKAEhcteitOIkxRgBaKUgSiDC0UxQhFEUIURQhGETQKhbwercg2MAukC7zKfSOxesRahAf6OcaSFBER2S0msBog3ljFEouoptG9fQNBKEFPZYr8vm4O5Or8yNW5EWCvpgMuwbhAyMdDiYRwfwT6eMDP013mVojHziq6Imc5ArzdEd3SV14XMx0tfDzk+WKZJS2/DGXllQjx90LXVoFIzStDiaoCkYE+SFLvh9/vr8KnJA2+uAI3hVaW2IqS3EL4Y4nnK02O7T8V9+BZ9+/l7ErtOEMhf0+vW9/E7V3aNuPZJSIiah6Xmhkx1GdEzC5s83pKLlcYascl3rQzEYLBVTkjzSG6kwoFZRX1qleMPUck2X7gqUuyNVSS+2nlzXjE47cmx/JU+TSo4K6rHKrZaTWwNTD6dYevHCIiIvvBmZFGiByIEV0iqzuwFl96Ba22PFH1XW29N/q5FQ80OxCpG2DU/dnGnKNPshVhY0MluXe4bzdqLKKpmcgHWa/qg77K45g1uAUSO3dyqp4qRETkWFwyGBHEssmA+JCqaxOBiIDqrqR6YkZEBCIN7VBrLUlG7BcjckFytQEIRrHBfBCN+H20IUjWdJLXw4N88eDY+5HI5FQiIrIxlw1G6hHLE53G1OrPMX6t2zXNiJiLsfvFLFcPwkNuqxvMByke+jLeadnLcfqlEBGRS2AwUpNYpogbIg/FXMjC8Pp72NiCsfvFrNf0wV5NJ8wW+SCo3UlWMfp1dOxyGzpabphERETNwmDEyNyS9Ucz8cn2szYZh1haEVUzNfeLMZRkK86LCPLFn7c+iij/VDnD42x77BARkfNxyWoac+72ay01W9bXDEj0hch7kuZD3XEsl1+IiMjh3r8ZjDRz91tRhSNyLy6VqvDyr8esEqCMVibrlmBYkktERA6Apb1WqcLRGdWtVXWAcja3DN8mpyGz6Gpw0tLXQxYMX0ufEeGPgOu5BENERE6HwYgFApRpNybUmj0RSydCzdt6x7TEvnOXTD5HtwTTxma/KxERkblxmYaIiIhs+v5t+yYaRERE5NIYjBAREZFNMRghIiIim2IwQkRERDbFYISIiIhsisEIEREROV4wsnDhQsTGxsLb2xv9+vVDcnJyo+d/99136NSpkzz/uuuuw6pVq5o7XiIiInL1YGTp0qWYPn065syZg/379yMxMRGjRo1Cdna2wfN37NiB8ePH4+GHH8aBAwdwxx13yMvhw4fNMX4iIiJytaZnYiakb9++WLBggbyu0WgQHR2NJ598EjNnzqx3/rhx41BaWopffvml+rb+/fujR48e+OCDD4x6TDY9IyIicjwW2ZumvLwc+/btw6xZs6pvUyqVGD58OHbu3GnwPuJ2MZNSk5hJWb58eYOPo1Kp5EVP/BL6X4qIiIgcg/59u6l5D5OCkdzcXKjVakRERNS6XVw/fvy4wftkZmYaPF/c3pB58+Zh7ty59W4XMzBERETkWIqLi+UMiUNtlCdmXmrOpoiloPz8fISEhEChEBvFmS9iEwHO+fPnufxjBD5fxuNzZTw+V8bjc2U8Plf28VyJGRERiERFRTV6nknBSGhoKNzc3JCVlVXrdnE9MjLS4H3E7aacL3h5eclLTS1atICliCeff6zG4/NlPD5XxuNzZTw+V8bjc2X756qxGZFmVdN4enqid+/e2LBhQ61ZC3F9wIABBu8jbq95vrBu3boGzyciIiLXYvIyjVg+mTx5Mvr06YOkpCTMnz9fVstMmTJFfn/SpElo3bq1zPsQnn76adxwww34z3/+gzFjxmDJkiXYu3cvPvzwQ/P/NkREROT8wYgo1c3JycHs2bNlEqoo0V29enV1kmpaWpqssNEbOHAgFi9ejBdeeAHPPfcc2rdvLytpunXrBlsTS0GiX0rdJSEyjM+X8fhcGY/PlfH4XBmPz5VjPVcm9xkhIiIiMifuTUNEREQ2xWCEiIiIbIrBCBEREdkUgxEiIiKyKQYjNdx2221o27YtvL290apVKzzwwAPIyMiw9bDsztmzZ+UuzHFxcfDx8UF8fLzMxBZ7F1F9r776qqwq8/X1tWjzPke0cOFCxMbGyv9zYhPO5ORkWw/JLm3duhVjx46VXSxFF+rG9vZydaKthNjMNSAgAOHh4XKX+BMnTth6WHZp0aJF6N69e3WzM9H/67fffrPJWBiM1DBs2DAsW7ZM/uH+8MMPOH36NO655x5bD8vuiH2IRLO7//3vfzhy5AjeeecduQOzKN2m+kSQ9pe//AVTp0619VDsytKlS2XfIhHI7t+/H4mJiXITzezsbFsPze6IXk7i+RHBGzVuy5YteOKJJ7Br1y7ZYLOiogIjR46UzyHV1qZNG7z++utyA1zR/+vGG2/E7bffLl/XrU6U9pJhK1as0CoUCm15ebmth2L33nzzTW1cXJyth2HXPvvsM21QUJCth2E3kpKStE888UT1dbVarY2KitLOmzfPpuOyd+Jl+6effrL1MBxGdna2fM62bNli66E4hJYtW2o//vhjqz8uZ0YaIDbm++abb+T0uoeHh62HY/cKCwsRHBxs62GQA80WiU9jw4cPr75NNEsU13fu3GnTsZHzvTYJfH1qnFqtlh3SxQySLbZrYTBSx4wZM+Dn5yd3CBbdZFesWGHrIdm9lJQUvPfee/jrX/9q66GQg8jNzZUvfvrOzXriuujsTGQOYjn5mWeewaBBg+yi67c9OnToEPz9/WX31ccffxw//fQTunTpYvVxOH0wMnPmTJnw1dhF5EDo/eMf/8CBAwewdu1auUOx2GvHVZrUmvpcCenp6Rg9erTMiXj00UfhKprzXBGRdYnckcOHD8tP/GRYx44dcfDgQezevVvmtYm9544ePQprc/p28GIfnby8vEbPadeundyRuK4LFy4gOjoaO3bscIldhk19rkSl0dChQ9G/f398/vnntfYkcnbN+bsSz5H4lFZQUABXJ5ZpRHXR999/L6sd9MQLoXh+OCPZMBHoik+vNZ83qm/atGny70hUIonKPzKOWCoVFZKiQMGuN8pzNGFhYfLS3Ck+QaVSwRWY8lyJGRFRfdS7d2989tlnLhWIXOvfFUEGaeJvZ8OGDdVvquL/m7gu3kSImkt8vn7yySdlwLZ582YGIiYS/w9t8Z7n9MGIscQU1Z49ezB48GC0bNlSlvW++OKLMkJ0hVkRU4hARMyIxMTE4K233pKzBHqRkZE2HZs9ErlHIiFafBV5EmJKVEhISJBrta5KlPWKmZA+ffogKSkJ8+fPl8lzU6ZMsfXQ7E5JSYnMzdJLTU2Vf0ciKVP0RqLaSzNip3gxKyJ6jehzkIKCgmRfJLpq1qxZuPnmm+XfUHFxsXzeRAC3Zs0aWJ3V63fs1J9//qkdNmyYNjg4WOvl5aWNjY3VPv7449oLFy7Yemh2WaIq/nQMXai+yZMnG3yuNm3apHV17733nrZt27ZaT09PWeq7a9cuWw/JLom/FUN/Q+Jvi2pr6LVJvG5RbQ899JA2JiZG/v8LCwvT3nTTTdq1a9dqbcHpc0aIiIjIvrnWQj8RERHZHQYjREREZFMMRoiIiMimGIwQERGRTTEYISIiIptiMEJEREQ2xWCEiIiIbIrBCBEREdkUgxEiIiKyKQYjREREZFMMRoiIiMimGIwQERERbOn/AVyTDpTVNbQPAAAAAElFTkSuQmCC",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(x, y)\n",
"plt.scatter(x_train, X_train_s @ beta + y_offset)\n",
"plt.scatter(x_test, X_test_s @ beta + y_offset)"
]
},
{
"cell_type": "markdown",
"id": "2993cfae",
"metadata": {},
"source": [
"## Exercise 4 - Testing multiple hyperparameters"
]
},
{
"cell_type": "markdown",
"id": "e6cdc89a",
"metadata": {},
"source": [
"**a)** Compute the MSE of your ridge model for polynomials of degrees 1 to 5 with lambda set to 0.01. Plot the MSE as a function of polynomial degree.\n",
"\n",
"**b)** Compute the MSE of your ridge model for a polynomial with degree 3, and with lambdas from $10^{-1}$ to $10^{-5}$ on a logarithmic scale. Plot the MSE as a function of lambda.\n",
"\n",
"**c)** Compute the MSE of your ridge model for polynomials of degrees 1 to 5, and with lambdas from $10^{-1}$ to $10^{-5}$ on a logarithmic scale. Plot the MSE as a function of polynomial degree and lambda using a [heatmap](https://matplotlib.org/stable/gallery/images_contours_and_fields/image_annotated_heatmap.html)."
]
}
],
"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.9.15"
}
},
"nbformat": 4,
"nbformat_minor": 5
}