From 73cda3ca72adc67449c1042d2d3b52f11a300e31 Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Mon, 11 Sep 2023 22:25:12 +0200 Subject: [PATCH] update --- .../week37/programs/codeexamplesscaling.ipynb | 223 +++++++----------- ...codeexamplesscaling.do.txt => dill.do.txt} | 157 ++---------- 2 files changed, 107 insertions(+), 273 deletions(-) rename doc/src/week37/programs/{codeexamplesscaling.do.txt => dill.do.txt} (59%) diff --git a/doc/src/week37/programs/codeexamplesscaling.ipynb b/doc/src/week37/programs/codeexamplesscaling.ipynb index e85829f91..a4b261ba5 100644 --- a/doc/src/week37/programs/codeexamplesscaling.ipynb +++ b/doc/src/week37/programs/codeexamplesscaling.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "a6f6db7d", + "id": "c385e1d6", "metadata": {}, "source": [ " 34\u001b[0m interceptRidge \u001b[38;5;241m=\u001b[39m y_scaler \u001b[38;5;241m-\u001b[39m \u001b[43mX_train_mean\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m@\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mbeta_Ridge\u001b[49m\n\u001b[1;32m 35\u001b[0m \u001b[38;5;28mprint\u001b[39m(interceptOLS)\n\u001b[1;32m 36\u001b[0m \u001b[38;5;28mprint\u001b[39m(interceptRidge)\n", + "\u001b[0;31mValueError\u001b[0m: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 3 is different from 2)" + ] } ], "source": [ - "\n", "np.random.seed(2018)\n", "n = 100\n", - "# we do not include the intercept\n", "d = 2\n", "Lambda = 0.01\n", - "true_beta = [2, 0.5, 3.7]\n", - "\n", "# Make data set.\n", "x = np.linspace(-3, 3, n)\n", - "y = 2 + 0.5*x + 3.7*x**2\n", + "y = 2.0 + 0.5*x + 5.0*(x**2)+ np.random.randn(n)\n", "\n", "#Design matrix X does not include the intercept. \n", "X = np.zeros((len(x), d))\n", "for p in range(d): \n", " X[:, p] = x ** (p+1)\n", "\n", - "\n", "#Split data in train and test\n", "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n", - "\n", "# Scale data by subtracting mean value,own implementation\n", "#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable\n", "X_train_mean = np.mean(X_train,axis=0)\n", @@ -501,7 +501,7 @@ "\n", "#Calculate beta\n", "beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)\n", - "beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d)\n", + "Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d)\n", "print(beta_OLS)\n", "print(beta_Ridge)\n", "# calculate intercepts and print them\n", @@ -513,8 +513,6 @@ "#predict value with intercept\n", "ytilde_test_OLS = X_test_scaled @ beta_OLS+y_scaler\n", "ytilde_test_Ridge = X_test_scaled @ beta_Ridge+y_scaler\n", - "\n", - "\n", "#Calculate MSE\n", "\n", "print(\" \")\n", @@ -535,11 +533,9 @@ }, { "cell_type": "markdown", - "id": "0686ad2c", + "id": "edf2436a", "metadata": {}, "source": [ - "We see that we get the same values for the parameters! As it should be. The MSE may however change (not the case here).\n", - "\n", "Finally, instead of using our own function we repeat the same example\n", "using the **standardscaler** functionality of the library\n", "**Scikit-Learn**. Here we limit ourselves to Ridge regression only." @@ -547,29 +543,25 @@ }, { "cell_type": "code", - "execution_count": 45, - "id": "a721dffd", + "execution_count": 41, + "id": "f1ec76ea", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "[0.5 3.7]\n", - "[0.49997445 3.69992817]\n", - "1.9999999999999947\n", - "2.00021345091295\n", " \n", "test MSE of OLS\n", - "1.7086727207087115e-28\n", + "1.1394311129039263\n", " \n", "test MSE of Ridge\n", - "6.420926689334773e-08\n" + "96.31278139904835\n" ] }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABQGklEQVR4nO3deXxU5dk//s+ZyWSSkAVCCDORECIgApFNQCIoSwkGNeLS/ixpLTxaK4i0wOMjVb+WUAoCbQUtNe3jgvhggFqLaJVAKJsskTVKCCBLIIFkEiBkIetk5vz+GM6QZSY5k8xyZubzfr146Zw5OblzZzJznXu5LkEURRFEREREbqLydAOIiIjIvzD4ICIiIrdi8EFERERuxeCDiIiI3IrBBxEREbkVgw8iIiJyKwYfRERE5FYMPoiIiMitAjzdgJbMZjOKiooQFhYGQRA83RwiIiKSQRRFVFVVISYmBipV22Mbigs+ioqKEBsb6+lmEBERUQcUFhaiV69ebZ6juOAjLCwMgKXx4eHhTr220WjE9u3bMWXKFGg0Gqde29ewr+RjX8nHvnIM+0s+9pV8ruqryspKxMbGWj/H26K44EOaagkPD3dJ8BESEoLw8HC+ONvBvpKPfSUf+8ox7C/52Ffyubqv5CyZ4IJTIiIicisGH0RERORWDD6IiIjIrRh8EBERkVsx+CAiIiK3YvBBREREbsXgg4iIiNyKwQcRERG5FYMPIiIiP2EyiziUXwYAOJRfBpNZ9Eg7GHwQERH5gczcYoxbsRPPrjsMAHh23WGMW7ETmbnFbm8Lgw8iIiIfl5lbjNnrj6G4oq7ZcUNFHWavP+b2AITBBxERkQ8zmUUs/jIPtiZYpGOLv8xz6xQMgw8iIiIfdii/rNmIh9EM7C0WIN6KNUQAxRV11rUg7sDgg4iIyIeVVjWfajGLArZcav3x3/I8V2LwQURE5MOiw4JaHRvUTUTLyve2znMVBh9EREQ+bHR8JPQRzQOLEVG313cIAPQRQRgdH+m2NjH4ICIi8mFqlYBFKYOaHBExuKsl+JAGPxalDIJaJbT6Wldh8EFEROTjkhP0mDigBwBAJQCBastxXUQQ0n8+AskJere2J8Ct342IiIjczmQWceJKJQBgweT+wM3T+HDGKIzpF+3WEQ8JRz6IiIh8XPaF67h2sx5dQzT4r/v7ALCsBfFE4AEw+CAiIvJ5X+QUAQCmJugRGOD5j37Pt4CIiIhcpqHRjK230qc/NjTGw62x4JoPIiIiHyRVsN15ugSVdY2IDtNidHwkzKZGTzeNwQcREZGvycwtxuIv85qlVa+ub0RWngE/GhDlwZZZcNqFiIjIh9irYFvdYMLs9cew41SJh1p2G4MPIiIiH9FWBVvJ8q2n3dYeexh8EBER+YiWFWxbEgEYKt1XQM4eh4KP9PR0DBkyBOHh4QgPD0diYiK2bt1qfX7mzJkQBKHZvzFjxji90URERNSaOyvTdoZDC0579eqF5cuXo1+/fgCAdevWYdq0aTh+/DgGDx4MAEhOTsbatWutXxMYGOjE5hIREZE97qxM2xkOBR8pKSnNHi9duhTp6enIzs62Bh9arRY6nc55LSQiIiJZpAq2hoo6m+s+BAC68CAA1W5uWXMd3mprMpnw6aeforq6GomJidbju3fvRnR0NLp27Yrx48dj6dKliI6Otnud+vp61NfXWx9XVlpyzxuNRhiNxo42zybpes6+ri9iX8nHvpKPfeUY9pd87KvbfvfIAMzbmIN6MwAICFSJEITbFWwXPtQfxks5LvuMlUMQRbGtRbGtnDhxAomJiairq0NoaCgyMjLw8MMPAwA2bdqE0NBQxMXFIT8/H2+88QYaGxtx9OhRaLVam9dLS0vD4sWLWx3PyMhASEiII00jIiIiAP+5IuCLAjX6hon4dYLJLd+zpqYGqampqKioQHh4eJvnOhx8NDQ0oKCgAOXl5fjss8/w/vvvY8+ePRg0aFCrc4uLixEXF4eNGzfiySeftHk9WyMfsbGxuHbtWruNd5TRaERWVhaSkpKg0Wicem1fw76Sj30lH/vKMewv+dhXzT3214M4ZajCc/fHYUivcESFanFvXDeoVYLL+qqyshJRUVGygg+Hp10CAwOtC05HjhyJw4cP4+2338bf//73Vufq9XrExcXh7Nmzdq+n1WptjopoNBqXvYBceW1fw76Sj30lH/vKMewv+dhXwNmSKpwyVEGjFjB38l3oGmJ744ez+8qRa3U6z4cois1GLpq6fv06CgsLodfrO/ttiIiISIbPc64AAMbfFW038PA0h0Y+XnvtNUydOhWxsbGoqqrCxo0bsXv3bmRmZuLmzZtIS0vDU089Bb1ej4sXL+K1115DVFQUnnjiCVe1n4iIiGDJbvrthevYcKgQAJAyVLk3/g4FHyUlJXjmmWdQXFyMiIgIDBkyBJmZmUhKSkJtbS1OnDiBjz/+GOXl5dDr9Zg4cSI2bdqEsLAwV7WfiIjI79kqJPfm16egDVAhOUF5QYhDwccHH3xg97ng4GBs27at0w1yFam0MGBJPzumXzTUKqGdryIiIlI2qZBcy90jJZX1mL3+GNJ/PkJxAYhf1HbJzC3GuBU78ey6wwCAZ9cdxrgVO5GZW+zhlhEREXVcW4XkpGOLv8yDyezQxlaX8/ngw15pYUNFHWavP8YAhIiIvJacQnLFFXXWkX+l8Ongw1sjQiIiIjnkFpJTWsE5nw4+WkaEjWZg3Q+3f2SlRoRERERyyC0kp7SCcz4dfLSM9EwicOy6Ci0HOpQWERIREckhFZKzRwCgjwjC6PhI9zVKBp8OPlpGetIP2zL4UFpESEREJIdaJWBRSuvyJsDtQnKLUgYpbnenTwcfUkQodbnU96ZbwYdSI0IiIiK5hsV2g63QQhcRpMhttkAHart4EykinL3+GARYgg+tSkS96favSYkRIRERkVxffHcFIoCRcd3w31MGoLSqDtFhlhtrpX6++fTIBwAkJ+iR/vMR0EUEQRCAoZGWYY+QQLViI0IiIiK5Nh8vAgA8MeIOJPbtjmnDLP9VauAB+PjIhyQ5QY+kQTpknyvFgexvcegqEKASMGFAtKebRkRE5DApa3dO4Q2cKq5EoFqFR+7xnptpvwg+AMsUzOj4SJTmidCFa2GorMeu06WY6kW/LCIiIlt1XAQByL5w3WtG831+2qUllQA8dqvS37+OX/Fwa4iIiOSzl7W7vtHsVVm7/S74AIDHh8YAAHafKcX2kwZsybmCg+evM9MpEREpVltZuyXekrXbb6ZdmurfMxSx3YJReKMWv/q/o9bj+oggLEoZ5DXDVkRE5D8cqeOS2Le7+xrWAX458rHjVAkKb9S2Os5ic0REpFTeWsfFFr8MPpZvPW3zOIvNERGRUnlrHRdb/DL4MFR6X/lhIiLyby2zdrfkTVm7/TL4kMMbhq2IiMh/SFm7bY3LK7mOiy0MPuzwhmErIiLyL8kJejwxPKbVcSXXcbHFL3e76MKDUHCj3m70qPOSYSsiIvIvoigi90olAOCXD8TjnjsiFF/HxRa/HPn47dS7AcDuvJm3DFsREZF/OVlUibOlNxEYoMLcSf29oo6LLX4ZfEwe2NNabK6prsEarxq2IiIi/2Ayizh4/jpWbrPs1kwa2BMRwRoPt6rj/HLaBbhdbO5Qfhk+OpCPbSdLcE+vCAYeRESkKLZquRw4fw2ZucVe+5nllyMfErVKQGLf7njt4YEAgP3nrsHQRvY4IiIid7JXy+VGjdGrk2L6dfAhieveBaP6dINZBN7+z1nWeiEiIo/zpVouLfnttEtLA3qG4fDFG9hwqAAbDhUAYK0XIiLyHF+q5dISRz5gGdZa/21Bq+Os9UJERJ7iS7VcWvL74EMa1rKFtV6IiMhTfKmWS0t+H3w4MqxFRETkLlItF3u8qZZLS34ffPjysBYREXkvqZaLLd5Wy6Ulvw8+fHlYi4iIvNvkgT0RHtR6b4i31XJpye93u0jDWoaKOtZ6ISIiRdl79ioq6xoR2SUQb/90GMqqG7yylktLfh98SMNas9cfgwDYDEC8dViLiIi826dHLgMAHh92Bx7o38PDrXEev592ASyp1m3VeglQCXj3Z8O9dliLiIi8143qBuw4VQIA+PG9vTzcGufy+5EPSdNaL5euV+N3W06iwWSGvmuIp5tGRER+xGQWcSi/DP88WgijScTgmHAMign3dLOcisFHE1Ktl8S+3ZF94To+zynCp0cKMSy2q6ebRkREfsBWEbnCshqvLiJnC6dd7PjxvbEAgC++K8KeM1dZ74WIiFzKXhG5yrpGn8u27VDwkZ6ejiFDhiA8PBzh4eFITEzE1q1brc+Looi0tDTExMQgODgYEyZMwMmTJ53eaHe4v293RHYJRFVdI2asPYTfbMzB9PeyMW7FTp96ARARkef5chE5WxwKPnr16oXly5fjyJEjOHLkCCZNmoRp06ZZA4yVK1firbfewpo1a3D48GHodDokJSWhqqrKJY13pe15BpRVN7Q6znovRETkbP6Wbduh4CMlJQUPP/ww7rrrLtx1111YunQpQkNDkZ2dDVEUsXr1arz++ut48sknkZCQgHXr1qGmpgYZGRmuar9LsN4LERG5k79l2+7wglOTyYRPP/0U1dXVSExMRH5+PgwGA6ZMmWI9R6vVYvz48Thw4ABeeOEFm9epr69HfX299XFlZSUAwGg0wmg0drR5NknXa++6h/LLUHazFlo10GACRAhQCyICmoRqZTdrkX2u1GeTj8ntK2JfOYJ95Rj2l3ze3ldRIQHQqi03tEYTYLbxuSOd19mf0VV95cj1BFEUHbp9P3HiBBITE1FXV4fQ0FBkZGTg4YcfxoEDBzB27FhcuXIFMTEx1vN/9atf4dKlS9i2bZvN66WlpWHx4sWtjmdkZCAkxPPbXA+VCvjkvBrdtSL+33ATmGuMiIhcpaIBWHRUDRECXhvWiJ7Bnm6RfDU1NUhNTUVFRQXCw9veGuxw8NHQ0ICCggKUl5fjs88+w/vvv489e/agvLwcY8eORVFREfT629uBnn/+eRQWFiIzM9Pm9WyNfMTGxuLatWvtNt5RRqMRWVlZSEpKgkajsXveofwyPLvuMABAFIEGMwAI0KjEZsHHhzNG+fTIh5y+IvaVI9hXjmF/yecLfbXjVAnmbshBoyhAgIhAteW49LGz6ulhmDywZ6e/j6v6qrKyElFRUbKCD4enXQIDA9GvXz8AwMiRI3H48GG8/fbbWLhwIQDAYDA0Cz5KS0vRs6f9ztJqtdBqta2OazQal72A2rv2mH7RiAwNblXvxWi2vASkei9j+kX7fNp1V/4efA37Sj72lWPYX/J5c18l33MHIr88jdKqeogQUG+yHNdHBGFRyiCn5/lwdl85cq1OJxkTRRH19fWIj4+HTqdDVlYWhg8fDsAySrJnzx6sWLGis9/GrVjvhYiI3O3wxRsorapHl0A11qSOQGWd0SeKyNniUPDx2muvYerUqYiNjUVVVRU2btyI3bt3IzMzE4IgYN68eVi2bBn69++P/v37Y9myZQgJCUFqaqqr2u8yUr2XlpnmwoMCsPLHQ3wq0xwREXmOlE59VdYZAMAjQ/SYeHe0h1vlWg4FHyUlJXjmmWdQXFyMiIgIDBkyBJmZmUhKSgIAvPLKK6itrcWLL76IGzdu4L777sP27dsRFhbmksa7WtN6L5uOFODz40WIjQxh4EFERE5hK536f06V+lw69ZYcCj4++OCDNp8XBAFpaWlIS0vrTJsURar3MkAXhq++L8bJokqcLKrA4JgITzeNiIi8mJROveXU/vXqBsxefwzpPx/hswEIa7vIFNklEEmDLAtn395xlrVeiIiow/wtnXpLrGrrgL49ugAAtueVYHteCQDXrUImIiLf5Ug69cS+3d3XMDfhyIdMmbnF+MvO862Os9YLERE5yt/SqbfE4EMG1nohIiJnig4Lcup53obBhwz+Vm2QiIhca3R8JPQR9gMLAZZpfV/Nos3gQwZ/Hx4jIiLnkpJZ2iKlE/PlZJZccCqDvw+PERGR8w3U265/ovODjQwMPmSQhsda1nqRSLVefHV4jIiInO8fRwoBAA/0j8KLE/qhtKrOZ9Opt8TgQwbWeiEiImdqNJnx6ZHLAIDpo3v75HbatnDNh0xSrRedjQVCy55I8OnhMSIicg6TWcTB89exfKulem1kl0BMHmi/8ruv4siHA5rWeimtqkP67vM4bahCaVWDp5tGREQKZ6uOS73RhJ2nS/zuBpYjHw6Sar1MG3YHZo3vC8Ayb8ccH0REZI9Ux6Vl2obqBpNfJqpk8NEJyQk6RARrcKW8Fv+79wLrvRARUSv+XsfFFk67dEKQRo0Rvbti15mrWJF52nqc9V6IiEji73VcbOHIRydk5hZj15mrrY6z3gsREUmYqLI1Bh8dxHovREQkBxNVtsbgo4NY74WIiOTw9zoutjD46CAOoxERkRz+XsfFFgYfHcRhNCIikmvywJ7oFqJpdVwXEYT0n4/wuw0K3O3SQaz3QkREcu0+cxU3aozoFqLB2z8djhs1DX5Tx8UWBh8dxHovRETUFpNZtGbE/mj/RQDAj+/thQfv6uHZhikAg49OkOq9tEyXqxKANdOH+90wGhERWdhKpQ4AvboFe6hFysLgo5Oa1nu5Ul6DtC/ycLO+EUGBak83jYiIPEBKpW5rRDztizz0DA/y+5tTBh9OINV7AbrjjKEK732TjzU7z6GqrtGv5/SIiPyN3FTqSYN0fv25wODDyaQhtWMF5ThWkAOA6daJiPwFU6nLw622TpSZW4y0L1pnPWW6dSIi/8AcUPIw+HCStobamG6diMg/MAeUPAw+nITp1omISMoBZW81hz+mUreFwYeTcKiNiIikHFD2kk8CzAEFMPhwGg61ERERYEnBMCMxrtVxf02lbgt3uzgJ060TEZHku8sVAICnR8Xi/r7dmXahBQYfTtJeunURHGojIvJlUjr1o5fKkFNYDo1awP88NABRoVpPN01xGHw4kb106wDw8/t6c6iNiMhH2UqnrlYJOHKxjO/9NjD4cLKm6dZLqyy7Wz75tgDHCsohiiIEgSMfRES+xF469TqjGbPXH+M6DxsYfLjA7XTrwPi7euCfRy8jr7gS/3fwEiJCNJz7IyLyEUyn3jEMPlysa0gg7o3rhgPnr+N3X5y0HmfKdSIi78d06h3DrbYulplbjAPnr7c6zpTrRETejzmeOsah4OPNN9/EqFGjEBYWhujoaDz++OM4c+ZMs3NmzpwJQRCa/RszZoxTG+0tpOE4W5hynYjI+zHHU8c4FHzs2bMHc+bMQXZ2NrKystDY2IgpU6agurq62XnJyckoLi62/vv666+d2mhvwZTrRES+TcrxZA/Tqdvm0JqPzMzMZo/Xrl2L6OhoHD16FA8++KD1uFarhU6nc04LvRiH44iIfJuU42nW+mOtnmM6dfs6teC0osKSwS0ysnlEt3v3bkRHR6Nr164YP348li5diujoaJvXqK+vR319vfVxZWUlAMBoNMJoNHamea1I13P2de2JCgmAVn17SqXRDJhEAQJEBKqbn+euNsnl7r7yZuwr+dhXjmF/yefJvhp3Zzd00apRXW9CgEqE+lacoQsPwm+n3o0fDYhS1O/QVX3lyPUEURQ7tOBAFEVMmzYNN27cwDfffGM9vmnTJoSGhiIuLg75+fl444030NjYiKNHj0KrbZ3lLS0tDYsXL251PCMjAyEhIR1pmmJdrwOWHFdDhIBXhzZC51s/HhGRXzpUKuCT82p0CxTxxgiTNfjwNzU1NUhNTUVFRQXCw8PbPLfDwcecOXPw1VdfYd++fejVq5fd84qLixEXF4eNGzfiySefbPW8rZGP2NhYXLt2rd3GO8poNCIrKwtJSUnQaDROvbY9O06VYP6mHACWNR5GE2CGALUgQqMCVj09DJMH9nRLWxzhib7yVuwr+dhXjmF/yefuvjKZRRy9dANXq+rw9s4LyL9eg5eT+uOFB+Nd/r07y1V9VVlZiaioKFnBR4emXebOnYsvvvgCe/fubTPwAAC9Xo+4uDicPXvW5vNardbmiIhGo3HZC8iV125p6pBeEFTqVml3zaKAN58aiqlD2u4/T3NnX3k79pV87CvHsL/kc0df2UqlDgA9I4K86vfk7L5y5FoOBR+iKGLu3LnYvHkzdu/ejfj49iO869evo7CwEHq9/ybTappyvaSyFn/c9gOulNeiolY5c4BERNQ+e6nUAeCVf55AWJCGySNlcGir7Zw5c7B+/XpkZGQgLCwMBoMBBoMBtbW1AICbN2/i5ZdfxsGDB3Hx4kXs3r0bKSkpiIqKwhNPPOGSH8BbSCnXHx/eC7PG3wkA+PveC9h8/DIOnr/OXB9ERAonN5U638/b59DIR3p6OgBgwoQJzY6vXbsWM2fOhFqtxokTJ/Dxxx+jvLwcer0eEydOxKZNmxAWFua0Rnu7sKAACLDk+Ji/6TsATLdORKR0TKXuPA5Pu7QlODgY27Zt61SDfF1mbjHmb/quVeQspVtn9UMiImVi7ibnYW0XN2pryI7p1omIlI2p1J2HwYcbMd06EZH3klKp20vjwVTq8jH4cCMO2REReS8plbqtsWmmUncMgw834pAdEZF3S07QY1Sfbq2O6yKCuGbPAZ2q7UKOkYbsDBV1drdqcciOiEi5SirrcLygHACw7Il70EWrRnSY5X2bIx7yMfhwI2nIbvb6YxAAmwEIh+yIiJTHZBZxKL8Ma/fno9EsYmRcN6Te19vTzfJanHZxs+QEPdJ/PgK6iNZTK3frwjhkR0SkMJm5xRi3Yiemv5eN7XklAIBzV28iM7fYwy3zXhz58ICm6dZLq+ogQMC8Tcdx2lCFfxwuhFaj4jAeEZEC2EunXl5jZG6mTmDw4SFSunXJ2gP5OF5Qjlc++956jFlPiYg8R2469aRBOt4oOojTLgqQmVtsXcDUlJT1lEN7RETux9xMrsPgw8OkyNoWZj0lIvIc5mZyHQYfHsbImohImZibyXUYfHgYI2siImWScjPZw3TqHcfgw8MYWRMRKZNaJeDlKXfZfI7p1DuHwYeHsVAREZFylVY1AAACWgQYTKfeOdxq62HMekpEpCxSNtPiilq8980FAMDSJ+5B78gQlFbVMQ+TEzD4UAAp6+niL/NaLT5Ne4x5PoiI3CUzt7jVe7FKAEICVc1yM1HnMPhQiJZZT9cduIhjBeX4oeQmDp6/zmibiMjF7GUzNYvArzfkQKNW8WbQSRh8KEjTrKfRYUGY/l42Pvm2AJ98W2A9h1lPiYicj9lM3YsLThWqvKbe5nFmPSUicj7mXHIvBh8KZDKL+P2/T9l8jllPiYicjzmX3IvBhwIxAicici/mXHIvBh8KxAiciMi9mM3UvRh8KBAjcCIi95JyLtnCbKbOx+BDgZj1lIjI/cbfFY0ugepWx5nN1Pm41VaBmPWUiMh9pIymnx0rRHWDCb26BWPFk0Nwrbqe+ZVchMGHQtnLehqkUWH108MYgRMROYGtjKaVtUZU1RsxbdgdHmyZb+O0i4IlJ+ixb+EkbHh+DF5JHgAAaGg0o85oxpacKzh4/jq32xIRdZCU0bTl7sLKukbmU3IxjnwonJT1NLFvd3z1fTFOFlVi3qYc6/PMeEpE5DhmNPUsjnx4icxcS+DREjOeEhE5jvmUPIvBhxeQInRbmPGUiMhxzKfkWQw+vAAjdCIi52I+Jc9i8OEFGKETETkXM5p6FoMPL8AInYjIudQqAXMm9rX5HDOauh6DDy/AjKdERM73Q8lNAEBgQPOPQmY0dT1utfUC7WU8FcEInYhIDimb6fmrVdhwqAAA8OGMUVCrBJRW1TGjqZsw+PAS9jKeAsCLE+5khE5E1A5b2Uw1agFVdQ2Yek+MB1vmfxyadnnzzTcxatQohIWFITo6Go8//jjOnDnT7BxRFJGWloaYmBgEBwdjwoQJOHnypFMb7a+aZjx9+6fD8NhQyx/LwQtlOHDuGrOeEhHZYS+bqdEk4sVPjjNXkps5FHzs2bMHc+bMQXZ2NrKystDY2IgpU6agurraes7KlSvx1ltvYc2aNTh8+DB0Oh2SkpJQVVXl9Mb7Iynj6bRhd+D/PToQASoBxwvKkfr+t/jNxhxMfy8b41bs5B8SEdEtcrOZ8sbNfRwKPjIzMzFz5kwMHjwYQ4cOxdq1a1FQUICjR48CsIx6rF69Gq+//jqefPJJJCQkYN26daipqUFGRoZLfgB/duzSDTTa+GNh1lMiotuYK0l5OrXbpaKiAgAQGWnZZZGfnw+DwYApU6ZYz9FqtRg/fjwOHDjQmW9FLTDrKRGRPMyVpDwdXnAqiiIWLFiAcePGISEhAQBgMBgAAD179mx2bs+ePXHp0iWb16mvr0d9fb31cWWlpX6J0WiE0WjsaPNskq7n7Ot6wqH8MpTdrIVWDRjNgFkUoBJEaJqEk2U3a5F9rrRDW3B9qa9cjX0lH/vKMewv+drqq6iQAGjVIkQRaDADgAC1IKLFDltEhQT4RV+76nXlyPUEURQ7dGs8Z84cfPXVV9i3bx969eoFADhw4ADGjh2LoqIi6PW3d188//zzKCwsRGZmZqvrpKWlYfHixa2OZ2RkICQkpCNN8ztF1cCK7wMgQMSrw0zoGezpFhERKc+hqwI+OadGaICIRSNMCFR7ukW+paamBqmpqaioqEB4eHib53Yo+Jg7dy4+//xz7N27F/Hx8dbjFy5cQN++fXHs2DEMHz7cenzatGno2rUr1q1b1+patkY+YmNjce3atXYb7yij0YisrCwkJSVBo9E49drudii/DM+uO2x9bDQBZrQe/fhwxqgOj3z4Sl+5GvtKPvaVY9hf8rXXV9vzDHhpw3cQW4x6SNk8Vj09DJMH9mz1db7IVa+ryspKREVFyQo+HJp2EUURc+fOxebNm7F79+5mgQcAxMfHQ6fTISsryxp8NDQ0YM+ePVixYoXNa2q1Wmi12lbHNRqNy/7YXHltdxnTLxqRocEwVNQ1W8FtFgXUmyx/ULqIIIzpF92pZDm+0Ffuwr6Sj33lGPaXfC37SkoqdqSgEiIECABMogCTyfK8PiIIi1IG+WWuJGe/rhy5lkPBx5w5c5CRkYEtW7YgLCzMusYjIiICwcHBEAQB8+bNw7Jly9C/f3/0798fy5YtQ0hICFJTUx37KahN7WU9BZj1lIj8m62kYiGBavzqwTvRJ6oLs5l6kEO7XdLT01FRUYEJEyZAr9db/23atMl6ziuvvIJ58+bhxRdfxMiRI3HlyhVs374dYWFhTm+8v5OynupsVGZ8blwf1DeamXSMiPySvaRi1Q0mrN5xFtoAFRL7dmfg4SEOT7u0RxAEpKWlIS0traNtIgckJ+iRNEiHQ/llKK2qw8rMM7hSXov39120nuPPw4pE5H/kJhVLGqRj8OEhrGrrA6Ssp9oAFa6U17Z6nknHiMifMKmY8jH48BFMOkZEZMGkYsrH4MNHMNInIrKIDmu9Dq4z55HzMfjwEYz0iYgsRsdHQm9jIb5EgGUtXEdyIJFzMPjwEYz0iYgs1CoBPx/T2+Zz0vJSpiLwrA7XdiFlkSL9lknHJFLSMUb6ROTLDuWX4VpNI7adLAEABGvUqDWarM/ruPtPERh8+Ij2ko6JYKRPRL5rxylLsPHsusOoN91+n3t16gD07xmO0qo6JhVTEE67+JC2ko49co8OEcGB2JJzhYnHiMinZOYWY/6mHJvPLfoiDxW1DZg27A4mFVMQjnz4mJZJx86XVuOdnWfx9QkDvjphsJ7HxGNE5AtaJhSzdV/FhGLKw5EPHyQlHZs27A7crQ8F0HoahonHiMgXtEwzYDI3f55pBpSJwYcPM5lFLPn3KZvPMfEYEfmCpukDCm4CZtge3WCaAWVh8OHDmHiMiHxd0/QBWwvtf6QxzYCyMPjwYUw8RkS+TkozIIpAXrkKLSeZmVBMmRh8+DAmHiMiXyelGWi8tdaj6ZpSJhRTLu528WFMPEZEvspkFq27+spuNsAMASpBRIAANNx6w2NCMeVi8OHD2ks8BvCOgIi8T2ZuMRZ/mddqTduYaBFzHhmNazWNTCimcJx28XH2Eo8JAJ5/IB71jWYmHSMir5GZW4zZ64/ZWEwvYsodZlTWGZlQzAtw5MMPNE08VlJZi+Vbz8BQWYf//Sbfeg6TjhGR0rVMKNaUSgC6aYHlW09jSsIdDDwUjiMffkJKPBakUcNQ2Xp3C5OOEZHStZU+IOBWrGGoZPoAb8Dgw49Idw22MOkYESldW2kBBEHeeaQMDD78CJOOEZE3Y/oA38Hgw48w6RgReTMpfUBbdOFMH+ANGHz4Ed41EJE3k9IH2CLNuvx26t1cbOoFGHz4Eemuwd6fJdMQE5FSmcwiDp6/jur6RoQHa1o93zPcctM0eWBPdzeNOoBbbf1Ie0nHRNxOOmY2eaCBREQ22EoqphKAFyf0Rf+eYYgOC8LwXmHYlrnVg60kR3Dkw8/YSzoGAHfrwhARHIgtOVe46JSIFMFeUjGzCPx113loA1RMKOaFOPLhh5omHSutqgNEYMGn3+G0oQrT38sGAGjVIlaOBnacKsHUIb083GIi8kdtJRWTLP4yD0mDdG5rEzkHRz78lJR0bNqwO6DVqOzm9pi/KYeJx4jII5gewHcx+PBzbSUekzDxGBF5AtMD+C4GH37O3p2FeCvW4J0FEXkK0wP4LgYffs7eHUPLgQ7eWRCRu7WXVIzpAbwXgw8/Z++OoVEEGkztn0dE5GxSTo9/f1+EJDt5O6S9LVJ6APIu3O3i56Q7C0NFXYsV5QL2GATeWRCRW9nK6QEAGrUAo+n2u5QuIgiLUgYhOUHv7iaSEzD48HNtJR7bcUUFswg8nGDZljs6PpJ3GETkMlJOD1vL240mEfMn90efqC6IDgvi+5GX47QL2Uw8JkBEnUmASQQ+2H8R09/LxrgVO7ntlohcor2cHgKAjYcL8eiQGCYV8wEMPgiAJQDZt3ASNjw/Bs+O7YOAW6+MJqOcMFTUYfb6YwxAiMjpmNPDvzD4ICu1SsDo+EhszTVAJQB3R5iBJmXopDiEeT+IyNmY08O/MPigZprefaTEmdGy/BzvPojIFZjTw784HHzs3bsXKSkpiImJgSAI+Pzzz5s9P3PmTAiC0OzfmDFjnNVecrGmdxW9ulgqR7Z3HhFRZzGnh39xOPiorq7G0KFDsWbNGrvnJCcno7i42Prv66+/7lQjyX1a3lUE2Ak+ePdBRM4k7byzhTk9fI/DW22nTp2KqVOntnmOVquFTscqg95Iuvu4cbMWACDY+Dvn3QcROZPJLOJQfhmKymuhFgSYxObTvczp4Xtckudj9+7diI6ORteuXTF+/HgsXboU0dHRNs+tr69HfX299XFlZSUAwGg0wmg0OrVd0vWcfV1f87tHBuC3nx4HAGhVIkQRaDADgAC1IGL6yBj8O6cQUaFa3BvXze/vRPi6ko995Rh/6K8dp0qwfOtpGCrrYDQDZlGANkDAiw/2QZ+o0GbvM231gz/0lbO4qq8cuZ4gimKHty0IgoDNmzfj8ccftx7btGkTQkNDERcXh/z8fLzxxhtobGzE0aNHodVqW10jLS0NixcvbnU8IyMDISEhHW0aOdnBEgEbL6gRrBbxxnATumg83SIi8iUXq4BVuQEQIOJ/hphwRxdPt4gcVVNTg9TUVFRUVCA8PLzNc50efLRUXFyMuLg4bNy4EU8++WSr522NfMTGxuLatWvtNt5RRqMRWVlZSEpKgkbDT8+2SH0VeddIlNWakH+tGm/vPA/x1uiHlAdEGvNY9fQwTLZTg8HX8XUlH/vKMb7cXyaziIdW74Whsg6iCBjNgAgBKkGERmV5b+kZHoRt8x6UNbrqy33lbK7qq8rKSkRFRckKPlyeXl2v1yMuLg5nz561+bxWq7U5IqLRaFz2AnLltX3NmH7RUKkDMG7FToi3Qg2TKMDUpOicAOD3X53BlIQ7/HoKhq8r+dhXjvHF/jpy/jou3ahH01xCgGXapf7W+8ulG/U4frkKiX27y76uL/aVqzi7rxy5lsvzfFy/fh2FhYXQ67lQyFsx8yARORuTivk3h0c+bt68iXPnzlkf5+fnIycnB5GRkYiMjERaWhqeeuop6PV6XLx4Ea+99hqioqLwxBNPOLXh5D58kyAiZ2NSMf/mcPBx5MgRTJw40fp4wYIFAIAZM2YgPT0dJ06cwMcff4zy8nLo9XpMnDgRmzZtQlhYmPNaTW7FNwkichZpW62hohZdgzUor7W9Q0KAZYstt/X7JoeDjwkTJqCtNarbtm3rVINIeaTcH4aKOpsVJ/kmQURyZOYWY/GXeW1O4wJMKuYPWNuF2tU086CttwERwE9HxeLf3xfh4PnrLDpHRK1k5hZj9vpj7QYegOVmJv3nI5hUzIe5fLcL+YbkBD3Sfz7C5l1LYIAKq3bc3s2kZzZCImrCZBax+Ms8myOnksguGrzx6GDowi2jqBzx8G0MPki25AQ9kgbpcCi/DKVVddh39ho+PXoZDY3mZucZKuowe/0x3rkQEYD2d8wBQFm1EbrwIIe21ZL34rQLOUStEpDYtzseHRKDfeeu2TxHurtZ/GUep2CIiDvmqBUGH9QhzP1BRHJxxxy1xOCDOoR3MkQkl7Rjzh4BrJbtbxh8UIfwToaI5JDyeiQPtl37idtq/RMXnFKHMPcHEbVHTl4PHXfH+SUGH9QhUu6P2euPQQBaBSAigIcTLDtjuG2OyP9IeT3sLTl/bmwfTB6k4/uDn+K0C3WYlPtDZ2cu94P9FzH9vWyMW7ETmbnFbm4dEXlKe3k9BABf5xoYePgxBh/UKckJeuxbOAkbnh+DZ8f2sXmOlPeDAQiRf+BuOGoPgw/qNLVKwOj4SGzNNdh8nnk/iPwLd8NRexh8kFPwToeIJNwNR+1h8EFOwTsdIjKZRRw8fx2Gilp0DdHYPY95PYi7XcgpeKdD5N/kbKsFmNeDLDjyQU4h5f1o660ksosGhso6HDx/nWs/iHyItK22vcADsOT1YNFJ4sgHOUV7eT8AS9XK+ZtyAFiGXJlYiMj7tbetFrDceLzx6GDowoO4vZYAcOSDnKi9vB9NcfstkW9ob7E5YLnx0IUHIbFvdwYeBIAjH+RkyQl6JA2yZDY1VNRiyVenUFbd0Oo8EZa538Vf5iFpkI5vSEReiovNqSM48kFOp1YJSOzbHbqIYJuBh4Tbb4m8HxebU0cw+CCX4R0Rke+TFpvbw221ZAuDD3IZ3hER+TaTWcSh/DJMuKuHzee5rZbs4ZoPchnpjshQUWd3JXzXYA3MogiTWeSbE5EXkZPXQ8ddbWQHgw9yGTnbb8trjfjZ+99y6y2RF5Hyeti7qXhubB9MHqTjtlqyi9Mu5FJyt99y6y2Rd2gvr4cA4OtcAwMPahODD3K55AQ99i2chE+euw9dg23Xe2DlWyLvwCKS5AwMPsgt1CoBKpWA8lqj3XP4pkWkfNzFRs7ANR/kNnzTIvJe0s6WsyVVss7nLjZqC4MPchtuvSXyTnIr1gKWNR865vWgdnDahdyGlW+JvI8jFWuZ14PkYvBBbiNtvQVgNwCRKt9Ofy8b41bs5O4XIg+SU7G2KV1EENJ/PoJb5qldDD7IrVj5lsh7yKlYCwAvTeyHDc+Pwb6Fkxh4kCxc80Fux8q3RN5B7uLv/j1Dkdi3u4tbQ76EIx/kEax8S6R8XCROrsLggzyK22+JlMlkFmE2i3YTAwKsWEsdx2kX8ijeWREpj5yttdzZQp3BkQ/yKDnbb5tWviUi15K7tZY7W6gzHA4+9u7di5SUFMTExEAQBHz++efNnhdFEWlpaYiJiUFwcDAmTJiAkydPOqu95GPkbL+VKt9y6y2Ra8nZWts1WINPfnkfd7ZQpzgcfFRXV2Po0KFYs2aNzedXrlyJt956C2vWrMHhw4eh0+mQlJSEqip5KXnJ/7DyLZEyyNlaW15rhEoQONVCneLwmo+pU6di6tSpNp8TRRGrV6/G66+/jieffBIAsG7dOvTs2RMZGRl44YUXOtda8lnS9tvs89cxJ+OYzQJ03HpL5FpcAE7u4tQFp/n5+TAYDJgyZYr1mFarxfjx43HgwAGbwUd9fT3q6+utjysrKwEARqMRRqP9CqgdIV3P2df1RR7rK9GE2oYGaNWAKAINZgAQoBJEaG6N05XdrEX2uVLFrLDn60o+9pVj3NVfJrOIo5du4LyhHIEqEUYzIEKAAMvfndAizo8KCVDc75CvLflc1VeOXE8QRbHDq/gEQcDmzZvx+OOPAwAOHDiAsWPH4sqVK4iJibGe96tf/QqXLl3Ctm3bWl0jLS0NixcvbnU8IyMDISEhHW0a+YhzFcCaPDVECJjR34QRUVx0SuRK/7qowp5iFUIDRCwcakJ4oKdbRN6ipqYGqampqKioQHh4eJvnumSrrdAiTBZFsdUxyauvvooFCxZYH1dWViI2NhZTpkxpt/GOMhqNyMrKQlJSEjQa+3vXyXN9dSi/DM+uO9zsmEoATCKw7qwKG85b7sJeeOBO3Hdnd9wb183j0y98XcnHvnKMq/trx6kSzN+UY11gahKBRrPl76nODPwhR209V/orW/X0MEwe2NPpbeksvrbkc1VfSTMXcjg1+NDpdAAAg8EAvf72KujS0lL07Gn7xarVaqHValsd12g0LnsBufLavsbdfTWmXzQiQ4NhqKizseJeuDUNA7yzOx/YnQ99RBAWpQxSxKp7vq7kY185xhX9ZTKL+P1XZ1Bnsh28N5oFNDZ5rKS/tbbwtSWfs/vKkWs5Nc9HfHw8dDodsrKyrMcaGhqwZ88e3H///c78VuSj5Gy9bYo7YIg6hkXjyJMcDj5u3ryJnJwc5OTkALAsMs3JyUFBQQEEQcC8efOwbNkybN68Gbm5uZg5cyZCQkKQmprq7LaTj3Kk8q00OrL4yzwmISNygKNF4zw9vUm+xeFplyNHjmDixInWx9J6jRkzZuCjjz7CK6+8gtraWrz44ou4ceMG7rvvPmzfvh1hYWHOazX5vKaVb/efu4o1u87bPbdp8TlW1iRqm8ks4lB+Gc6WyMu9xNIG5AoOBx8TJkxAWxtkBEFAWloa0tLSOtMuImvlW+YeIHIOOTVbJAIsKdSVsqWdfAsLy5Hiyb3zulZVD5NZ5PAwkQ1SzRY5k5MsGkeuxsJypHhyis8BwJKvTrH+C5ENcmq2NMWiceRqDD5I8RzZAcPdL0StcWcLKQ2DD/IKcnfAcPcLUWvc2UJKw+CDvEZygh77Fk7CG48MbPO8prtfiPydySziWlV9+yeCO1vIfbjglLyKWiUgKqx1RlxbuPuF/J3c3S3c2ULuxpEP8jpy787OltzEwfPXOf1Cfkna3SIn8AC4s4Xci8EHeR25u1/W7DqH6e9lcwcM+R1HdrdwZwt5AoMP8jqs/0LUNrm7W954ZCB3tpBHMPggr8T6L0T2yV3vFBWm5VQLeQQXnJLXYv0XouZYt4W8BYMP8mqO1n/ZemvqZXR8JO/4yKewbgt5EwYf5BPk3sF9fPASPj54CfqIICxKGcS5bvIJrNtC3oZrPsgnyN0BI+EiVPIVrNtC3ogjH+QTpB0ws9cfgwC0+0YswnIHuPjLPCQN0vEOkLyWI3VbxvaL4pQjKQJHPshnOLIDBmAadvJ+JrOI/eeuyTqXdVtISTjyQT6l6Q6YrbnF+PjgpXa/Zv+5q7wbJK/jyAJTgDtbSFk48kE+R9oBM1XmnPaaXeeZBZW8itzU6YBlelHPnS2kMAw+yGc5sgiVC1DJWziywJQ7W0ipGHyQz3IkDTuzoJK3kLvAFODOFlIuBh/k0xxNw84FqKRUJrOIg+evWxPlteeliX1Zt4UUiwtOyedJi1BXZf2ANbvOtXs+s6CS0uw4VYLff3VG9ogHAIzt14OvX1IsjnyQX1CrBIztFyXr3I8PXsL097K5CJUUY/6mHNmBBxeYkjdg8EF+g1lQydtI64/krkLiAlPyFgw+yG84sgAV4CJU8ryjl244dD4XmJK3YPBBfoVZUMlbmMwivr1wHQAgthP7/iIxDhueH8MFpuQ1GHyQ30lO0GPfwknY8PwY/CIxTtbX7D93laMf5DaZucUYt2In/v7NBZhFoLGdl97UBD1Tp5NXYfBBfolZUEmpmmYvFUVg43kVzKLtoIKLS8lbMfggv8YsqKQkLbOXNorAt1dVsLXklItLyZsx+CC/xiyopARSArFVWc1zeZhFAQJEBNh4cXJxKXkzJhkjvyctQpVTIVRagLoq6weM7RfF4W7qtPaq0/5/d5rx+SUVGk23j700sS/mJw3giAd5LY58EOH2ItSXJvaTdf6aXeesich2nCpxcevIV7VXnTZAEHF/z9ajbMxeSt6OwQfRLY5kQZUYKuowf1OOaxpEPk1OdVp1i3doLjAlX+G10y4mkwlGo9GhrzEajQgICEBdXR1MJlP7X+DHXNVXGo0GarXaaddzNmkBqqGiTlZWyabnmMwiNK5qGPkcR6rTAlxgSr7F64IPURRhMBhQXl7eoa/V6XQoLCyEIPCPty2u7KuuXbtCp9Mp8ncgLUCdvf4YBMhLay2dk/HtJcwY148fDNQuk1nE/nPXHPoaXUQQFqUM4gJT8gleF3xIgUd0dDRCQkIc+gAzm824efMmQkNDoVJxxqktrugrURRRU1OD0tJSAIBer8w3UUcWoDa1YtsZvLe/gB8Q1Kb2Fpja8uGMURjTL5qBLfkMrwo+TCaTNfDo3r27w19vNpvR0NCAoKAgBh/tcFVfBQcHAwBKS0sRHR2t2CmY5AQ9kgbpcCi/DPvPXcWaXeftnts09bWUC4RbIMkWaYGpnBE1AYAuPAhANUbHRzLwIJ/i9E/gtLQ0CILQ7J9Op3PKtaU1HiEhIU65HnmG9PtzdM2Ou0lZUOcnDWgzEVmjGWi4tSyGuUCoJSmHx+Zjl/Ha5lzZgQcA/Hbq3a5sGpHHuGTkY/DgwdixY4f1sbPvbpW4VoDk87bfX3vrQMwQ8O4ptXUExFYuEN61+qeOTLEAt9d3/GhAFL7Od1HjiDzIJcFHQECA00Y7iJSg7XUgIvKrLJkom1qz6xzW7DoHPRcK+iVHpliaappATOmjg0Qd5ZLg4+zZs4iJiYFWq8V9992HZcuW4c4777R5bn19Perr662PKysrAViG5Fv+4RmNRoiiCLPZDLPZ7HC7xFu3ptI1yD5X9pXZbIYoijAajYpd82HLjwZEYUL/B5Dx7SWs2HbGejxAEKFRAeUNlnERjQpoOtBx42Yt5m04ilVPD8PkgT3d33AFkf6mff1D1WQW8eZXJxGoFpsck6rTWgJVjQqwNQiYGN8NZlMjzCb/6S9nYF/J56q+cuR6giiKTp2Y3rp1K2pqanDXXXehpKQEf/jDH3D69GmcPHnS5iLRtLQ0LF68uNXxjIyMVms7pBGV2NhYBAYGOrPZLvfiiy9iw4YNACw/R7du3TB48GA89dRTSE1Nlb2oMyMjA6+++iouXbrkyua6VENDAwoLC2EwGNDY2Ojp5jhFRQPw3mk1CqsFqAURP+tnxr1RXPNBgFkEMi+rsO2y5W/83igzUvuaEcA17+RjampqkJqaioqKCoSHh7d5rtODj5aqq6vRt29fvPLKK1iwYEGr522NfMTGxuLatWutGl9XV4fCwkL06dMHQUFBDrdFFEVUVVUhLCwMZhE4fLEMpVX1iA7TYlQf187L/9d//RdKSkrw4YcfwmQyoaSkBNu2bcPy5csxbtw4bNmyBQEB7Q9EffTRR1iwYAHKyspc1lageV85e41GXV0dLl68iNjY2A79HpVgx6kSa2bTQJWIJSPN+H+HVahuFGC+tVxQLYhQC83vbhc+NACp98X57RoQo9GIrKwsJCUlQaPxrZRsO06VYPnW0zBU3p6WE0XLguS2XhPA7QWmLUfHfLm/nI19JZ+r+qqyshJRUVGygg+Xb7Xt0qUL7rnnHpw9e9bm81qtFlqtttVxjUbTqlNMJhMEQYBKperQ9k9p+mDbyRIs+epUs7l7V8/LC4KAoKAgxMTEAABiY2MxcuRIJCYm4kc/+hE+/vhj/PKXv8Rbb72FtWvX4sKFC4iMjERKSgpWrlyJ0NBQ7N69G8899xyA24t4Fy1ahLS0NKxfvx6rV6/GmTNn0KVLF0yaNAmrV69GdHR0h9or9ZXU386kUqkgCILN37G3mDqkFwSVGou/zEPZzVoAQIN4O/AAAJMowNQitP/91z8wFwhs/317s8zcYryY8d2t9R32A0tbrwmg/fcfX+svV2JfyefsvnLkWi4f+Kuvr8epU6cUk1DqP2euY07G8VaLBqX8DJm5xW5tz6RJkzB06FD861//AmD5YH7nnXeQm5uLdevWYefOnXjllVcAAPfffz9Wr16N8PBwFBcXo7i4GC+//DIAy1TGkiVL8N133+Hzzz9Hfn4+Zs6c6dafxd9Ixeg+nDHKoa/z1GuNXENOjRZbIrtosOrpYdjw/BjsWzjJr4NR8j9OH/l4+eWXkZKSgt69e6O0tBR/+MMfUFlZiRkzZjj7WznMZBaxcscFm28SIiz3K4u/zEPSIJ1bh8XvvvtufP/99wCAefPmWY/Hx8djyZIlmD17Nt59910EBgYiIiLCZu6UZ5991vr/d955J9555x2MHj3amqWUXEOtEjA6PhJfn7IkhCq4Ud/uh5D0/GubT6DWaIYuPIjbcb2UySzio/35HarRsuyJexhwkN9yevBx+fJlTJ8+HdeuXUOPHj0wZswYZGdnIy4uztnfymGHL5ahpKrB7vNSfoZD+WVI7Ot4BtWOEkXRuq5i165dWLZsGfLy8lBZWYnGxkbU1dWhuroaXbp0sXuN48ePIy0tDTk5OSgrK7NOmxQUFGDQoEFu+Tn83W+n3o0XM76TXROmrNpoXTfC7bjep7M5PPi7Jn/m9OBj48aNzr6k05RW1bd/EoDSKsfeTDrr1KlTiI+Px6VLl/Dwww9j1qxZWLJkCSIjI7Fv3z4899xzbW5hqq6uxpQpUzBlyhSsX78ePXr0QEFBAR566CE0NNgPtsi5Jg/s2aGaMADTsnubjuTweGliPyadI7rFq2q7dFZ0WOuFrbbPc98OjJ07d+LEiROYP38+jhw5gsbGRvz5z3+2LvL8xz/+0ez8wMDAViXuT58+jWvXrmH58uWIjY0FABw5csQ9PwA1I9WE+Wh/PpZ8dUr213EqRvlMZhGH8stgqKjFkq9OyQ48BFhGO+Yn3cXfJ9EtfhV8jOoTiZ5hgSitarD5xiG9SYyOj3TJ96+vr4fBYLButc3MzMSbb76JRx99FL/4xS9w4sQJNDY24i9/+QtSUlKwf/9+/O1vf2t2jT59+uDmzZv4z3/+g6FDhyIkJAS9e/dGYGAg/vKXv2DWrFnIzc3FkiVLXPIzUPvUKgEzx8bj/X35MFTUOXR3zKkYZeroFIsUaixKGcTAg6gJv0pzo1YJeGWyJdNqy7cBd7xJZGZmQq/Xo0+fPkhOTsauXbvwzjvvYMuWLVCr1Rg2bBjeeustrFixAgkJCfjkk0/w5ptvNrvG/fffj1mzZuHpp59Gjx49sHLlSvTo0QMfffQRPv30UwwaNAjLly/Hn/70J5f8DCSPVA8GaGvjZdu4K0YZpCkWRwMPwHIzw6k0otb8auQDAH40oDv+mjq8VZ4PVy8C++ijj/DRRx+1e978+fMxf/78ZseeeeaZZo/T09ORnp7e7Nj06dMxffr0ZsdcnD+O2tF2PZj2Sb+93352AmFBGoy5szvvnt3IZBaRff46fvvZCYe30QLAG48MxMyx8fydEdngd8EHACQn6PBQgh6H8stQWlWH6DDOr5NrSGtAmq4VuFFte9rPnvJaI372/rechnGjjk6zALenbxl4ENnnl8EHYBkWd+d2WvJfTV9rwYFqzF5/TPZ23Ka4I8Z1pMWkpVV1uHitBqt3/NCh0Q6u8SCSx2+DDyJP6MxUDHfEuEZnRjlaYg4PInkYfBC5WWenYrgjxnk6kq+jpcguGrzx6GAGg0QOYPBB5AHOmooprqjDrPXH8NzYPpg8SMcPP5k6u5gUYJp0os7wq622REokTcXoIjqe3O6D/Rcx/b1sjFuxk1tz25GZW4xxK3biZx98i/Ja+5mD28NttEQdx5EPIgWQpmKyz1/HnIxjHf5QNNwaCZk/uT/6RHXhTi44bzGppGuwBn/92QhufSbqBAYfRAqhVgkY2z8Ky5+6B7PXHwPg+DSMdP6qHWetx/x5XYgzF5NKYcbyp+7B2H5Rnb4ekT/jtAuRwjhjGqYpaV3Iki9P4uD56zCZfTv5nMks4uD56/j9lycxq4OZSW3hNAuR83DkQ0EEQcDmzZvx+OOP23z+4sWLiI+Px/HjxzFs2DC3ts2WPn36YN68eZg3b56nm+JznJGcrKUP9l/EB/sv+tRISNMpleiwINyobsCSrzo/0iEt/uX0FZFrMPhwk5kzZ2LdunUAALVajZiYGDzyyCNYtmwZunXrBgAoLi62/r9STJo0CXv27Gl13Gg04vDhw+jSpYv1WHvBEznGWTtiWvKVdSHOnFJpifk6iFyLwYcbJScnY+3atWhsbEReXh6effZZlJeXY8OGDQAAnU7n4Rba9vzzz+P3v/99s2MBAQHo0aOHh1rkfzpbJ6YpW+tCdOFaTB/dW9HBiLMXjtrCxaRE7sE1H26k1Wqh0+nQq1cvTJkyBU8//TS2b99ufV4QBHz++efWx4cOHcLw4cMRFBSEkSNH4vjx462u+cUXX6B///4IDg7GxIkTsW7dOgiCgPLycus5Bw4cwIMPPojg4GDExsbi17/+Naqrq2W3OyQkBDqdrtk/wDLtsnr1auv/A8ATTzwBQRCsj8l5khP02LdwEjY8PwbPju0DoOMVc1syVNZj1Y6z+M3GHOuW3a+/L8LB89exJeeKR9aKSGs3tuRcwds7zmLs8p2Y/l42frMxB6ucHHgIt/5Ji0kZeBC5ltePfIiiiFqjSda5ZrMZtQ0mBDQ0QqXqfNwVrFFDEDr2JnXhwgVkZmZCo9HYfL66uhqPPvooJk2ahPXr1yM/Px+/+c1vmp1z8eJF/PjHP8ZvfvMb/PKXv8Tx48fx8ssvNzvnxIkTeOihh7BkyRJ88MEHuHr1Kl566SW89NJLWLt2bYfabsvhw4cRHR2NtWvXIjk5GWq12mnXptukqZjEvt0xOj7SZdMOxRV1eDGjebCrjwjCG48MRLcuWusai3vjuuHopRvNCjQCaLYOo+k5USGWtxyTWcSR89ftXsdZazfk4jQLkXt5ffBRazRh0O+2eeR75/3+IYQEyu/Cf//73wgNDYXJZEJdneVN9a233rJ57ieffAKTyYQPP/wQISEhGDx4MC5fvozZs2dbz/nb3/6GAQMG4I9//CMAYMCAAcjNzcXSpUut5/zxj39EamqqdVFo//798c4772D8+PFIT09HUFD7OyreffddvP/++9bHL7zwAv785z83O0eagunatatip498TdNFqU2nIoDOrwuxxVZAohKApgMiXUMswXR5jdHmOVq1iJWjgfF/3IWSm412r+NKXExK5HleH3x4k4kTJyI9PR01NTV4//338cMPP2Du3Lk2zz116hSGDh2KkJAQ67HExMRm55w5cwajRo1qdmz06NHNHh89ehTnzp3DJ598Yj0miiLMZjPy8/MxcODAdtv9s5/9DK+//rr1cdeuXdv9GnKPltWZB+hCXTYaYkvLgKFp0GHvHAC3kqgJbZ7jKhzlIPI8rw8+gjVq5P3+IVnnms1mVFVWISw8zGnTLo7o0qUL+vXrBwB45513MHHiRCxevBhLlixpda4otv9uLIpiq2mfll9nNpvxwgsv4Ne//nWrr+/du7esdkdERFjbTcrWdDQkK8+AD/dfdMoOGV/A+jdEyuH1wYcgCLKnPsxmMxoD1QgJDHBK8NFZixYtwtSpUzF79mzExMQ0e27QoEH4v//7P9TW1iI4OBgAkJ2d3eycu+++G19//XWzY0eOHGn2eMSIETh58qRbggeNRgOTSd76G3Idd60L6QijCfjLSTUa3Pgy8aW8JkS+wvOfwH5swoQJGDx4MJYtW9bqudTUVKhUKjz33HPIy8vD119/jT/96U/NznnhhRdw+vRpLFy4ED/88AP+8Y9/4KOPPgIA64jIwoULcfDgQcyZMwc5OTk4e/YsvvjiC7vTPZ3Rp08f/Oc//4HBYMCNGzecfn1yXNMdMm//dBjmT77LurPDE8wQcK5SgOiiFugjgvBu6nDrz7vh+THYt3ASAw8ihWHw4WELFizAe++9h8LCwmbHQ0ND8eWXXyIvLw/Dhw/H66+/jhUrVjQ7Jz4+Hv/85z/xr3/9C0OGDEF6erp1bYZWqwUADBkyBHv27MHZs2fxwAMPYPjw4XjjjTeg1zv/zfjPf/4zsrKyEBsbi+HDhzv9+tQx0kjItGF34DeT+zs1dbujAgQRM/qboFF1fiJICl/mT+7fLNB4eEiM9edN7Mt8HURK5PXTLt5CGpFoKTU1FampqQBar9cYM2YMcnJymh1rec5jjz2Gxx57zPp46dKl6NWrV7NdLKNGjWqWT8QRO3futDtFdfHixWaPU1JSkJKS0qHvQ+5ja5fMhkMFMFS6fmpGrQJGRInYeKHz1+LCUSLvxeDDy7377rsYNWoUunfvjv379+OPf/wjXnrpJU83ixSu5S6Zlyb1c0mNFGfyhiysRCQPgw8vd/bsWfzhD39AWVkZevfujf/+7//Gq6++Kutrv/nmG0ydOtXu85cvX3ZWM0nhWgYjAPBQgq7dgMTRPB/W84I1beb5sJXQjMEGke9g8OHlVq1ahVWrVnXoa0eOHNlqWkdiNps70SryBXICko5kOL12Kht7/mcijl+uavM6DDSIfBeDDz8WHBxsdwuu2WxGZWWlm1tESmcrIGn52NYx6bHRaMTXp+Rfh4h8E3e7EBERkVt5ZfDBKQHvxt8fEZF/86ppl8DAQKhUKhQVFaFHjx4IDAx0qKqs2WxGQ0MD6urqFJHhVMlc0VeiKKKhoQFXr16FSqVCYGCgU65LRETexauCD5VKhfj4eBQXF6OoqMjhrxdF0Zqu3JGgxR+5sq9CQkLQu3dvBoBERH7Kq4IPwDL60bt3bzQ2NjpcR8RoNGLv3r148MEHodFoXNRC3+CqvlKr1QgICGDwR0Tkx7wu+AAsdUs0Go3DH4pqtRqNjY0ICgpi8NEO9hUREbkKx72JiIjIrRh8EBERkVsx+CAiIiK3UtyaD6lqqyuyaxqNRtTU1KCyspLrGNrBvpKPfSUf+8ox7C/52FfyuaqvpM/tltXXbVFc8FFVVQUAiI2N9XBLiIiIyFFVVVWIiIho8xxBlBOiuJHZbEZRURHCwsKcvh2zsrISsbGxKCwsRHh4uFOv7WvYV/Kxr+RjXzmG/SUf+0o+V/WVKIqoqqpCTExMu3mcFDfyoVKp0KtXL5d+j/DwcL44ZWJfyce+ko995Rj2l3zsK/lc0VftjXhIuOCUiIiI3IrBBxEREbmVXwUfWq0WixYtglar9XRTFI99JR/7Sj72lWPYX/Kxr+RTQl8pbsEpERER+Ta/GvkgIiIiz2PwQURERG7F4IOIiIjcisEHERERuZXfBh+PPfYYevfujaCgIOj1ejzzzDMoKirydLMU5+LFi3juuecQHx+P4OBg9O3bF4sWLUJDQ4Onm6ZYS5cuxf3334+QkBB07drV081RlHfffRfx8fEICgrCvffei2+++cbTTVKkvXv3IiUlBTExMRAEAZ9//rmnm6RIb775JkaNGoWwsDBER0fj8ccfx5kzZzzdLMVKT0/HkCFDrMnFEhMTsXXrVo+0xW+Dj4kTJ+If//gHzpw5g88++wznz5/Hj3/8Y083S3FOnz4Ns9mMv//97zh58iRWrVqFv/3tb3jttdc83TTFamhowE9+8hPMnj3b001RlE2bNmHevHl4/fXXcfz4cTzwwAOYOnUqCgoKPN00xamursbQoUOxZs0aTzdF0fbs2YM5c+YgOzsbWVlZaGxsxJQpU1BdXe3ppilSr169sHz5chw5cgRHjhzBpEmTMG3aNJw8edL9jRFJFEVR3LJliygIgtjQ0ODppijeypUrxfj4eE83Q/HWrl0rRkREeLoZijF69Ghx1qxZzY7dfffd4m9/+1sPtcg7ABA3b97s6WZ4hdLSUhGAuGfPHk83xWt069ZNfP/9993+ff125KOpsrIyfPLJJ7j//vtZilmGiooKREZGeroZ5EUaGhpw9OhRTJkypdnxKVOm4MCBAx5qFfmaiooKAOD7kwwmkwkbN25EdXU1EhMT3f79/Tr4WLhwIbp06YLu3bujoKAAW7Zs8XSTFO/8+fP4y1/+glmzZnm6KeRFrl27BpPJhJ49ezY73rNnTxgMBg+1inyJKIpYsGABxo0bh4SEBE83R7FOnDiB0NBQaLVazJo1C5s3b8agQYPc3g6fCj7S0tIgCEKb/44cOWI9/3/+539w/PhxbN++HWq1Gr/4xS8g+knCV0f7CgCKioqQnJyMn/zkJ/jlL3/poZZ7Rkf6i1oTBKHZY1EUWx0j6oiXXnoJ33//PTZs2ODppijagAEDkJOTg+zsbMyePRszZsxAXl6e29sR4Pbv6EIvvfQSfvrTn7Z5Tp8+faz/HxUVhaioKNx1110YOHAgYmNjkZ2d7ZEhKHdztK+KioowceJEJCYm4n//939d3DrlcbS/qLmoqCio1epWoxylpaWtRkOIHDV37lx88cUX2Lt3L3r16uXp5ihaYGAg+vXrBwAYOXIkDh8+jLfffht///vf3doOnwo+pGCiI6QRj/r6emc2SbEc6asrV65g4sSJuPfee7F27VqoVD41YCZLZ15bZHnDu/fee5GVlYUnnnjCejwrKwvTpk3zYMvIm4miiLlz52Lz5s3YvXs34uPjPd0kryOKokc+93wq+JDr0KFDOHToEMaNG4du3brhwoUL+N3vfoe+ffv6xaiHI4qKijBhwgT07t0bf/rTn3D16lXrczqdzoMtU66CggKUlZWhoKAAJpMJOTk5AIB+/fohNDTUs43zoAULFuCZZ57ByJEjrSNoBQUFXD9kw82bN3Hu3Dnr4/z8fOTk5CAyMhK9e/f2YMuUZc6cOcjIyMCWLVsQFhZmHVmLiIhAcHCwh1unPK+99hqmTp2K2NhYVFVVYePGjdi9ezcyMzPd3xi3769RgO+//16cOHGiGBkZKWq1WrFPnz7irFmzxMuXL3u6aYqzdu1aEYDNf2TbjBkzbPbXrl27PN00j/vrX/8qxsXFiYGBgeKIESO4JdKOXbt22XwNzZgxw9NNUxR7701r1671dNMU6dlnn7X+/fXo0UP80Y9+JG7fvt0jbRFE0U9WWBIREZEi+N/kPREREXkUgw8iIiJyKwYfRERE5FYMPoiIiMitGHwQERGRWzH4ICIiIrdi8EFERERuxeCDiIiI3IrBBxEREbkVgw8iIiJyKwYfRERE5FYMPoiIiMit/n+U184AxTcB0AAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABeG0lEQVR4nO3deVzUdf4H8Nd3BpjhVkQYUFBEPBDzPlBLLSGtTHNr28yy7fKqzdrWDmuFbDWtzG1dra3N3DXT9ldmbkVSXnmtJyliHognICrKfQwz398f+B055vgOzD2v5+Pho/jOd758+DDMvL+f4/0WRFEUQUREROQgCmc3gIiIiLwLgw8iIiJyKAYfRERE5FAMPoiIiMihGHwQERGRQzH4ICIiIodi8EFEREQOxeCDiIiIHMrH2Q1oSq/XIz8/H8HBwRAEwdnNISIiIhlEUURZWRmio6OhUJgf23C54CM/Px8xMTHObgYRERG1wPnz59GxY0ez57hc8BEcHAygvvEhISE2vbZWq8WmTZuQmpoKX19fm17b07Cv5GNfyce+sg77Sz72lXz26qvS0lLExMQYPsfNcbngQ5pqCQkJsUvwERAQgJCQEL44LWBfyce+ko99ZR32l3zsK/ns3VdylkxwwSkRERE5lFXBR1paGgRBaPRPo9EYHhdFEWlpaYiOjoa/vz9GjRqFo0eP2rzRRERE5L6sHvno1asXCgoKDP+OHDlieGzx4sVYsmQJli1bhn379kGj0SAlJQVlZWU2bTQRERG5L6uDDx8fH2g0GsO/9u3bA6gf9Vi6dCnmzp2LSZMmISkpCatWrUJlZSXWrFlj84YTERGRe7J6wenJkycRHR0NlUqFIUOGYMGCBejSpQvy8vJQWFiI1NRUw7kqlQojR47Erl27MG3aNKPXq6mpQU1NjeHr0tJSAPULYrRarbXNM0u6nq2v64nYV/Kxr+RjX1mH/SUf+0o+e/WVNdcTRFEU5Z78/fffo7KyEt26dcOlS5fw5ptv4tdff8XRo0dx/PhxDB8+HBcvXkR0dLThOU8//TTOnj2LH374weg109LSkJ6e3uz4mjVrEBAQIPsHISIiIueprKzE5MmTUVJSYnG3qlXBR1MVFRWIj4/HnDlzMHToUAwfPhz5+fmIiooynPPUU0/h/PnzyMjIMHoNYyMfMTExuHLlil222mZmZiIlJYVbsSxgX8nHvpKPfWUd9pd87Cv57NVXpaWlCA8PlxV8tCrPR2BgIHr37o2TJ09i4sSJAIDCwsJGwUdRUREiIyNNXkOlUkGlUjU77uvra7cXkD2v7WnYV/Kxr+RjX1mH/SUf+0o+W/eVNddqVZ6PmpoaHDt2DFFRUYiLi4NGo0FmZqbh8draWmzbtg3Dhg1rzbchIiIiG9DpRezNKwYA7M0rhk7f4smPVrFq5OPFF1/E+PHjERsbi6KiIrz55psoLS3F1KlTIQgCZs+ejQULFiAhIQEJCQlYsGABAgICMHnyZHu1n4iIiGTIyC5A+sYcFJdXYfFg4PFV+xAW5I954xMxNinK8gVsyKrg48KFC3jooYdw5coVtG/fHkOHDsWePXvQqVMnAMCcOXNQVVWFmTNn4tq1axgyZAg2bdokK887ERER2UdGdgFmrD4IEYBKefN4YUk1Zqw+iBVT+js0ALEq+Fi7dq3ZxwVBQFpaGtLS0lrTJiIiIrIRnV5E+sYcGJtgEQEIANI35iAlUQOlwnJdFltgbRciIiIPtjevGAUl1SYfFwEUlFQb1oI4AoMPIiIiD1ZUZjrwaMl5tsDgg4iIyINFBKttep4tMPggIiLyYIPjwhAVqoa0mkOrB3YUCpBSjAoAokLVGBwX5rA2MfggIiLyYEqFgHnjEw1f60UBX52p//iXApJ54xMdttgUYPBBRETk8cYmRWHFlP4I9KvfZ5vUVoQgAJpQtcO32QKtTK9ORERE7iElUYNA1VFU1OowsL2IORMHY2jXCIeOeEg48kFEROQFduVeQVFZDdr4+yKxjYjBcWFOCTwABh9EREReYf2hiwCAu3pHwsfJn/4MPoiIiDxcZW0dMrILAQAT+kQ7uTVc80FEROSRpAq2RWXVOHmpHJW1OsSGBaBfTCi+z3Zu2xh8EBEReRipgm3TtOq9O4RAEJyzzqMhBh9EREQepGEF26a+PVKIe3pHOrxNTXHNBxERkYcwV8EWqE8q9tb3vzqySUYx+CAiIvIQcirYFpY6roCcKQw+iIiIPIQjK9O2BoMPIiIiD+HIyrStweCDiIjIQzStYNuUAEAT4vwAhcEHERGRh2hawbYhKSB5eVwPxzXIBAYfREREHkSqYOvv2/gjXqpgO6an87faMs8HERGRh7k1ob0hmdjsMQkYEtfOUEhOq9U6uXVeNPIhpZkF6rci6fSmdkETERG5tx+OFqKyVofO7QLw3B0JSI5v57QKtsZ4xciHlGb2alkVnu0FzFm1D2FB/pg3PhFjk6Kc3TwiIiKb+vLgBQDApP4dXSKdelMeP/IhpZktKKlGrR5YnqOEKAKFJdWYsfogMrILnN1EIiIim8m/XoVduVcBAPf16+Dk1hjn0cGHsTSzVToBesBwLH1jDqdgiIjIY6w/dBGiWL/tNiYswNnNMcqjg4+maWaVN0ae9Pr6/4oACkqqDWtBiIiI3JkoivjqxpTL/f07Ork1pnl08NE0zay01kZv4TwiIiJ3dPhCCXIvV0Dtq8C43hpnN8ckjw4+mqaZVQhA5yARaJL77eSlcuzOvcrpFyIicls6vYi/bzkFABgQ2xYBfq67p8R1W2YDUprZwpJqwxqPIRF6nClXNjpv2ZZTWLblFDQhKjw0OBadwwMREaw27IkmIiJyZRnZBUj75igKS2sAADtzr2LEos0uu6vTo0c+GqaZlUKIvu1EAMZHOApLa/Dejyfx3NosPPTRHoxYtJm7YYiIyKVJuzqlwEPiyrs6PTr4AG6mmdWE1k/BBPjcXPthiSv/4oiIiIzt6pS48q5Ojw8+gPoAZMdLt+OTqYMA3Nz1Yokr/+KIiIia7upsylV3dXpF8AHUT8EMjgsD0HS5qXmu+osjIiKSu1vT1XZ1ek3w0VBLMs262i+OiIio6a7O1p7nKF4ZfGhCrP8luNovjoiIaHBcGMKD/Ew+LgCIClUbRv5dhVcGHy+P6yF76sVVf3FERERKhYBe0aFGH5M+5+aNT3S5tBFeGXyM6RmJFVP6I0RtPs2JK//iiIiIaup0yDp/HQAQFuDb6DFNqBorpvR3yTwfHp1kzJyxSVEYHNcOQxb8CK1OxOTBsdj8axEKS2+u7dCEql02QQsREdFPx4pQUqWFJkSN7XNG48DZaygqq3b5RJleG3wAQFigH8YmRWHjL/nwUQrY+fLt2JtX7Ba/OCIiov/sPw8AmNS/A/x8FEiOb+fkFsnjldMuDT0woL7q34asfGh1eiTHt8OEvh2QHN+OgQcREbmsotJqbDtxGQBw/wDXrWBrjNcHH8O7hiM6VI2SKi0ycy45uzlERESyfHXoIvQiMKBTW3RpH+Ts5ljF64MPpULAb25EjP85cMHJrSEiIrJMFEXDlMtv+rvXqAfA4APAzeGqn09eRv71Kie3hoiIyLyD564j93IF1L4KjO/jfpsiGHwA6NQuEEO7hEEUga8OcvSDiIhcmzTqcVfvKASrfS2c7XoYfNzwwIAYAPVTL6LIInJEROSaKmrqsPGXfADAgwNjnNyalmHwccO43hoEqXxw9moli8gREZHL+u5IASpqdejcLsBts28z+LghwM8H99xSP2/2xX5OvRARkWv64saUywMDYyC0pFKqC2Dw0cADN4avvjtSgPKaOie3hoiIqLHTl8ux78w1KASgc7tAbMi6iN25V6HTu9dyAa/OcNpU/9g2iG8fiNzLFfjvL/n43eBYZzeJiIjIQBqZ91UqMGvNQcPxKDcrB8KRjwYEQcCDg+pHP9buO+/k1hAREd1Up9Njzf/OAgBq6vSNHissqcaM1QeRkV3gjKZZjcFHE5P6d4SPQkDW+es4Xljm7OYQERFBpxexYmsuSquNLwmQJl3SN+a4xRQMg48mwoNUGNMzEgCwjqMfRETkZBnZBRixaDPezTxh9jwRQEFJtVvs2GTwYcSDg+unXr46dAE1dTont4aIiLxVRnYBZqw+iIKSatnPKSqTf66zMPgw4raE9ogKVeN6pRabjrLYHBEROZ5OLyJ9Yw6snUSJCFbbpT22xODDCKVCwAM36r1I+6mJiIgcaW9esVUjHgLqd724Q+IxBh8mSDk/fj55BeeLK53cGiIi8jbWTJ9IqcbmjU+EUuH6iccYfJgQExaAEV3DAdws4ENEROQo1kyfaELVWDGlP/N8eAIp58d/Dlxwi61LRETkOQbHhSEyRGX2nDb+vvjsySHY8dLtbhN4AK0MPhYuXAhBEDB79mzDMVEUkZaWhujoaPj7+2PUqFE4evRoa9vpFKm9ItEmwBcFJdXYfvKys5tDREReRKkQMKp7e6OPCTf+vfWb3hjeNdwtploaanHwsW/fPvzjH//ALbfc0uj44sWLsWTJEixbtgz79u2DRqNBSkoKysrcL2GXykeJ+/p1AACs3XvOya0hIiJvIooifjlfAgAIUTeuhuJu0yxNtai2S3l5OR5++GF89NFHePPNNw3HRVHE0qVLMXfuXEyaNAkAsGrVKkRGRmLNmjWYNm2abVrtQL8bFIuVO8/gp2NFKCqtRkSI629hIiIi9/fLhRL8WlgGPx8Ftv5pNI4XlqGorBoRwfU7WtxttKOhFgUfs2bNwt13340xY8Y0Cj7y8vJQWFiI1NRUwzGVSoWRI0di165dRoOPmpoa1NTUGL4uLS0FAGi1Wmi12pY0zyTpetZct0s7NfrFhOLQ+RKs23sW00d2sWmbXFVL+spbsa/kY19Zh/0lnyf1lU4v4sDZa1i2NRcAcGdiBIL9BAyMDQEQAgDQ6+qgb2EOTHv1lTXXszr4WLt2LQ4ePIh9+/Y1e6ywsBAAEBkZ2eh4ZGQkzp49a/R6CxcuRHp6erPjmzZtQkBAgLXNkyUzM9Oq83v6CTgEJT79+SQ6lv8KNw42rWZtX3kz9pV87CvrsL/k85S+qtEBB88oAQjoXHcB3313webfw9Z9VVkpPy2FVcHH+fPn8dxzz2HTpk1Qq01PPwhC409nURSbHZO88soreOGFFwxfl5aWIiYmBqmpqQgJCbGmeRZptVpkZmYiJSUFvr6+sp83ulaHjW9vw9XqOoR2GwyVjwJXymsQHqTCgE5t3Xroy5SW9pU3Yl/Jx76yDvtLPk/oqx+PXcLz67IgAtDpgTpRgAAR//hVCUEA3nuwr6H2WGvYq6+kmQs5rAo+Dhw4gKKiIgwYMMBwTKfTYfv27Vi2bBmOHz8OoH4EJCrq5iKYoqKiZqMhEpVKBZWq+VYiX19fu72ArLm2Ti/iaGEp+sa0wc8nr2DmmixUNyhlHBWqxrzxiW676McSe/4ePA37Sj72lXXYX/K5a1/p9CLe+PY4qnVNbt4hoEZfv7PljW+PIzWpg81ueG3dV9Zcy6rdLnfccQeOHDmCrKwsw7+BAwfi4YcfRlZWFrp06QKNRtNoKKe2thbbtm3DsGHDrPlWLkGqJPjQR3vw88krANAo8ACAwpJqzFh9EBnZBc5oIhEReQBLqdTdqWKtHFaNfAQHByMpKanRscDAQLRr185wfPbs2ViwYAESEhKQkJCABQsWICAgAJMnT7Zdqx1AqiRoKbWYiPqINH1jDlISNR45BUNERPYlN5W6O1SslaNFu13MmTNnDqqqqjBz5kxcu3YNQ4YMwaZNmxAcHGzrb2U31lYSbBiRJse3s2fTiIjIA8lNpe4OFWvlaHXwsXXr1kZfC4KAtLQ0pKWltfbSTmNtJUGJp0SkRETkWIPjwhAVqjb52SOgPrGYO1SslYO1XYxoaRDhKREpERE5llIhYN74RKOPuVvFWjlsPu3iCawNIjwtIiUiIsfr2NZ4biuNB+6qZPBhhDT8VVhSbXHdhydGpERE5HhrbtQQu7t3FKYM7eQxqdSN4bSLEQ2Hvyz9uiNDVPj75H4I9ffDhqyL2J17FTq93KWqREREQHlNHTYcuggAeHhoLJLj22FC3w5Ijm/ncYEHwJEPk8YmRWHFlP5I35jTaAGQJkSF+wd0xMc78lCt1eO3A2Mw/9tjjc7x9MRjRERkOzq9iKWZJ1BRq0NUqBqDO3v+FD6DDzPGJkUhJVGDvXnFzYa/yqrrsGr3Wby/+VSz50mJx9y53DEREdlfRnZBo5vcgpJq3Lp4i8ffwHLaxQKlQjA6/PXgoFiTz5EmXdI35nAKhoiIjJKSWTbdXusNmbMZfLRQSZX50sGelgqXiIhaT6cXsTv3KtYfvIBX12cb3dTgDTewnHZpIW9LhUtERK3TdIrFHE/PnM3go4W8LRUuERG1nNx6YU156g0sp11aSMoFYoqA+l0vTDxGROTdrK0X1pCn3sAy+Gghb0uFS0RELdOSemGefgPL4KMVxiZF4YMp/eHn07gbNaFqbrMlIiIA1k+deMMNLNd8tNLYpCgs/S0wc81BBKl8sPzh/hjeNdxjXzBERGQda6dOPLGWS1MMPmwgpVckIkNUuFRag2uVtQw8iIjIQE69sLBAX7x+Ty9oQjyzlktTnHaxAV+lApMHdwIA/Hv3WSe3hoiIXImlemECgAX39cZ9/Ty3lktTDD5s5KHBMfBRCNh/9hpy8kud3RwiInIhUr0wTZNdkpoQlVeuEWTwYSMRIWrcmaQBAPx7D0c/iIiosbFJUdjx0u3o0zEUAHDPLVHY+fIdXhd4AAw+bOqRofVTL18fuojSavPp14mIyPucL67E4YslAIAXU7t7xRSLMQw+bGhIXBi6RQahSqvDlwcuOLs5RETkYlbvOQtRBEZ1b4/O4YHObo7TMPiwIUEQ8EhyZwDAR9tP4+tDF7A796rHFgYiIiL5Kmvr8MX+8wCAR5M7Obk1zsWttjYWrFJCAJBfUo3Z634BUJ+lztP3bBMRUWM6vYi9ecUoKqtGRLAap6+Uo7S6DrFhARjZLcLZzXMqBh82lJFdgOfX/dJsH3dhSTVmrD7olSuaiYi8kbEKtj431nc8MrST1671kHDaxUbMFQ6SjqVvzOEUDBGRh5Mq2Dat51J34/0/LNDXGc1yKQw+bMRS4SARQEFJNfbmFTuuUURE5FByKti+s+mE19+IMviwEbmFg6wtMERERO5DTgVb3ogy+LAZuYWDrC0wRERE7oM3ovIw+LARqXCQuSVEUaH1BYOIiMgz8UZUHgYfNmKpcBAAzBuf6PUrnImIPJmlG1EBvBEFGHzYlKnCQQAwqFNbhPr7YUPWRSYeIyLyUOZuRKWveSPKPB82NzYpCimJGkNimYoaHV5dfwT7zl7DQx/tMZzHxGNERJ5JuhF97etsXCmvNRzX8H3fgCMfdqBUCEiOb4cJfTuY3M8tJR7LyC5wcOuIiMjexiZF4bZu7QEAAzq1xedPDcWOl25n4HEDgw87kvZ7G8PEY0REnutKeQ3++0v9zeVrd/dEcnw7r59qaYjBhx0x8RgRkXf6/H/nUKvTo09MG/SLbevs5rgcBh92xP3eRETep7ZOj3/vOQsAeHx4Z+c2xkUx+LAj7vcmIvI+32cXoKisBu2DVRjHNR5GMfiwI+73JiLyPp/uOgMAmDKkE/x8+DFrDHvFjph4jIjIu2Sdv45D567DVylg8pBYZzfHZTH4sDNTiccCVUqsmNKf266IiDzIqhujHuNviUb7YJVzG+PCmGTMARomHss8VohPdpyBn1KBUd0jnN00IiKykYLrVfjml3wA9bk9dHqRI9smcOTDQaTEY3PvSkTHtv64VqnF14cuOrtZRERkAxnZBUh5b5shb9Pcr7MxYtFmJpI0gcGHgykVAh4b1hkA8MnOPIgiE4wREbmzjOwCTF99EOU1ukbHmcnaNAYfTvDAwBgE+Clx4lI5duVedXZziIiohZjJumUYfDhBqL8vHhjQEQDwyY48J7eGiIhaipmsW4bBh5M8NjwOAPDTr0XIu1Lh5NYQEZE5Or2I3blXsSHrInbnXjWMZDCTdctwt4uTxIYFoF9MGxw6fx0LvzuGFVMGcFU0EZELysguQPrGnEYjHFGhaswbn8hM1i3EkQ8nyMguwIhFm3Ho/HUAwKacSxi28CcuSiIicjEZ2QWYsfpgs6kVaTFpTv51s89nJmvjOPLhYNILuenSo0tlNZi++iCeH5OAzuGBiAiuf7FyNISIyDmkxaTGloqKqA8s3tl00nBMABqdK717M5N1cww+HMjcC1ny3o83X8jSsB6zoBIROZ6cxaRV2vrttX9M7YY1/zvX6HwN38NNYvDhQJZeyE1Jw3pMw05E5HhyF4nGtPXHM6O7YuaortibV4yismqOXlvA4MOBrF3tLA3rpW/MQUqihi9iIiIHkrtI9N6+HSAIApQCkBzfzs6t8gxccOpALVntzD3iRESOJW2rLSypQlign8mq5ACgEIBnb+/qsLZ5Co58ONDguDBEhapRWFJtdt2HMdwjTkRkf8a21ZpzX78OUPsq7dwqz8ORDwdSKgTMG58IAGYjaWO4R5yIyL5Mbas1Re2rwLx7e9m5VZ6JwYeDjU2Kwoop/aEJlRdMcI84EZH9ydmNGBboi/ce7IvEqBAAwCNDOyFE7euYBnoYTrs4wdikKKQkagyros9cqcTSH080e9FzjzgRkWPI2Y1YXKFFVa0OOQWl8FEI6NOxDTZkXeTOlhZg8OEkSoXQaFV0d00Q0r45isLSGsMx7hEnInIMuevqvjx4AQDgoxTwzOeHDMeZl8k6Vk27rFixArfccgtCQkIQEhKC5ORkfP/994bHRVFEWloaoqOj4e/vj1GjRuHo0aM2b7QnGpsUhZ0v34HZYxIAAEEqH2x6/jakJGqMFjMiIiLbkbuu7uC5awCAaq2+0XEpLxPLZMhj1chHx44d8dZbb6Fr1/ptRatWrcKECRNw6NAh9OrVC4sXL8aSJUvw6aefolu3bnjzzTeRkpKC48ePIzg42C4/gCdRKgQ8e3sCvj50EWeuVuLN/+Zg+8krRosZMbomIrIdS7sRBQD+fkpU1uqMPp95maxj1cjH+PHjcdddd6Fbt27o1q0b/vKXvyAoKAh79uyBKIpYunQp5s6di0mTJiEpKQmrVq1CZWUl1qxZY6/2exylQsCTt3YBAKzbf8FkMSNG10REtmNuN6JUs6VOZ37kmXmZ5GvxbhedToe1a9eioqICycnJyMvLQ2FhIVJTUw3nqFQqjBw5Ert27bJJY73Fff06wFTQLL300zfmcAqGiMiGTO1G1ISqcV+/aNTq9Cae2RjzMllm9YLTI0eOIDk5GdXV1QgKCsL69euRmJhoCDAiIyMbnR8ZGYmzZ8+avF5NTQ1qam4usiwtLQUAaLVaaLVaa5tnlnQ9W1/X1rLOFkO4MYgnQISvAhCaBCPF5VXYc6rIbltw3aWvXAH7Sj72lXXYX/LZqq/u6B6OUQm34sDZa7hSXoPwIBWSokMwZukOAICPIEJp4bY9PMDHpX9n9npdWXM9QRRFq26fa2trce7cOVy/fh1ffvklPv74Y2zbtg3Xr1/H8OHDkZ+fj6iom+sRnnrqKZw/fx4ZGRlGr5eWlob09PRmx9esWYOAgABrmuZRKrRA2kElavUCZibq0D2UoxxERM6w+5KAtaeVaOMn4vV+OvgwQ5ZRlZWVmDx5MkpKShASEmL2XKuDj6bGjBmD+Ph4vPTSS4iPj8fBgwfRr18/w+MTJkxAmzZtsGrVKqPPNzbyERMTgytXrlhsvLW0Wi0yMzORkpICX1/XTQyzN68Yj6/ahzo9oBPrRz/8jGTv/WTqILuOfLhDX7kC9pV87CvrsL/ks1df6fQixr2/E3lXK/HK2G6IDfPH8+uyAKDRwlRpcPq9B/tiTM/IppdxKfbqq9LSUoSHh8sKPlqd50MURdTU1CAuLg4ajQaZmZmG4KO2thbbtm3DokWLTD5fpVJBpVI1O+7r62u3PzZ7XtsWhnaNQFiQv2GxqQgBNQ0WWAuon4Mc2jXC7iuqXb2vXAn7Sj72lXXYX/LZuq9+yi5A3tVKhKh98HByHIJUPhAUymb1X9xxJ6Kt+8qaa1kVfLz66qsYN24cYmJiUFZWhrVr12Lr1q3IyMiAIAiYPXs2FixYgISEBCQkJGDBggUICAjA5MmTrf4hvJm06nrG6oPMekpE5CSiKGLFttMAgEeTOyNIVf+R2TRLNTOcWs+q4OPSpUt45JFHUFBQgNDQUNxyyy3IyMhASkoKAGDOnDmoqqrCzJkzce3aNQwZMgSbNm1ijo8WkFZdv/Z1Nq6U1xqOM+spEZFj7DldjF/OX4efjwKPDe/c6LGmWarJOlYFH//85z/NPi4IAtLS0pCWltaaNtENUnR93/KdOHyhBCmJkfhgygBG10REDvDh9lwAwAMDOiI8qPnyAGo5rtl1cUqFgFfG9QQAbD9xGcUVtRaeQURE1tDpxWZlLI4VlGLr8ctQCMDTt3VxdhM9DgvLuYGhXcLQL7YNDp27jpU78zBnbA9nN4mIyCNkZBcYXTwa09YfADCudxQ6tQt0VvM8Fkc+3IAgCJgxMh4A8O/dZ1Fa7brJa4iI3EVGdgFmrD7YrIxFQUk19p6pLyA3/bZ4ZzTN4zH4cBNjekYiISIIZTV1WPjtMVa5JSJqBZ1eRPrGHKNF5CQqHwUSo22bb4rqcdrFTSgUAoZ3bYeTReX4fN95fL7vPAD33FtOROQsOr2IvXnF2HnqcrMRj6Zq6vTYm1fMXS12wODDTWRkF+DTXc1r5EhVbldM6c8AhIjIDGPrOyxhkTj74LSLG5CGB41hlVsiIstMre+wJCJYbfkkshqDDzewN6/Y7B+MiBsLpPKKHdcoIiI3IWd9hzFRoWq71c/ydgw+3IDcYT8ODxIRNWfpBs4YASxjYU8MPtyA3GE/Dg8SETVn7Y2ZJkTFdXR2xgWnbmBwXBiiQtUoLKk2OWzI4UEiIuOsuTF7+tYueGlcD4542BlHPtyAVOUWuFnVtikODxIRGSfdwFl6h9SEqvHind35XuoADD7chFTlVhPaPIKPCw/Enb00TmgVEZHrk3MDBwDTbusCPx9+LDoCp13ciFTldm9eMYrKquGnVOC5dVnIu1KB3blXMaxruLObSETkkqQbuKZ5PkL9fVBSVYfwIBUeGhzrxBZ6FwYfbkapEBpl29tz+ipW7T6LZVtOMfggIjKj6Q1cu0AVXl1/BCVVdZh2WxeofZXObqLX4PiSm3t6ZDx8FAJ25V7FgbPM80FEZI50AzehbwcUlVXjXHElwgL98PBQjno4EoMPN9ehjT9+078jAOBvm085uTVERO5BpxexbEv9e+YTI+IQ4MeJAEdi8OEBZoyKh1IhYOvxy8g6fx06vYjduVdZ+ZaIyITvjhTg9OUKhPr74tHkTs5ujtdhqOcBOocHYkLfaHx18CJe//oIrpTXNlpQxcq3REQ36fUilt0YKX58eByC1b5ObpH34ciHh3j29gQIAnDkYmmzNMJS5duM7AIntY6IyHVsyrmE45fKEKzywWPDOzu7OV6JwYeHiA0LgNrE/nRWviUiqieKIv62+SQAYOqwzgj156iHMzD48AA6vYhPd+ahSqs3eY5U+fa9zBNcB0JEXisz5xKO5pci0E+JJ0bEObs5XotrPtxcRnZBs6Q55izbcgrLtpziOhAi8jqiKGLpj/WjHo8N74y2gX5ObpH34siHG8vILsCM1QetLhUNcB0IEXmfH45eQk5BKYJUPnhyRBdnN8erMfhwUzq9iPSNOSar3FrCdSBE5E30ehFLfzwBAHhsGEc9nI3Bh5vam1fcohGPhqR1IHvzmBmViDzbppxC/FpYv8PlyVu51sPZGHy4qaKy1gUe9roWEZGrqR/1qF/r8fvhndEmgKMezsbgw01FBKtd8lpERK4m4+jNUY8nuNbDJTD4cFOD48IQFaqGYOJxAYDa1/yvV0B99tPBcWG2bh4RkUvQ60X8VRr1GBGH0ADm9XAFDD7clFIhYN74RABoFoBIX796Vw8oTEQn0uF54xOhNHUSEZGb+y67oD6bqdqHeT1cCIMPNzY2KQorpvSHJrTxtIkmVI0VU/rj0eQ4PDAgBgDgp1QYPYd5PojIU+n0It7LrN/h8sSIOGYzdSFMMubmxiZFISVRg715xSgqq0ZEcP00ijSa8ewdXfHVoQuo1enx+j2JCA/yQ0SwGgM6tcWBs9ewIetis+cQEXmCrw9dRO7lCrQJ8OWoh4th8OEBlAoByfHtjD7WsW0AHhwUg9V7zuGHo4VY9/RQ/HC0ECPf3sLKt0TktnR6sdFNV7+OwY0er63TY+lP9aMe00fGs3Kti2Hw4QVmje6KL/ZfwN68Yiz98QTe/+lUs+RkUsZTTsUQkaszVlaiU1sVXuhx85z/HDiP88VVCA9S4dHkTk5oJZnDNR9eICrUH5MHxwIAlm/NNZoVlRlPicgdmCorcam0/usfj11CtVaHv/10CgDwzOh4BPjxPtvVMPjwEjNHx8NPqYBWZzqwYMZTInJl5spKSMfe+v5XrN5zFoWl1YgOVeOhIbGObCLJxODDS0QEq3Fbt3BZ5zLjKRG5IjllJQpKqvH+T/V5PZ69IwEqH6UjmkZWYvDhRR4cJO8OgBlPicgVybkx0olAaXUdOrULwP0DOjqgVdQSDD68yO09IhCsMj33yYynROTKLN0YVdbVBx8AcHfvKCgEpg9wVQw+vIhSIeDN+3oZfYwZT4nI1VkqK/HTRQWkd7PlW3MxYtFmZGQXOKx9JB+DDy8zoW9HPDiw+VAkM54SkaszV1ZCFIFthY2PSikEGIC4Hu4/8kJvTEzCjlNXcfF6Fe7tE42HBscywykRuQWprETTPB91IqAXBQgQId4ITaQdMK+uP4IqrR6aEGZzdhUc+fBCKh8lZo9JAABsP3kZvTqE8I+RiNzG2KQo7Hjpdnz+1FA8MzoeACClJ/Ix8qlWXKHF8+uy8NBHezgV4yIYfHipSf07omtEEK5XavHhtlxnN4eIyCpSWYmESCmtuoDENnqTlbwlnIpxDQw+vJRSIeBPd3YHAPxzR54hOyARkTupqNHd+D8R98TqLZ7PbM6ugcGHF0tNjMSATm1RrdVj6Y8nnd0cIiKrfX9jBEMhAB0C5T2H2Zydj8GHFxMEAS+Pq6/E9MX+8zhVVO7kFhERGafTi9idexUbsi5id+5V6PQidp26gp9PXoFSIcCnBcvWmM3ZebjbxcsN6hyGMT0j8eOxS3j7h1/x4SMDnd0kIqJGjFWx1YSoDKnTHx4Si+S4NqjNO2DVdZnN2Xk48kF4aWx3KATgh6OXcOAshyGJyHWYqmJbWFqDs8WVUPsq8Ic7EjCmZyQA4JOpg/Deb/sgLNDPZDIyZnN2PgYfhITIYPx2YAwAYOF3v0IUuQiLiJzPXBVbiY9CgbYBfoavB8eF4b7+HbHgviQAzZORMZuza2DwQQCA2WO6Qe2rwP6z1/DjsSJnN4eISFYV2/KaOqMLR6VkZJrQxlMrzObsGrjmgwDU/0E+PjwOy7fmYlHGrxjRZaizm0REXk7ugtD680KaHR+bFIWURA325hWjqKwaEcHMcOoqGHyQwfRR8fh87zmcKirHFwcuoo2zG0REXk3uglBz50nJyMi1cNqFDELUvpg9phsA4K+bT6G6zskNIiKvZqmKLReOui8GH9TI5CGx6BIeiOIKLTLz+fIgIseTcnr893A+fjcoFkDzhaMSLhx1T5x2IQOdXsT+M9dwe48InN6Rh235AvKvV6FTe19nN42IvISxnB5tAnxRp9Oj3JBKvX7EY974RC4cdVMMPghA8z94ASK0ooA/fXkEX0wf7uTWEZE3kHJ6NN1ae71Sa/j/QZ3b4oWU7lw46uY4rk5Gk/hIZan3nrmOP6w5aEhnTERkD3JyegDAew/2RXJ8OwYebs6q4GPhwoUYNGgQgoODERERgYkTJ+L48eONzhFFEWlpaYiOjoa/vz9GjRqFo0eP2rTRZDum/uAVAjAovL5C5DeHC/DQR3swYtFmlqEmIruQk9MDAM4XVzmgNWRvVgUf27Ztw6xZs7Bnzx5kZmairq4OqampqKioMJyzePFiLFmyBMuWLcO+ffug0WiQkpKCsrIymzeeWs/cH/zdsXqgQVhSWFKNGasPMgAhIpuzLqcHuTur1nxkZGQ0+nrlypWIiIjAgQMHcNttt0EURSxduhRz587FpEmTAACrVq1CZGQk1qxZg2nTptmu5WQT5v6Q26oApQDobsQfIupXnKdvzEFKoobDnkRkM7bI6UHuo1ULTktKSgAAYWH1e6zz8vJQWFiI1NRUwzkqlQojR47Erl27jAYfNTU1qKmpMXxdWloKANBqtdBqtc3Obw3pera+rjsLD/CBStl8llWlqD8WoBRRVgcAApSCCB8FUFxehT2niri3/ga+ruRjX1nHm/qrX8dgdGqrwqXSasN4q1YP6EUBAkT4KeozMffrGGy0P7ypr1rLXn1lzfUEsYVVxERRxIQJE3Dt2jX8/PPPAIBdu3Zh+PDhuHjxIqKjow3nPv300zh79ix++OGHZtdJS0tDenp6s+Nr1qxBQEBAS5pGNrb3soDPTimhUoh4rZ8OIX6Wn0NE1Br5FcDiw0qIEPCHXnWIb549nVxMZWUlJk+ejJKSEoSEmP+FtXjk45lnnsHhw4exY8eOZo8JQuPheFEUmx2TvPLKK3jhhRcMX5eWliImJgapqakWG28trVaLzMxMpKSkwNeXuSskPx67hOfXZQG4ucJDpRAxf6Aer+9XoFpXf+dRoxcw76ASvor6stUc+ajH15V87CvreGN//XjsEhZ+dwznr9dAhAAFRGzID8TLfXpgTM9Ik8/zxr5qKXv1lTRzIUeLgo9nn30W33zzDbZv346OHTsajms0GgBAYWEhoqJuJn4pKipCZKTxF41KpYJKpWp23NfX124vIHte2x2Nu6UjBIWyWWIfAKjRC6jV3wwc9aKAILUfhnaN4JqPJvi6ko99ZR1v6q9xt3SEKCgx87OD8FEIeOeBvhjfJ1r2+4039VVr2bqvrLmWVbtdRFHEM888g6+++gqbN29GXFxco8fj4uKg0WiQmZlpOFZbW4tt27Zh2LBh1nwrcrCxSVHY8dLt+PypoXh8eGez57YJ8AXjDiKyh2qtDgu+OwYAmD4yHhP7deCNjgeyauRj1qxZWLNmDTZs2IDg4GAUFhYCAEJDQ+Hv7w9BEDB79mwsWLAACQkJSEhIwIIFCxAQEIDJkyfb5Qcg25GqPybHt8OgTqGozTvQ6PGIYBWuV2qRe7kCGw8X4N4+0SauRETUMh9tP40L16qgCVFj5uh4ZzeH7MSq4GPFihUAgFGjRjU6vnLlSjz22GMAgDlz5qCqqgozZ87EtWvXMGTIEGzatAnBwcE2aTA5xpiekfgur35tx5XKOkQE11eOXL7lFN7NPIGF3x3DmJ4RCPDzgU4vYm9eMYrKqg3n8U6FiKyVf70Kf996CgDwyl09EODHCiCeyqrfrJyNMYIgIC0tDWlpaS1tE7mQwXFhjebxnrqtC9btP48L16qwfEsukjqENFsrwoJPRNQSC747hmqtHoM6t+XIqodjbReyitpXidfvSQQAfLAtF9Ob1IQBmAmViKz3v9NX8d/DBVAIQNq9vUzukCTPwOCDrKLTiwhW+aCHJhh1JgrNSUfTN+awGB0RWaTTi0jbmAMA+N3gWPSKDnVyi8jeOKFGsmVkFxjdjmuMCKCgpBp784qRHN/O/o0jIrf1r91ncKygFP6+Sozs1h46vch1Yx6OIx8kS0Z2AWYYmWKxhEWgiMicL/adwxs3Rj2qtDpM+/cBVtD2Ahz5IIt0ehHpG3PQkgkUFoEiooYa7o47c6US7/14otk50rqxFVP6c+G6h2LwQRbtzSu2esRDQH0RKKZgJyKJ3KlbVtD2fJx2IYtaMnUiArgrSYO9ecVcdEpEVk/dNlw3Rp6HIx9kkbVTJwoB0IvAP3eewT93nkFUqBqv390TbQNVTERG5IVaM3XLdWOeicEHWTQ4LgxRoWoUllSbfPNQ+ypQrdUDqA88GiooqcbMNYcaHWMiMiLv0ZKpWwnXjXkmTruQRUqFgHnj6xOLNR2rEG78e+u+3lBakRSIiciIvEdLRi8E1N+kcN2YZ2LwQbKMTYrCiin9oQltfBeiCVVjxZT+iAz1h05G+n0JE5EReQ9rRy+k25h54xM5PeuhOO1Cso1NikJKosZoEbkNWRetvh4TkRF5B2nqVu7Ui4bTsh6PwQdZRakQjAYKrZmX5YIyIs8mTd1OX32w2WMC6m9Enh+TgM7hgVyQ7iUYfJBNyFmUagoXlBF5vrjwIMNOuIY4yuGdGHyQTUh3NjOM3NmYwkRkRN5Brxfx6voj0ItASmIkHh8ex233Xo7BB9mMtChVTgZDLigj8h5r953HgbPXEOinRPq9vRDdxt/ZTSInY/BBNtVwUeoH23Kx7cRltAv0g49SwKXSGsN5HGol8g5FZdV46/tjAIA/pnZn4EEAGHyQHUiLUm/pGIrU97bj4vUqPDEiDmN6RnKolcjLpH+Tg9LqOvTuEIqpwzo7uznkIpjng+wmUOWDN+9LAgCs3JmHQJUSE/p2QHJ8OwYeRF4gI7sQ3x4pgFIh4K3f9ObfPRkw+CC7Gt09AhP6RkMvAi99eQRand7ZTSIiByip1OL1DdkAgOkju6BXdKiTW0SuhMEH2d3r9ySiTYAvjhWU4qOfTzu7OUTkAAu+O4bLZTXo0j4Qz96e4OzmkIth8EF2Fx6kwut319eGWfrjSZwqKndyi4jInnaeuoJ1+89DEIDFv7kFal+ls5tELobBBznEpP4dMKp7e9TW6THn/35hPRciD6TTi9jyaxH+8Hl9FespQzphYGfm8aHmGHyQQwiCgAX39UaQygcHz13Hyp15zm4SEdlQRnYBRizajN9/ug9XK2oBAJtyClm5moxi8EEOE93GH3Pv7gkAeGfTcZy5UuHkFhGRLWRkF2DG6oPNkgsWldZgxuqDDECoGQYf5FC/GxSDEV3DUa3VY87/HYae0y9Ebk2nF5G+McdoTSfpWPrGHE61UiMMPsihBEHAwkm9EeCnxN4zxfjX7jPObhIRtcLevGKz5RREAAUl1dibV+y4RpHLY/BBDhcTFoBXxvUAACzKOI6zVzn9QuSOdHoRO09dkXVuUZn5ek/kXRh8kFM8PKQThnYJQ5VWhxf/w90vRO5GWmC6bMspWeefvFSO3blX+bdOABh8kJMoFALevr8PAv2U2HfmGj5m8jEit2Fqgak5y7acwkMf7cGIRZu5AJUYfJBz6PQiLlyrwr19owEA7246geOFZU5uFRFZYm6BqRyFJdXcAUMMPsjxpOHahz7ag8/3ngcA1Or0eHLVPtTWsfYLkSuztMDUEu6AIYDBBzmYueHa89eq8NzaQ05oFRHJZYuFo9wBQww+yGHkDNd+n12I/Wf4hkTkanR6Ebtzr+LkJdtNj3IHjPfycXYDyHvIHa59Zs0h/PTHkQhU8eVJ5AoysguQvjFH1t+vAKBtoC+KK7QWz40IVtugdeSOOPJBDiP3LqewtBrpG4/auTVEZI400vHGxqOYLnNni3Djv29OSEJUqNrwtbHzokLVGBzHonPeireW5DBy73IEAF/sv4DbEtqjXZAKRWXViAiuf6NSKky9nRGRrVgz0tGQJlSNeeMTMTYpCgqFgBmrD0IAGk21Sn/B88Yn8u/ZizH4IIcZHBeGqFA1Ckuqja77EFD/5jWxXwes2JqLZ9cegtjgxKgGb2xEZB/SonBr9qE8M7orhncNb3SDMDYpCium9G8WxGj4d0xg8EEOpFQImDc+0eLdkLT9Tmzy7iflB1gxpT/fuIjsoKU5PBIig5Ac367Z8bFJUUhJ1GBvXjFHMKkRrvkgh5LuhjShjadgNKFqrJjSHymJGrz57TGjz2V+ACL7amkOD3NTqkqFgOT4dpjQtwOS49sx8CAAHPkgJzB3N7Q796rsCpnG7rSIqOWs3foqTZVy4ShZi8EHOYV0N9SQNRUyv7+RmplDuES2Y83WVy4cpdZg8EEuwdrV9f/afRb/2n2Wi1CJbMjSovCGuHCUWoNrPsjpWlIhU8IiVUS2Iy0KN+eJ4Z3x+VNDseOl2xl4UIsx+CCnam2FTC5CJbKtsUlReOb2+GbHo0LV+GBKf7w+vhcXjlKrcdqFnKq1FTKBm4tQP92Zh8eGx/FNkagVCkqqsHrPOQDAmJ6RGN8niltkyeYYfJBT2bKw1Pxvj+HjHXmchyZqIa1Oj2fXHMK1Si16RYdg2eR+UPsqnd0s8kCcdiGnsnVhKa4BIbKeVMdl2r8PYP/ZawhS+WD5w/0ZeJDdcOSDnEpOyvXIEBUAAZdKLa/AF288J31jDlISNRwmJrLA2E4zH6WAYwWl6NQu0IktI0/GkQ9yqoar65uGCdLXaff2Qtq9xs8xpmEiMiIyzdROs5JKLUcQya4YfJDTWUq5PjYpyuQ55thyPQmRpzG304y7yMjeOO1CLkFOASrpnE935mG+ifovDdl6PQmRJ7G004ylDMieGHyQyzCWct3YOY8Nj8PHO/LMrhNhvQki8+SODHIEkeyB0y7kdsytE5Gw3gSReVfLa2WdxxFEsgcGH+SWTK0BUSoELPltH+b5IDKjsKQaf99yyuw5AuqzmnIEkeyB0y7kthquEzl1uQzvZZ5EcUUtNh4uwL19O3Dkg8iIaq0O0/69H1cratGhjT8uXq+CADSawmTFWrI3jnyQW5PWiTwytDM+/f0gqHwU2PxrEZZkHnd204hcjiiKmPN/h/HLhRIE+Cnxpzu7Y/lk8zvNiOyBIx/kMW7p2AaLfnMLZq/Lwt+35CIxKhR338I3T/IeOr1odsfY8+uy8M0v+QCAylodZq/LQlSoGq/f3RNtA1Umn0dka1aPfGzfvh3jx49HdHQ0BEHA119/3ehxURSRlpaG6Oho+Pv7Y9SoUTh69Kit2ktk1sR+HfDUrXEAgD/+JwtHLpQ4uUVEjpGRXYARizbjoY/24Lm1WXjooz0YsWizIVHYW98fw9dZ+c2eV1hSjVlrDqGkqhYT+nZgxVpyCKuDj4qKCvTp0wfLli0z+vjixYuxZMkSLFu2DPv27YNGo0FKSgrKyspa3VgiOV4a2wO3dWuPaq0eT/5rHwpbWTWXyNWZylRaWFKN6asPYu76w/hg22mjz2VCMXIGq4OPcePG4c0338SkSZOaPSaKIpYuXYq5c+di0qRJSEpKwqpVq1BZWYk1a9bYpMFElvgoFVg2uR8SIoJwqbQGT6zah4qaOmc3i8gu5GQq/ex/581egyUJyNFsuuYjLy8PhYWFSE1NNRxTqVQYOXIkdu3ahWnTpjV7Tk1NDWpqagxfl5aWAgC0Wi20Wq0tm2e4nq2v64ncva/8lcCHU/ri/g//h6P5pXju84NY9lBfw3CyTi/iwNlruFJeg/AgFQZ0atvioWZ37ytHYl9ZR05/7c0rRnF5FVRGCtCKIqDVA+KN/Sx+CkAw8zIvKqmAVhvSylY7B19b8tmrr6y5niCKYovH2QRBwPr16zFx4kQAwK5duzB8+HBcvHgR0dHRhvOefvppnD17Fj/88EOza6SlpSE9Pb3Z8TVr1iAgIKClTSMCAOSVAcuOKlEnCrg9So8JnfXObhKRQ+hFYOUJBQ4XKxCgFDG7tw6R/s5uFXmyyspKTJ48GSUlJQgJMR/E2mW3i9AktBZFsdkxySuvvIIXXnjB8HVpaSliYmKQmppqsfHW0mq1yMzMREpKCnx9fW16bU/jSX0Ve7gAL/znCDYXKKCJicVXB843G6KWXp3vPdgXY3pGWnV9T+ore2NfWcdcf/147BLe+v5XFJYaX9NUpwd0Yv2Ih1YE3j1iZGjkBgFAZIgaP8y+zW0Xm/K1JZ+9+kqauZDDpsGHRqMBABQWFiIq6uYWx6KiIkRGGn9DV6lUUKlUzY77+vra7QVkz2t7Gnfsq6bbDSf0i8G5a9VY+uNJrNl7AaaSsgsA3vj2OFKTWpagzB37ylnYV9Zp2l8Z2QWYueaXG0G0pdeqAK2ZAT/p2a/c3QtqlV/rGuoC+NqSz9Z9Zc21bBp8xMXFQaPRIDMzE/369QMA1NbWYtu2bVi0aJEtvxWRURnZBUjfmNNo1X9UqBp/vqcnRnZrj20nLpt8Lqt4kjswt8C0JTShaswbn8iEYuRQVgcf5eXlOHXqZk2AvLw8ZGVlISwsDLGxsZg9ezYWLFiAhIQEJCQkYMGCBQgICMDkyZNt2nCipqTthk3flAtLqjHzs0N4LLkTtsm4Dqt4kivbm1fcbEttS7Tx98XfH+6PoV2Y14Mcz+rgY//+/Rg9erTha2m9xtSpU/Hpp59izpw5qKqqwsyZM3Ht2jUMGTIEmzZtQnBwsO1aTdSEpe2GAoANh5snWDKGVTzJlckNju/sFYFNR4sAGK/b8tZvemN413DbNo5IJquDj1GjRsHcBhlBEJCWloa0tLTWtIvIKpbuBkUAxRVahAX6objCeClxAfVD0KziSa5MbnD82LAuuK9fx2bTkJxmIVfA2i7kEeTeDU7sG42VO8+YnC9nFU9ydYPjwhAVqkZhSbXJ13FU6M36LFLlZ9ZtIVfCqrbkEeTeDaYkarBiSn9ENaniqfJR4P3f9UWovx82ZF3E7tyrTDVNLkmpEDBvfKLJxwU0DqKlys+s20KuhCMf5BEs3Q02nFJpeDe4K/cKPtiWi5o6PeZ8eRhVDfYkRnF4mlzU2KQovP9QX/zxP4dRW8fXLLkfBh/kEaS7wRmrD95IJH2TdJ9n7G4wOb4d9HoRf9+a2yjwAOp3ycxYfRArpvTnmzm5lNo6Pb4+lI/aOj3UvgrMGBmPwXHtOKVCboPTLuQxxiZFYcWU/tA0mVLRhKrx98n9jE6p6PQivjp00ej1WO2TXFGdTo/n12Xhp1+LoPJR4NPfD8ZzY7pxSoXcCkc+yKOMTYpqtsDuWkUt5n/bPPHYvPGJCPX3s7hLpmHisabZU7kzhhxJrxcx58vD+PZIAXyVAj58ZACGdmFCPHI/DD7I40hTKkB94rFZa4wnHpux+iAeH95Z1jWLyqpNZ0+9u7uNWk5kmrZOj+mrD+CnX4ugEID3H+qHUd0jnN0sohbhtAt5LEuJxwBgfZbxKZemzlypxIzVB5uNkhSUVGP2uiwA9blGOD1DtiKNsgHA8i2ncEv6Jvz0a33SML0IvLExBxnZBc5sIlGLMfggj2VN4jFzM+UqHwVW7TadG0Ty+Kp9GLFoMz8QqNUysgswYtFmPL5qH/Qi8NctuajS6hqdI43e8fVG7ojBB3ksaxKPAaZrg9bU6U1mRW2KHwjUEjq9iN25V7Eh6yL++uNJwyibKAKf5yqgF5u/OrkgmtwZ13yQx7Im8djguLBm6zlaQqojk74xBymJGu4+IIuMrSWS1InA3ssK3HxlNcZKzOSuGHyQx2pJ4rE9uVcxa81BXK/Stvj78gOBGjK2Q0oKSk1VYpboRQEKQYRCAOr0Jk4CKzGT+2HwQR6rJYnHFAqhVYFHQ/xAIFM7pOaNT0RKosbkguibRDyWoMdnuQrUmTmLlZjJ3XDNB3k0c4nHjGUutWXAwA8E7yaNajSdTpHWBS3bfNLiNJ+vAujTzkwVcdwsIkfkTjjyQR7PWOIxU2mobREwNJzOIe9kaZu3AGDlzjMWr2NuyZCx0Tsid8Hgg7xCw8Rj5sgpVx6sVqKqVo+6GzsMxAYn8gPBu5hazyFnm3drp/c0LCJHbozBB1EDctaJvH1/H3RoE4DJH+9BWXUdavVAUVX9Y5pQNV6/u6ehjkxEsBoDOrXFgbPXLI66kHsxt56jxtzq0Aba+PvKDkI0ISo8NDgWncMD+Toit8fgg6gJaZ1I0w+Wpnea/312BB75516cK67EX7OVeHVcT7QPCWhWR0Yh1GeklLDsufsztUtFWs8xe0yCrOvcmhCOjYeb54SRAt9Zo+KByhP4ZOogDO0awWCDPAaDDyIj5KwT6dQuEF/OGIbHVv4PR/PLMP+744apmIaaHpI+oIwteCXXJ2c9x+d7z0ETosalUtPbvIPVPvjvkfrAQ+WjaDRaIgW6d3QPx3ffneAoB3kcBh9EJshZJ9I+WIXPHh+Eh//+I45ek7d5jInI3Juc9RyFpTV4fkw3LP3xRLPpO+mc0ur6zbNThsbi9bsTcfDc9WaBrlZrm23fRK6GW22JWilQ5YMnu+uhFOSnuG6YiIzci9zt2J3DA4xu8/ZT3nzbfWlsD8yfkASVrxLJ8e0woW8HJMe3Y0BKHo8jH0Q2oBAAHwWg01k+tyEmInM/crdjXymrwWPD4wzTd4cvXMcnO/NwqbQGgX5KTB8Zj+g2auw5XcxpFfI6DD6InOjkpXLszr3KDx83Imc7NgDM//YYPt6Rh3njE+GrVOBvm0+hvKYO7QL9IAjAu5knDOdyETJ5GwYfRDaiCVHj3LUaC+myG1u25RSWbTmFqBtbdNsGqrgl18WZ247dVEFJNaavPmj4OiEiCCeLypudx0XI5G0YfBDZyMvjemDmml8sfiAZU1BSjZlrDjU6Zuxu2FyRMnIcU9uxzZk8JBabj10y+hgXIZO3cdvgQ6fTWb0SXKvVwsfHB9XV1dBZOznvZezVV76+vlAqlTa7nisZ0zPS6AdSS4IRoPndsLmkVrxbdjxpO/anO/Mw/9tjFs+PDw/EmtIak4+zGjJ5E7cLPkRRRGFhIa5fv96i52o0Gpw/fx6CwDsLc+zZV23atIFGo/HI34Gx/CADOrXFv3afwZLME6islR/INbwb1uuBWWtMJ7XicL1zKBUCwoNVss49W1wp6zwuQiZv4HbBhxR4REREICAgwKoPML1ej/LycgQFBUGh4C5jc+zRV6IoorKyEkVFRQCAqCjP/LA0lh/kyVu7QOWrwOtfH7XqWtLd8Gsbss0mteJwvf2ZmvKSu/ulU1iArPNYDZm8gVsFHzqdzhB4tGtn/bCkXq9HbW0t1Go1gw8L7NVX/v7+AICioiJERER47BSMMV3bB7f4ucUVtSYf43C9/Zmb8kpJ1Jit0SJVOX4kuTM+3pFncpcMqyGTN3GrT2BpjUdAgLw7CHJN0u/P27I3Sls07TU2weF6+5DquDRdWFp4YyfLwx/vMRt4APVVjv18FJg3PrHRcWPncfSKvIFbBR8ST1wr4E289fcnbdEEmn/42AKH623PUh0XANhzuhgKAbgrSQNNSOP1H5pQdaP1ONIumaZZT5ueR+Tp3GrahcjdtWSLphxt/H2hF0Xo9CLvnFvA1HoOS3VcJK/fk4jfD4+TtRVaTtFCIk/H4IPIwYx9+Fy8VonXNmSjWqu3fAEjrldp8fDH/+PW2xYwt56jYaVZc8IC/QDIK0ZozXlEnsotp13c0WOPPQZBECAIAnx9fREZGYmUlBR88skn0Ovlf+B8+umnaNOmjf0aSg4hffhIhcTuHxiDo+lj8buBHVt1XWnrbUZ2gY1a6n50ehG7c69iQ9ZF7M69Cp3edJYVc+s5Zqw+iDNXKmR9T055EVnHa0c+nJEpcuzYsVi5ciV0Oh0uXbqEjIwMPPfcc/i///s/fPPNN/Dx8dpfB+FGQNI1HGv3X7B4boCf0mjOEG/femtNIjZL6zkEAJ/976zJvga4Q4Wopbxy5CMjuxAjFm3GQx/twXNrs/DQR3swYtFmu98tqlQqaDQadOjQAf3798err76KDRs24Pvvv8enn34KAFiyZAl69+6NwMBAxMTEYObMmSgvr68FsXXrVvz+979HSUmJYRQlLS0NALB69WoMHDgQwcHB0Gg0mDx5siGfBrkPuXfQ5pKVNdx6K5c1owWuytIoRsO/b51exKc788yu5xABFJXVmg08AO5QIWoJrws+fjp+FbPWHJL1BuUIt99+O/r06YOvvvoKAKBQKPD+++8jOzsbq1atwubNmzFnzhwAwLBhw7B06VKEhISgoKAABQUFePHFFwEAtbW1mD9/Pn755Rd8/fXXyMvLw2OPPebQn4VaT852XLkfdHK33mZkFzglGLclObtS0jfmQKcXDT+vnJToABCs9sF9fTtAE8IdKkS24lXj/Dq9iMU/nna5TJE9evTA4cOHAQCzZ882HI+Li8P8+fMxY8YMLF++HH5+fggNDYUgCNBoNI2u8fjjjxv+v0uXLnj//fcxePBgQ5ZScg9yKqbKHZWQM4oijRa4e9p2S7tSpNGgZZtPYemPJ6yqtVNWXYf1WRehCVHh+TEJ6BweyB0qRK3kVSMf+84U41KZvEyRjiSKoiH3xZYtW5CSkoIOHTogODgYjz76KK5evYqKCvML3w4dOoQJEyagU6dOCA4OxqhRowAA586ds3fzycZM5YKIClVj4aQkPDDA/KJU4ca5ltYhWDNa0BoNp3Ts9bcld5Rn5c68FhX5A4BLpTVY+uNJqHwUSI5vx8CDqBW8auSjqMx0RcnG5zk2U+SxY8cQFxeHs2fP4q677sL06dMxf/58hIWFYceOHXjiiSfMZgOtqKhAamoqUlNTsXr1arRv3x7nzp3DnXfeidpa08EWuS5zuSAeGtwJvTuG4M8bcow+V4TldQhy1zy0Nm170wWgKqWIxYOBH49dwrhbWrezpyG5a2VMZSKVw9sX8xLZklcFHxEyq086ctvc5s2bceTIETz//PPYv38/6urq8O677xrqqXzxxReNzvfz82tW4v7XX3/FlStX8NZbbyEmJgYAsH//fsf8AGQ35nJBPJoch4hgNV77OhtXyhsHmJ3aBUCnh8mEY8Z2hJhjLhg3t2vM1JQOADy/LguCQmmzKR1prYy5uimhZuqvyMU6OkS24VXBx6DOYYgM9kNRWa1TCjvV1NSgsLCw0VbbhQsX4p577sGjjz6KI0eOoK6uDn/7298wfvx47Ny5Ex988EGja3Tu3Bnl5eX46aef0KdPHwQEBCA2NhZ+fn7429/+hunTpyM7Oxvz58+3y89ArqPh6Eh2fgn+d/oqfj55BWevVmLWmoNoE+CLMT0j8fyYbujQtr6gn7mAwJSTl8qxO/dqszUOloqtmZrSkdhyBMHSWhkRQFhg64MPCevoELWOV635UCoEzBnTBYBzCjtlZGQgKioKnTt3xtixY7Flyxa8//772LBhA5RKJfr27YslS5Zg0aJFSEpKwmeffYaFCxc2usawYcMwffp0PPjgg2jfvj0WL16M9u3b49NPP8V//vMfJCYm4q233sI777xjl5+BXIs0OvLUrV3w8dRBmD+xF4JU9ZWCr1dq8X8HLmD4os24d9kO/JhzCWnfHLV6zcOyLaea7YCxtK112eaTsqZ03ss8YbOtvabWykhOX6k0+3wB9QGKHEwqRtQ6giiKLrWhv7S0FKGhoSgpKUFISEijx6qrq5GXl4e4uDio1db/8ev1epSWlmLXuUrM//aYrERE3krqq5CQEMMUkK209vfoarRaLb777jvcdddd8PWV9+FlDy0Z1bCGFJL/fXK/Zn8/TQWqlKioaZ4fo37Nhw5z9ipRo7sZ5Nvq70+r02PXqSv4dNcZ7Dh1BVqdvN5o+rNZKnu/46XbHbLmw1VeW+6AfSWfvfrK3Od3U1417SIZm6TBnUlRLOxEHsPczhVbkRZcvrYhG8UV5qcvjAUe5rRma2+1VoefT17B99kF+OlYEUpaMLWiaRD8KBSC0ekbJhUjsh2vDD4AFnYizyK3+mpriYDFwKOl1226k8TUYladXkROfil2nLqCHacuY29ecaMRjnaBfugT0wabf7Wc4feZ0V0xvGt4o5sPU5WHNRwdJbIZrw0+iDyJJyyAbLiTpKSqttmHf5DKB53bBeDC9SpcrzQeAIUF+mH+hF7Q6kVZwUdCZJDRmxCWvSeyLwYfRB7AVgsgO7bxx4XrVRbPC/RTosJMfZnWeP+nE9h9unkysvKaOmTnlwIA1L4KVGubV4O+VlGLWWsOYfaYBFnfy1y/cXSUyH68arcLkaeSUxPGnKhQNT6Y0h/b5oyGJsRyPpyWBh5aPbDyuAJaM083Fng0FOinhL+v0uhj0uTL53vPQRNiuj/kZoElIvtg8EHkAaQ8F0DzbeTmPDO6Kz5/aih2vHQ7xiZFQakQkHZvLwhmrnNbQjhGdmsvK0hpSi8KyCpWQG/i6nJmNSpqdbhmYtoFqA9ACktr8NDgWADO2VZPROZx2oXIQ5haKGmMtGX0+ZRuzT6ATV3H2HbYOp0emTmXkHu5HB9uP42y6jqT3zNIpcQdPSLgU3IBG84qUCc2/+D//bDO+OfOM/J+YAs6hwdw4SiRi2LwQeRBGi6UzMwpxCc7z7Roy6jcBZc+SgXG9a7/EO8aEYQZqw8CJr7fOw/0wR3dw/Hdd+exv0SNs9du1lqSAptQfz+bBR8RwWokx7fjwlEiF8Tgg8jDSAslk+PbYXBcWIvv/K1dcClni6pUIPGH2bfh0IUyo9tozdVokaNpmQQuHCVyPQw+XIggCFi/fj0mTpxo9PEzZ84gLi4Ohw4dQt++fR3aNmM6d+6M2bNnY/bs2c5uCpng6C2jcr9f04BApxexO/cqisqq8btBsVj64wmjNVos4XoOIvfA4MNBHnvsMaxatQoAoFQqER0djbvvvhsLFixA27ZtAQAFBQWG/3cVt99+O7Zt29bsuFarxb59+xAYGGg4Zil4Iudw9J2/td/PWIG6NgH1KZ9N5fMwhes5iNwDgw8HGjt2LFauXIm6ujrk5OTg8ccfx/Xr1/H5558DADQajZNbaNxTTz2FN954o9ExHx8ftG/f3kktIk9hqh5NSaUWIoDn7uiKVbvOmq1GGxboi9fv6QVNCNdzELkLbrV1IJVKBY1Gg44dOyI1NRUPPvggNm3aZHhcEAR8/fXXhq/37t2Lfv36Qa1WY+DAgTh06FCza37zzTdISEiAv78/Ro8ejVWrVkEQBFy/ft1wzq5du3DbbbfB398fMTEx+MMf/oCKigrZ7Q4ICIBGo2n0D6ifdlm6dKnh/wHgvvvugyAIhq+JTDFXj0ZKt/7F/gtYcF9vo1t/pWML7uuN+/p1QHJ8OwYeRG7C7YMPURRRWVsn+19Vrc6q8839a01B4NOnTyMjI8NkRcGKigrcc8896N69Ow4cOIC0tDS8+OKLjc45c+YM7r//fkycOBFZWVmYNm0a5s6d2+icI0eO4M4778SkSZNw+PBhrFu3Djt27MAzzzzT4rYbs2/fPgDAypUrUVBQYPiaPIu0NmND1kXszr0Knb7lfwOW6tFI6dbbBvphxZT+0IQ2zkaqCVW3qBAdETmf3aZdli9fjrfffhsFBQXo1asXli5diltvvdXm36dKq0Pin3+w+XXlyHnjTgT4ye/C//73vwgKCoJOp0N1df2b7pIlS4ye+9lnn0Gn0+GTTz5BQEAAevXqhQsXLmDGjBmGcz744AN0794db7/9NgCge/fuyM7Oxl/+8hfDOW+//TYmT55sWBSakJCA999/HyNHjsSKFStklbRfvnw5Pv74Y8PX06ZNw7vvvtvoHGkKpk2bNi47fUStY2xthrHcH3LJrUdTVFaNCX07cMsskQexS/Cxbt06zJ49G8uXL8fw4cPx4YcfYty4ccjJyUFsbKw9vqVbGD16NFasWIHKykp8/PHHOHHiBJ599lmj5x47dgx9+vRBQECA4VhycnKjc44fP45BgwY1OjZ48OBGXx84cACnTp3CZ599ZjgmiiL0ej3y8vLQs2dPi+1++OGHG42otGnTxuJzyLOYWptRWFKNGasPtmgEQm49Guk8bpkl8hx2CT6WLFmCJ554Ak8++SQAYOnSpfjhhx+wYsUKLFy40Kbfy99XiZw37pR1rl6vR1lpGYJDgqFQtH7GyVR9CVMCAwPRtWtXAMD777+P0aNHIz09HfPnz292rpwpHVEUIQhCs2MN6fV6TJs2DX/4wx+aPV9uIBgaGmpoN3kfOWsz0jfmICVRY9VIhFSPxlROj6b5OojIc9g8+KitrcWBAwfw8ssvNzqempqKXbt2NTu/pqYGNTU3Mx2WltZXrdRqtYaERBKtVmu4a9frb1a0VPvICyREUUDdjaJUTT+0W0IURdnrPqRzG7b79ddfx913341p06YhOjoaAAw/W48ePfDvf/8bFRUV8Pf3BwBD/0nndO/eHd9//32ja0prLaRz+vXrh6NHj6JLly5G29XwuU3ba6rdTc+THvP19YVWqzV5bsPvKYoitFotlErrAjhXJL1Om75ePcXevGIUl1dBZeZXVVxehT2niiwGCg37yhfAn+/ujufXZQEwnhX1z3d3h15XB719Cui6PE9/bdkS+0o+e/WVNdezefBx5coV6HQ6REZGNjoeGRmJwsLCZucvXLgQ6enpzY5v2rSp0ZQDUL+9U6PRoLy8HLW1tS1uY1lZWYuf21JarRZ1dXWG4AoA+vfvjx49eiA9Pd2wbqOqqgqlpaW455578Nprr2Hq1Kl48cUXce7cObzzzjsA6hejlpaWYvLkyXjvvffw/PPP45FHHsGRI0ewcuVKw8+oUCgwc+ZMpKam4umnn8bUqVMREBCA48ePY+vWrVi8eLHFdut0OtTW1jZqt0Sv16O6utrwWGxsLDIyMnDLLbdApVKZnJ6pra1FVVUVtm/fjro607VA3E1mZqazm2A3iwdbPufKsT347pi86zXsq0Vmrl2bdwDf5cm7pifz5NeWrbGv5LN1X1VWVso+124LTo1NBxgbbXjllVfwwgsvGL4uLS1FTEwMUlNTERIS0ujc6upqnD9/HkFBQbIWSjYliiLKysoQHBxsk5EPa/j6+sLHx6fZz/THP/4RTzzxBF577TUAgL+/P0JCQhASEoJvvvkGM2fOxMiRI5GYmIhFixbhgQceQGBgIEJCQtC7d2988cUX+NOf/oQPP/wQycnJmDt3LmbNmoX27dtDrVZj2LBh2LJlC1577TXcddddEEUR8fHx+O1vf9usLQ1JfaVUKuHn52f0XIVCAbVabXjs3XffxYsvvoh//etf6NChA06fPm302tXV1fD398dtt93Wot+jq9FqtcjMzERKSorJ3UvubG9eMR5fZXn30idTB8ka+TDWVzq9iANnr+FKeQ3Cg1QY0KktF5PC819btsS+ks9efWXsJtUUmwcf4eHhUCqVzUY5ioqKmo2GAPW5L1Sq5qW5fX19m3WKTqeDIAhQKBQtWrMhTQdI13AkKbtpU1OmTMGUKVMANF+vMWzYMGRlZTU61vSciRMnNsoo+pe//AUdO3ZsNGo0ZMgQqyNcqa+2bNlisq/OnDnT6OsJEyZgwoQJFq+tUCggCILR37E787SfRzK0awTCgvwtrs0Y2jVCdsDQtK98AQzv1vz9gep56mvLHthX8tm6r6y5ls0/gf38/DBgwIBmH3aZmZkYNmyYrb+d11u+fDn27duH06dP49///jfefvttTJ061dnNIg+iVAiYNz4RgPFEXwBrqRCRdewy7fLCCy/gkUcewcCBA5GcnIx//OMfOHfuHKZPn26Pb+fVTp48iTfffBPFxcWIjY3FH//4R7zyyiuynvvzzz9j3LhxJh+/cOGCrZpJbk5OxVoiIrnsEnw8+OCDuHr1Kt544w0UFBQgKSkJ3333HTp16mSPb+fV3nvvPbz33nsteu7AgQObTetILO1YIe/j6Aq5ROS57LbgdObMmZg5c6a9Lk824O/vbzJ/h16vt2rxEHkHJvoiIltw+9ouRERE5F7cMvjglIB74++PiMi72W3axR78/PygUCiQn5+P9u3bw8/Pz6p8HXq9HrW1taiurnb4Vlt3Y4++EkURtbW1uHz5MhQKBfz8/GxyXSIici9uFXwoFArExcWhoKAA+fn5Vj9fFEVUVVXB39/f4UnG3I09+yogIACxsbEMAImIvJRbBR9A/ehHbGws6urqoNNZV/BBq9Vi+/btuO2225iExgJ79ZVSqYSPjw+DPyIiL+Z2wQeAFmfHVCqVqKurg1qtZvBhAfuKiIjshePeRERE5FAMPoiIiMihGHwQERGRQ7ncmg+paqs9smtqtVpUVlaitLSU6xgsYF/Jx76Sj31lHfaXfOwr+ezVV9LndtPq68a4XPBRVlYGAIiJiXFyS4iIiMhaZWVlCA0NNXuOIMoJURxIr9cjPz8fwcHBNt+OWVpaipiYGJw/fx4hISE2vbanYV/Jx76Sj31lHfaXfOwr+ezVV6IooqysDNHR0RbzOLncyIdCoUDHjh3t+j1CQkL44pSJfSUf+0o+9pV12F/ysa/ks0dfWRrxkHDBKRERETkUgw8iIiJyKK8KPlQqFebNmweVSuXsprg89pV87Cv52FfWYX/Jx76SzxX6yuUWnBIREZFn86qRDyIiInI+Bh9ERETkUAw+iIiIyKEYfBAREZFDeW3wce+99yI2NhZqtRpRUVF45JFHkJ+f7+xmuZwzZ87giSeeQFxcHPz9/REfH4958+ahtrbW2U1zWX/5y18wbNgwBAQEoE2bNs5ujktZvnw54uLioFarMWDAAPz888/ObpJL2r59O8aPH4/o6GgIgoCvv/7a2U1ySQsXLsSgQYMQHByMiIgITJw4EcePH3d2s1zWihUrcMsttxiSiyUnJ+P77793Slu8NvgYPXo0vvjiCxw/fhxffvklcnNzcf/99zu7WS7n119/hV6vx4cffoijR4/ivffewwcffIBXX33V2U1zWbW1tXjggQcwY8YMZzfFpaxbtw6zZ8/G3LlzcejQIdx6660YN24czp075+ymuZyKigr06dMHy5Ytc3ZTXNq2bdswa9Ys7NmzB5mZmairq0NqaioqKiqc3TSX1LFjR7z11lvYv38/9u/fj9tvvx0TJkzA0aNHHd8YkURRFMUNGzaIgiCItbW1zm6Ky1u8eLEYFxfn7Ga4vJUrV4qhoaHObobLGDx4sDh9+vRGx3r06CG+/PLLTmqRewAgrl+/3tnNcAtFRUUiAHHbtm3OborbaNu2rfjxxx87/Pt67chHQ8XFxfjss88wbNgwlmKWoaSkBGFhYc5uBrmR2tpaHDhwAKmpqY2Op6amYteuXU5qFXmakpISAOD7kww6nQ5r165FRUUFkpOTHf79vTr4eOmllxAYGIh27drh3Llz2LBhg7Ob5PJyc3Pxt7/9DdOnT3d2U8iNXLlyBTqdDpGRkY2OR0ZGorCw0EmtIk8iiiJeeOEFjBgxAklJSc5ujss6cuQIgoKCoFKpMH36dKxfvx6JiYkOb4dHBR9paWkQBMHsv/379xvO/9Of/oRDhw5h06ZNUCqVePTRRyF6ScJXa/sKAPLz8zF27Fg88MADePLJJ53UcudoSX9Rc4IgNPpaFMVmx4ha4plnnsHhw4fx+eefO7spLq179+7IysrCnj17MGPGDEydOhU5OTkOb4ePw7+jHT3zzDP43e9+Z/aczp07G/4/PDwc4eHh6NatG3r27ImYmBjs2bPHKUNQjmZtX+Xn52P06NFITk7GP/7xDzu3zvVY21/UWHh4OJRKZbNRjqKiomajIUTWevbZZ/HNN99g+/bt6Nixo7Ob49L8/PzQtWtXAMDAgQOxb98+/PWvf8WHH37o0HZ4VPAhBRMtIY141NTU2LJJLsuavrp48SJGjx6NAQMGYOXKlVAoPGrATJbWvLao/g1vwIAByMzMxH333Wc4npmZiQkTJjixZeTORFHEs88+i/Xr12Pr1q2Ii4tzdpPcjiiKTvnc86jgQ669e/di7969GDFiBNq2bYvTp0/jz3/+M+Lj471i1MMa+fn5GDVqFGJjY/HOO+/g8uXLhsc0Go0TW+a6zp07h+LiYpw7dw46nQ5ZWVkAgK5duyIoKMi5jXOiF154AY888ggGDhxoGEE7d+4c1w8ZUV5ejlOnThm+zsvLQ1ZWFsLCwhAbG+vElrmWWbNmYc2aNdiwYQOCg4MNI2uhoaHw9/d3cutcz6uvvopx48YhJiYGZWVlWLt2LbZu3YqMjAzHN8bh+2tcwOHDh8XRo0eLYWFhokqlEjt37ixOnz5dvHDhgrOb5nJWrlwpAjD6j4ybOnWq0f7asmWLs5vmdH//+9/FTp06iX5+fmL//v25JdKELVu2GH0NTZ061dlNcymm3ptWrlzp7Ka5pMcff9zw99e+fXvxjjvuEDdt2uSUtgii6CUrLImIiMgleN/kPRERETkVgw8iIiJyKAYfRERE5FAMPoiIiMihGHwQERGRQzH4ICIiIodi8EFEREQOxeCDiIiIHIrBBxERETkUgw8iIiJyKAYfRERE5FAMPoiIiMih/h/tNjXjSqd+OAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -579,52 +571,7 @@ } ], "source": [ - "from sklearn import linear_model\n", - "np.random.seed(2018)\n", - "n = 100\n", - "d = 2\n", - "Lambda = 0.01\n", - "true_beta = [2, 0.5, 3.7]\n", - "\n", - "# Make data set.\n", - "x = np.linspace(-3, 3, n)\n", - "y = (2 + 0.5*x + 3.7*x**2)\n", - "\n", - "\n", - "\n", - "#Design matrix X does include the intercept. \n", - "X = np.zeros((n, d))\n", - "for p in range(d): \n", - " X[:, p] = x ** (p+1)\n", - "\n", - "\n", - "#Split data in train and test\n", - "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n", - "\n", - "\n", - "# Scale data by subtracting mean value using scikit-learn\n", - "from sklearn.preprocessing import StandardScaler\n", - "scaler = StandardScaler()\n", - "#scaler.fit(X_train)\n", - "#scaler.fit(y_train)\n", - "#X_train_scaled = scaler.transform(X_train)\n", - "#X_test_scaled = scaler.transform(X_test)\n", - "#y_train_scaled = scaler.transform(y_train)\n", - "\n", - "#Calculate beta\n", - "OLS = LinearRegression(fit_intercept=True)\n", - "OLS.fit(X_train,y_train)\n", - "ypredictOLS = OLS.predict(X_test)\n", - "RegRidge = linear_model.Ridge(Lambda,fit_intercept=True)\n", - "RegRidge.fit(X_train,y_train)\n", - "ypredictRidge = RegRidge.predict(X_test)\n", - "print(OLS.coef_)\n", - "print(RegRidge.coef_)\n", - "print(OLS.intercept_)\n", - "interceptRidge = RegRidge.intercept_\n", - "print(RegRidge.intercept_)\n", - "#predict value without intercept\n", - "ytilde_test_Ridge = X_test @ RegRidge.coef_+ RegRidge.intercept_\n", + "X_test @ RegRidge.coef_+ RegRidge.intercept_\n", "ytilde_test_OLS = X_test @ OLS.coef_+ OLS.intercept_\n", "\n", "#Calculate MSE\n", @@ -644,7 +591,7 @@ { "cell_type": "code", "execution_count": null, - "id": "cf021a10", + "id": "de73251c", "metadata": {}, "outputs": [], "source": [] diff --git a/doc/src/week37/programs/codeexamplesscaling.do.txt b/doc/src/week37/programs/dill.do.txt similarity index 59% rename from doc/src/week37/programs/codeexamplesscaling.do.txt rename to doc/src/week37/programs/dill.do.txt index 41d94b7e2..4a7448043 100644 --- a/doc/src/week37/programs/codeexamplesscaling.do.txt +++ b/doc/src/week37/programs/dill.do.txt @@ -7,9 +7,22 @@ DATE: today ===== This note contains code examples with a simple scaling ===== -The programs here use both ordinrary least squares and Ridge regression with one value only for -the hyperparameter $\lambda$. The first example has no scaling and includes the intercept as well and we are trying to fit a second-order -polynomial. +The programs here use both ordinrary least squares (OLS) and Ridge +regression with one value only for the hyperparameter $\lambda$. The +first example has no scaling and includes the intercept as well and we +are trying to fit a second-order polynomial. The second code takes out +the intercept and subtracts the mean values of each column of the +design matrix and the mean value of the outputs. + +The third and final code uses _Scikit-Learn_ as library in order to +calculate the optimal parameters for OLS and Ridge regression. Note +that it is highly recommended to not include the intercept in Ridge +and Lasso regression, in order to avoid penalizing the optimization by +the intercept. The second and third codes do thus not include the +intercept. In the second code we do the scaling ourselves while the +last code uses the standard scaler option included in _Scikit-Learn_, known as centering (where +we subtract the mean values). + !bc pycod import matplotlib.pyplot as plt @@ -34,14 +47,14 @@ def Ridge_fit_beta(X, y,L,d): np.random.seed(2018) n = 100 d = 3 +# hyperparameter lambda Lambda = 0.01 -true_beta = [2, 0.5, 3.7] -# Make data set. +# Make data set, simple second-order polynomial x = np.linspace(-3, 3, n) -y = 2 + 0.5*x + 3.7*x**2 +y = 2.0 + 0.5*x + 5.0*(x**2)+ np.random.randn(n) -#Design matrix X includes the intercept and scaling is made +# The design matrix X includes the intercept and no scaling is made X = np.zeros((len(x), d)) for p in range(d): X[:, p] = x ** (p) @@ -74,17 +87,17 @@ plt.plot(x, X @ beta_Ridge, label="Ridge_Fit") plt.grid() plt.legend() plt.show() - !ec In this example we do not include the intercept and we scale the data by subtracting the mean values. This follows the discussion in the "lecture material":"https://compphysics.github.io/MachineLearning/doc/LectureNotes/_build/html/chapter3.html#more-on-rescaling-data". see also the weekly slides "for week 36":"https://compphysics.github.io/MachineLearning/doc/pub/week36/html/._week36-bs029.html". +It is recommended whrn we use Ridge and Lasso regression to not include the intercept in the optimization process. Before we discuss the code, we repeat some of the basic math from the slides of week 36. Let us try to understand what this may imply mathematically when we -subtract the mean values, also known as *zero centering*. For +subtract the mean values, also known as *zero centering* or simply *centering*. For simplicity, we will focus on ordinary regression, as done in the above example. The cost/loss function for regression is @@ -208,137 +221,14 @@ Now we try to implement this. !bc pycod -np.random.seed(2018) -n = 100 -# we do not include the intercept -d = 2 -Lambda = 0.01 -true_beta = [2, 0.5, 3.7] - -# Make data set. -x = np.linspace(-3, 3, n) -y = 2 + 0.5*x + 3.7*x**2 - -#Design matrix X does not include the intercept. -X = np.zeros((len(x), d)) -for p in range(d): - X[:, p] = x ** (p+1) - - -#Split data in train and test -X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) - -# Scale data by subtracting mean value,own implementation -#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable -X_train_mean = np.mean(X_train,axis=0) -#Center by removing mean from each feature -X_train_scaled = X_train - X_train_mean -X_test_scaled = X_test - X_train_mean -#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered, note) -y_scaler = np.mean(y_train) -y_train_scaled = y_train - y_scaler - - -#Calculate beta -beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled) -beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d) -print(beta_OLS) -print(beta_Ridge) -# calculate intercepts and print them -interceptOLS = y_scaler - X_train_mean @ beta_OLS -interceptRidge = y_scaler - X_train_mean @ beta_Ridge -print(interceptOLS) -print(interceptRidge) - -#predict value with intercept -ytilde_test_OLS = X_test_scaled @ beta_OLS+y_scaler -ytilde_test_Ridge = X_test_scaled @ beta_Ridge+y_scaler - - -#Calculate MSE - -print(" ") -print("test MSE of OLS:") -print(MSE(y_test,ytilde_test_OLS)) -print(" ") -print("test MSE of Ridge") -print(MSE(y_test,ytilde_test_Ridge)) - - -plt.scatter(x,y,label='Data') -plt.plot(x, X @ beta_OLS+interceptOLS,'*', label="OLS_Fit") -plt.plot(x, X @ beta_Ridge+interceptRidge, label="Ridge_Fit") -plt.grid() -plt.legend() -plt.show() - !ec -We see that we get the same values for the parameters! As it should be. The MSE may however change (not the case here). Finally, instead of using our own function we repeat the same example using the _standardscaler_ functionality of the library _Scikit-Learn_. Here we limit ourselves to Ridge regression only. !bc pycod - -from sklearn import linear_model -np.random.seed(2018) -n = 100 -d = 2 -Lambda = 0.01 -true_beta = [2, 0.5, 3.7] - -# Make data set. -x = np.linspace(-3, 3, n) -y = (2 + 0.5*x + 3.7*x**2) - -#Design matrix X does not include the intercept. -X = np.zeros((n, d)) -for p in range(d): - X[:, p] = x ** (p+1) - -#Split data in train and test -X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) - - -# Scale data by subtracting mean value using scikit-learn -from sklearn.preprocessing import StandardScaler -scaler = StandardScaler() -#scaler.fit(X_train) -#scaler.fit(y_train) -#X_train_scaled = scaler.transform(X_train) -#X_test_scaled = scaler.transform(X_test) -#y_train_scaled = scaler.transform(y_train) - -#Calculate beta -OLS = LinearRegression() -OLS.fit(X_train,y_train) -ypredictOLS = OLS.predict(X_test) -RegRidge = linear_model.Ridge(Lambda) -RegRidge.fit(X_train,y_train) -ypredictRidge = RegRidge.predict(X_test) -print(OLS.coef_) -print(RegRidge.coef_) -print(OLS.intercept_) -interceptRidge = RegRidge.intercept_ -print(RegRidge.intercept_) -#predict value without intercept -ytilde_test_Ridge = X_test @ RegRidge.coef_+ RegRidge.intercept_ -ytilde_test_OLS = X_test @ OLS.coef_+ OLS.intercept_ - -#Calculate MSE -print(" ") -print("test MSE of OLS") -print(MSE(y_test,ytilde_test_OLS)) -print(" ") -print("test MSE of Ridge") -print(MSE(y_test,ytilde_test_Ridge)) -plt.scatter(x,y,label='Data') -plt.plot(x, X @ RegRidge.coef_ + RegRidge.intercept_ , label="Ridge_Fit") -plt.grid() -plt.legend() -plt.show() !ec @@ -346,6 +236,3 @@ plt.show() - - -