diff --git a/doc/src/week37/programs/dill.do.txt b/doc/src/week37/programs/codeexamplesscaling.do.txt similarity index 68% rename from doc/src/week37/programs/dill.do.txt rename to doc/src/week37/programs/codeexamplesscaling.do.txt index 4a7448043..abd6c4c70 100644 --- a/doc/src/week37/programs/dill.do.txt +++ b/doc/src/week37/programs/codeexamplesscaling.do.txt @@ -221,14 +221,130 @@ 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 + +# Make data set. +x = np.linspace(-3, 3, n) +y = 2.0 + 0.5*x + 5.0*(x**2)+ np.random.randn(n) + +#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 + 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 + +# Make data set. +x = np.linspace(-3, 3, n) +y = 2.0 + 0.5*x + 5.0*(x**2)+ np.random.randn(n) + +# 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) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) + +#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 @@ -236,3 +352,11 @@ _Scikit-Learn_. Here we limit ourselves to Ridge regression only. + + + + + + + + diff --git a/doc/src/week37/programs/codeexamplesscaling.ipynb b/doc/src/week37/programs/codeexamplesscaling.ipynb index a4b261ba5..407aea948 100644 --- a/doc/src/week37/programs/codeexamplesscaling.ipynb +++ b/doc/src/week37/programs/codeexamplesscaling.ipynb @@ -2,8 +2,10 @@ "cells": [ { "cell_type": "markdown", - "id": "c385e1d6", - "metadata": {}, + "id": "11ec219f", + "metadata": { + "editable": true + }, "source": [ "\n", @@ -12,8 +14,10 @@ }, { "cell_type": "markdown", - "id": "b1745002", - "metadata": {}, + "id": "f5df87ed", + "metadata": { + "editable": true + }, "source": [ "# Scaling examples with own code and the library Scikit-Learn\n", "**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University\n", @@ -25,8 +29,10 @@ }, { "cell_type": "markdown", - "id": "52e4b0f1", - "metadata": {}, + "id": "95d166ea", + "metadata": { + "editable": true + }, "source": [ "## This note contains code examples with a simple scaling\n", "\n", @@ -49,35 +55,13 @@ }, { "cell_type": "code", - "execution_count": 37, - "id": "fc969dee", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[1.79934087 0.47179152 5.01549939]\n", - "[1.79909592 0.47176716 5.01550546]\n", - " \n", - "test MSE of OLS:\n", - "1.13943111290393\n", - " \n", - "test MSE of Ridge\n", - "1.1395235273363686\n" - ] - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAACDY0lEQVR4nO3dd3gUVRfA4d/MbElCCiQBEkILTYSggEhXQARBaVYEROxIUbFhQSVYQLChINg+EUWKjaZSpYmA1AABQUooCYmUkEJCkt2d+f4Iu5BKAkk25bzPwwM7c3f25mbZPXPLuYphGAZCCCGEECVEdXcFhBBCCFGxSPAhhBBCiBIlwYcQQgghSpQEH0IIIYQoURJ8CCGEEKJESfAhhBBCiBIlwYcQQgghSpQEH0IIIYQoUSZ3VyA7Xdc5ceIEPj4+KIri7uoIIYQQogAMwyA5OZkaNWqgqvn3bZS64OPEiRPUqlXL3dUQQgghxBU4fvw4NWvWzLdMqQs+fHx8gMzK+/r6Fum1bTYby5cvp3v37pjN5iK9dnkjbVVw0lYFJ21VONJeBSdtVXDF1VZJSUnUqlXL9T2en1IXfDiHWnx9fYsl+PDy8sLX11fenJchbVVw0lYFJ21VONJeBSdtVXDF3VYFmTIhE06FEEIIUaIKFXyEh4ejKEqWP0FBQa7zhmEQHh5OjRo18PT0pHPnzuzZs6fIKy2EEEKIsqvQPR9NmzYlNjbW9Wf37t2uc5MmTeLDDz9k6tSpbNmyhaCgILp160ZycnKRVloIIYQQZVehgw+TyURQUJDrT9WqVYHMXo/JkyczZswY7rrrLsLCwpg5cyapqanMnj27yCsuhBBCiLKp0BNODxw4QI0aNbBarbRp04bx48dTr149oqKiiIuLo3v37q6yVquVTp06sWHDBoYOHZrr9dLT00lPT3c9TkpKAjInxNhstsJWL1/O6xX1dcsjaauCk7YqOGmrwpH2Kjhpq4IrrrYqzPUUwzCMghZesmQJqampNGrUiP/++4+3336bffv2sWfPHvbv30+HDh2IiYmhRo0aruc88cQTHD16lGXLluV6zfDwcMaNG5fj+OzZs/Hy8irwDyKEEEII90lNTWXgwIEkJiZedrVqoYKP7FJSUqhfvz6jR4+mbdu2dOjQgRMnThAcHOwq8/jjj3P8+HGWLl2a6zVy6/moVasWp0+fLpaltitWrKBbt26yFOsypK0KTtqq4KStCkfaq+CkrQquuNoqKSmJwMDAAgUfV5Xno1KlSjRr1owDBw7Qr18/AOLi4rIEHydPnqR69ep5XsNqtWK1WnMcN5vNxfYGKs5rlzfSVgUnbVVw0laFI+1VcNJWBVfUbVWYa11Vno/09HT++ecfgoODCQ0NJSgoiBUrVrjOZ2RksHbtWtq3b381LyOEEEKIIuDQDTZHxQOwOSoeh37Fgx9XpVA9Hy+88AK9e/emdu3anDx5krfffpukpCSGDBmCoiiMGjWK8ePH07BhQxo2bMj48ePx8vJi4MCBxVV/IYQQQhTA0shYxi3eS/y580xqDY/M3IK/tydjezehR1jw5S9QhAoVfERHRzNgwABOnz5N1apVadu2LZs2baJOnToAjB49mvPnzzN8+HDOnj1LmzZtWL58eYHyvAshhBCieCyNjGXYrO0YgFW7eDwuMY1hs7Yz/YGWJRqAFCr4mDt3br7nFUUhPDyc8PDwq6mTEEIIIYqIQzcYt3gvuQ2wGIACjFu8l25NgtDUy+/LUhRkbxchhBCiHNscFU9sYprrsW7A1D0qzukeBhCbmOaaC1ISJPgQQgghyrGTyWlZHusGHEi6GHzkVa44SfAhhBBClGPVfDyyPHYYWf/Oq1xxkuBDCCGEKMdah/pf8sjgffNnDNJW4kVaPuWKlwQfQgghRDmmqQoPd6gLwLXKMe7R/uQN07eYcLjKPNyhbolNNoWrzHAqhBBCiNJvbO+mBPt6wIrvAfhDb0kS3gC82rMxT3SqX6L1keBDCCGEqADa16tCVe0vABY4OqCQudKlfYPAEq+LBB9CCCFEBRB8djMBSgJJig8Bda4jLMOXuMR0ArwtJV4XCT6EEEKICiDg0AIAKrW4h7aGxriebTBUDatJy/+JxUAmnAohhBDlXUYK7F0EgNHsXiAzK7k7Ag+Qng8hhBCiXHLuYHsyOY0mp5fR0JYCVepihNwIu5a4tW4SfAghhBDljHMHW2da9W/M39BQg4NBt1NHKbkltXmRYRchhBCiHHHuYOsMPAJJ5CZ1FwCPR9Rn5T//ubN6gAQfQgghRLmR2w62fbQNaIrBDr0BR4xg3l2yz231c5LgQwghhCgnsu9gC9BPWw/AL46OGEBcUsltIJcXCT6EEEKIciL7zrT1lRiuU6OwGRq/OtoCYBi5PbNkSfAhhBBClBPZd6a980Kvxxr9es7iC0CGLhNOhRBCCFFEWof6U9nTDICC7go+Fjg6usoEeLl/oasEH0IIIUQ5oakK797dDIC26j+EKGdIMrxYqbd0lXm7X1N3Vc9Fgg8hhBCiHOkRFsyrPRtzt/YnAL862pJO5v4tr/ZszK3XVndn9YAKFnzsjklk6h6V3TGJ7q6KEEIIUWz6NPHjdm0zAGcb3U29wEoEelvp3byGm2uWyf0DPyXEoRt8ti6KA0kqn6+LYlrtADTV/ZNuhBBCiKIWdGIlkIZRJZQRgwcxHMhw6FhNGjabzd3VK//BR/TZVH7bFcuXfx4m+VwyzZQYlu01aDN+JY/fVI87rgumZhUvd1dTCCGEKDoRswFQrh8AioICbttELjflPvjoOHE1kLnWeZX1DQwUWqd/yulzChOW7GPCkn0cefcON9dSCCGEKCKJ0RC1LvPf1/d3b13yUK7nfDh0w7XkKMoI5hye+CmpdFV3uMpU9jTj0EtBxhUhhBCiKOyaBxhQpwNUqevu2uSqXAcfm6PiSTifObalo/KL4yYA7tHWucoknLexOSreLfUTQgghitKu42eJWTsj88H1A9xbmXyU6+Aje5rZXy4kWblZ3UVVzuZZTgghhCiLNq1fSYj9ODbFCk36urs6eSrXwUf2NLNRRg2OmBtgUnT6aX+5jh/47xwbD52R4RchhBBlTvTZVHZHJ7LzeAI++38EYLlxIztP6eyOTiT6bKqba5hTuZ5w2jrUn2A/D+IS0zAAswoJQTfB8YPcrf3Jl447AIWpqw8ydfVBgnytDGhdm7qBlajm40HrUH9ZjiuEEKJUcy6sMGNns/UvUGBeRgfWfXrxJru0Lawo1z0fmqowtncTABRAUSCmcmvSDTON1eOEKVFZysclpfPRygM8MzeCAV9uouPEVSyNjHVDzYUQQoiCebhDXQBuUbdTRTnHf0Zl/tLDcpwvTcp18AGZaWanP9CSIL/MIRi7qRIrjFYArtSzeYlLTGPYrO0SgAghhCiVHLrB0sg44OJiivmOm3BwMafH0si4UjetoNwHH5AZgKx/6Ra+HnIjAAv0zFUvfbW/MGPP83nOX9W4xXtL3S9OCCGE2BwVT2xiGoEk0kWNAOBHx81ZysQmppW6VZ0VIviAzCGY1qH+APxlNCPOqIK/co5bLsn5kRuD0vmLE0IIIZyrNftp6zEpOjv0BhwyQvIsV1pUmODjUjoqCy4su70050d+StsvTgghhMhc1Wm4vst+ytbrkbVc6VEhg4/ASlZ+upBwrLMaQQCX3+W2tP3ihBBCiNah/nT2OUFj9TjphpnFjnZZzitAsJ+Hq+e/tKiQwUdMUgYHjZpE6PUwKw76ahvyLFtaf3FCCCGEpiqMq505fWCZ3ookKrnOORNFjO3dpNSljaiQwcf79zRDVeAnRycA7tXWcnF66UWl+RcnhBBCYE+nTszvAKzyuDXLqSA/D6Y/0JIeYcHuqFm+ynWSsbz0vT6YxsF+DJxyjtdNs7hWPUZT5Qh7jNAs5YL8PBjbu0mp/MUJIYQQ7F8CaQngU4MPnhlF/6OJnExOK/WJMitk8OGUhDfL9Fb00TbS37SGN2yhjL+zGZWsWqn/xQkhhKjYdkUn4Fj0KS0Arr8fzWSiXf0Ad1erQCrksAtAgLeFqt5Wtla5HYA7TRsJqaTQpXFV+jYPoV39AAk8hBBClFrL/97JdWlbMh80H+TeyhRShe35CPbzZP3LXbAoneHjz/FJimZNnxTMfp7urpoQQgiRq+izqZxNsaEoYI78AU0xiOAaTOlVMaITqVLJTM0qXu6u5mVV2OADwGq6kH62+QBY9x7mXXPg+nvdWykhhBAiD85N5MBgpWUVqPCDrSOzp6x3lSltm8jlpsIOu2TRfGDm34dWQWK0e+sihBBC5GFy/+aYVIWWygEaqCc4b1hYdCG3h0lVmNy/uXsrWEASfAD414M6HQEDds5xd22EEEKIXPVrEcKCER0upIiA3/U2nCNzmGXBiA70a5EztXppJMGHU4sHAEjf+h0DPt/IrugE99ZHCCGEyIVqS6G3thGAefbOKGVwbYQEH05N+oDFB2vSURxH/uKX7THurpEQQgiRQ40Ty/BW0jih1qBv33toFuJHVW8rAd4Wd1etwCr0hFOnzNnDdkJC78B//1zuM61lws7rueeGmhgGZWb2sBBCiPKv8r55AAR3fpRBbesysE0dMhz6xUUUZYAEH1ycPdxSacwvVrhd/ZvwlER6lbHZw0IIIcq50wfh2EYMRWV7lZ5ER8S4kmKWJTLswsXZw9uNhhzUa+ClpHOHtgkoW7OHhRBClHM7vgNgg9Kcu78/wjNzIxjw5SY6TlzF0shYN1eu4CT44OLsYVCY5+gMwP1aZm9IWZo9LIQQovzadew0ZzbMBODbtJuznItLTGPYrO1lJgCR4COb+fpN2AyNFupBGinH3V0dIYQQAodusPa3OQQYZzlt+LJKb5nlvHNf9nGL9+LQc+7SXtpI8HGBc6+XGiG1iQ3qDMBDHn+WqdnDQgghypfos6l8vvYQbcavpNGJBQDMd3TElsuUTQOITUxjc1R8yVbyCsiE0wtce71oKsqBVJj9BwM8/kKpJPGZEEII93AuiAgkkVusOwD44cL0gLycTE4r7mpdNflmvYTVpKEoCjToCj41UM6fhX2/urtaQgghKiCHblDZ0wzAXdo6zIqDCL0+B4ya+T6vmo9HSVTvqkjwkRtVc2U8Zft37q2LEEKICmlzVDwJ522AQX9tDQBzHLfkWV4Bgv3KxrJbCT7y0mJQ5t+HV8PZI26tihBCiIrHOXzSWtlHfTWWFMPKr462uZZ1Zlgf27sJmlr6861L8JGXKnWhXmcA/lv3PwZ8sUn2exFCCFFinMMn95sy530scrQnBc9cywb5eTD9gZb0CAsusfpdDZlwmp8Wg+HwGjwi5/J3cmt+2e7DdTUru7tWQgghKoDWof408rVze/rfAMx1dMlRprKnmU8HtaRtvYAy0ePhdFU9HxMmTEBRFEaNGuU6ZhgG4eHh1KhRA09PTzp37syePXuutp4lLvpsKpE+N2G3VsbPdpKb1V0s3nmCyJhEdkcnEn021d1VFEIIUY5pqsLkJv/iodjYp9dip1HfdU658Ofdu5vRoUFgmQo84CqCjy1btvDFF19w3XXXZTk+adIkPvzwQ6ZOncqWLVsICgqiW7duJCcnX3VlS1LHiavp9dlWvk3JHF+7X1tNfEoGvaasp/fU9a7lT0IIIUSxMAyaxC4A4DdLdy7O7Ch7wyzZXVHwce7cOQYNGsSXX35JlSpVXMcNw2Dy5MmMGTOGu+66i7CwMGbOnElqaiqzZ88uskqXBOd+L85urq7qdgI5C8h+L0IIIUpAzHb4LxI0K6OefZ05j7fl4/ubM+fxtqx/6ZYyG3jAFc75GDFiBHfccQe33norb7/9tut4VFQUcXFxdO/e3XXMarXSqVMnNmzYwNChQ3NcKz09nfT0dNfjpKQkAGw2Gzab7Uqqlyfn9Qpy3TvCqlHXvw39phts0xtyg3qAe7V1THP05aehbWhaw7fI61eaFKatKjppq4KTtiocaa+CK09t5dANftwWjf/qd7kDcDTujW7xplVtAF8AdIcd3XFl1y+utirM9QodfMydO5ft27ezZcuWHOfi4uIAqF69epbj1atX5+jRo7leb8KECYwbNy7H8eXLl+Pl5VXY6hXIihUrClTu+DkAE3McXbhBPUB/bTXTHb1Zv349R72LpWqlTkHbSkhbFYa0VeFIexVceWmrdYcz+ChjHSiwMb0RZ37/vchfo6jbKjW14HMhCxV8HD9+nGeeeYbly5fj4ZF3BjVFyTrxxTCMHMecXnnlFZ577jnX46SkJGrVqkX37t3x9fUtTPUuy2azsWLFCrp164bZbL5s+djENL6J2sQR3+5kJM2mjv0kPb3+pXaT+zAMg0BvKzfUqVLmJvoURGHbqiKTtio4aavCkfYquPLQVvO2HufNxXsxgLuUTXib0zisB/HI3iaAwhu9m9C/Va2rfp3iaivnyEVBFCr42LZtGydPnuSGG25wHXM4HKxbt46pU6eyf/9+ILMHJDj44ljUyZMnc/SGOFmtVqxWa47jZrO52N5ABb127UAz60Z3IeJYAjF/9SY0ag53GSt4bF5jV5lgPw/G9m5Spsfe8lOcv4fyRtqq4KStCkfaq+DKals5dIPXFv6Dc1LpAEvmooZ5ji5k6JnTM19b+A8D2oQW2Q1vUbdVYa5VqAmnXbt2Zffu3URERLj+tGrVikGDBhEREUG9evUICgrK0pWTkZHB2rVrad++fWFeqlRYGhlL1w/WMvCrvxm+L3NVz836ZgJIdJWJS0xj2KztLI2MdVc1hRBClHGX7kR7rXKU5uohbIbGz46b8yxXlhWq58PHx4ewsLAsxypVqkRAQIDr+KhRoxg/fjwNGzakYcOGjB8/Hi8vLwYOHFh0tS4BSyNjGTZrO8aFx/8YdYjQ69NcPcTd2jq+cPQGMrcwVoBxi/fSrUlQuRyCEUIIUbwu3Yl2oPYHAMv0VpzGL89yZVmRp1cfPXo0o0aNYvjw4bRq1YqYmBiWL1+Oj49PUb9UsXHoBuMujLtdavaFDX0GaKvgkrMGmfNDyktEKoQQomQ5U6l7kUY/7S8AZju65lmurLvq9Opr1qzJ8lhRFMLDwwkPD7/aS7vN5qh4YhNzRpe/OtrxumkWoep/tFP3slFvmuV8eYlIhRBClKzWof4E+3lw07nV+CjnidKrs1FvkqVMWdmxtiBkY7lc5BVEpOLBIkfm3JXM3o+syktEKoQQomRpqsLY3k0YdGHIZY7jFoxLvqIVys6OtQUhwUcu8gsinEMvt6lb8CfrsqLyEpEKIYQoeT0CTnK9ehgbpiwTTYPLeCr13Miutrlwdn/FJablmPexxwhllx7KdWoUd2vr+NLRC4CHO9QtNxGpEEIIN9g6AwCtaR+mtryNk8lpVPPJHGopb98v0vORC2f3F1y6jc9F3ztuBTJnJCvovNzjGro3CWJhRAwbD53BoWcPWYQQQoi8RR6O5vz2uQCorR6mXf0A+jYPoV39gHIXeID0fOSpR1gw0x9oybjFe3NMPl3saMeYCxNP26t7+Gq9B6fP7XedL++Jx4QQQhQdh24QseR/hBnniTPXpGrtjmjurlQxk56PfPQIC2b9S7e4dhKcOqAFgd5WGtSsTmydvgAM0v7g9LmMLM+TxGNCCCEuJ/psKp+vPUSb8Su57r/5AHyV2ok2E/7g87WHiD5b8L1Syhrp+bgMTVVoVz/A9bhb0+pYNBX9vyrw2Vy6q1upxllOUsVVRhKPCSGEuJyOEzNTqIcph7nOGkW6YeJnx02cPZfBhCX7mLBkH0fevcPNtSwe0vNRSFaThqIobE4JYrN+DSZFp7+2Okc5STwmhBAiO4dusPHQGeZvj6aSNXNwxbm8dqnemrNc3FC1sqe53M4hlJ6PK3QyOY1V9q60tuxngGkV0xx9ceQySieJx4QQQkDmth3Z5xH6kuLKaDrLfmuW8gnnbWyOis/S+15eSM/HFarm48FSvTXxhjc1lHg6qxF5lhNCCFGxOfcLy76A4S7tTzyVDPbrNdliXJPjeeX1BlaCjyvUOtQffz9ffnJ0Ai52mzkplK9UuEIIIa5MXvuFgeH67pjluJXckjuU1xtYCT6ukKYqjOxS35XxtLO6k5rKKdd5AxjZpb5MNhVCiAour/3C2qr/0FCNIcWwMt/RMcu58n4DK8HHVRizYA9HjGD+dIShKgb3Z9vvZcyCPW6qmRBCiNIir6GTB7SVACxwdOQcXq7jzlvW8rSXS3YSfFyFyf2bY1KVC91l0F9bjQUbACZVYXL/5m6snRBCiNIgt6GTqpzlNnULgOs7xCmoHO7lkp2sdrkK/VqE0KCaN/2m2IgzqhCknKWHuoVFensWjOhAWIifu6sohBDCzXLbL6y/tgaz4mCb3pB/jDr4VzLzeq+mBPmWz71cspOejyJgx8RsR1cABpuWu7k2QgghSpPs+4Wp6AwwZQ7Tf2+/FQUYf2cz7mxRfvdyyU6Cj6sU4G2hqreVndX6oismblT/pX2lWAK8Le6umhBCiFIiLMSPl3s2JsDbwi3qDkKUM8Qb3mzyvImXezaucD3lMuxylYL9PFn/chcsmoryU2/YM5/vrtuF5veYu6smhBCilHCmUgd4wJw50fRHRydOpFPuU6nnRno+ioAz5To3ZgYcWuSPkJbo5loJIYQoLZwLFOoocdys7gJwDddXxAUKEnwUpTodoOq1YEuFiDnuro0QQohSol+LEBaM6MAD2kpUxWC143qOGkEALBjRgX4tQtxcw5IlwUdRUhS48VEAEv+cTrcP1jBr09FyuzGQEEKIglPs57lPWwPAt47uKOV/XmmeJPgoYivMnUnBA7+UI1Q98zevLYik48RVLI2MdXfVhBBClCDnDrYLI2LYeOgMwccW46ekEqsG0a3PQJqF+FHV21ohFyjIhNMi9P2mI4xZ8C9vmm7iQdMKBmsr2KCHEZuYxpOztvNOv6YMalvX3dUUQghRzHLuYGuw3HMq/kBQ1xEMbFuPAW1CyXDoWE05d0Qv76Tno4g4dMOVTv07RzcAuqnbCOKMq8yYBXtkCEYIIcq53HawbaXsp5FxhDTDzCrPzO8IRVEqZOABEnwUmc1R8a5/HzBqstHRBJOiM8j0R57lhBBClC957WD7oGkFAAsdHXht2YkKfyMqwUcRyb5x0ExHdwAGaKuwkpFnOSGEEOVHbjvYVuUsPdXNQOZE09jEtAp/IyrBRxHJvnHQCv0GYowAApUk7lA35VlOCCFE+ZHbDeZAbRVmxcFWvRF7jLp5lqtIJPgoIs6Ng5wcaHxvz9ypcIhpOWAQ7Je5YZAQQojyKfsNpgk7Ay8Mv39r755nuYpGgo8i4tw46NJl23McXUg3zFyvHqaFcpCxvZtUiA2DhBCiosp+g3mbupXqSgKnDD+W6K3zLFfRSPBRhHqEBTP9gZauHpCz+LLI0Q6AqQ224Odpca33ruiTjYQQojzSVIWHO9R1PR5iWgbAbMct2C5kt3i4Q90KfyMqeT6KWI+wYLo1CWJzVDwnk9MIzXgBfl9HteNL6fflUk5RGYBgPw/G9m5Cj7Bg91ZYCCFEkRrbuynBvh4sXLqE1up+bMbFYfhXezbmiU713VxD95Oej2KgqQrt6gfQt3kIJ7wasVVvhBk7A7WLy27jEtMYNmu7ZD4VQohyqH2DQB6+0Ovxu96GU1RxHRcSfBQr53rvmRcmGQ0y/YEZO4BrDfi4xXtlCEYIIcqZqmoifbQNAFg7DqdZzYqbSj03MuxSjJzrvZfQmjijCkFK5lrvRXp7IDMAca73blc/wL2VFUIIUWSq/zsPsGOE3ECP23pzW3ejwqZSz430fBQj5zpuOya+t3cF4CHT0jzLCSGEKAccNtjyFQBKmycz/67AqdRzI8FHMbp0HfccR1fSDRMt1YNcrxzMs5wQQogybu9COBcH3tWhST9316ZUkuCjGF2aeOw0fiy+MNzy8CW9H5J4TAghyo9d0Qn8u+j9zAetHgWTzPHIjQQfxciZeMxphr0HAHeof1OdzLz+knhMCCHKj7//XEEj2z84FBO0etjd1Sm1JPgoZj3CgplwZxiqAnuMuvytN8asOHjQtIIJd4ZJng8hhCjjos+msjs6kciYRGr8+y0AS2hPZKKV3dGJRJ9NdXMNSx9Z7VICBrSpQ98WIUQcS0D5dxhsfobhPn+itKzm7qoJIYS4Sh0nrgYyd6/9y/oXKPBZWncip6x3lTny7h3uql6pJD0fJcTLYqJ9g0Ba9xgMlWujnI+HXT+4u1pCCCGu0uT+zVEVGGxageXC7rWRRj0AVCXzvMhKgo+SpmrQemjmvzdNB0MSjAkhRFnmYVYxGxkMupDF+n/2nq5zupF5XmQlLeIOLR4AcyU49Q9vT/2MXdEJ7q6REEKIK+DMZH2ntp4AJZloI5DleqssZSSTdU4SfLiDZ2VoMQiAtifn8cv2GPfWRwghxBXJzGR9nke0JQDMsN+Gg6zJxJyZrMVFEnyUMOes6H/rDATgVm0HERHbiIxJlFnRQghRSjl0g42HzrAwIoaNh864ejJOJqdxs7qLRmoM5wwPfnB0yfX5ksk6K1ntUsKcs6IB/mduQVdtB3dmLKLXlIt7u8isaCGEKD2WRsYybvFeYhMvBhDBfh6M7d2Eaj4ePHqh1+MHR2eS8cr1GpLJOivp+ShhzlnRAP9zZE5Kuldbhx/nZFa0EEKUMksjYxk2a3uWwAMgLjGNJ2dt53+//E4nbRe6oTDDcVuO5ytIJuvcSPBRwjzMKs55Rxv0pvyj18ZLSWeQ9ge6AUfPpOTo1hNCCFHynJNJc/skdh67JfFnAJbrrYg2qmcp48xdLZmsc5LgowQ538gXKXxpvx2AIaZlmLHz0coDPDM3ggFfbqLjxFUsjYx1T2WFEKKCy5xMmvdcjSokcZeWmUgso9VQAryz7uMS5OfB9AdaSibrXEjwUYJyeyMv1tvzn1GZ6koCvdUNWc7FJaYxbNZ2CUCEEMINLjdJdJD2Bx6KjV16KE9v9OT0uQwAPr6/OXMeb8v6l26RwCMPEnyUoNzeyDZMzLRnjhM+ZloCl3TwOf8la8SFEKLk5TdJ1IKNIablAHxt7wkomFSFyf2b07d5CO3qB8hQSz4k+ChBeb2Rv3d0JdWw0kQ9Sjt1b5ZzBrJGXAghSpJzWW1c4nn8K1nILYToq/1FVSWRE4Y/v+ptAVgwogP9WoSUbGXLKFlqW4Jah/oT7OdBXGJalglMiXjzo+NmhphW8Lj2Gxv1pjmeK2vEhRCi+OW2rDY7BZ0ntN8AmGHvgUMxkeusVJEn6fkoQZqqMLZ3E4AckfTXjp7ohsItWgT1lZwZT2WNuBBCFK+8ltVm10ndSUM1hhQ8adhzBM1C/Kjqbc0x4VTkTYKPEtYjLJjpD7QkyC9rMHHUCGKFfgMAj2q/ZzlX1ccqa8SFEKIY5bes1sm/kpmP+jfn41qZK1y82j3KfR3DWDiiA+tf7kKwn2fJVLYckODDDXqEBbP+pVuY83hbPr6/ueu4c9nt3dp6Akl0HT+VnC4Tl4QQohhdblktQHyKjfq2g/j9txFDNbE9+H4WRsSw6XA8JlW+TgtD5ny4iaYqtKufmVLdMOC5HyLYalxDhF6f5uohHjQt40P7fagKfHhfc/dWVgghyrmCzqs7tfwDAJbTjqGzjwHHgIvp1mVpbcEUKlSbPn061113Hb6+vvj6+tKuXTuWLFniOm8YBuHh4dSoUQNPT086d+7Mnj17irzS5U2/FiEsGtkRUPjM3huAwdpKPElj/vAOVPf1kKynQghRjAoyr64Gp+lk+xOAT873yHJO8jIVTqF6PmrWrMm7775LgwYNAJg5cyZ9+/Zlx44dNG3alEmTJvHhhx/yzTff0KhRI95++226devG/v378fHxKZYfoLxZYbQiSq9OqPof92lreXSmrytxDUh0LYQQxSGv1YiXesS0FJOi85ejKXuM0CznDDIXEoxbvJduTYJkqPwyCtXz0bt3b26//XYaNWpEo0aNeOedd/D29mbTpk0YhsHkyZMZM2YMd911F2FhYcycOZPU1FRmz55dXPUvNwK8LVT1ttI0pAqnmj0OwGPa75w9dz5LOYmuhRCi6OW3GhHAh1T6a5m7kn/pyH3nccnLVHBXPEPG4XAwd+5cUlJSaNeuHVFRUcTFxdG9e3dXGavVSqdOndiwYUM+VxIAwX6erH+5CwtHdOCGPiOIx5da6iluV//OUk6yngohRPHIazUiwCBtJT7Kef7VQ1ijX5/vdSQv0+UVesLp7t27adeuHWlpaXh7ezN//nyaNGniCjCqV8+6q1/16tU5evRontdLT08nPT3d9TgpKQkAm82GzWYrbPXy5bxeUV+3qKiA3a6z+WgyWx3deUb7iSfNv7Lc3pbssXj8ufNsOniy2Jbglva2Kk2krQpO2qpwpL0Krqjaqus1gXRueBPbjp7l9Ll0Ar2teKp2an43HIDP7b1xfh6bVYPcRlcCvUyl+ndWXO+rwlxPMQyjULfPGRkZHDt2jISEBH7++We++uor1q5dS0JCAh06dODEiRMEB1+cj/D4449z/Phxli5dmuv1wsPDGTduXI7js2fPxsvLqzBVK1fM9mS673kWk57B/JCXmHaqGX3q6NT2dnfNhBCiYvGJWcMtJ7/mhOFPp/TJ2NEwUHihmZ1a8pnskpqaysCBA0lMTMTX1zffsoUOPrK79dZbqV+/Pi+99BL169dn+/bttGjRwnW+b9++VK5cmZkzZ+b6/Nx6PmrVqsXp06cvW/nCstlsrFixgm7dumE2m4v02kVpc1Q8j8zcwhh1Jg9qy1mnX8eDGS+jKQamSwbKvh5yY7H2fJSFtioNpK0KTtqqcKS9Cq7Y2kp3wPR2mBMO85XXoxyuP4Q5W45jABYVlAs9H84OkI/6N+fWa6vndbVSobjaKikpicDAwAIFH1ed58MwDNLT0wkNDSUoKIgVK1a4go+MjAzWrl3LxIkT83y+1WrFarXmOG42m4vtP1txXrsotG1QDV8vD744dzuD1BXcrO7iWuUo/xh1cDgyy1T1sdK2QbVin1Fd2tuqNJG2Kjhpq8KR9iq4Im+rvUsg4TCGhx+PPh2O4uHLTY2qMm7xXuKSLt44l8WViEXdVoW5VqGCj1dffZWePXtSq1YtkpOTmTt3LmvWrGHp0qUoisKoUaMYP348DRs2pGHDhowfPx4vLy8GDhxY6B+iItNUhVPJ6UA1ftfb0FvbxFDTYkbZRrrKSNZTIYQoZoYBf00GQLnxcfDIvJvv2awG3ZsGszkqnpPJaVTz8aB1qL98JhdCoYKP//77j8GDBxMbG4ufnx/XXXcdS5cupVu3bgCMHj2a8+fPM3z4cM6ePUubNm1Yvny55Pi4ApP7N+e5HyL4zN6H3tomeqmbeF+5j2ijmmQ9FUKIYrYrOoEF83/gjTPbQLNCm6FZzl+apVoUXqGCj//973/5nlcUhfDwcMLDw6+mToLMrKcNqnnTawqsdVxHJ20XT2i/8Yb9YRaN7EhYiJ+7qyiEEOXWL9tj6HRyFmhAi0HgXc3dVSpXZCecMmC6ow8A92lrsmw4J4QQ4uo5dIONh84w468ovt90lJ3HE/gnYgNdtJ04UPm3/kPsjk4k+myqu6tabsjGcqWYM+tpql9bTjmuo2rCLkZ4riDA+053V00IIcqFpZGxjFu8N8eOth+Z54MGSxytGTkzGogG4Mi7uWc3FYUjwUcp5sx6atFUlP0vw9yBPGRZgWLNADzdXT0hhCjTlkbGMmzW9hx7udRS/qO3uhGAz+y9ADCpCu/fm39mU1FwMuxSyllNGoqiQKOeULUxSnoy23/5kG4frmXWpqOSYl0IIa6AQzcYt3hvrpvIDdV+xaTorHVcR6RRD4AFIzrQr0VIyVayHJPgo6xQVXbVfQiAWvtncOxkPK8tiKTjxFWyyZwQQhSQc37HRyv25xhqAajKWe7V1gLwqb1vrpvMiasnwUcZ8f2mI9z1Zw2ijUCqKonco60DMndQfHLWdr7fdMS9FRRCiFJuaWQsHSeuYsCXm5i6+lCuZR4z/Y5VsbNVb8RmozE1/T2p6m0lwNtSwrUt3yT4KAMcusGYBXuwY+JLe+Zkp6HaYjQcrjJjFuyRIRghhMiDc35Hbr0dTn6c4wFtJZDZ6wEKE++6jvUvdyHYT+bZFSUJPsqAzVHxrn/Pc3TmjOFDbfUUvS5MiMqtnBBCiEz5ze+41EPaMiop6fyj12aN3pxgPw/a1AvAatJKpJ4ViQQfZcDJ5IuRehpW/mfvCcBI00IU9FzLCSGEyLQ5Kj7fHg8AL9J42JS5+7qz12Ns7yaSMr2YSPBRBlTz8cjy+DtHd5IMLxqqMdymbs2znBBCiILdmA3Q/qCykkKUXp0d3jcz/YGWZWqTuLJGgo8yoHWoP8F+FwOLZLyY4bgNgJGmBYBBsF/mxkZCCCGyutyNmQUbT5h+A+Ab9U6mDW5NSGUvyWhajCT4KAM0NbP779LOvxn2HqQYVsLUI3RRI6R7UAgh8uC8gcvrE/JebS3VlQROGP7MSWtP30//ovfU9XScuLpE61mRSPBRRvQIC2b6Ay1dPSAJ+DDLcSsAHwUtp0fTIHdWTwghSi3nDRyQIwAxYWeYaREAX9h7kXEh8bdJVZjcv3kJ1rJikfTqZUiPsGC6NQlic1Q8J5PTCDHVx5i/ksrxOzm05Xdeiwjgldsbc13Nyu6uqhBClCrOG7js+7g86LWRmvppThl+zHHc4jq+YEQH2T28GEnwUcZoqkK7+gEXDxwbApu/QPnzAzaeep5ftsdI8CGEELnIfgNXrZKJlotegWT4wtGLdCwoChiSMqnYybBLGRZ9NpV99R7GUEzUS95GS+VfFu88QWRMomz/LIQQuXDewPVtHkK71DVYk49yFl8iqt/FO3eG0SzETzKalgDp+SjDnJOh3jV15H7TGp4yzeeRlEb0mrLeVUa2fxZCiFzoDlj3HgA+XZ7hh5u7oigKA1vXJsOhS2KxYiY9H2XY5P7NMakK0x19sBsqXbSdXKccBEBTFAa3rc3GQ2ck7boQQlxiV3QCH095D84cAI/KmNo8kbl7OKAoigQeJUCCjzKsX4sQFozowFEjiIV6BwCeNs0HwGEYfLfpGAO+3CQ73wohxCXmbzvObWdmZT5oOxw8fN1boQpIgo9yYqqjHw5Doau2g2bK4Szn4hLTGDZruwQgQogKK/psKrujE4mMSSR550Iaq8c5hxd7aw+QOXJuIHM+yrgAbwtVva14+zXm9zMd6c2fPG36hcdtL7jKGGSubR+3eC/dmgRJMjIhRIVzMWGYwa+WH0GFb+zdef+L3a4yMkeu5EjPRxkX7OfJ2tGd6Xt9MB+l98VhKHTTttNUicpSzgBiE9P4aMW/Mg9ECFHhOOfIdVO3EaYe4ZzhwVcXNumUhGIlT3o+yrilkbGXJM2pwSK9PXdqfzHK9AuP257PUX7q6oNMXX2QYD8PxvZuIhsnCSEqhH4tQmhQtRLKFy8D8I3jNhLwASShmDtIz0cZtjQylmGztmfJ1jfV3g/dUOimbaOpciTP58o8ECFEReNzZBlN1aMkG558Zb8dRUag3UaCjzLKoRuMW7yX7IMnh4wQFuvtAHja9Euez3c+b9zivTIEI4Qo/3SdGhGTAfjdqw8v3tlOEoq5kQy7lFGbo+Kz9Hhc6hP7nfRWN3KbtpUm9iPsNermWs45D2RzVHzWlO1CCFHe7PsV8+m9GFYf7nvqXRQvf0ko5kbS81FGnUzOPfCAzN6PX/W2ADxr+vmqriWEEGXdruPxHPvlDQCUNsNQvPwz/y0JxdxGgo8yqpqPR77nP7bfdWHly7YceT8Key0hhCjL9q36ntr2KNLUStBuuLurI5Dgo8xqHepPsJ8Hec2XOmSEsOBC1tPnTD/mWkYBgv08aB3qXzyVFEIIN3ElFYs+S8uoLwD41ridyHhVkoqVAjLno4zSVIWxvZswbNZ2FMgx8RQgqskI7Ps30EXbSUv7v2w3GrnOOYOWsb2bSNIxIUS540wqdoe6iU8tx0gyvJia1p3xsvFmqSA9H2VYj7Bgpj/QkiC/rMMmQb5Wpg9qwQsDbkdrMQiAlz2yrnwJ8vNg+gMtJc+HEKJcmty/ORbV4FnTTwB87ehBEpUASSpWGkjPRxnXIyyYbk2C2BwVz8nkNKr5ZA6jOHszlE6jYedcWuu7+OLmNN7bF8iQ9nW5r1Utth09y8KImBzPEUKIsq5fixBaJSyl5toTnDW8+Z/9dtc5SSrmfhJ8lAOaquS9VLZybWj5IGz9H432fsKB0y/yxz//8enqg1mW6krGUyFEWeLQjSw3XS1q+mQtYM+g2vaPAPjc0YtkvFAUMCStUakgwUc5F302leSGT3DNjlnUPRdBBzWS1fub5SjnzHgqQzFCiNIu67YSmepUsfJc40sKRczCknycM1RmW/V7ead1Q+ZtOU5sQpokFSsFJPgo55yTrsaaOvOwaRnPm37kr4wwyLZOJvvOt0IIURo5t5XI3oHxX1JmILLyn//oeW0grH0PAL/uL/FDu64oiiJJxUoRmXBazjl3cpxm78t5w0JL9SBd1e25lr0046kQQpQ2eW0rAVm3jJg55Q1IPgG+NTG1fhTlwiYuklSs9JDgo5zr1yKEBSM6cIrKzHD0AOAF0w8o6Hk+RzKeCiFKo/y2lXBKOZfEHYlzMh90ehFM1hKomSgsCT4qkM/tvUgyvLhWPU4fdUOe5STjqRCiNMrvxkg34Pg5GKwuI1BJ4jhBRFbrJQnFSikJPiqAAG8LVb2t1A4J4RulDwDPmX7ChD1LOcl4KoQozfK7MbLpCl/sTmeo6VcAPsi4i16f/k3vqetdc99E6SHBRwUQ7OfJ+pe7sGhkB67tN5pThi911JP019a4ykjGUyFEaZf/thIGT5oW46eksl+vySK9PQCqgiQUK4Uk+KggrCYNRVHo1rw+p1o8BcDTpl/wIB2QjKdCiNLPua0EZF+vByHqWR7WlgLwnr0/+oWvN90AD7N81ZU28hupgJr0egbDrxbVlQQe1JbTo2kQ61+6RQIPIUSpl9e2EiO0+XgqGWzVG7JSb5nl3KvzdzN/RwwbD53BoUuWsdJA8nxUMNFnUzmbYqNK81HUXPs8w0yL6RfVg39iG2AYUKWSmZpVvNxdTSGEyNOl20r8dfAUv65Zz31q5ryOSbb7yd4vEp9i49l5EYBkcy4tJPioYJwTr1Sqs8wSQkM1hvsyfqHXFLOrjOz0KIQo7TRVoZJVY9me/3je9CMmRSfW53oiMhrnvs33BZLNuXSQYZcKxpl0TEflPft9ADyiLaEaZ2WnRyFEmfLL9hgsp3bTW9uEbijsq3EvymXmy1+ajEyGYNxHgo8Kxpl0DGC53oqteiM8lQxGmX5mwYgO9GsR4uYaCiFE3qLPprI7OpHImEQW7zzBaNM8ABbr7dmj1y7QxnGSzdn9ZNilAlMUhXdt9/OT9U3u09ZwOOEghNzg7moJIUQOzl1sB3y5yXWsnbqHTpZd2AyN9+33cHx34b7SJJuz+0jPRwXkTDrWLMSPO/vdw9+WNpgUndo73nd31YQQIoelkbF0nLgqS+ChoPOyKTON+mzHLRw3qqMqBial4EMpks3ZfaTnowJyJh2zaCqKomCEfowxvT0eB3/n4LaVvL7Nmxe7N3B3NYUQIs9dbHupm7hePcw5w4Mp9rsAeC7MQYsb23LynI23fvuHsykZuc49VcjMbSTZnN1Hej4qKGfSMQCl2rUoLR4AwLJqHBsPn2ZBRKw7qyeEEHnuYmvB5prr8Zm9N2cUP9e51qH+3NmyJuPvDANyJiOTbM6lgwQfguizqfxzzUh0zYPaKbu4Vd3Ob7vjOH4OImOSZFMmIYRb5LWL7WBtBbXUU8QZVfjKcTuhAZUI9LbgczFjQJ7JyCSbc+kgwy7ClfvjRVN3RpgW8ZJpLj1SmvP+bhPv784cY5XcH0KIkpbbhFBfzvGUaT4AH9rvIQ0rT3dtwK2Nq/LH8qVZyl6ajOxkchrVfDKHWqTHw/2k50O4cn98Zu9DvOFNQzWG+7XMgEST3B9CCDfJbULocNMiKisp7Ndr8pOjEwDVfT2xmnL/OtNUhXb1A+jbPIR29QMk8CglJPgQrtwfyXgx2X43AM+afsKbVH4e2kZyfwgh3CL7LrYhnOJhbRkAE+wDMFAJlomjZZIEHyKLOXpXDunBBCpJDDMtcnd1hBAVkEM32HjoDL/uOsHtYcGuCafPm3/Eqtj4y9GUNXpzDGBkl/rSm1EGyZwPAWTm/gj0tuLrUYk//UdS/9gYHtWWcEYfAwS4u3pCiApiaWQs4xbvzTHRtJlymLu09UBmr4dz3cqYBXsY1LZuCddSXC3p+RAA7DyegEmFw6dTCP+3Ln/r1+Kh2EhZMpYBX2xiV3SCu6sohCjnnDk9cq5wMXjd/B0Avzg6EmnUA5D9qMowCT6E6z98XFL6hSMKEx2DAGh0cikpUZuZtvqQbMIkhCg2eeX0AOipbqa1up/zhoVJtv6u47IfVdlVqOBjwoQJ3Hjjjfj4+FCtWjX69evH/v37s5QxDIPw8HBq1KiBp6cnnTt3Zs+ePUVaaVF08voPv9sIZZ9P5gZ0Y8zfs3RPLG3Gr+TztYck74cQosjlldPDSgavmGYD8IWjF3EEXHbnWlH6FSr4WLt2LSNGjGDTpk2sWLECu91O9+7dSUlJcZWZNGkSH374IVOnTmXLli0EBQXRrVs3kpOTi7zy4url9R/epis8fOp+0gwzbdR93KZu5fS5DCYs2efKCyKEEEUlr03eHtKWUftCQrHP7L24r1VNmoX4UdXbSoC3pYRrKYpKoSacLl2aNYHLjBkzqFatGtu2bePmm2/GMAwmT57MmDFjuOuuzFz7M2fOpHr16syePZuhQ4cWXc1Fkch7V0eDWAL40nEHT5kW8Krpe1ZnNCcDM5U9zTh0Q2aYCyGKTG45PQJIZKRpAQDv2fpzHg/ubFGTtvX8yXDoWE1aCddSFJWrWu2SmJgIgL9/5hrrqKgo4uLi6N69u6uM1WqlU6dObNiwIdfgIz09nfT0dNfjpKQkAGw2Gzab7Wqql4PzekV93bIs0MuEVcs5ympSDFLsCtPtfbhPW0Md9SSPaEv4n9Gb8xkZbDp4UtbWXyDvq4KTtiqcitReLWr6UKeKlf+S0jAAuw7PaT/ho5wn0qjL70oH6lSx0qKmD3a7HRWw2XTX8ytSW12t4mqrwlxPMQzjimYRGoZB3759OXv2LH/++ScAGzZsoEOHDsTExFCjRg1X2SeeeIKjR4+ybNmyHNcJDw9n3LhxOY7Pnj0bLy+vK6maKALHz8H7u00oGNyp/smHls84Z3jwQ71JBFSu7O7qCSHKofh0SLnw/fXHPzH8pL2Kphj8UvNVTng2ppIZ/K3uraPIW2pqKgMHDiQxMRFfX998y15xz8fIkSPZtWsX69evz3FOyTYbyDCMHMecXnnlFZ577jnX46SkJGrVqkX37t0vW/nCstlsrFixgm7dumE2my//hApi5T//8ey8CADXxFOzYuBjNkixwWKjAw/qK2iuHsL34M+M1p/g6yE3Ss/HBfK+Kjhpq8KpaO3V8PXlF/5l8J35ezTFYInjRp47GOYqc+Ct7rk+t6K11dUorrZyjlwUxBUFH0899RSLFi1i3bp11KxZ03U8KCgIgLi4OIKDL+4YePLkSapXr57rtaxWK1ZrzlDWbDYX2xuoOK9dFvW8riaKqmVN7KPB+BscvLJFI0PXeNM2mF+s4dylrmOBuSdT1gTwqpeV62pWdmvdSxN5XxWctFXhVJT2mty/OS/8uJOubOEmLZJ0w8x4+0AgM6fH+/def9l2qChtVRSKuq0Kc61CrXYxDIORI0fyyy+/sGrVKkJDQ7OcDw0NJSgoiBUrVriOZWRksHbtWtq3b1+YlxIlrEdYMOtfuoU5j7flkQ51ATCpuJa0bTcascDRAVUxeMP8HZsOn+GX7THuq7AQotzp1yKEhUNv4DXTLAC+cNzBcSPzxlVyepQvher5GDFiBLNnz2bhwoX4+PgQFxcHgJ+fH56eniiKwqhRoxg/fjwNGzakYcOGjB8/Hi8vLwYOHFgsP4AoOs7dH9vVD+DGOn5kRG3Lcv4Ly4PcYWyjUfoeeqsbWbzTyj031MQwoEolMzWryBwdIcTVCdz1BdXVU8Qa/kyz90FR4MpmJorSrFDBx/Tp0wHo3LlzluMzZszgoYceAmD06NGcP3+e4cOHc/bsWdq0acPy5cvx8fEpkgqLknHrtdX5PQq+HnIjp1PtPDM3gr0pPnyi9eJ580+8Yp7NrSkt6TXl4pyfj+9vTjWfzB0mZRmuEKIwdkUn8PnidXxyeioA3/s+ymsdWzFvy3FiE9Ikp0c5U6jgoyALYxRFITw8nPDw8CutkyhFWof6YzabMQx44cedfOHoRX/TGmoqpxlmWsQH9vtcZZ+ZGwFAsJ8HY3s3oUdYcB5XFUKIrH7ZHsNtJ6ahaWnotdry/MOvoqgqA1vXlpwe5ZDs7SIKpF+LEBaM6EA6Ft6yDQbgCe1X6ihxOcrGJaYxbNZ2lkbGlnQ1hRBlSPTZVHZHJxIZk0h0xEr6aBvRUTjc6g0iTyQTfTYVRVEk8CiHrirJmKiYluutWOu4jk7aLsaavuUR24s4t7eGzOW6CjBu8V66NQmSIRghRK6cWzWo6Pxq+R+oMMd+C2PmJACZQ7pH3r3DfRUUxUZ6PkSB7Y5OQFXAQCHcPoQMQ+MWLYKu6vYcZQ0gNjGNzVHxJV9RIUSZMLl/c0yqwiBtJU3UoyQaXrxvvxfIXFo7uX9z91ZQFBsJPkSBLI2M5dX5kegXpv1EGcF85ci8Ixlr+hYrGbk+L++9Y4QQFV2/FiGM7RLIi6YfAHjP3p+zZCaXHH3bNbK0thyT4ENclkM3GLd4L9mnG0+19yPW8Ke2eoqh2q+5Pje3zaKEEBWXQzfYeOgMCyNi+HjlAbzWvYmvkspuvS6zHV1d5cYv2SfzxsoxmfMhLmtzVPzFzKeXSMWDd2yDmGqZwnDTQn7RbyLaqOo6X9XHKinYhRAuSyNjs2RSvlHZx4/WP9ENhddtj6Bnux+WeWPll/R8iMvKb+jkV70tGxxN8FBsjDV9m+XcqeR0NkfF49AlQ5AQFd3SyFiGzdruCjw0HLxpngHAXEdnIowGOZ4j88bKLwk+xGXlP3Si8Ib9IWyGRjdtG93UrVnODvhyEx0nruL3XSdcXa0bD52RgESICiS3oduHtGVcqx7nrOHNe/b+eT5X5o2VTzLsIi6rdag/wX4exCWm5Zj3AXDQqMmXjjsYblrEWPO3rE8P4zwXA5bYxDSGz96R5TmSiEyIiiP70G01zjLK9DMAE+33uyaZ5kbmjZVP0vMhLktTFcb2bgJcms0jqyn2OzmuV6WmcppnTPMve01JRCZExZG99+I18yx8lPNE6PWZ5+ic63MUMm9SZN5Y+STBhyiQHmHBTH+gJUF+We9CqvlY8fUwERzoz1j7EAAe1X6nkXI83+s5e1DGLd4rQzBClHOX9l7crO6kj7YRh6Hwmu1hjFy+hpw3OWN7N5HJpuWUDLuIAusRFky3JkFsjornZHKaaxM5u66zZHcso+alsNRxIz20Lbxt/pr+Ga/n+sHidGkisnb1A0ruBxFClCjn0O3ZxETeNn0NwEzHbUQa9XItHyTDsuWeBB+iUDRVyREoaKpGdV9PAMbZHuQmdRet1f3co63jxzy6VC8lE8qEKN+cQ7eH5o6mtnqKWMOfDy5kMnV69taG1A2sJDtjVxAy7CKKhHNcNpYAPrLfA8Crptn4k3TZ58qEMiHKt13RCfyxdi1Pmn8DYKxtCClk3rCoCky4M4xnbm1E3+YhtKsfIIFHBSA9H6JIaKrCwx3qMuOvI3zjuI27tT+5Vj3Ga+ZZPGcbnutzFDK7V2VCmRDl2/xtx+n/3wdoqgN7o5481Hokd5xLp6q3lea1K+Nlka+iikZ6PkSRGdu7Ka/2bIwdEy/bHkM3FO7S1nOzujNHWZlQJkT5Fn02ld3RiUTGJKJFfEcr9V9S8eDQDWPx8TBzQ50qtG8QKIFHBSW/dVGk2jcIBGCX0YBvHLfxiGkp75i+ZoD5I6JTLsa6MqFMiPKt48TVAFQlgZXW70CB92338vWMw8BhAI68e4cbayjcSXo+RJEK8LZQ1dtKs5p+VOoZzkm1KrXUUyy57k/e7hdGw2revN0vjPUv3SKBhxDl2OT+zTGpCmPNM/G7sHHcTEd3AEyqwuT+zd1bQeFW0vMhilSwnyfrX+6CRVNRFAWj6jSYfS8+EV+S2qg1B056cvDkORlqEaKc69cihBYp66mz8m/shsrLtidwoAGwYEQHwkL83FxD4U7S8yGKnNWkoSiZwUVM1Y4k1O8Lhs7N+9/ChJ3FO08QGZPI7uhEos+murm2Qohicf4sNf4aA8Dnjl7sMeqiyD2HuEB6PkSx6jhxNQH0YKX1DxorR3hc+53PUvrQa8p6VxkZ9xWi/NgVncCE3/fxqffX+J8/xVFqsKr6w7zTuj7zthwnNiGNAG+Lu6sp3EyCD1GsJvdvzgs/7uRt2wN8YPmMUaafWa7fwCEjBJOq8P6917u7ikKIIvTL9hiUI2vxt8wDIPjBr/gptD2KojCwdW0yHDpWk+bmWgp3k+BDFKt+LUJoUM2bXlN0+jg20EnbxSTzF9ybMZYFI26WcV8hyoHos6mcTbGhGwZLth/kR9OXAJy+9kFiLWFUSThPzSpeKIoigYcAJPgQJUbhFfvjLFNHc4N6gIe1pcDN7q6UEKIIOJfVArxmmkNt0ylijAC67+hEyo7MIVYZXhWXkgmnotg5l98GhtRjb9iLALxo/oHkE/sZ8MUmdkUnuLeCQoir8nCHugDcoOznEW0pAGNsj7pSqDvPC+EkwYcods7ltwtHdKDNPc9h1OuMBxnUWPsCmw6f4pftMe6uohDiCjl0g6WRcXiQzvvmz1AVgx/tN7NGb+4qszQyDoduuK+SotSR4EOUCOfy2+iE8+y/8R0cJi/qnNvJg9oKWXorRBm2OSqe2MQ0RpvmEar+R6zhz1v2wVnKxCamsTkq3k01FKWRzPkQJco5NvyA1p+3zTN4yTSXNanN6TUlw1VGxoaFKBscusEv26NprfzDI6bM4ZaXbY+TRKUcZU8mp5V09UQpJj0fokQ5Uy5/7+jKRkcTvJR03jN/hoouKZeFKEOWRsbSceIqft12kPfMnwMwx96FtXruy+cP/HeOjYfOyPCLACT4ECWsX4sQFozogIHKi/ahnDM8aK3u5zHtNxaM6EC/FiHurqIQ4jK+33SEJ2dtJzYxjZdMc6mjniTGCOAd+6A8nzN19UEGfLmJjhNXsTQytgRrK0ojCT6E28QYVXnzwtjw86Yfscbvd3ONhBCX49ANxizYA0A7dQ8PmZYDMNr2BOfwuuzz4xLTGDZruwQgFZwEH6LE7Y5OQFXAAH5wdGalowVWxU7g8qfYffSkLL8VohRzThz1JYUPzNMBmGXvyl96swI93znoMm7xXhmCqcAk+BAlamlkLK/Oj+TiZ47CK7bHiTe8qZK8jyM/j2Xj4TOy/FaIUso5cfRN8wxqKPFE6dUZn89wS24MZAVMRSfBhygxDt1g3OK9ZL/XOUVlxtgeBeD2xDm0UA7I8lshShmHbrDx0BkO/JdMb3UD/bQN2A2VZ20jSMXjiq4pK2AqLllqK0qMMx9AbpbobZjv6MCd2l98aJ5Gr5QJsvOtEKXE0shYxi3eS2xiGsGcYan1awA+dfQjwmiQo7x/JTPxKbbLXreaz5UFLaLsk54PUWIud5cz1jaEE4Y/oep/vG76FkCW3wrhJs6ejjcX73GtbFHQec/8GX5KKhF6PabY++V4ngK83TeMYD8PlDyurQDBfh60DvUvxp9AlGbS8yFKzOXucpLw5tmMEcyxvM39pjWs0Ztzc99HURTYeOgMrUP90dS8Ps6EEEXl0p6OSz2sLaOjtofzhoVnbSOwZ/sKCfbzYGzvJvQIC0ZVFYbN2o4CWYZanf+Dx/ZuIv+fKzAJPkSJaR3qT7CfB3GJaTnmfTj9bVzLdEdvRpgW8a75S3rOr08sAUDWDzYhRPFYGhnLsFnbc/wfvVY5ykumuQC8Yx9ElHHx/+HILg3o0CAwyw1Cj7Bgpj/QMkcQEyT/jwUSfIgSpKkKY3s3yfVuyKlOgBdxoaPYuSuS69XDfGieziDbq+iorvwA0x9oKR9cQhSDvCaFe5LGFPMUrIqNPxwtmOW4Ncv5htW9aVc/IMf1eoQF061JEJuj4jmZnEY1Hw/pwRSAzPkQJcx5NxTkl3UIJtjPgykDmvPHc51Yuv8sz9hGkGJYaaft5QntVyAzWDGA1xfukfwAQhSDvCaFv276jgbqCf4zKvOibShkm82R35Cqpiq0qx9A3+YhtKsfIIGHAKTnQ7hBfndDGw+d4VRyOqcIJtw+hPfMX/C86Uc26E3ZZdQH4FRyOpuj4nO90xJCXLncJoXfrm5ioGk1uqHwrG048fhmOS8TR8WVkJ4P4Ra53Q05dIO/Dp52lfnR0YnfHK0xKw6mmKfgw8V8H0siY2WTKiGKWPYejBBO8a75KwCmO3qzQQ/L8RyZOCquhAQfolRw7pA5dfXBS45mZj89rleljnqSCeavcM4U+XbjUdmkSogi1jrUn6o+VgA0HHxs+RRfJZUdegM+st+TpWyQr5XPZP6VuEISfAi3c86uz22sOYlKPGV7Cpuh0UvbxEBtVZbzskmVEEVHUxVOJacDMMr0M63Uf0kyPHnKNjLLsto5j7flr5e7SuAhrpgEH8Kt8ppdf6kIowET7fcDMNb0LdcqR13nZJMqIYrW5P7N6azu5CnTAgDG2B4l2qgGgKpknpeJo+JqSfAh3Cq/lOuX+p+jJ384WmBVbEw1f4IXF5/j3KTqm7+iJAAR4ir1qwdf+XwBZO5Wu1hv7zq3aGRH+rUIcVfVRDkiwYdwq4JuLGWg8oJtKLGGP/XVWN42f032TCFv/faPzAER4grtik7ggc/XkzL7QUzpZ4nU6/KWYzAAinRyiCImwYdwq8JsLHUWX57KGIndULlLW59j/gfIHBAhroRDN/h09SE6Hp9Opf+2olt8eM38IteEVOWdO8NoFuJHVW8rAd4Wd1dVlBOS50O41eVSritAdV8roBCXlMZWozGT7P151TyHsaaZROp1Xfk/ILMvRCFzDki3JkEyLi1EPqLPpvLbrli+/PMw16du5ElLZkK/l+1D6daxDX2a16CWfyUGtq5NhkPHatLcXGNRXkjPh3ArZ8p1yJ4z8eLj8D5NCe/TxHX8C0cvljlaYVXsTLN8TGWSszzPOQdkc1R88VVciHKg48TVTFiyD8+U43xg/gyAr+09+CG1Je8t/5ebJq0BQFEUCTxEkZLgQ7hdXinXg/w8XPu49AgLZsKdYWR2ZCi8YHuSKL06NZXTfGz+FBU9x3ULOp9EiIrIoRtU9jTjQTqfmydTWUlhh96ACfaBrjKVPc0yiVsUCxl2EaVCQTagGtCmDn1bhDB701He/n0fw2zPMt/yBp20XTxt/MLkbEmQCjOfRIiKZnNUPAnnM/jI/BVN1KOcMnwZlvEMtku+FhLO22QrA1EspOdDlBoF2YDKy2Li4Y71ANhn1OYV22MAPK3Np7O6I0tZ2W9CiNztik7gjYWRPKQt407tL+yGysiMZ4gjZ5AhPYiiOEjwIcocTVV4uENdABboHfnOfiuqYvCJ+VNClcxVLg93qCuTTYXIwy/bY6hyagtjTN8DMN4+iL+Na3MtKz2IojjIsIsok8b2bkqwrwfjl+zjTfuDNFaPcaP6L1+aP+DPTnN4uGtTd1dRiFIl+mwqZ1NsKApsitjNd5ZPMCsOFjja87WjR47yCpnzrqQHURQHCT5EmdW+QSAAdkwMy3iWRdYxNFBPUO3om+w6/hUTlvzLK7c35rqald1bUSFKgY4TVwNgJYN5lklUVRP5R3cOXWbtJXQ+kh1rRXGRYRdRZgV4W6jqbaVZTT+evbMD71d+gzQs+B77g4Tfwtl4+Ay/bI9xdzWFKBUm92+OSYWJ5i9orh7irOHNUNuznCfnsMqlK82EKA7S8yHKrGA/T9a/3AWLpqIoCjc37E9shE7ouue4OW4mt6t+LN5p4Z4bamIYUKWSmZpVvNxdbSGKjUM3cl0xtis6gXlbjvNZ3XXcemIDNkNjuO0ZjhnVAXi5xzVcX6tKnivNhChqhe75WLduHb1796ZGjRooisKCBQuynDcMg/DwcGrUqIGnpyedO3dmz549RVVfIbKwmjSUCxtP3DRpNV2WB/Gl/XYAPjB/Ro3UffSasp7eU9e7up2FKI+WRsbSceIqBny5iWfmRjDgy02uvY5+2R5DpSPLuSXmcwDG2R9ko35xXtS7S/eTeD4j35VmQhSlQgcfKSkpXH/99UydOjXX85MmTeLDDz9k6tSpbNmyhaCgILp160ZycnKu5YUoKpndygrv2gew1nEdnkoGX1o+oDrxmFSFyf2bu7uKQhSLpZGxDJu1PccO0bGJaTw5azvbt/zFZPOnqIrBt/ZuzHJ0y3GNcYv3SkIxUWIKHXz07NmTt99+m7vuuivHOcMwmDx5MmPGjOGuu+4iLCyMmTNnkpqayuzZs4ukwkLkpV+LEBaM6IADjZG2p/lXDyFIOcv/LO+z6Inm1KtaiQFfbGJXdIK7qypEkXHoBuMW7811bySAABKZqkzCW0njL0dT3rQPzrWcbEkgSlKRzvmIiooiLi6O7t27u45ZrVY6derEhg0bGDp0aI7npKenk56e7nqclJQEgM1mw2azFWX1XNcr6uuWR2W1rex2OwDnFC8esb3IAssbhKlHSFz1NB9UeY2Nh8/w09bjNKrqxbajZzl9Lp1Abys31KlyxV3NZbWt3EHaqnAK0l6bo+KJP3cea7atVxw6aEYGX1k+oLZ6iiN6dYbbnsGOCTAwKaBlu/08mZiCzeZbxD9FyZD3VsEVV1sV5nqKYRhX3M+mKArz58+nX79+AGzYsIEOHToQExNDjRo1XOWeeOIJjh49yrJly3JcIzw8nHHjxuU4Pnv2bLy8ZHKgKJyEdHh/t0ZlC7SrrnP2xCGmGu9gVWzM0O9gXMYgvE0GT17rAKCSGfytbq60EEXo2DlYdFSlT207PeM+pcG5LSQYlbgrYxyHjczP5Rea2anl7eaKinInNTWVgQMHkpiYiK9v/kFssax2cU4AdDIMI8cxp1deeYXnnnvO9TgpKYlatWrRvXv3y1a+sGw2GytWrKBbt26YzeYivXZ5U5bb6u7eOhZNQVEUGr6+nBfVoXximcrD6m/8qwUxx34L7++++Na3aoYrr8FH/Ztz67XVC/V6ZbmtSpq0VeHk114r//mPd5fsIy4p6zwPuw4OQ8Fj3480MG0h3TDxRMazFwIPA1D4eI/GpR19ClDd14Nlo24us5NN5b1VcMXVVs6Ri4Io0uAjKCgIgLi4OIKDL64PP3nyJNWr5/6BbrVasVpz3nqazeZiewMV57XLm7LYVqp2cbnhiC4N+HytQqgtlmfNP/O26Wv+M6qwSm/pKp/uyPywVYA3f9tP97CQK/oALott5S7SVoWTvb2WRsYyfPbOC/M8cr5XH9BW8LjpNwBetA1ls9HkwpnMsjb94nOc/3rljqZ4WC1FX/kSJu+tgivqtirMtYo0yVhoaChBQUGsWLHCdSwjI4O1a9fSvn37onwpIXKVfbnhp6sPUtnLzMeOu/jB3glNMfjU/AnNlYM5nmsgk+5E6Xe5CaZd1B2MM30DwHu2+1ikd8j3epJQTLhDoXs+zp07x8GDFz+4o6KiiIiIwN/fn9q1azNq1CjGjx9Pw4YNadiwIePHj8fLy4uBAwcWacWFyM653DD7h/LpcxmAwqv2R6mmJNBZ28n/LO9xd0Y4R4ycH7iyi6cozTZHxedYUuvUQjnANPPHaIrBPHtnPnX0zfM6lT3NfDqoJW3rSV4PUfIK3fOxdetWWrRoQYsWLQB47rnnaNGiBW+88QYAo0ePZtSoUQwfPpxWrVoRExPD8uXL8fHxKdqaC3GJy90NAhiqieG2Z9ilhxKgJDPTPJFAEnOUk108RWmWV3BcX4nha8t7eCoZrHFczxj7I+Q2JKNc+PPu3c3o0CBQAg/hFoUOPjp37oxhGDn+fPPNN0DmZNPw8HBiY2NJS0tj7dq1hIWFFXW9hcgiv7tBJ4cOqXjwSMZojurVqKOe5GvLJCpxPks52cVTlGa5BcfBnOE7ywSqKOfYoTdgmO0ZHBc6tgO9s87jkGEWURrI3i6iXCjoUEmnhoGsPQBDbC/xsyWc69QovjR/wMO20aRj4eEOdeVOUJRqrUP9qepj5VRyZn6kyiTzreVdaijxHNRr8HDGi2hWb8ICKxGXmMb8Ee05Hn9e9m0RpYoEH6JcKOhQyZOdG9ChQSDjl8BDGS8x2/IO7bW9TGUKUV0+xcvLi24frmVI+7oMaF1bPqRFqaOpiivw8CSNry3v0VCNIdbw58GMl0nAB9LtLBrZgQyHjtWkyYaKotQp0tUuQrhL61B/gv08chnhzqQAwX6Zd33tGwQCEGnU47GMF0gzzHTTtlHzz9G8vmAXB06e47UFka5NuYQobSb3b46HksEX5g9pqR4kwajEgxkvc4JAVCXzvKIoWE3a5S8mhBtI8CHKBU1VGNs7M5dB9gDE+Xhs7yZoqkKAt4Wq3laa1fSjT797ed3yIjZD43ZjHeGmmXBh2qpzU67vNx0pqR9DiHztik5gwBebqO9vZnPD77hJi+Sc4cHDGaM5YNQEYNHIjvRrEeLmmgqRvzI77OJwOAqdl95ms2EymUhLS8PhcBRTzcqH4mors9mMphXP3ViPsGCmP9CScYv3Zpl8GuTnwet3XIufp4WFETFU8/Fg7ejOeJo1dAPGzA8jXR3GZPOnDDGtIAVPJtn74wxbxizYw/2t68gQjHC7X7bHsPnwSdT5E/FNWEWaYeYx24vsMBqiKHDlm2UIUbLKXPBhGAZxcXEkJCRc0XODgoI4fvx4nuneRabibKvKlSsTFBRULL+DHmHBdGsS5MpwWs3Hg7MpGbz1W9aAJNjPg7G9m+DnmbkSYJHeHh97Ku+Yv2a4aRE2ND6y3+sqvzkqnnb1A3DoRpZry8oYUdzi0yEyJgmz2cSvEdFMMn9B04Q/0VUzzxkvklCtNe+0q8O8LceJTUgjwLvsZykV5V+ZCz6cgUe1atXw8vIq1BeYruucO3cOb29vVFVGnPJTHG1lGAapqamcPHkSIEsK/qKkqQrt6gcAmYnHRszOmXgsLjGNYbO280iHuq5j3ztuxYKNsebveMY0H91Q+dhxNwBbj8aTeD4jR69KsJ8Hb9xxTbH8HEIAjNtugu2bUNB52zSDu01/YjdUhqeNZLkeBueTGdSmDgNb13ZNMBWitCtTwYfD4XAFHgEBAYV+vq7rZGRk4OHhIcHHZRRXW3l6egKZ+/1Uq1at2IZgIP/EY5nba8H8iJgsx2c4eqKi87r5e541/4wDlamOO/njn5N8uPzfHNeKTUxj1LwIJrXO7B1p26CaDM+IIuHsZRvcwMHsQwpvaf9joGk1uqHwnG04y/UbURX48L7mADLBVJQpZSr4cM7x8PKSZWNlmfP3Z7PZijX4uFziMQOIT7HhX8lCfEqG6/j/HHdgQucV8xxeMP+Ijsr0433zzZ4K8MjMLfh7ezK2dxNJ4CSuytLIWMYt3svp5POEVjKYYPof92lrcBgKz9uGsUjP3CtLN8DDLDdSouwpk+9ama9RtpXU76+gicf6Na+RY4XM547eTLL1B2C0eR4jtV/gsuHHxeEcWaIrCsOhG2w8dIaFETF8vPIAw2ZtzwycDZ3Hz3/lCjyetQ1ngd4xy3PHLd6LQ5eZpqJsKVM9H0IURkETj3VrEkTrUP8c8zmmOfoCBqPNP/C8+Se8lHQm2u8nt/0ynJzDOeMW76VbkyAZghGX5ezlyN5Lp6Lzrulz7tbWYzdURtlG8KveLsfznTsxO+c5CVEWSPAhyi1n4rG4xLRc+ywUMpfhOtNNd2sSxKZDZxgxezsJ5zOH+KY5+pGGlTfM3zHMtBhP0hlnfxAjn05DA/lCEBfltkLKGZTmtROzBRsfmqfRS/sbu6HytG0kv+tt83wN2YlZlDVlctilLHrooYdQFAVFUTCbzVSvXp1u3brx9ddfo+t6ga/zzTffULly5eKraDlSmMRjzvKqqrgCD6evHT151fYouqHwkGk575q+QuXi70w3MpdDZidfCGJpZCwdJ65iwJebeGZuBAO+3OTKnJvXhGhP0vjK/D69tL/JMDRG2J7ON/AA2YlZlD0VNvi4dIx146EzJTJm2qNHD2JjYzly5AhLliyhS5cuPPPMM/Tq1Qu73V7sr18ROROPBfll/XDOa2fPvAKG2Y6uPG97Eoeh0N+0hk/MU7GQGaTYdCVzOWQ28oVQsTl7NbIPpzjnBU1ddSDHOV/O8Z3lXW7WdpNqWHnENppleus8X+PSbQOEKEsq5LDL0sg43vrtn1yTThXnKgWr1UpQUBAAISEhtGzZkrZt29K1a1e++eYbHnvsMT788ENmzJjB4cOH8ff3p3fv3kyaNAlvb2/WrFnDww8/DFyctDl27FjCw8OZNWsWkydPZv/+/VSqVIlbbrmFyZMnU61atWL7ecqK3BKP5bWzZ34Bw3z9Js7brHxinkIvbRNVSGao7VnO4Ul1z4u9H5cO54iKqSDLvGf8dSTL8UAS+c4ygWvVYyQaXjycMZrtRsNsz7oot947IcqKCtfz8cf+M4yYvSPPu5GSXqVwyy23cP311/PLL78AoKoqn3zyCZGRkcycOZNVq1YxevRoANq3b8/kyZPx9fUlNjaW2NhYXnjhBQAyMjJ466232LlzJwsWLCAqKoqHHnqoRH+W0syZeKxv8xDa1Q/I88P6chvULdVb84htNOcMDzpoe5hreZtqJPDfeQXdkC+EiiavHtSCLPO+dHivjhLHT5ZwrlWPccrwo3/GG2w3GgEK99Vz5Pp+zKv3ToiyoEL1fDh0g0krD+d7N+KOVQqNGzdm165dAIwaNcp1PDQ0lLfeeothw4Yxbdo0LBYLfn5+KIri6kFxeuSRR1z/rlevHp988gmtW7d2ZSkVBeOcJzJs1nYUcl9c+5fejPszXmOGZRJh6hF+soTzoO1ljhhBBHpbeKxjKBl2w7WPzA11qrDt6NnL9rqIsiW3VSrOHtR0e8HmcSkKtOBfvrR8QICSzHG9Kg/YXuGokfn/W8OgQ3WDBUegSiULA9vUoW5gJXkfiTKvQgUfW47E819yRp7n3bVKwTAM1zDK6tWrGT9+PHv37iUpKQm73U5aWhopKSlUqlQpz2vs2LGD8PBwIiIiiI+Pd01iPXbsGE2aNCmRn6O8yGuDumo+VtJsDuoGVmJXdD3uzhjHd+YJ1FFP8pNlHI9lvEDEuQa8u3R/luupSuakVKeSGOITxSuvVSrOHtRRtzbM9XnZdVe28LF5Kh6KjV16KI9mvMgpKrvOP9O1PqT+y4yHbpTsuaJcqVDDLieTc1mSkGu5kl2l8M8//xAaGsrRo0e5/fbbCQsL4+eff2bbtm18+umnAPnu4JuSkkL37t3x9vZm1qxZbNmyhfnz5wOZwzGi8HqEBbP+pVuY83hbPr6/OXMeb8vGV7qy5bVbWTiiA5P7N+eEEsTdGePYrdclUElinuUteqkbc1wr+1xmSURWtl1uPgfAnM3HCPLNe/gO4GHTUqabJ+Oh2FjpaMH9Ga9nCTxe7dmYYZ0bAEgvhyh3KlTwUc3HWsByJbdKYdWqVezevZu7776brVu3Yrfb+eCDD2jbti2NGjXixIkTWcpbLJYcW9zv27eP06dP8+6773LTTTfRuHFj1+Zt4srlNk/EatJQFIV+LUJYMKIDp8kcn1/paIFVsTHVMoWnLpMN1XlGMlOWTQWZzxGXlM6A1rWBnMu8TdgJN33DWNO3qIrBLHtXhtqeI5XMz516gZUI9LbSu3mNYvoJhHC/ChV83FjXn+o+ljzvRop72Vp6ejpxcXHExMSwfft2xo8fT9++fenVqxcPPvgg9evXx263M2XKFA4fPsx3333HZ599luUadevW5dy5c/zxxx+cPn2a1NRUateujcVicT1v0aJFvPXWW8XyM4iczisePGF7ni/tPQF43vwTH5qnu5bi5ubSIT5RthS0Z7RuoFeOZd5+nGOmeSIPmZYDMNF2P6/ZH8GBhqbA5P7N+eP5Tvz1cheC/TyLpf5ClAYVKvjQVIXRt9YDCpZ0qqgtXbqU4OBg6tatS48ePVi9ejWffPIJCxcuRNM0mjdvzocffsjEiRMJCwvj+++/Z8KECVmu0b59e5588kn69+9P1apVmTRpElWrVuWbb77hxx9/pEmTJrz77ru8//77xfIziIsCvC1U9bYSVsOXe+oZvGMfzCu2R7EbKndp65ljeZtqnM33GpKIrOwpaM/o6eR0ujUJYv1Lt3Bb0yAaKNGs8A6ng7aHFMPKExnPMt3RB+enz8KRHenXIkR2pxUVgmIYRqnq901KSsLPz4/ExER8fX2znEtLSyMqKorQ0FA8PAo/NKLrOklJSWw4luqWPB9libOtfH19UdWijVGv9vdYmqTbHSi6gyVLlvDcJg2HodBB3c0088f4KamcNCozPONpthqNc33+yC4N6NAgsMKM6dtsNn7//Xduv/12zGazu6tzRRy6QceJq/JM23+pQG8LT9xUj71rf+Atx8f4KOc5awnm/uRn2G/UzlL21Z6NeaJT/SzHykN7lRRpq4IrrrbK7/s7uwq12sWpR1gQt4UFFyjplBD5sZo0bLbMlUU1K3tw9Gw6f+nN6JPxNp+bP6Kxepw5lnd40z6Y7xzdyN7nNnX1QaauPkiwnwev33EtVSpZ5T1ZyhVkObZT/Lk00la8yWTTAlBgk34tw5Ke4Sw5P5jHL9lH7QAvuQESFUKFDD7g4mRCIYrKyz0bM2z2TgCOGkHclTGOieYv6K1t4i3zNzRXDzHG9ghp5Jz4HJuYxvDZO7Icy603Lr9NykTJyWs59qUCSORj81Q6ansA+Nbejbfsg7Hl8bEruyGLiqRCzfkQojjdem11JtwZhvN7IxUPnrI9xVu2QdgNlbu1P1loeZ2GSnSBrpd9SW5+m5SJkudcjv36HdfmONdK2cdv1lfpqO0h1bDydMYI3rA/nGfgATIJWVQsEnwIUYQGtKlD5LjbmP1YGwa3rY2mKPzPcQeDba9wyvDjGjWaRZbXuF9bRf4d9lmX5P6+K/9NyiQAcQ9NVUi1XVz6rqLzpLaIuZa3CVLOclCvQZ+Mt1ikdyjwNWUSsqgIJPgQooh5WUy0bxDIW/2asXBkRwA26k3pmf4u6xzN8FQyeNf8FVPNU/AhNd9rOe+GX1sYmW9SK8kZUvzy2sclMiYJgGDO8L15PC+b52JSdBY52tEn420OGjUBuD0sKM9rX0p2QxYVQYWd8yFESVIUOG34McT2Ek/ov/GC6Qd6aZu4XjnEC7Yn+dvI2XV/qfiU0rctQEWS2z4uzpUsW4/Ec7u6iQnmr/BTUkkxrIyzP8gPjs44JxgH+VqZfH8Ltk9aRVxS7pmWZTdkUZFIz4cQxciZC6RZiB9v9Q1DVVQ+d/Tmvow3OK5XpZZ6innWt3jd9B1Wri4VvnTXFw/nPi7Zh7xOn8tg6pJtvJLxCdMsn+CnpBKh1+OOjPH84OjCpSubwvs0xWJSCe/TFAX35BkSojSR4EOIYhTs58n6l7uwcEQHBrerw9SBLQDYYTSkZ8YEZtu7APCoaQm/W16huXLwil9LuuuLXn77uHRVt7HcOpp7tHU4DIUp9n7ckxHOEePi6iRVgQl3hrlWLDlXyVya9RQyezymP9BSltmKCkOGXYQoZpdmq+zZrAav9jzP+CX7OIcXr9ofZ7l+IxPNX1BfjeVny1i+ctzOZPvdnKfgwURlTzO6YeDQDblzvgJ5LWHObR+XKiQRbv6WvtoGAA7rQYy2PeFKJPdOvzC8PUxU9bbSvHZlvCxZP2Z7hAXTrUmQLJkWFZoEH8Llm2++YdSoUSQkJLi7KuVa7+Y1+OLPKHw9TLSqW4XV+6x0P9eAceZv6KdtYKjpN+7Q/uYN20Os0lsW6JoJ520M+upvydR7BXKbz+FsxwMnz11S0qCPuoE3zN8RqCThMBS+dPTiI/vdpGNBUcAw4PpalQkL8cv3NSXPkKjoZNilhB0/fpxHH32UGjVqYLFYqFOnDs888wxnzpxxlencuTOjRo3K8xqrV6+mS5cu+Pv74+XlRcOGDRkyZAh2u/2yr79mzRoURcnx57XXXqN///78+++/rrLvvvsuLVsW7MtPFFywnyd/vdyFP57vxKR7rmfzmK7MGN6NUbaRPJzxItFGIDWV03xteZ9p5slUp+B5H2Tpbd6rUnKT13wOZzv+8c9/ADRSjjPX8jafWD4lUElin16LOzPe5F37ANKx8GjHUJqF+FHV20qAt6VYfz4hyoMK3fOxKzqBCb/v45XbG3NdzcrF/nqHDx+mXbt2NGrUiDlz5hAaGsqePXt48cUXWbJkCZs2bcLfP/+Z7nv27KFnz548/fTTTJkyBU9PTw4cOMBPP/2ErusFrsv+/fuz5N739vbG09MTT0/ZSbMkXDoUoygKFi3z8Rq9Bd3Sr+UZ0y88pv3O7dpmblZ3Mc3el/85epJO5hebqkBu36kGFTtTZn69GNl7g/Kbz+E8dig6ltdMP/OQtgyTonPesDDV3o8vHL1cCcOCfK28evu1qApkOHTZFE6IAqjQPR+/bI9h4+Ez/LI9pkReb8SIEVgsFpYvX06nTp2oXbs2PXv2ZOXKlcTExDBmzJjLXmPFihUEBwczadIkwsLCqF+/Pj169OCrr77CYin4HVe1atUICgpy/fH29uabb76hcuXKQOYQzMSJE9m5c6erd+Sbb765wp9cXI5zVUxoYCXO48G79oH0zniH7XoDvJU0RpvnsdLyIneomwAj18DD6UoyZRamt6C0ulwvxqW9QQ7d4Ju/ovJMja7h4H5tFX9YXuAx0xJMis5Sx43cmv4enzr6ZclUGt6nKZqqyG60QhRChev5OJGYxtFkA01VWbzzBACLd57gnhtqYhhQpZKZmlW8ivx14+PjWbZsGe+8806O3oWgoCAGDRrEvHnzmDZtWr7XCQoKIjY2lnXr1nHzzTcXeT2d+vfvz44dO1i9ejUrV64EwM8v/3FsceWcq2I0ReGmSauJTUzjH6MOd2eE00fdwMvmudRST/Gp5RMe0hsx3jaIHUbDfK9Z0KW3hektKK0u14txaW/Qir1x+ezJYtBN3cZo0zwaqpk3JVF6dcLtD7FWvz5LSVXJnFxaVtpIiNKkwgUft0/f5vq3s0M6PiWDXlPWu44fefeOIn/dAwcOYBgG116bezKpa6+9lrNnz3Lq1Kl8r3PvvfeybNkyOnXqRFBQEG3btqVr1648+OCDl93C+FI1a9bM8vjo0aNZHnt6elKpUiVMJhNBQQXLzCiujvOueWzvJjw5azsABioL9Y4sT2/FE9pvDDX9yo3qv8y3jmWVozkf2u8h0qiX6/UKsvTW2VuQ/Uvb2VtQVpZ/5rYq5VLO3qCpqw4yeeW/uQYprZR9vGSey41q5ryneMObKfY7+d5xKxlkbjs+4MaatK0fmOdKFiFEwVS4YZd3ejfEdGEc3PkB5PzbpCpM7t/cHdXCMDJroSj5j9FrmsaMGTOIjo5m0qRJ1KhRg3feeYemTZsSG1vwSYZ//vknERERrj9VqlS5qvqLotMjLDjLBnUA5/HgY8fddEn/gHn2ztgNlVu0CH61vsbn5g9prBzLco3KXmY+/uNfdkUn5Pk6BZnzUBRp2y8d0imuTdMK2ssz46+obD+vQTt1D3PMb/OT9U1uVP8lzTDzqb0PndInM8PR0xV4AMzZEo3VpNK+QaAEHkJchQr3v+eOptVoVrsqfT7dkOPcghEdLrtE7ko1aNAARVHYu3cv/fr1y3F+3759VKlShcDAwAJdLyQkhMGDBzN48GDefvttGjVqxGeffca4ceMK9PzQ0FDX/A5R+gxoU4e+LUKIOJbAqXPpVPW2YjGp3PPZRl6yP8F0R2+eNs2nn/oXt2lbuU3byh+OFnxh78XfRmMSUm1sOhzPL9tjcp1Mfbk5D1A0aduzD+lYNYNJrWHlP//R87qal3l2wRU0wVrCeduFfxl0UnfxlGk+rS70dGQYGj85buYT+13EkffPW1En8wpRlCpc8HEp57p859/FKSAggG7dujFt2jSeffbZLPM+4uLi+P7773nwwQcv2/ORmypVqhAcHExKSkpRVhmz2YzD4bh8QVEsnBvUOUXGJAKZ79cjRjDP2YYzTenDKNPP3K5upqu2g67aDnbq9Zip9GVB+g2u+Uz/xiUz6++jhPdpyomE8/nMecgpv16FvJJzQd5DOgDPzotAUbUiG9JpHepPsJ8HcYlp+e4V7EE6d2rreUhbxjVqNADphpm5js58bu/NCS4f/Ms+OkJcvQoZfDhXFgRX9qD/jbWYt+U4sQlpxb4+f+rUqbRv357bbruNt99+O8tS25CQEN555x1X2VOnThEREZHl+UFBQSxcuJCIiAjuvPNO6tevT1paGt9++y179uxhypQpRVrf2rVrExUVRUREBDVr1sTHxwer1VqkryEKLvv7dsZfRzh4siYjbc8QqsTymPY792jruF49zId8xPPWAOaldeHRKXH8R+YS7o9W/Mua/afy/YLO7sB/59h46EyOLJz5TVTt1iQozyEdp6LsQdBUhbG9mzBs1nYUyPG6IZxisGkl92urqKxkBukphpU5jlv43N6LUxRu2FH20RHi6lTI4MO5ssCiqSiKwsDWtUtkfX7Dhg3ZunUr4eHh9O/fnzNnzhAUFES/fv0YO3Zslhwfs2fPZvbs2VmeP3bsWPr27cv69et58sknOXHiBN7e3jRt2pQFCxbQqVOnIq1vnz59WLp0KV26dCEhIYEZM2bw0EMPFelriILL/r69LsSP3lP/QlEgyghmjP1RPrTfw4Om5QzWVhCinOE58088bfqFVXoLZjtuYf2/12FQuPf51NUHmbr6YJYVMJebqDrq1oYFGtL5aMW/dGgQeFXpxS/N1zP9gZaugMiDdHqoW7hXW0sHbY+r/DG9KjMdt/GDozPJZF3Z5mVRSc24fL4c2UdHiKtTIYMPyJnkqaTW59epU4cZM2bkW2bNmjX5nv/uu++u+PU7d+7smtya3UMPPZQluLBarfz444+oaoWbl1xqXfo+DfSxZukJ+WLdYY6egY/s9zLN3pee6mYGmFbRRt1Hd20b3bVtnDZ8+d3RhsWOdmw1GmEUYs65M7D4dGAL3vrtn3wnqn7x5+ECXTO3wCY/uSUGdObr+XbDUYa0DWFWl/NELP0f3Y2N+CjnXc9d72jKTMdt/KG3RM/j537/nut589e9su29EMWswgYfQpR1l/aE6AZM+eOA61w6FhboHVmQ0ZEGSjQDtVX01f4iUEniQdMKHjSt4IThzzLHjazSW/C3fm2WVR25cebLeG1hJPEptnzLpqQXbq5QXkt7swcbWQKN9gqKAssijtBV3Ubb3Z9Ra882Kisp1Cezssf0qvzk6MQv+k1EG1XzfH1nzo7br6uBqioMcy11vki2vRei6EjwUc707NmTP//8M9dzr776Kq+++moJ10gUJ2dPyObDZ/K8Wz9o1ORN+4O8Yx9EBzWSPtpGuqtbqKHE87BpGQ+zjBTDygY9jDX69WzSr+WQUYOLX7cXGXDZwCM/ea3adR5+es4O5g31oEXtKjh0g09XH2Lj4TNM+H0fL/dszOKdJ1DQ+WfHn/y6czod1EhWq/vwsFys02nDl+WOG1jo6Mhm45p8e3dGdmlA+/oBWXJ2OLe9zz6fJaiMJV4TojST4KOc+eqrrzh//nyu5y63b4wouwoyAdKBxjr9etbp12PlEW5Wd9FV3U4XLYLqSgLdtG100zKT8J0xfNimN2KLfg0RegP2GbVzzI+4EpdLGZLhMJi84gDtGwTw5Z+HOX0ug8okYzmykz8++5JJ6iFaWA/gr5zL8jxnL84SR2u2GtfkOaySXcPq3llWFDnJtvdCFC8JPsqZkJAQd1dBuEFhJ0CmY2GF3ooVeiuwGzRVjtJZjeBmbRfXK4cIUJJd80Scoo1A/tFrs9+oxVGjOtFGVWKMQGKNgCx7neTHYcCWU+DINqfTkzSqKQnUUM5Q6/AJzFGxfKCcoL71BDWV0zmuk2x4skm/lr/0MNbrYRw0Qsitp+Zy8ms32fZeiOIjwYcQ5UBB81zkTmGPUZc9jrp86uiHGTvNlMO0Uvdzo/ovTdQjhChnqKmcpqZ2mm5sz/Js3VA4gy+JRiUSqUSC4U0ilbAZJhwoONBwoKJi4KWk43k0ja6mdLyUdAJJpKqSgK+Se2+d0yE9mJ1GfXbp9dip12e3EYo9n48vk6rgZdFISrPn8RPLxFEh3EmCDyHKgcvlucjLyC4NCPbzYMyCSFeyPTsmthuN2O5oxBcX5o36co7GynEaq8dopERnBiLKKWoqp/BQbFQlkapK4lX9DKmGlTijClFGMIeMGhw2gjmk12C/UZMkvAt1Lbtu8GjHeny08t8c52TiqBDuJ8GHEOVEXhMlc+O883+2WyNOJqdRdWXWpHuHTp4jJePiipUkvNlsXMtmR/aNEQ0CSaKqkoCfkoIf5y78nYIZOxo6mqKjomOgcN6wkoqV81hJNTyIx4eTRmVOGlVIxpMrGTrJS91ALz6TiaNClEoSfAhRjlw6UXLF3ji+/utIjp6Q7Hf+uSXdOxafyt3TN+LrYaJV3SpsPXKW+JQMEs7bsm1HoHAGP04bfijG5XtcSmIrA6dqPh60qx8gE0eFKIUk+BCinHFOlGxXP4DWof4FuvPPnnSvTkAl/rokIDEMg2PxqdwzfaOrh+S7jUfZ/18yjYN8eKBtHeZtOc7hU5k9JtkDjJDKngzvUp+5m49x+L9ErBYLlb0stKpbhU2HznA84TyNg3wY2LoOb/66B5vjyiOU7PM5ZOKoEKWPBB+liKIozJ8/P9ddbwGOHDlCaGgoO3bsoHnz5iVat9zUrVuXUaNGMWrUKHdXReThapaM5haQZO8hOZdux9tqyrJNgUM3suzG26SGL36eZhRF4d4WwSz+bQm39ehEJQ+LK7BJPG9j74kkTp1LZ0TnBkz+40Ch5q646nnhb5nPIUTpJsFHCXnooYeYOXMmAJqmUaNGDe644w7Gjx9PlSqZm1rFxsa6/l1a3HLLLaxduzbHcZvNxpYtW6hUqZLr2OWCJ+EeRXnnnz0g8fEwZ3nsPJ9b7gxnGZMKVpPq2sF52Z64HL0zlb0yr5uQWriEZjKfQ4iyQYKPEtSjRw9mzJiB3W5n7969PPLIIyQkJDBnzhwgc9fa0ujxxx/nzTffzHLMZDJRtWre6aqFKIi8NqhLTLVhAM90bcDMDUdJOJ93EOJfyczrvZoS5CvzOYQoK2THsBJktVoJCgqiZs2adO/enf79+7N8+XLXeUVRWLBggevx5s2badGiBR4eHrRq1YodO3bkuOaiRYto2LAhnp6edOnShZkzZ6IoCgkJCa4yGzZs4Oabb8bT05NatWrx9NNPk5KSUuB6e3l5ERQUlOUPZA67TJ482fVvgDvvvBNFUVyPhciLQzcYt3hvnhvUKcAPW6MZf2czFHKug3EeG39nM+5sEUK7+gESeAhRRpT94MMwICOl4H9sqYUrn9+fq5i2f/jwYZYuXYrZnPtmXikpKfTq1YtrrrmGbdu2ER4ezgsvvJClzJEjR7jnnnvo168fERERDB06lDFjxmQps3v3bm677Tbuuusudu3axbx581i/fj0jR4684rrnZsuWLQDMmDGD2NhY12NRvjh0g42HzrAwIoaNh87guFy+9HxsjorPd0mwAcQmplGlkoXpD7QkyC9rNtIgP48cG9EJIcqGYht2mTZtGu+99x6xsbE0bdqUyZMnc9NNNxX9C9lSYXyNAhVVgcpF+dqvngBLpcuXu+DXX3/F29sbh8NBWlrmh+6HH36Ya9nvv/8eh8PB119/jZeXF02bNiU6Opphw4a5ynz22Wdcc801vPfeewBcc801REZG8s4777jKvPfeewwcONA1KbRhw4Z88skndOrUienTp+Phcfm03NOmTeOrr75yPR46dCgffPBBljLOIZjKlSuX2uEjcXWWRsbmmJsRfBVzLAqyH42zXN/mIbJkVohypFiCj3nz5jFq1CimTZtGhw4d+Pzzz+nZsyd79+6ldu3axfGSZUKXLl2YPn06qampfPXVV/z777889dRTuZb9559/uP766/HyuriZV7t27bKU2b9/PzfeeGOWY61bt87yeNu2bRw8eJDvv//edcwwDHRdJyoqimuvzZ40KqdBgwZl6VGpXLnyZZ8jype85mbEJaYxbNb2K+qBKOh+NM5ysmRWiPKjWIKPDz/8kEcffZTHHnsMgMmTJ7Ns2TKmT5/OhAkTivbFzF6ZPRAFoOs6ScnJ+Pr4oKpFMOJkLtwun5UqVaJBgwYAfPLJJ3Tp0oVx48bx1ltv5ShrFGBIxzAM14qBvJ6n6zpDhw7l6aefzvH8ggaCfn5+rnqLiqcgczPGLd5LtyZBheqJuNx+NLL/ihDlV5EHHxkZGWzbto2XX345y/Hu3buzYcOGHOXT09NJT093PU5KSgIyl3LabFlnuNtsNtddu65fsi2mybNAdTMMA8wODLMXulIE3bWGUeB5H4ZhuOru9Prrr3PHHXcwdOhQatTIHDpy/myNGzfmu+++IyUlBU/PzJ/P2X7OMtdccw1LlizJck3nXAtnmRYtWrBnzx7q1auXa72ytGO2+uZV7+zlnOfMZjM2my3Pspe+pmEY2Gw2NE3Lt2xZ4HyfZn+/lhebo+KJP3ceaz6/qvhz59l08ORlA4VL28oMvHHHNTw7LwLIPQvrG3dcg+6wozuokMr7e6soSVsVXHG1VWGuV+TBx+nTp3E4HFSvXj3L8erVqxMXF5ej/IQJExg3blyO48uXL88y5ACZyzuDgoI4d+4cGRkZV1zH5OTkK37ulbLZbNjtdldwBdCyZUsaN27MuHHjXPM2zp8/T1JSEr169eK1115jyJAhvPDCCxw7doz3338fyJyMmpSUxMCBA/noo4949tlnGTx4MLt372bGjBmun1FVVYYPH0737t154oknGDJkCF5eXuzfv581a9YwadKky9bb4XCQkZGRpd5Ouq6TlpbmOle7dm2WLl3Kddddh9VqzXN4JiMjg/Pnz7Nu3Trs9tx3HS2LVqxY4e4qFJtJrS9f5vQ/m/j9n4Jd79K2mpjPtTOitvF7VMGuWZ6V5/dWUZO2KriibqvU1NQCly22Cae5DQdkPwbwyiuv8Nxzz7keJyUlUatWLbp3746vr2+WsmlpaRw/fhxvb+8CTZTMzjAMkpOT8fHxybUuxclsNmMymXL8TM8//zyPPvoor732GgCenp74+vri6+vLokWLGD58OJ06daJJkyZMnDiRe++9l0qVKuHr60uzZs344YcfePHFF/n8889p164dY8aMYcSIEVStWhUPDw/at2/P6tWree2117j99tsxDIP69etz33335ajLpZxtpWkaFosl17KqquLh4eE698EHH/DCCy/w7bffEhISwuHDh3O9dlpaGp6entx8881X9HssbWw2GytWrKBbt255rl4qyzZHxfPIzMuvXvp6yI0F6vnIra0cusG2o2c5fS6dQG8rN9SpIpNJKf/vraIkbVVwxdVWud2k5qXIg4/AwEA0TcvRy3Hy5MkcvSGQmfvCarXmOG42m3M0isPhQFEUVFW9ojkbzuEA5zVKkjO7aXYPPPAADzzwAJBzvkb79u2JiIjIcix7mX79+mXJKPrOO+9Qs2bNLL1Gbdq0KXSE62yr1atX59lWR44cyfK4b9++9O3b97LXVtXM7Ja5/Y7LsvL28zi1bVANf2/Py87NaNugWoEDhuxtZQY6NMr5+SAyldf3VnGQtiq4om6rwlyryL+BLRYLN9xwQ44vuxUrVtC+ffuifrkKb9q0aWzZsoXDhw/z3Xff8d577zFkyBB3V0uUI5qqMLZ3EyD3RF8ge6kIIQqnWIZdnnvuOQYPHkyrVq1o164dX3zxBceOHePJJ58sjper0A4cOMDbb79NfHw8tWvX5vnnn+eVV14p0HP//PNPevbsmef56OjooqqmKON6hAUz/YGWBdohVwghLqdYgo/+/ftz5swZ3nzzTWJjYwkLC+P333+nTp06xfFyFdpHH33ERx99dEXPbdWqVY5hHafLrVgRFc/V7JArhBCXKrYJp8OHD2f48OHFdXlRBDw9PfPM36HreqEmD4mKQRJ9CSGKQtnf20UIIYQQZUqZDD5kSKBsk9+fEEJUbMU27FIcLBYLqqpy4sQJqlatisViKVS+Dl3XycjIIC0trcSX2pY1xdFWhmGQkZHBqVOnUFUVi8VSJNcVQghRtpSp4ENVVUJDQ4mNjeXEiYLt53IpwzA4f/48np6eJZ5krKwpzrby8vKidu3aEgAKIUQFVaaCD8js/ahduzZ2ux2Ho3AbPthsNtatW8fNN98sSWguo7jaStM0TCaTBH9CCFGBlbngA7ji7JiapmG32/Hw8JDg4zKkrYQQQhQX6fcWQgghRImS4EMIIYQQJUqCDyGEEEKUqFI358O5a2txZNe02WykpqaSlJQk8xguQ9qq4KStCk7aqnCkvQpO2qrgiqutnN/b2Xdfz02pCz6Sk5MBqFWrlptrIoQQQojCSk5Oxs/PL98yilGQEKUE6brOiRMn8PHxKfLlmElJSdSqVYvjx4/j6+tbpNcub6StCk7aquCkrQpH2qvgpK0KrrjayjAMkpOTqVGjxmXzOJW6ng9VValZs2axvoavr6+8OQtI2qrgpK0KTtqqcKS9Ck7aquCKo60u1+PhJBNOhRBCCFGiJPgQQgghRImqUMGH1Wpl7NixWK1Wd1el1JO2Kjhpq4KTtiocaa+Ck7YquNLQVqVuwqkQQgghyrcK1fMhhBBCCPeT4EMIIYQQJUqCDyGEEEKUKAk+hBBCCFGiKmzw0adPH2rXro2HhwfBwcEMHjyYEydOuLtapc6RI0d49NFHCQ0NxdPTk/r16zN27FgyMjLcXbVS65133qF9+/Z4eXlRuXJld1enVJk2bRqhoaF4eHhwww038Oeff7q7SqXSunXr6N27NzVq1EBRFBYsWODuKpVKEyZM4MYbb8THx4dq1arRr18/9u/f7+5qlVrTp0/nuuuucyUXa9euHUuWLHFLXSps8NGlSxd++OEH9u/fz88//8yhQ4e455573F2tUmffvn3ous7nn3/Onj17+Oijj/jss8949dVX3V21UisjI4N7772XYcOGubsqpcq8efMYNWoUY8aMYceOHdx000307NmTY8eOubtqpU5KSgrXX389U6dOdXdVSrW1a9cyYsQINm3axIoVK7Db7XTv3p2UlBR3V61UqlmzJu+++y5bt25l69at3HLLLfTt25c9e/aUfGUMYRiGYSxcuNBQFMXIyMhwd1VKvUmTJhmhoaHurkapN2PGDMPPz8/d1Sg1WrdubTz55JNZjjVu3Nh4+eWX3VSjsgEw5s+f7+5qlAknT540AGPt2rXurkqZUaVKFeOrr74q8detsD0fl4qPj+f777+nffv2shVzASQmJuLv7+/uaogyJCMjg23bttG9e/csx7t3786GDRvcVCtR3iQmJgLI51MBOBwO5s6dS0pKCu3atSvx16/QwcdLL71EpUqVCAgI4NixYyxcuNDdVSr1Dh06xJQpU3jyySfdXRVRhpw+fRqHw0H16tWzHK9evTpxcXFuqpUoTwzD4LnnnqNjx46EhYW5uzql1u7du/H29sZqtfLkk08yf/58mjRpUuL1KFfBR3h4OIqi5Ptn69atrvIvvvgiO3bsYPny5WiaxoMPPohRQRK+FratAE6cOEGPHj249957eeyxx9xUc/e4kvYSOSmKkuWxYRg5jglxJUaOHMmuXbuYM2eOu6tSql1zzTVERESwadMmhg0bxpAhQ9i7d2+J18NU4q9YjEaOHMn999+fb5m6deu6/h0YGEhgYCCNGjXi2muvpVatWmzatMktXVAlrbBtdeLECbp06UK7du344osvirl2pU9h20tkFRgYiKZpOXo5Tp48maM3RIjCeuqpp1i0aBHr1q2jZs2a7q5OqWaxWGjQoAEArVq1YsuWLXz88cd8/vnnJVqPchV8OIOJK+Hs8UhPTy/KKpVahWmrmJgYunTpwg033MCMGTNQ1XLVYVYgV/PeEpkfeDfccAMrVqzgzjvvdB1fsWIFffv2dWPNRFlmGAZPPfUU8+fPZ82aNYSGhrq7SmWOYRhu+d4rV8FHQW3evJnNmzfTsWNHqlSpwuHDh3njjTeoX79+hej1KIwTJ07QuXNnateuzfvvv8+pU6dc54KCgtxYs9Lr2LFjxMfHc+zYMRwOBxEREQA0aNAAb29v91bOjZ577jkGDx5Mq1atXD1ox44dk/lDuTh37hwHDx50PY6KiiIiIgJ/f39q167txpqVLiNGjGD27NksXLgQHx8fV8+an58fnp6ebq5d6fPqq6/Ss2dPatWqRXJyMnPnzmXNmjUsXbq05CtT4utrSoFdu3YZXbp0Mfz9/Q2r1WrUrVvXePLJJ43o6Gh3V63UmTFjhgHk+kfkbsiQIbm21+rVq91dNbf79NNPjTp16hgWi8Vo2bKlLInMw+rVq3N9Dw0ZMsTdVStV8vpsmjFjhrurVio98sgjrv9/VatWNbp27WosX77cLXVRDKOCzLAUQgghRKlQ8QbvhRBCCOFWEnwIIYQQokRJ8CGEEEKIEiXBhxBCCCFKlAQfQgghhChREnwIIYQQokRJ8CGEEEKIEiXBhxBCCCFKlAQfQgghhChREnwIIYQQokRJ8CGEEEKIEiXBhxBCCCFK1P8BtqD5pxfNaEwAAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "execution_count": 1, + "id": "eef9ae09", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], "source": [ "%matplotlib inline\n", "\n", @@ -147,8 +131,10 @@ }, { "cell_type": "markdown", - "id": "50154095", - "metadata": {}, + "id": "c368bb8d", + "metadata": { + "editable": true + }, "source": [ "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).\n", "see also the weekly slides [for week 36](https://compphysics.github.io/MachineLearning/doc/pub/week36/html/._week36-bs029.html).\n", @@ -165,8 +151,10 @@ }, { "cell_type": "markdown", - "id": "25dbedec", - "metadata": {}, + "id": "9ba5c048", + "metadata": { + "editable": true + }, "source": [ "$$\n", "C(\\beta_0, \\beta_1, ... , \\beta_{p-1}) = \\frac{1}{n}\\sum_{i=0}^{n} \\left(y_i - \\beta_0 - \\sum_{j=1}^{p-1} X_{ij}\\beta_j\\right)^2,.\n", @@ -175,8 +163,10 @@ }, { "cell_type": "markdown", - "id": "2d0c55cc", - "metadata": {}, + "id": "be61f447", + "metadata": { + "editable": true + }, "source": [ "Recall also that we use the squared value. This expression can lead to an\n", "increased penalty for higher differences between predicted and\n", @@ -190,8 +180,10 @@ }, { "cell_type": "markdown", - "id": "69cee15d", - "metadata": {}, + "id": "2efa5b58", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\frac{\\partial C}{\\partial \\beta_j} = 0,\n", @@ -200,16 +192,20 @@ }, { "cell_type": "markdown", - "id": "a4049e25", - "metadata": {}, + "id": "cd3b54c0", + "metadata": { + "editable": true + }, "source": [ "for all $j$. For $\\beta_0$ we have" ] }, { "cell_type": "markdown", - "id": "f3d52473", - "metadata": {}, + "id": "1a141a26", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\frac{\\partial C}{\\partial \\beta_0} = -\\frac{2}{n}\\sum_{i=0}^{n-1} \\left(y_i - \\beta_0 - \\sum_{j=1}^{p-1} X_{ij} \\beta_j\\right).\n", @@ -218,16 +214,20 @@ }, { "cell_type": "markdown", - "id": "d69193fd", - "metadata": {}, + "id": "5c2a9db1", + "metadata": { + "editable": true + }, "source": [ "Multiplying away the constant $2/n$, we obtain" ] }, { "cell_type": "markdown", - "id": "87f92158", - "metadata": {}, + "id": "f279815d", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\sum_{i=0}^{n-1} \\beta_0 = \\sum_{i=0}^{n-1}y_i - \\sum_{i=0}^{n-1} \\sum_{j=1}^{p-1} X_{ij} \\beta_j.\n", @@ -236,8 +236,10 @@ }, { "cell_type": "markdown", - "id": "d5ad854b", - "metadata": {}, + "id": "f77462b8", + "metadata": { + "editable": true + }, "source": [ "Let us specialize first to the case where we have only two parameters $\\beta_0$ and $\\beta_1$.\n", "Our result for $\\beta_0$ simplifies then to" @@ -245,8 +247,10 @@ }, { "cell_type": "markdown", - "id": "6cd6700d", - "metadata": {}, + "id": "8ef3d129", + "metadata": { + "editable": true + }, "source": [ "$$\n", "n\\beta_0 = \\sum_{i=0}^{n-1}y_i - \\sum_{i=0}^{n-1} X_{i1} \\beta_1.\n", @@ -255,16 +259,20 @@ }, { "cell_type": "markdown", - "id": "83d9895f", - "metadata": {}, + "id": "d7ca1b4a", + "metadata": { + "editable": true + }, "source": [ "We obtain then" ] }, { "cell_type": "markdown", - "id": "c551c1d4", - "metadata": {}, + "id": "6859a8d9", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\beta_0 = \\frac{1}{n}\\sum_{i=0}^{n-1}y_i - \\beta_1\\frac{1}{n}\\sum_{i=0}^{n-1} X_{i1}.\n", @@ -273,16 +281,20 @@ }, { "cell_type": "markdown", - "id": "2e70abb5", - "metadata": {}, + "id": "a018b414", + "metadata": { + "editable": true + }, "source": [ "If we define" ] }, { "cell_type": "markdown", - "id": "b1d65476", - "metadata": {}, + "id": "9f3fd8e8", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mu_{\\boldsymbol{x}_1}=\\frac{1}{n}\\sum_{i=0}^{n-1} X_{i1},\n", @@ -291,16 +303,20 @@ }, { "cell_type": "markdown", - "id": "c8d07c9e", - "metadata": {}, + "id": "df1038ff", + "metadata": { + "editable": true + }, "source": [ "and the mean value of the outputs as" ] }, { "cell_type": "markdown", - "id": "87a136ed", - "metadata": {}, + "id": "c979cecf", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mu_y=\\frac{1}{n}\\sum_{i=0}^{n-1}y_i,\n", @@ -309,16 +325,20 @@ }, { "cell_type": "markdown", - "id": "75a7a5e6", - "metadata": {}, + "id": "1380fd5d", + "metadata": { + "editable": true + }, "source": [ "we have" ] }, { "cell_type": "markdown", - "id": "95fb12ac", - "metadata": {}, + "id": "40e7054e", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\beta_0 = \\mu_y - \\beta_1\\mu_{\\boldsymbol{x}_1}.\n", @@ -327,16 +347,20 @@ }, { "cell_type": "markdown", - "id": "cf52ab75", - "metadata": {}, + "id": "f785b8c5", + "metadata": { + "editable": true + }, "source": [ "In the general case with more parameters than $\\beta_0$ and $\\beta_1$, we have" ] }, { "cell_type": "markdown", - "id": "e3b71b32", - "metadata": {}, + "id": "1b94b3e7", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\beta_0 = \\frac{1}{n}\\sum_{i=0}^{n-1}y_i - \\frac{1}{n}\\sum_{i=0}^{n-1}\\sum_{j=1}^{p-1} X_{ij}\\beta_j.\n", @@ -345,16 +369,20 @@ }, { "cell_type": "markdown", - "id": "ee47da8f", - "metadata": {}, + "id": "1346bbb5", + "metadata": { + "editable": true + }, "source": [ "We can rewrite the latter equation as" ] }, { "cell_type": "markdown", - "id": "b35c78dd", - "metadata": {}, + "id": "044dc319", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\beta_0 = \\frac{1}{n}\\sum_{i=0}^{n-1}y_i - \\sum_{j=1}^{p-1} \\mu_{\\boldsymbol{x}_j}\\beta_j,\n", @@ -363,16 +391,20 @@ }, { "cell_type": "markdown", - "id": "296179e8", - "metadata": {}, + "id": "1ae565fb", + "metadata": { + "editable": true + }, "source": [ "where we have defined" ] }, { "cell_type": "markdown", - "id": "4035eda9", - "metadata": {}, + "id": "2fe619e9", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mu_{\\boldsymbol{x}_j}=\\frac{1}{n}\\sum_{i=0}^{n-1} X_{ij},\n", @@ -381,8 +413,10 @@ }, { "cell_type": "markdown", - "id": "099d1c83", - "metadata": {}, + "id": "025a9267", + "metadata": { + "editable": true + }, "source": [ "the mean value for all elements of the column vector $\\boldsymbol{x}_j$.\n", "\n", @@ -391,8 +425,10 @@ }, { "cell_type": "markdown", - "id": "27ce7366", - "metadata": {}, + "id": "22131558", + "metadata": { + "editable": true + }, "source": [ "$$\n", "C(\\boldsymbol{\\beta}) = (\\boldsymbol{\\tilde{y}} - \\tilde{X}\\boldsymbol{\\beta})^T(\\boldsymbol{\\tilde{y}} - \\tilde{X}\\boldsymbol{\\beta}).\n", @@ -401,16 +437,20 @@ }, { "cell_type": "markdown", - "id": "4464c2ce", - "metadata": {}, + "id": "ab3ab42d", + "metadata": { + "editable": true + }, "source": [ "If we minimize with respect to $\\boldsymbol{\\beta}$ we have then" ] }, { "cell_type": "markdown", - "id": "ea48afd1", - "metadata": {}, + "id": "ca424c8e", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\hat{\\boldsymbol{\\beta}} = (\\tilde{X}^T\\tilde{X})^{-1}\\tilde{X}^T\\boldsymbol{\\tilde{y}},\n", @@ -419,8 +459,10 @@ }, { "cell_type": "markdown", - "id": "c4ed0a05", - "metadata": {}, + "id": "2367454f", + "metadata": { + "editable": true + }, "source": [ "where $\\boldsymbol{\\tilde{y}} = \\boldsymbol{y} - \\overline{\\boldsymbol{y}}$\n", "and $\\tilde{X}_{ij} = X_{ij} - \\frac{1}{n}\\sum_{k=0}^{n-1}X_{kj}$.\n", @@ -430,8 +472,10 @@ }, { "cell_type": "markdown", - "id": "011f1e31", - "metadata": {}, + "id": "2cf181e1", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\hat{\\boldsymbol{\\beta}} = (\\tilde{X}^T\\tilde{X} + \\lambda I)^{-1}\\tilde{X}^T\\boldsymbol{\\tilde{y}}.\n", @@ -440,43 +484,31 @@ }, { "cell_type": "markdown", - "id": "92267a4c", - "metadata": {}, + "id": "b20efe04", + "metadata": { + "editable": true + }, "source": [ "Now we try to implement this." ] }, { "cell_type": "code", - "execution_count": 40, - "id": "715fe165", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[0.47179152 5.01549939]\n", - "[1.79909592 0.47176716 5.01550546]\n" - ] - }, - { - "ename": "ValueError", - "evalue": "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)", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", - "Input \u001b[0;32mIn [40]\u001b[0m, in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 32\u001b[0m \u001b[38;5;66;03m# calculate intercepts and print them\u001b[39;00m\n\u001b[1;32m 33\u001b[0m interceptOLS \u001b[38;5;241m=\u001b[39m y_scaler \u001b[38;5;241m-\u001b[39m X_train_mean \u001b[38;5;241m@\u001b[39m beta_OLS\n\u001b[0;32m---> 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)" - ] - } - ], + "execution_count": 2, + "id": "5b3cfc82", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], "source": [ + "\n", "np.random.seed(2018)\n", "n = 100\n", + "# we do not include the intercept\n", "d = 2\n", "Lambda = 0.01\n", + "\n", "# Make data set.\n", "x = np.linspace(-3, 3, n)\n", "y = 2.0 + 0.5*x + 5.0*(x**2)+ np.random.randn(n)\n", @@ -486,8 +518,10 @@ "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 +535,7 @@ "\n", "#Calculate beta\n", "beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)\n", - "Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d)\n", + "beta_Ridge = 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,6 +547,8 @@ "#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", @@ -522,7 +558,6 @@ "print(\"test MSE of Ridge\")\n", "print(MSE(y_test,ytilde_test_Ridge))\n", "\n", - "\n", "plt.scatter(x,y,label='Data')\n", "plt.plot(x, X @ beta_OLS+interceptOLS,'*', label=\"OLS_Fit\")\n", "plt.plot(x, X @ beta_Ridge+interceptRidge, label=\"Ridge_Fit\")\n", @@ -533,8 +568,10 @@ }, { "cell_type": "markdown", - "id": "edf2436a", - "metadata": {}, + "id": "9f177d5f", + "metadata": { + "editable": true + }, "source": [ "Finally, instead of using our own function we repeat the same example\n", "using the **standardscaler** functionality of the library\n", @@ -543,35 +580,53 @@ }, { "cell_type": "code", - "execution_count": 41, - "id": "f1ec76ea", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - " \n", - "test MSE of OLS\n", - "1.1394311129039263\n", - " \n", - "test MSE of Ridge\n", - "96.31278139904835\n" - ] - }, - { - "data": { - "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": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "execution_count": 3, + "id": "9130dfd9", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], "source": [ - "X_test @ RegRidge.coef_+ RegRidge.intercept_\n", + "from sklearn import linear_model\n", + "np.random.seed(2018)\n", + "n = 100\n", + "d = 2\n", + "Lambda = 0.01\n", + "\n", + "# Make data set.\n", + "x = np.linspace(-3, 3, n)\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((n, d))\n", + "for p in range(d): \n", + " X[:, p] = x ** (p+1)\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 using scikit-learn\n", + "from sklearn.preprocessing import StandardScaler\n", + "scaler = StandardScaler()\n", + "scaler.fit(X_train)\n", + "X_train_scaled = scaler.transform(X_train)\n", + "X_test_scaled = scaler.transform(X_test)\n", + "\n", + "#Calculate beta\n", + "OLS = LinearRegression()\n", + "OLS.fit(X_train,y_train)\n", + "ypredictOLS = OLS.predict(X_test)\n", + "RegRidge = linear_model.Ridge(Lambda)\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", "ytilde_test_OLS = X_test @ OLS.coef_+ OLS.intercept_\n", "\n", "#Calculate MSE\n", @@ -587,35 +642,9 @@ "plt.legend()\n", "plt.show()" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "de73251c", - "metadata": {}, - "outputs": [], - "source": [] } ], - "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.10" - } - }, + "metadata": {}, "nbformat": 4, "nbformat_minor": 5 } diff --git a/doc/src/week37/programs/scaling.py b/doc/src/week37/programs/scaling.py deleted file mode 100644 index c1f0c364f..000000000 --- a/doc/src/week37/programs/scaling.py +++ /dev/null @@ -1,86 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -from sklearn.linear_model import LinearRegression -from sklearn.preprocessing import PolynomialFeatures -from sklearn.model_selection import train_test_split -from sklearn.preprocessing import StandardScaler - -def MSE(y_data,y_model): - n = np.size(y_model) - return np.sum((y_data-y_model)**2)/n - -def OLS_fit_beta(X, y): - return np.linalg.pinv(X.T @ X) @ X.T @ y - -def Ridge_fit_beta(X, y,L,d): - I = np.eye(d,d) - return np.linalg.pinv(X.T @ X + L*I) @ X.T @ y - - -np.random.seed(2018) -n = 1000 -d = 3 -L = 0.001 -true_beta = [2, 0.5, 3.7] - -# Make data set. -x = np.linspace(-3, 3, n) -y_real = 2 + 0.5*x + 3.7*x**2 - -y = np.sum( - np.asarray([x ** p * b for p, b in enumerate(true_beta)]), - axis=0) + 0.1 * np.random.normal(size=len(x)) - - -#Design matrix X including the intercept -X = np.zeros((len(x), d)) -for p in range(d): # (d-1) - X[:, p] = x ** (p+1) # (p+1 if not intercept included) - - -#Split datamatrix -X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) - -# Scale data by subtracting mean value,own function -#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,L,d) -print(beta_OLS) -print(beta_Ridge) - -interceptOLS = y_scaler - X_train_mean @ beta_OLS -interceptRidge = y_scaler - X_train_mean @ beta_Ridge -print(interceptOLS) -print(interceptRidge) -#predict value -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,y_real,label='no noise') -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()