diff --git a/doc/src/week37/programs/.ipynb_checkpoints/codeexamplesscaling-checkpoint.ipynb b/doc/src/week37/programs/.ipynb_checkpoints/codeexamplesscaling-checkpoint.ipynb new file mode 100644 index 000000000..265ef3fb7 --- /dev/null +++ b/doc/src/week37/programs/.ipynb_checkpoints/codeexamplesscaling-checkpoint.ipynb @@ -0,0 +1,278 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a85a3d6f", + "metadata": {}, + "source": [ + "\n", + "" + ] + }, + { + "cell_type": "markdown", + "id": "e7c10e6b", + "metadata": {}, + "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", + "\n", + "Date: **Sep 11, 2023**\n", + "\n", + "Copyright 1999-2023, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license" + ] + }, + { + "cell_type": "markdown", + "id": "e7fed16e", + "metadata": {}, + "source": [ + "## This note contains code examples with a simple scaling" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "6fdaaf05", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[1.97993409 0.49717915 3.70154994]\n", + "[1.97958122 0.49715476 3.70158419]\n", + " \n", + "test MSE of OLS:\n", + "0.011394311129039391\n", + " \n", + "test MSE of Ridge\n", + "0.011401577230117069\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB76klEQVR4nO3deVhUZfvA8e85Z4YBBNxQQUUlzRVNzQ21XErDzNSW19QW2zXtzZZfm5n6lprti6WtWqmpZW6VuORWprlFCmqZYm7ghgKCwMw55/cHzjhsCgoMDPfnurh0znnmzMPDLPc8y/0opmmaCCGEEEKUEtXTFRBCCCFExSLBhxBCCCFKlQQfQgghhChVEnwIIYQQolRJ8CGEEEKIUiXBhxBCCCFKlQQfQgghhChVEnwIIYQQolRZPF2B3AzD4OjRowQGBqIoiqerI4QQQohCME2T1NRUateujapevG+jzAUfR48eJSwszNPVEEIIIcRlOHToEHXr1r1omTIXfAQGBgLZlQ8KCirWa9vtdlasWEHv3r2xWq3Fem1vI21VeNJWhSdtVTTSXoUnbVV4JdVWKSkphIWFuT7HL6bMBR/OoZagoKASCT78/f0JCgqSJ+clSFsVnrRV4UlbFY20V+FJWxVeSbdVYaZMyIRTIYQQQpQqCT6EEEIIUaok+BBCCCFEqZLgQwghhBClSoIPIYQQQpQqCT6EEEIIUaok+BBCCCFEqZLgQwghhBClSoIPIYQQooLQDZPN8UkAbI5PQjdMj9RDgg8hhBCiAoiOTaDrlNU88OUWAB74cgtdp6wmOjah1OsiwYcQQgjh5aJjExgxazsJyRkYJkyNUzFMSEzOYMSs7aUegEjwIYQQQngx3TCZsHQXzgEWw4S9KdnBh/PYhKW7SnUIRoIPIYQQwottjk8iITnDdbu9sgt/MtDPxxomkJCc4ZoLUhok+BBCCCG82PHUC4FHEGf5yvoa22zDqUVSgeVKmgQfQgghhBerGejr+n8fbQs+is4BsxbHqF5guZImwYcQQgjhxTqEVyO0cnZgcav6GwBL9c45yoRW9qVDeLVSq5MEH0IIIYQX01SFcf2aU4PTRKq7AFhiRHJhuimM69ccTVVKrU4SfAghhBBeLioilOlt/kVVTLYbjegcXh0FUBWYPDCCqIjQUq2PpVQfTQghhBAecW3KagAqtxtEF9Pk1m7taBsejL9P6YcC0vMhhBBCeLukeDiyFRSVetcNAaDjVdU9EniABB9CCCGE94tdkP1vg+sgoJZn64IEH0IIIYT3cwYfLe/wbD3OkzkfQgghhBdy7mB77shOeh7fhan5oDTr5+lqARJ8CCGEEF4nOjaBCUt3kZCcwTOWefS0wHqzNef2ZXBDkwBPV0+GXYQQQghv4r6DLZj0UzcC8G1GR0bM2s6q3cc8W0Ek+BBCCCG8hm6YjF0c50of1lrZR331OGmmjVVGW0zg1R92e7KKgAQfQgghhNfYHJ/EidRM1+1btex06iuNa8nABsCRlCyP1M1dkYKPadOm0apVK4KCgggKCiIyMpJly5a5zg8bNgxFUXL8dOrUqdgrLYQQQoi83Hem1dDpp2UPuSxx28vFoph57lfaijThtG7durz22ms0atQIgC+//JL+/fvzxx9/0KJFCwCioqKYMWOG6z4+Pj7FWF0hhBBCFMR9Z9rOahw1lGSSzADWG61cx7UyMOZRpOCjX7+cS3QmTpzItGnT2LRpkyv4sNlshISEFF8NhRBCCFEozh1sE5IzGKBtAOBHvRMOt4/7kCBfIM1DNcx22UttdV3n22+/JS0tjcjISNfxtWvXUrNmTapUqUK3bt2YOHEiNWvWLPA6mZmZZGZeGJ9KSUkBwG63Y7fbL7d6+XJer7iv642krQpP2qrwpK2KRtqr8KStLni5bxOen/s7N6lbAFhqdEbBxARsKjx309XY/40psc/YwlBM0yzS4M/OnTuJjIwkIyODgIAA5syZw8033wzAvHnzCAgIoH79+sTHxzN27FgcDgfbtm3DZrPle73x48czYcKEPMfnzJmDv79/UaomhBBCCKD26U20P/ARaT7BrGr+FiYKugmWEhxySU9PZ8iQISQnJxMUFHTRskUOPrKysjh48CBnzpxhwYIFfPbZZ6xbt47mzZvnKZuQkED9+vWZO3cut912W77Xy6/nIywsjJMnT16y8kVlt9tZuXIlvXr1wmq1Fuu1vY20VeFJWxWetFXRSHsVnrRVTtr8oah7l3OoxXC2NxxJcICNa+tXRVOVEmurlJQUgoODCxV8FHnYxcfHxzXhtF27dmzZsoX33nuPjz/+OE/Z0NBQ6tevz969ewu8ns1my7dXxGq1ltgTqCSv7W2krQpP2qrwpK2KRtqr8KStgPQk2PczAGHd7iesZr18ixV3WxXlWlfcAWOaZo6eC3enTp3i0KFDhIaGXunDCCGEEKIw4haC4YCQllCzqadrk68i9Xy8+OKL9OnTh7CwMFJTU5k7dy5r164lOjqas2fPMn78eG6//XZCQ0M5cOAAL774IsHBwQwcOLCk6i+EEEIILmwk13jTLKoDRsR/ymwm0SIFH8eOHeOee+4hISGBypUr06pVK6Kjo+nVqxfnzp1j586dfPXVV5w5c4bQ0FB69OjBvHnzCAwMLKn6CyGEEBWecyM5NfkQG3y3Y5gKA9aH8FiVBKIiyt7oQ5GCj88//7zAc35+fixfvvyKKySEEEKIwnNuJGcCj51Pp77JaMbOlEqMmLWdaXe3LXMBSFntkSl2zu4oyM59rxueTy8rhBBCXImcG8mZ9D+fWGyR0QUTMIGxi+PK3GfeZScZK0+c3VEnU8/RIEDlwO9bCA70Y1y/5mUuGhRCCCEKy30juWbKQZqoh8k0LUTrHVxlTqRmsjk+iciG1T1VzTy8vufD2R2VkJyBYcLeFBXDhMTkDEbM2k50bIKnqyiEEEJcFveN5AZqvwKw2mhDCpUKLFcWeHXw4d4d5UcGt6nrGGOZhW5SprujhBBCiMJwbiSnYriGXBbqXQssV1Z4dfDh3h1VlbO8Yf2Ehy0/UYcTrjLO7ighhBCivHFuJNdZjaOWcobTZgBrjDY5yoRW9qVDeDUP1TB/Xh18uHczHSWYjXp2Cvj+52cD51dOCCGEKC80VWFcv+YM1H4B4Ae9E/Zc0znH9WuOpiqeqF6BvDr4yN3N9L2R3RV1m/YrYBZYTgghhCgvohoH0d9nO5BzyCUkyMb0MrjMFrx8tYuzOyohObtnY5nenlcsM2ikHqWlEs9O86oy2R0lhBBCFNruH7Do6ZhVw/m/vvdw/GwmNQOzP9vKWo+Hk1f3fDi7o5zS8GeP/7UA3Ha+i6osdkcJIYQQhbZjHgBKq0FENgqmf+s6RDasXqY/27w6+ACIighl+t1tCQmyYVXhXGgXAPpbNvLxkJZlsjtKCCGEKJTURNi/Jvv/rf7j2boUgVcPuzhFRYTSq3kIm/45zoldWWT5Vqdaxilusu0C8t9qWAghhCjLdhw+w5/z3+Ie04C6HaB6Q09XqdC8vufDSVMVOoRXw1Q0tJZ3ZB/cMdezlRJCCCEug26YfLhmH23PZO+pZrQa5OEaFU2F6PnIzWj5H7QtH2Pu+Yktu+NJyPQp85NzhBBCiMOn0/lxRwKf/rKfamn7aGH7lyxTo8/yavzn3D76tgqlblV/T1fzkipk8EFIK1KDGhGY8g/fzfqQ+XoPIDsRi+z3IoQQoqzqOmWN6/8PWrLTqa81WrMvzcbkZXuYvGwPB17r66nqFVqFGXZxt2rPcT461Q6A28+vegHZ70UIIUTZpRsmVfysQHY69QHn06l/r1/nKlPFz1outgypkMHHqz/sZpHeBcNU6Kjuoa6SnW5d9nsRQghRVm2OT+LMOTsAndU4QpUkzpiVWO2WTv3MOXu52DKkQgYfR1KySKA6m4xmAAxQf81xXvZ7EUIIUda4bwVyu7YegKV6JFlYCyxXVlXI4MOiZPdqfG9kd1VlJxzL2dNRHv54QgghKg7nViCVOEeUugWABfr1BZYryypk8KGd/62X6R1IN21cpSbSVtmbo0x5+OMJIYSoOJxbhtys/Y6fksU+I5QY80JuD4WyuYNtfipk8BESlB1YpOHHMqM9kHPiaXn54wkhhKg4nFuGOD+vsns9stNDOJNElJctQypk8PF8n6au/zu7rPppG7GRBZSfP54QQoiKJapOJp3U3RgoOXewrezLtDK6g21+KmSejxub1WL63W0ZvySOjSnNOWJWp45yijsr7aDrwEfKzR9PCCFEBfPn+U3kwq/n7a59OZ6aUS6TZFbI4AMu7PeyOT6Js7/fAX9/zP8a7OB3Px8Wxxwpl39MIYQQ3kk3TDbvP0XLLV8TAJjXDCayYXVPV+uyVdjgA7LHzyIbVofKj8DfH2PuW81/45ZxgqqAZDwVQgjhedGxCUxYuovaKX+ywHaIs6Yv/X4K4jlLQrn9fKqQcz5yi06sxFajMZpbxjiQjKdCCCE8Kzo2gRGztpOQnOHK7bFM78CBFMr151OFDz50w2TC0l0sOJ+e9g5tPc6cH87MHxOW7pKMp0IIIUqVbpiMXRyHCdjI4hZtEwALjOvLfUbuCh98bI5PIiE5gx/1TmSaVpqoh4lQ4l3nTSAhOUMyngohhChVm+OTOJGaCUBvdStByjkOm8H8blxYsVleM3JX+ODDmck0hUqsMK4Fcub8yF1OCCGEKA3unzt3ausAWKBfh5nro7s8fj5V+ODDPZPpd3o3APprG7DiKLCcEEIIUdKcnzuhnKKrGgvAd+U0nXpuFT746BBezbVF8S9GS46ZVaimnKWnut1VpoqfVTKeCiGEKFXOdOq3ab+gKiabjGYcMmvlKFNeM3JX+OBDUxVeu70lAAYq35+feOrs4gJ47faWku9DCCFEqdJUhXG3NOOO859H3zq65SlTXjNyV/jgA7ITjr14PuX6t+eHXrqrf1KD07zYp2m5XUcthBCifIsK+pdw9Rhp+PKT0cF1PCTIxvRylE49twqdZMxdv9a1+eSXePC9mgNGBA3SY7nHfxP9Wt/u6aoJIYSoqGJmAeDX+na+aNm93KZTz02Cj/NCK/ux4fke+GgqyvZHYOl/GVV1E2pQ+ZvII4QQwgtkpUHcIgDUNncTWb/8plPPTYZd3NgsGoqiQIuBYPFDPbUXDm/1dLWEEEJUILphsnHfKbZHz4Sss5jVroJ6kZ6uVrGS4CM/vkHQvH/2/893eQkhhBAlLTo2ga5TVjP4001kbf0agI+TOxEdl+jhmhUvCT4K0ubu7H93LoCsdM/WRQghhNdz38clTDlGJ3U3hqnwVVpkud7HJT9FCj6mTZtGq1atCAoKIigoiMjISJYtW+Y6b5om48ePp3bt2vj5+dG9e3fi4uKKvdKlon4XqFIfslI5uGEugz/ZxI7DZzxdKyGEEF7IfR8XcO4zBr8aERylernexyU/RQo+6taty2uvvcbWrVvZunUrPXv2pH///q4A4/XXX+ftt99m6tSpbNmyhZCQEHr16kVqamqJVL5EqSq0HgpA+uav2Lj/FB+t2ec1f3ghhBBlh/s+LiqGa5sPZ/oHKL/7uOSnSMFHv379uPnmm2ncuDGNGzdm4sSJBAQEsGnTJkzT5N1332XMmDHcdtttRERE8OWXX5Kens6cOXNKqv4l5vDpdGZndsZAoem5P6irnCA6LpGOk1bx8bp9HD4tQzFCCCGKh/v+LF3UWOoqJ0k2/VlhtCuwXHl22UttdV3n22+/JS0tjcjISOLj40lMTKR3796uMjabjW7duvHbb7/x6KOP5nudzMxMMjMzXbdTUlIAsNvt2O32y61evpzXK8x1u05ZA0A9awuu02K5Q1vHu47bOXk2i8nL9jB52R72vtL7Elcpv4rSVhWdtFXhSVsVjbRX4ZX3tgr2t2DTTAwTBmlrAVikdyETHwCsqomqZJe70t+xpNqqKNdTTNMs0jjCzp07iYyMJCMjg4CAAObMmcPNN9/Mb7/9RpcuXThy5Ai1a9d2lX/kkUf4999/Wb58eb7XGz9+PBMmTMhzfM6cOfj7+xelasVq6wmF2ftUblE28r7PVA6bwVyf+S4GKqpiMrShQbsaMgQjhBCi+JxITuWefU9gUxz0zZzELrM+JgrPtHQQFuDp2l1ceno6Q4YMITk5maCgoIuWLXLPR5MmTYiJieHMmTMsWLCA++67j3XrLuyDoig5M66ZppnnmLsXXniBp556ynU7JSWFsLAwevfufcnKF5XdbmflypX06tULq9VaYLnN8UnM37wFTYHlRjvOmJWoq5ykixrLL0YrNAXmx6tEXd++XG7oUxiFbSshbVUU0lZFI+1VeN7QVqt2H+PY/NewWRzEGg34i/rnz5h8EKfx7l2tubFZrYteozBKqq2cIxeFUeTgw8fHh0aNGgHQrl07tmzZwnvvvcdzzz0HQGJiIqGhF3LNHz9+nFq1Cm4sm82GzWbLc9xqtZbYE+hS1z6Z7iBTdwZMPizSuzDMsoJB2lp+MVphNxRXufL6JC+skvw7eBtpq8KTtioaaa/CK89t1adlHbr+/Dskwzy9Ow4z+7MmJMjG+FtbFPs+LsXdVkW51hWnVzdNk8zMTMLDwwkJCWHlypW0adMGgKysLNatW8eUKVOu9GFKVc3AnCnV5+vdGWZZQW91C1VJ4TRB+ZYTQgghLtvR7QQm/42p2ej3nydol2Xzin1c8lOk4OPFF1+kT58+hIWFkZqayty5c1m7di3R0dEoisLo0aOZNGkSV199NVdffTWTJk3C39+fIUOGlFT9S0SH8GqEVvYlMTkDE9hlNmCHEU4rNZ6B2gZm6H0IqezrtUMuQgghPGB7dkZTpfmtdGh+lYcrU7KKtNT22LFj3HPPPTRp0oQbbriB33//nejoaHr16gXAs88+y+jRo3nsscdo164dR44cYcWKFQQGBpZI5UuKpiqM69ccAGesOV/vDsAgbQ1gMq5fc6+LRIUQQpQ+3TD5ZsMe0rbNy77d+m4P16jkFSn4+Pzzzzlw4ACZmZkcP36cVatWuQIPyJ5sOn78eBISEsjIyGDdunVEREQUe6VLQ1REKNPubktI5eyhlSV6Z86ZPjRRDzOnj6XYx96EEEJUPM69XLb8OINKpHPQqMF18+xelUo9P1c858ObRUWE0qt5CJvjkziemsHZHX3x27+Qqw59z+KYll47FieEEKLkzd50gDGLsjOED/JZC2T3sh9NyWL4rO1MHNCCoZ0aeK6CJUiCj0vQVIXIhtUB+P30bdTYv5BKexfzQuxNpONLaGVfxvVrLj0hQgghCk03TFfg0UBJoKO6B91U+E6/3lVmzKI47upQ3yu/4MqutoUUHZvAXctV9hshBCgZ9NU2AZCYnOF1uw0KIYQoWe57tPxHy86Vtd5oRSLVCyznTST4KIQLuw0qfHt+4uldWnb6dfP8jzftNiiEEKJkOfdoseDgzvPBx1y9R4HlvI0EH4Xgvtvgd/p1OEyVa9W9NFIOu8p4026DQgghSpYzT9QN6nZqKMmcMCvzs9G2wHLeRoKPQnCPPE9QldVGdhI1Z+9HfuWEEEKIgjjzSQ0+/znyrd4NR65pmKFenE9Kgo9CyB15fqP3BOB27RdsZBVYTgghhMiPpipM7lmZ69UdQN4hFwW8Op+UBB+F4IxQndYZ13DUrEZV5Sw3qVsA745QhRBCFL/uactRFZOtSksOmhf2QAut7Mu0u9t69SpKWWpbCM6Mp8NnbQfAQGWeowdPWhcwxLKaJVldvDpCFUIIUTx0w8yeR5h8lpu2foUNaHvbk3zj34njqRkVJn+UBB+FFBURyvS72zJ+SRyJKZnM17vzX8v3dFJ381X/qlzvxRGqEEKIKxcdm8CEpbtISM6gp7qdW30SOUMgvxvtuKlh9UtfwItI8FEEuTOepmzpTtUja7g+9Segs6erJ4QQooyKjk1gxKztOBMyOCeafue4jonfxDJN8/HqYZbcZM5HETkznvZvXYeq1z0CgBkzh01/H2VxzBE27jsl+T6EEEK46IbJhKW7XIFHLZLoqWYP439zfqLphKW7KtRnh/R8XImre5PhWxPf9ON8PXMaPxqdACTluhBCCJfN8UkkJF9IxXCntg5NMfndaMo+sw4ACckZbI5Pcm3n4e2k5+MKRO8+wSdnuwAwWPvZdVxSrgshhHByzwGlYDBIWwvAXEePAst5Owk+LpMz5fo8R3cMU6GrFkd9JRGQlOtCCCEucM8Bdb26kzD1BMmmPz8ZHQss5+0k+LhMzpTrR6jBeqMVkDfjqaRcF0II4Z4raqi2CoAF+vVk4uMqU9FyRUnwcZncu8ecGU/v0NZhxVFgOSGEEBWPM1dUCKfoqf4BwGz9hhxlKlquKAk+LpN799jPRhuOmVWooaTQW91aYDkhhBAVU1REKJ9G7MKiGDkmmqoKTB4YUeEWKMhql8vk7EZLSM7AgYW5eg+esCxkqLYqx6qXitSNJoQQogC6g5bHlwBQ9bpHea9Ga2oE2Ghdrwr+PhXvo1h6Pi6TsxvNaa6jJ7qp0FnbRUPlCFDxutGEEEIUYO8KSDkC/tVp3GMo/VvXoXOj4AoZeIAEH1fEmXI9JMhGAtVZbbQF4CG/dUz38k2BhBBCXJxumGzcd4rFMUc4/cvH2QdbDwWLzbMVKwMqZshVjNxTrvPPw7BxG3f5/ILSpIqnqyaEEMJD3PdxqaucoJ/POlBgfVBfrvd05coA6fkoBs6U65G9/gNV6qFkJEPcQk9XSwghhAc493FxZjW9S1uNqpj8qkdw36JTkoASCT6Kl6rBtcMASP31E9nrRQghKhhnAkrnu74Fhyuj6Sz9RklAeZ4EH8VstV9vHGgEnvyDj+ctYvCnm+g6ZbVEukIIUQE4E1A69VK3UUNJ5rhZhVXn5wVKAkoJPopVdGwCD353kGi9PQBDzu/3Inu9CCFExZA7seTd5zOaztO743CbZlnRE1BK8FFM3LvaZuk3AjBA20AlzsleL0IIUUG4J5ZsqByhixaHbip84+hZYLmKSIKPYuLe1bbJaMY+I5QAJYOB2q+uMtLVJoQQ3s19Hxdnr8fPRluOEuwqIwkoJfgoNjm70BRX70f2k88soJwQQghv4kxA6UcGt2vrAfha75WjjCSglOCj2OTuQlugX0+6aaOpeoj2yl8FlhNCCOFdoiJC+bbzYYKUc8QbtfjViAAgJMgmCSjPkyRjxcR9rxeAFCqxSO/MEMsa7rWsYIu9qXS1CSFERWCaRBz9FgCl/YO8G9aWmoHZ7/8VvcfDSXo+iknuvV4AZp3vaotSt1CDM9LVJoQQFcHhLZC4Eyy+NLjhYfq3rkNkw+ry/u9Ggo9i5L7XC8AuswFbjcZYFZ05bfdIV5sQQngx3TCZtelfVn01EQAj4nbwl97u/MiwSzFz3+vleGoGlU+OgF+fpNGh79i4dxTH0xzS/SaEEF7GuZdLVvIxfrP9Cgo8GHcNgxolyBfPfEjwUQKce70A4BhK5pZXsKUeZeaMj1hudACyl1qN69dcnpRCCFHOzd50gDGL4gAYoa3DpjiIMa5iTWpd1szazsQBLRjaqYFnK1nGyLBLCYvek8RnadcBcK+20nVcsp4KIUT5pxumK/BQMRhqyc7tMcttee2YRZJgMjcJPkqQbphMWLqL2Y4b0E2FLlocDZUjwIXMHxOW7pInpRBClFPuiSN7qn9QVznJGbMSS/XIAsuJIgYfkydPpn379gQGBlKzZk0GDBjAX3/9laPMsGHDUBQlx0+nTp2KtdLlxeb4JBKSMzhKMD+f31DImfEOsgOQhOQMeVIKIUQ55Z448l5tBQBz9R5k4lNgOVHE4GPdunWMHDmSTZs2sXLlShwOB7179yYtLS1HuaioKBISElw/P/30U7FWurxwf7J9pfcG4HZtPZU4V2A5IYQQ5YczcWRD5QjXazvRTSXHkEvuciJbkSacRkdH57g9Y8YMatasybZt27j++utdx202GyEhIcVTw3LM/cm2wWjBPiOUhmoCA7Vfczw55UkphBDlkzPB5L1p2b0ePxttOWzWyFFGEkzmdUWrXZKTkwGoVi1no65du5aaNWtSpUoVunXrxsSJE6lZs2a+18jMzCQzM9N1OyUlBQC73Y7dbr+S6uXhvF5xX7cgbeoGUr+qjYTkDLIMlS/13vxP/ZJh2nJm6Tfio2Y/KdvUDSy1OhVWabdVeSZtVXjSVkUj7VV4nmyr//WuQ+clvwAwx+iNTTMxTVAARYGX+zbB0B0YeqlXLV8l1VZFuZ5imuZlzXY0TZP+/ftz+vRpfvnlF9fxefPmERAQQP369YmPj2fs2LE4HA62bduGzWbLc53x48czYcKEPMfnzJmDv7//5VStzDl0Ft7caSGAdDbaHidQOcfdWS/QqUUzwgI8XTshhBBXIvzEClodnkWqb21WN52cHXFUQOnp6QwZMoTk5GSCgoIuWvayg4+RI0fy448/8uuvv1K3bt0CyyUkJFC/fn3mzp3Lbbfdlud8fj0fYWFhnDx58pKVLyq73c7KlSvp1asXVqu1WK99Md9uO8RLi3djmibjLV9yn2UFq/S2HO/7OXdeG1Zq9SgKT7VVeSRtVXjSVkUj7VV4Hmsr08AyPRIlaR/23lPYHDyQk2czCQ6wcW39qmUymWRJtVVKSgrBwcGFCj4ua9jl8ccfZ8mSJaxfv/6igQdAaGgo9evXZ+/evfmet9ls+faIWK3WEnsCleS18zOk01UMaFuPmINnSE+sCqtWcIP2B0pjBcr4G0ppt1V5Jm1VeNJWRSPtVXil1Va6YfLN5oPErl/Ia+n7MH0CsV47lC62wBJ/7OJS3G1VlGsVabWLaZqMGjWK77//ntWrVxMeHn7J+5w6dYpDhw4RGlqxM3n6+1jo3CiYG7t2gYY3oGDCls88XS0hhBBFFB2bQNcpq3lpUSw3pi4CYL5+PdF7z3q2YuVIkYKPkSNHMmvWLObMmUNgYCCJiYkkJiZy7lz20tGzZ8/yzDPPsHHjRg4cOMDatWvp168fwcHBDBw4sER+gXKp43AAMrd8Sb+3ljNr07+SaEwIIcqB2ZsOMHzWdhKSMwhTjtFTjQFgenpPhs/azuxNBzxav/KiSMHHtGnTSE5Opnv37oSGhrp+5s2bB4CmaezcuZP+/fvTuHFj7rvvPho3bszGjRsJDCw/XVElLTqzBYeUUGyOVFomLeelRbF0nbJaUq0LIUQZ5p5KHeAebRWqYrJOb0W8md27L6nUC6dIcz4uNTfVz8+P5cuXX1GFvJ1zA6IHtBt52fo192nLmaP3JCE5g+GyAZEQQpRZ7tmo/cngLm0NADP1m/KUc20uKvIle7uUIveo+Vu9G2mmjSbqYSLVXa4yEjULIUTZ5J6N+nZtPUFKOvuNENYa1xRYTuRPgo9S5B41p+LPAj07K+z9WnSB5YQQQpQNzmzUCgbDtOxe/pn6TZi5Pkola/WlSfBRinJHw1+e3+/lRnU79ZRjBZYTQgjhec5U6t3UHTRUE0gx/VxfIp0klXrhSPBRinJHw/vMOqzRr0FVTFcUnV85IYQQnqUbJpvjk+gTEcID2jIA5uvdScPPVUYBxvVrXiYTi5U1EnyUImfU7O4LvQ8Ad2rrCCBdomYhhChjnHk9Bn+6iXW/bXDtXus+0TS0si/T7m5LVETFzmlVWFe0sZwoGk1VGNevOSNmbcc5pfQXoyV7jTpcrR7hP9o6OvQbI1GzEEKUEdGxCTnes51z9FYZ13LYrMmDXRpwY/MQOoRXk/fuIpCej1IWFRHKtLvbuvWAKHyhRwHwTOU1VLZpLI45wsZ9p2TVixBCeJBumExYussVeFTmLLdr2RupztCjUICfYhMl8LgM0vPhAVERofRqHsLm+CSOp2YQ4ncNWQu+wz/9MDNmfMQKoz2Q3Y03rl9z6cYTQggP2ByfRELyhQUAd2lr8FOy2GXUZ5PRDICE5AzJ63EZpOfDQzRVIbJhdfq3rsNpu8an6d0AuN9t4mlicgYjZm2XzKdCCOEB7isPNXTutawAYIZ+E9nTS/OWE4UjwYeHObv1vnb0wm5qRGq7aK4cAHB19U1YukuGYIQQopS5rzyMUrdQRznFSTOIJXrnAsuJwpHgw8Oc3XqJVOcnoyOQM+mYyYVuPSGEEKXnwgpFk4ctPwIwS7+RTHxcZWSF4uWR4MPD3LvrvnBkTzy9VfuNGpwpsJwQQoiS51yh2FbZS2t1H5mmlVmOXjnKSF6PyyPBh4e5d9f9aTZiq9EYm+LgnvNji/mVE0IIUTqiIkJ5t96vACzUu3CSygCoCkweGCELAi6TrHbxMGe3nnNG9WeOm2nn8zd3a6v4yNGfDGzSrSeEEKXImc30eGoGYRynzfHVADQZ8BzvWepTI8BG63pV8PeRj9DLJS3nYc5uveGztgOwwmjHv0ZN6qvHuUNbzyy9l3TrCSFEKYmOTWDC0l2uL4TjLF/S1mJwotZ1tGnXmTYerp+3kGGXMiAqIpTJAyNQFTBQXSnXH9SWMXmA5PkQQojS4Mxm6gw8gkjjP9paAJ461FXSHhQj6fkoIwZ3rE//NnWIOXiGpDNX41ixkPCsROpV3sXGfUEcT82gZqCvZNITQogSkDubKcBd2moqKZnsMcL41Yjgn6W76NU8RN6Di4EEH2WIv4+Fzo2CgWA48yD8+g475r/K4HMvucpI1lMhhCh+ubOZWnAwzJKd9PFzvQ8mimQzLUYy7FJGrak8ALup0cbcRUtlv+u4ZD0VQojilzudwc3qZmorSZwwK7NY71JgOXF5JPgog3TD5MVVp1hiRAK4ktuAZD0VQoiSkDOdgcnDlh8A+MrRiyysBZQTl0uCjzLI2f33maMvADerv1Obk67zkvVUCCGKV4fwalTxyw4yItVdtFQPcM70YZZ+o6tMFT+rpD0oJhJ8lEHObr3dZn1+1VtgUQzut0QXWE4IIcSV0VSF125vCcCjWnavx3y9G6cJcpV57faWMtm0mEjwUQa5d+t9qt8CwGBtNUGkFVhOCCHElYmKCOWN6zS6a3+imwqf6Te7zr3Yp6lM9C9GEnyUQRc2M4J1Rit2G2EEKBkM1X52lZGsp0IIUfz6nv0OgF+snYm89lquCq5EcICNfq1re7hm3kWW2pZBObOeKnziuIV3fKZxvyWaL/QoMvGRrKdCCFFMnOnUU44foNeehQB0G/YK3eteg2maZOkGNovm4Vp6F+n5KKPcs54uNSI5YlanpnKGgdoG2cxICCGKSXRsAl2nrGbwp5s48ONbqKaD7UoLlp/J7ulQFEUCjxIgwUcZNrhjfWIn3MRXD3UhqeVDAEyqtZrB7cM8XDMhhCj/3NOpB5LOEC17A7kPMm+WfEolTIKPMs6Z9bRlv8fBtzJq0j72rJvH4pgjbNx3SnJ9CCHEZcidTn2w9jOByjn+MuqyxmgNSD6lkiRzPsoLWyD76g+i4V+fcHbN2zyRlb38S9KtCyFE0bmnU7fi4IHz6Qw+1fsCSo58SpJOvfhJz0c5ER2bwF1/tibTtNBO/Ztrlb8ASbcuhBCXwz1PUn9tAyHKaRLNqjlSqecuJ4qPBB/lgG6YjF0cxwmqsEC/DoDh51P/mud/xi6Ok+5BIYQoJGeeJAWD4dpSAL5wRGHPNSAg+ZRKhgQf5cDm+CROpGYC8JneF8NU6KVto5Fy2FXmRGqmpFsXQohCcuZT6q1uo5F6lGTTn9luqdRB8imVJAk+ygH3br/9Zm1WGO0AGGFZWmA5IYQQBdNUhXG3NGOEZQkAX+q9ScMvRxnJp1RyJPgoB3J3+33kuBWAW9XfqMOJAssJIYTISTdMNu47xeKYI4Qlb6O1uo9zpg8zHVGuMqqC5FMqYbLapRxwdg86Z2bvMBvyq96CrlocD1t+ZLxjmHQPCiHEJUTHJjBh6S7Xe+lX1imgQcJVdzC1Sy9OnM2kRoCN1vWq4O8jH48lSXo+ygFnunV3H+n9AbhLW0N1kqV7UAghLsI9oRhAhLKf67WdOEyVe3d3JCXDTv/WdejcKFgCj1IgwUc5ERURyvS72xISZAPgN6MFMUZDfBU781vHSPegEEIUIHdCMcA112Ox0Zkj1JCEYqWsSMHH5MmTad++PYGBgdSsWZMBAwbw119/5Shjmibjx4+ndu3a+Pn50b17d+Li4oq10hVVVEQoG56/gW8e7sR7d7XBt8f/AdAwfi6x+w8x+JNN7Dh8xrOVFEKIMsY9oRhAuJJAH3ULANMdt+ZIKCZKR5GCj3Xr1jFy5Eg2bdrEypUrcTgc9O7dm7S0NFeZ119/nbfffpupU6eyZcsWQkJC6NWrF6mpqcVe+YpIUxUiG1anf+s6NO02CGo0hcxk/v7hPTbuP8VHa/ZJ9C6EEG5yrwR8VFuKqpis1K9lr1m3wHKi5BQp+IiOjmbYsGG0aNGCa665hhkzZnDw4EG2bdsGZPd6vPvuu4wZM4bbbruNiIgIvvzyS9LT05kzZ06J/AIV2eHkDFYHDwHgulPzsZFFdFwiHSet4uN1+zh8Ot3DNRRCCM9zXwkYyilu034BYJqjX4HlRMm6olk1ycnJAFSrlr3KIj4+nsTERHr37u0qY7PZ6NatG7/99huPPvponmtkZmaSmZnpup2SkgKA3W7HbrdfSfXycF6vuK/rKV2nrMFCA9bagqmrnORObR2z9Bs5eTaLycv2MHnZHva+0vvSF8qHt7VVSZK2Kjxpq6KR9iq8i7VVm7qB1K9qIyE5g0fUH/BRdH7Tm7PdbAyAj2oSWtmXNnUDK0Rbl9TzqijXU0zTvKw+etM06d+/P6dPn+aXX7KjyN9++40uXbpw5MgRateu7Sr7yCOP8O+//7J8+fI81xk/fjwTJkzIc3zOnDn4+/tfTtUqjK0nFGbvU7lbXcH/rF9y2Ayme+bbOLCgKiZDGxq0qyFDMEIIAXDyTDJD9z+Fr2JnSNaLbDRaYKLwTEsHYQGerl35l56ezpAhQ0hOTiYoKOiiZS+752PUqFHs2LGDX3/9Nc85Rcm55NM0zTzHnF544QWeeuop1+2UlBTCwsLo3bv3JStfVHa7nZUrV9KrVy+sVmuxXru0bY5PYv7mLWgKzNN78LhlEXWVkwzUfuVbvTuaAvPjVaKub39Z+T+8qa1KmrRV4UlbFY20V+EVpq32znkGX8XOdqMRm83s9AWqAtUatuTma8NKs7oeVVLPK+fIRWFcVvDx+OOPs2TJEtavX0/duhcm64SEhACQmJhIaOiFpZ/Hjx+nVq1a+V7LZrNhs9nyHLdarSX2YivJa5eWk+kOMnVnQOfDJ46+jLHO4TFtMQv067Ebqqvclfyu3tBWpUXaqvCkrYpG2qvwCmyr9CSaH/kOAL+ez/FWlbYVPqFYcT+vinKtIk04NU2TUaNG8f3337N69WrCw8NznA8PDyckJISVK1e6jmVlZbFu3To6d+5clIcSl5B7YtRs/UZOmwGEq8e4Rd1YYDkhhKiQfv8Yss5CrZY063anJBTzsCIFHyNHjmTWrFnMmTOHwMBAEhMTSUxM5Ny5c0D2cMvo0aOZNGkSCxcuJDY2lmHDhuHv78+QIUNK5BeoqJwp1519H+n48rmjDwCjLItQMKjiZ8UwTVl6K4So0GL3H+Ls+qnZN65/BgqYBiBKT5GCj2nTppGcnEz37t0JDQ11/cybN89V5tlnn2X06NE89thjtGvXjiNHjrBixQoCAwOLvfIVmXvKdefL6Ev9JlJMfxqrR7hJ3cqZc3aGfvY7XaesJjo2wXOVFUIID9ENk79/eJcA8ywJ1nroTftd+k6ixBV52CW/n2HDhrnKKIrC+PHjSUhIICMjg3Xr1hEREVHc9RZkZzyddndbQipnD62k4s9MPXtp7eOWhXA+mXBicgYjZm2XAEQIUWEcPp3Ox+v2cf3EH7n+1HwAXk/rS8fJqyUPUhkge7uUc1ERofz6XE9mP9iRKn5WvnD0Ic200UL9l57qHwCu/Qxk7wIhREXRdcoaJi/bw00Z0QQrKRw0arDE6OzKg9R1yhpPV7FCk+DDC2iqgqoqnDln5wyBfK33AuBxyyKcoYfsXSCEqCh0w6SKnxUbWQy3LAWydwLX0VxlqvhZ5cuYB0nw4SXc9yT43HEzGaaVNuo/XKfuLLCcEEJ4o83xSZw5Z2eI9jM1lTMcNoNZoF+fo8yZc3b5MuZBEnx4CfcltSeowhz9BgCesHwPbhtJy9JbIYS3O56akbPXw9Efez5preTLmOdI8OElnEtvnaY5+pFhWmmn/k1XNRaA0Mq+l5XtVAghypOagb7cpa2hlnKGI2Z1vtW7FVhOeIYEH17CfektwAmquno/RlsWACbj+jVHU2V9uxDCe22OT+J40hlGWgvu9VCQL2OeJsGHF4mKCGXywAic8YV778fnXVOp7OfD4pgjbNx3SiZaCSG8yqrdxwAYNnML2xa9R02SOGpWy9Pr4fz6JV/GPEvyynqZwR3r079NHWIOnuHE2UyS9gyl9p6ZBG97l/6/BuJ86YVW9mVcv+ZERYRe/IJCCFHGRccm8OS8GKZ0AItpZ4Rbr0cWOfcbCZH3vjJBgg8v5O9joXOjYADW2B+g2u7ZXMMeuqixbDBaAhcSj027u628CIUQ5ZZumIxdHIduwqGzcLu6llAliQSzGvP17gBU9rMw/tYIQoKyh1qkx8PzZNjFi+mGyYurTuaZ+wGSeEwI4R02xydxIjUTu6HwwU6DxyxLgOxhZ2evR/I5ByFBvkQ2rC6BRxkhwYcX2xyfREJyhmvuR3v1b7qcX/kCknhMCFH+XVguazJIW0Pt870e8/QeBZQTZYEEH17M+WJzX/nypFvvR+5yQghR3jiXywaoWYy0LAbgQ0d/MvHJt5woGyT48GLuLzb3lS/d1B0FlhNCiPJAN0w27jtFYvI5qlXyYai6ilrns5nm7vWQZbVlj0w49WLOxGMJyRmcoCpf6r151PIjT1vmsy6rFaDIi1IIUe5ExyYwYekuEpKze239yeARW/YKlw8cA/Pk9ZBltWWP9Hx4sdyJxz529CPNtNFKjae3uhWQF6UQonyJjk1gxKztrsAD4D5tBdWVFM761GSp2dV1PCTIxnRZ0VcmSfDh5dwTjyURxBd6HwCetnzHQ13qkekwJOmYEKJc0A2TCUt35Zi1FkA6j1h+AOCv0IEEVfLjnUGt+ebhTmx4/gYJPMooGXapANwTj51JCidz+SqacIhjG+fy2YbOgCQdE0KUfc4VfO4e0KKpqpxln1mbw1UjOf2P3bWsVpRd0vNRQTgTj6n+VZmakd37MdqyAA0duJB0LDo2wZPVFEKIAuVemRfEWR6y/ATA+/ptHEpTMU1ZwVceSPBRgTi7LL/Qo0gyA2ioJjBA3QBI0jEhRNmXe2Xew5afCFLS2W2EscTRiTd3WsgyFFnBVw5I8FGBOLss0/BjuqMfAE9YFmDFAUjSMSFE2dYhvBpV/LKzlgaTzAPaMgDeddyBiYqqmFT3t8gKvnJAgo8KxL0r8iu9NyfMytRTTzBIW1NgOSGEKCs0VeG127P3pxppWUQlJZMYoyHLjXYAPBWh8+qAFrKCrxyQ4KMCce+KzMDGB44BAPzXshBfMvMtJ4QQZUlURCiTugcxVFsFwOuOQTh36wa4sVktD9VMFIUEHxWIM+mY0zf6DRwyalBTOcP92nJAMgEKIcomZ0bTxTFHuO7oZ/goOtu0VtRtG8VVwZUIDvAh0OrpWorCkqW2FYgz6djwWdsBsGPhbccdvOMzjeGWJczWezKuX1s0VcHQPVxZIYQ4zz2jaSPlMMt9FoMCmdeP4fVu12CaJmkZWfy8ItrTVRWFJD0fFYx70jGAxUYX9hhhVFbSmdv8d8nzIYQoU3JnNH3G8i2aYhKtt2foMp3o2AQURcFmkY+z8kR6Piog96RjJ85mYp4ZC2sfodnBOSxYO4Tp29O4t1MYQZ6uqBCiQtMNk7GL41ypAK5R/iFK24JhKrzluBMTGLs4jl7NQzxZTXEZJFSsoJxJx/q3rkOzbv/hdPU2KI5zpK+azN7jZxm/dBcAq3Yf83BNhRAV1eb4JE6kXpgM/3+WeQAsNLqy16wLwInUTEkPUA5J8CGY/fu/PHI0O+/HXdoa6iuJ6CYcOgtPzI1h9qYDnq2gEKJCcl/231XdSVctjixT4x3H7QWWE+WDBB8VnG6YjFkUxxazKWv0a7AqOk9bvgVwZQscsyhOsp4KIUqdc9m/gsHzlm8AmK3fyGGzZr7lRPkhwUcF595d+YZjEAC3ahtpoRw4f9TMU04IIUqDMz3ArepvRKgHSDX9+MAxMEcZSQ9QPknwUcG5d1fuMhvwvd4VgBctcwATq5q3nBBClCRnTo8fdhxl6LUhPHO+N3a6ox9JuabCj+vXXDKalkOy2qWCy91d+Zb9Tvqqm+iixdFN38FvtMq3nBBClAT3nB4AD2o/EWY9wXGq8rnex1UuJMjG+FtbSHqAckqCjwrO2a3pfKEfoQYz9Zt41PIjY6xz6G+PILhyJenWFEKUOGdOD+cMsyDSGGVZBMBb9jsYcWMEDYIrUTMwe6hFejzKLxl2qeCcWU/dX8IfOvpzxqxEY+UQ/dVfuDkihM3xSTLpVAhRYnTDZMLSXbi/ywy3LKWqcpa9Rh0W6Nczd8shbmlVm8iG1SXwKOck+BBERYQy7e62rn1fUghgutEfgNHat8za8DeDP91E1ymriY5N8GRVhRBeanN8kqsHFiCEUzygLQNgiuMuHGgkJGfI5HcvIcGHALIDkF+f68k3D3figS4NmG30It0nmBDlNA9o2fslJCZnMGLWdglAhBDFLvek9qcs3+Gr2NlsNGGV0bbAcqJ8kuBDuGiqQofwaizdkUCG6cOvVe8AYIRlMVVJwQRXOmMZghFCFCf3Se3NlQPcoa0HYLJ9CLgNDMvkd+9Q5OBj/fr19OvXj9q1a6MoCosWLcpxftiwYSiKkuOnU6dOxVVfUcKc6YzthsKj/3Yl1mhAkHKOJyzfu8pIOmMhRHFzTn4HkzGW2aiKyRI9kj/Mq11lJKeH9yhy8JGWlsY111zD1KlTCywTFRVFQkKC6+enn366okqK0nOhS9PERGWiYygAd2uraKgcyaecEEJcOefk957qH3TR4sg0rbzuuCtHGcnp4T2KvNS2T58+9OnT56JlbDYbISGyy2B55OzStKpgN2Cj0YKVelt6adt5wTKHh+z/l6OcEEIUl6hmwURWXQBp8IUexWGzBgCqAhMHREhODy9SInk+1q5dS82aNalSpQrdunVj4sSJ1KxZM9+ymZmZZGZe2LUwJSUFALvdjt1uL9Z6Oa9X3Nf1Jm3qBlK/qo1TqeewGwpgMskxlO7qn9yo/UF3YyfxQe1oUzdQ2vE8eV4VnrRV0VSU9tINk2+3HebMumk8nhmP4VeN1jeP4z27jeqVfGhVtwp+PtpF26GitFVxKKm2Ksr1FNM0L3vmoKIoLFy4kAEDBriOzZs3j4CAAOrXr098fDxjx47F4XCwbds2bDZbnmuMHz+eCRMm5Dk+Z84c/P39L7dq4gqdyYQ3d2pU8YHIWgbtEmbzHzOaJFsYvzR7BRSZqyyEKD7L9p1jYvIzVFdS+bPufRyocYOnqySKKD09nSFDhpCcnExQUNBFyxZ78JFbQkIC9evXZ+7cudx22215zufX8xEWFsbJkycvWfmistvtrFy5kl69emG1Wov12t5mZdxRzsXHMH6bSpapEGSe5WefpwniLBubj+XlQ+24N7Ied15bt8KPwcrzqvCkrYrG29tr3tZD/O98YrGn1TkMt/zIXqM2tzpeQ0fj5X7NGdQurFDX8va2Kk4l1VYpKSkEBwcXKvgo8fTqoaGh1K9fn7179+Z73maz5dsjYrVaS+wJVJLX9ha9WtTmp39jmH5vB06mO6gZ6MvhP/9L8x2TaBT3HgmZbzN2SRofrYtnXL/mMhaLPK+KQtqqaLyxvXTD5KXFuwGFMOUY92vLAZjoGMo5I/uj6aXFuxncMbxIX3C8sa1KSnG3VVGuVeJ956dOneLQoUOEhsqHU3nUIbwa/VvXYf+JVPpvbsp+I4QaSgqPWZYAkJCcwfBZ25m96YBnKyqEKFfcl+uPsczBpjhYr7dkrdG6wHLCexQ5+Dh79iwxMTHExMQAEB8fT0xMDAcPHuTs2bM888wzbNy4kQMHDrB27Vr69etHcHAwAwcOLO66i1KiGyZjFsVhx8Kk80tvH9SWEaYcc5UZs0gSjwkhCs+5XL+zGkuUtgWHqfKK4x5Aybec8C5FDj62bt1KmzZtaNOmDQBPPfUUbdq04eWXX0bTNHbu3En//v1p3Lgx9913H40bN2bjxo0EBgYWe+VF6XD/5rHKaMuvegtsip0xljkFlhNCiIupGeiLhs7Llq8B+FrvxV6zbr7lhPcp8pyP7t27c7E5qsuXL7+iComyJ+c3D4UJjvtYpj5PlLaFLvpONhgt8yknhBB56YbJ5vgkEpPP8bDfWpqah0gyA3jHcXuespLR1HvJeklxSbm/eew16/K13guAcZavsODIPn7sLBv3nZLhFyFEvqJjE+g6ZTWDP93EhPm/MtyYC8DbjjtJISBHWQXJaOrNJPgQl3Rhz4UL3nHcTpIZQGP1CPdoKwGYuuYfBn+6ia5TVsvOt0KIHKJjExgxazsJydk9pE9avqOKksZuI4xv9J45yoZW9mXa3W1lFZ0Xk+BDXJJzzwX37x8pBPCmYxAAT1oWUI0U17nE5AxGzNouAYgQAsgeahm7OA5nn2gT5SB3a6sA+J/jXnQ0KvtZeGdQa755uBO/PtdTAg8vJ8GHKJSoiFCm3d02Rw/IXL0HcUZ9gpR0nrHMdx03z/+MXSwrYIQQF3bLzmbysuVrNMXkJ70DG40WACSfcxAS5Etkw+oy1FIBSPAhCi0qIpRfn+vJNw93YlSPhhiojLffB8Bd2hpaKPE5yp9IzZQVMEKIHJPRb1E30UWLI8O0MskxpMBywrtJ8CGKRFMVIhtW5+pa2Uunt5hNWaJHoiom/7PORMHIUV7eTIQQzknrlTjHS9ZZAHzk6M9hs2a+5YT3k+BDXBb3N4mJ9qGkmTauVfdyh7a+wHJCiIrJOWn9cctCQpTT/GvU5GP9lhxlZFltxSLBh7gs7itgjlGNd8+v0X/e8g2VOQvIm4kQIlvc0WSaWxN4UFsGwHjHfWTik6OMLKutWCT4EJfFuQLGaYYexd9GHaorqa7JpzdHhLA5PkkmnQpRwX2/7TAPJH+EVdH5RWnHGqON61xIkI3psqy2winxXW2F94qKCGX63W0ZvySOxBR42XE/c31eZaj2M/P17ny+AT7fcIDQyr6y860QFczh0+mcTrNjmCZn//g2e5IpVqrd8Q4TUyuTpRs0DQmiQ3g16fGogCT4EFckKiKUXs2zezhW7mrAos2rGaD9xivWLxiY9T9MVFfeD0kaJETF0XXKGgD8yeBn21egwDT7rbz39SHgEAAHXuvrwRoKT5JhF3HFNFWhQ3g1lu5IYKJ9KKmmH63V/QzS1gKS90OIiuj+Lg0AGG1ZQKiSxEGjBtP1fnnOi4pJgg9RLJxJhE5QlXccdwDwnGVujsynkvdDiIpBN0yiYxNppvzLA+cnmb7suD/HJNPo2ET5MlKBSfAhioV7Po8v9d7sNupRVTnLi9Y5BZYTQninzfFJHEtOZ5L1cyyKwQ96R9YarXOUSUjOkC8jFZgEH6JYuOfz0NF40f4ghqlwh7aeSDXOdU52vhXCe+mGycZ9p1gWm8AQ7WfaqP+QavrxP/u9+ZaXLyMVlwQfoljk3vn2D/NqZuk3AvCq5QtsZAGy860Q3io6NoGuU1Yz+NNNLNv4J89a5gHwpuM/HKdqvveRJIQVlwQfoljkzvsB8IZjEMfNKjRUExhhWZLjnOx8K4T3iI5NYMSs7SQkZ/dkjLV+TZCSzp/GVXyt98pTXkGSEFZ0EnyIYuPM+xESZAMgFX/Gn+9uHaEtoaFyxFXWOegyYekuGYIRohzTDZOxi+Ncr+nr1B3cqm1EN5Xs4ddcHzPOjB6S0bRik+BDFKuoiFA2PH+Da+fbn4yOrNZbY1McTLR+wYWwI/t/MulMiPLNudINwJdMXrHMAOBL/SbizPA85UMq+0rOHyFJxkTxc+58mz2ZTOFlx/1Eqv9HJ3U3d2rr+FbvnqO8TDoTovxyf/2OtiyggXqMo2Y13j6/5N7p3sj69IkIlYymApCeD1GCnJPJDps1eOf8xnNjLLOpwZl8ywkhyh/n6zdC2c/D2o8AvGR/gLP45yjXJyKUyIbVJfAQgAQfogS5r4D5XL+ZnUYDqihpjLfOdJWRSWdClF+6YWIYJsG+Cq9bP0VTTJbokaw22uYoJ69zkZsEH6LEuK+A0dF4zv4IDlOlr7aZm9TNgEw6E6K8ci6tHfr579xpX0xz9V9OmwFMyCenh7zORW4SfIgSFRURyuSBEagK7DIbuPZ2eMU6k7siAnhrxd/M2vSvrHgRohxxX1obriQw2rIAgFfsd3OKyq5yIUE2psvkUpEPmXAqStzgjvXp36YOMQfPcCq5GUmr/qTmuQO03fMWcx2P8tKiWD5c8w/j+jWXNykhyjj3pbUKBq9ZP8Wm2Fmvt+R74zoAgmwWpt1zLZ2ukjkeIn/S8yFKhb+Phc6Ngkmxqzx85j4MU+E/lnV0VXcC2Utuh8/azuxNBzxbUSHERbkvrR2iraajuoc008aLjodwZvFIyXSgKooEHqJAEnyIUqMbJmMWxbHNbMKXem8AJls+w58LS/XGLIqTIRghyjDn0tq6yglesGRvHPmm4z8cNmvkW06I/EjwIUqNezKxNxyDOGwGE6ae4HnLNwWWE0KULTUDfVEwmGL5hAAlg9+NpszUb8q3nBAFkeBDlBr3b0Lp+PKc/WEA7rWspLMam285IUTZ4NyxNjH5HI/4raWLFsc504dn7Y9g5vookaW14lJkwqkoNbm/CW0wWvK140busazidesnRGW+xln85RuTEGVMdGwCE5buIiE5g7rKcZb7fA0KTHHcxb9mSJ7ysrRWXIr0fIhS4550zGmyYwgHjRrUVU7yomU21Sv58N7Pf7Pj8BnPVFIIkYP7sloFgzcsn1BJyeR3o6lr7paTLK0VhSXBhyg1zqRj7t+H0vHl/+zDARhiWUO/gF1s2p/ER2v2ycRTITws9461d2uriNR2kW7a+D/7o5ioVPaz8M6g1nzzcCc2PH+DBB6iUCT4EKUqKiKUaXe3zdED8rvZjFn0AeDRM+8SRBrRcYl0nLSKj9ft4/DpdE9VV4gKzX1ZbX0lkRfOTw6f4riLg2YtAJLPOQgJ8pV9W0SRSPAhSl1URCi/PteTbx7uxHt3tQZgYsadxBu1CFWSeNn6NQAnz2Yxedkeuk5Z48HaClFxOSd/a+i8Y/0IfyWTjXpzvtJ75VtOiMKS4EN4hKYqRDaszi2talPFz8o5fHnGPhzDVLhDW0/U+b1fAKr4WWUIRohSphsmJ8/3eozQltBW/YcU04+n7cPzrG6RSeKiqCT4EB61OT6JM+fsAGwzmzDt/N4vk62fUZPTAJw5Z5fcH0KUIuemca/8uJuWyn6esHwPwMv2+zlKcI6ysqxWXA4JPoRH5e6ufddxBzuNBlRVzvKmdToKBgAb/jkhvR9ClAL31S2+ZPKO9SOsis4PekcWGV3ylJdlteJyFDn4WL9+Pf369aN27dooisKiRYtynDdNk/Hjx1O7dm38/Pzo3r07cXFxxVVf4WVyd9fasTDaPpIM08r12k7u1VYCMHXNPrpOWU10bIInqilEhZB7dctzlrk0Uo9yzKzCS/YHwG2tmiyrFVeiyMFHWloa11xzDVOnTs33/Ouvv87bb7/N1KlT2bJlCyEhIfTq1YvU1NQrrqzwPs7cH+7fm/aZdZjkGALAC5Y5NFIOA5CYnMGIWdslABGihLivbrlO3cH9luUA/J/9Uc4Q6Co3tm8zWVYrrkiRg48+ffrw6quvctttt+U5Z5om7777LmPGjOG2224jIiKCL7/8kvT0dObMmVMsFRbexZn7A8gRgHyl92ad3gpfxc571g+x4sAETGDsYtl8ToiS4BwGrU4yb1mnA/CVoxfrjWtylAsOtMlQi7gixZpePT4+nsTERHr3vpD1zmaz0a1bN3777TceffTRPPfJzMwkMzPTdTslJQUAu92O3W4vzuq5rlfc1/VGpdlWNzQJ5qMh1/Dasj0kpmSQqSuAwv/ZHyVafY4W6r88Y5nHZMdQAE6kZrLpn+NlZpKbPK8KT9qqaEq7vYL9Ldg0g7e0j6mpnOFv40IvJIBVNVGV7HJl7W8oz63CK6m2Ksr1FNM0L/srpKIoLFy4kAEDBgDw22+/0aVLF44cOULt2rVd5R555BH+/fdfli9fnuca48ePZ8KECXmOz5kzB39//8utmijHtp5QmL1PxTAVeqtb+MTnHQDuy3qOX8xWDG1o0K6G9HwIURKuOr6clkdmk2la6Z/1P/aY9VEwMVF4pqWDsABP11CUVenp6QwZMoTk5GSCgoIuWrZENpZTlJzdcaZp5jnm9MILL/DUU0+5bqekpBAWFkbv3r0vWfmistvtrFy5kl69emG1Wov12t7GU221OT6J+Zu3oClgmLDCaM9Xjl7ca1nJm9Zp3GqfzPz4KlStexUdr6rOtfWrerz7V55XhSdtVTSl0V6rdh/jtWV7OJqcwdXmARb5zAMFXnUM4S+zHhbFRDcBTD6I03j3rtbc2KxWidTlSshzq/BKqq2cIxeFUazBR0hI9u6GiYmJhIZemIh0/PhxatXK/8lqs9mw2Wx5jlut1hJ7ApXktb1NabdVp0Y1qRbgR0LyhSW4Ex1Daa/uoZl6iNe1adxrf57318bD2nhCK/syrl/zMjHxTZ5XhSdtVTQl1V7RsQk8NudPTMCPTN73+RAfxcFK/Vq+1nsDCo7znYwhQTbG39qiTLzWLkaeW4VX3G1VlGsVa56P8PBwQkJCWLlypetYVlYW69ato3PnzsX5UMJLuU9AdcrEh1H2/3LO9OE6LZZHtR9c52QFjBCXJ/ey2rGWr2mkHiXRrMqz9ocBRTaNEyWmyMHH2bNniYmJISYmBsieZBoTE8PBgwdRFIXRo0czadIkFi5cSGxsLMOGDcPf358hQ4Zc/MJCnBcVEcr0u9sSEnShR2yfWYfxjvsAeNryLa2VfwBkBYwQl8l9We2t6m8MsazBMBWetD/GabKHvGXTOFFSijzssnXrVnr06OG67Zyvcd999zFz5kyeffZZzp07x2OPPcbp06fp2LEjK1asIDAwsKBLCpFHVEQovZqHsDk+iQ3/nGDqmn3M07tznbqTW7RNfGD9gL5ZE0khe/bbidRMNscnEdmwuodrLkT54FxWe5VylEnWzwD4UO/PRqNFvuWEKE5FDj66d+/OxRbIKIrC+PHjGT9+/JXUSwjX5nMX3vwUXrA/RCtlH/XUE7xl/ZiH7U/hzBCy4Z8TdAivJt/QhLgI3TDZHJ/E3mOp+JLJh9b3CFAy2GQ0413H7XnKy6ZxoiTI3i6izHN/80vFn8fsT5BpWuilbeMRt/kfkoJdiItzbhg3+NNNTF2zj/GWL2mmHuKEGcTjWaPQ0XKUl03jREmR4EOUec4U7E6x5lX8z3EvAM9a5tFe2eM6JxNQhcif+4ZxALep67nLshbDVHjCPooTVM1zH9k0TpQUCT5EmZffCpjZ+g0s0jtjUQw+8PmA6iQDuGbuT1i6SyagCnFe7pUtjZTDvGqdAcB7jtv4zYjIUV42jRMlTYIPUS7kXQGj8KL9IfYadQhRTvOedSoqBpAdgCQkZ7A5Pslj9RWiLHFf2RJAOh9b38FfyeQXPYIP9IGucqN6NJJltaJUSPAhyo2oiFA2PH8Do3o0AiAdX0bYnyDdtNFVi+Mpy7c5ysssfSGyOV8LCgZvWafTUE0gwazGaPtIDLePgatrBciyWlEqJPgQ5YqmKnRpFOy6/Y9Zl+ftDwEwyrKYm9TNrnMyS1+I7CGXHYfOADBCW8pN2lYyTQsjskZziso5ysprRpQWCT5EuZN7AuoSowufO/oA8JZ1Oo2UwzJLXwgurG75fMMBrlN38IxlPgDjHMOIMRvlKCuvGVGaJPgQ5U5+E1AnOYawUW9OgJLBJ9a3uamRH1HvrmfWpn9l4qmokGZvOsDw86tb6irHed86FVUx+cbRg7l6zzzlZWWLKE0SfIhyKSoilMkDI3C+V+pojLT/lyNmda5SE+myYwz/HE/hpUWxkvtDVDi6YTJmURwAvmQy3fouVZWzxBhXubYpcJKVLcITinVXWyFK0+CO9enfpg4xB89w4mwmfx48w6Mbn2SBzwR6adt5wvyedx13kJCcwfBZ25k4oAVDOzXwdLWFKDHO7KUb/jlx/ojJG9aPiVAPcMoMZETWk2Ti4yo/tm8zhnUJlx4PUeok+BDlmr+Phc6NgtENkyfmxgBX8YL9Qd72mc5oy/fsMeoRbXQAYMyiOO7qUF/eaIVXio5NYMLSXa4kYgCjtEX00zZhNzVGZI0mgZx7HwUH2uT1IDxChl2EV3DP6fG9cb1rAuo71o9oocTnW04Ib5E7eynATepmnrFmLz9/yfEAm81mee4nq1uEp0jwIbxC7pwekxxDWKtfg5+SxWc+b1GT0wC8vDiWHYfPeKCGQpQM3TCZsHQX7tOqmysHeMc6DYAZjpuYp/fIcz9Z3SI8qdwOu+i6jt1uL9J97HY7FouFjIwMdF0voZp5h5JqK6vViqZply5YRLm/weloPG5/nO+VcVytHuFTn7cYlDWWvcfhozX7+HBoW+luFl5hc3xSjh6PYJL51Oct/JVM1ustedVxd577KMjqFuFZ5S74ME2TxMREzpw5c1n3DQkJ4dChQyiKvOgupiTbqkqVKoSEhBTrdZ25P9zfhFPx50H7MyzyGcs16n7etH7MKPvjRMcl0nHSKh6+7ir6tgqlblX/YquHEKVJN0w2/HPSddtGFh/7vE0d5RT7jFBG2R/Pd6facf2ay+oW4VHlLvhwBh41a9bE39+/SB9ghmFw9uxZAgICUFUZcbqYkmgr0zRJT0/n+PHjAISGFt+bnzP3x4hZ23N0Px80azE860lm+UziFm0T+8xQ3nHcycmzWUxetofJy/Zw4LW+xVYPIUpL7gmmCgZvWz/iWnUvyaY/D9ufJoUAV/lRPRrRpVEwHcKrSY+H8LhyFXzouu4KPKpXr37pO+RiGAZZWVn4+vpK8HEJJdVWfn5+ABw/fpyaNWsW6xBMVEQo0+5um2fG/2azGWMcD/KG9ROesCzksFmDb/XuAPhbNTb8c5JOV8l+FqL8cE4wdQ+0n7PMpa+2mSxT45Gsp9lv1nadC63sy5O9GstzXJQZ5eoT2DnHw99fusnLM+ffr6hzdgojKiKUX5/ryTcPd2JUj4au49/q3fnAMQCAyZbPuF79E4B0u87Qz36XRGSi3Mhvgund2kqGW34A4P/sj/J7rpUtMr9DlDXlKvhwkvka5VtJ//00VSGyYXWurhWY4/hbjjv5Xu+KRTH4yPoezZUDrnOJyRmMmLVdAhBR5uWeYNpd/YMJlpkAvGm/k8VGV9c5VYHJAyNkfococ8pl8CFEYeTNYaDwnP0Rfju/B8wMn9epTfZkPfP8z9jFcbIXjCjT3JeVRyj7+dD6PppiMt/Rjan6ANe5kd2vInbCTQzuWN8DtRTi4iT4EF7LuQLGvZ/FjoXh9if5y6hLLeUMM3xeJ4izrvMnUjMlEZko05xBdbiSwEyf16mkZPKLHsGLjgfB7dne9eqa+PuUq2l9ogKR4KOUDBs2DEVRUBQFq9VKrVq16NWrF1988QWGYRT6OjNnzqRKlSolV1Ev4r77rXsAkkIl7s96lkSzKk3Uw3zh8yZ+XPg2uSw2gY37TkkPiChTdMNk1qZ/Gbsoloa+KXztM5lgJYWdRgOG25/E4bZ+QBKIibKuwgYfumGycd8pFsccKbUPmqioKBISEjhw4ADLli2jR48ePPHEE9xyyy04HI4Sf/yKyLkCJqRyziGYowRzX9ZzJJv+tFP/Zpr1Paxk/w2+2vgvgz/dJJNQRZmxavcxuk5ZzUuLYjl+IpEPjYnUVU4Sb9RiWNZzpOGXo7xMMBVlXYUMPqJjE+k6ZTWDP93EE3NjSu2DxmazERISQp06dWjbti0vvvgiixcvZtmyZcycOROAt99+m5YtW1KpUiXCwsJ47LHHOHs2e1hg7dq13H///SQnJ7t6UcaPHw/ArFmzaNeuHYGBgYSEhDBkyBBXPo2KzrkCZvaDHaniZ3Ud/8usx/1Zz5Ju2uiu/cmb1ukoXOiFkkmooixIyoQn5saQkJyBL5l87vMmTdVDHDercI/9BU5R2VVWJpiK8qLCBR8//3WKkXP+yDFbHDz3QdOzZ0+uueYavv/+ewBUVeX9998nNjaWL7/8ktWrV/Pss88C0LlzZ959912CgoJISEggISGBZ555BoCsrCxeeeUV/vzzTxYtWkR8fDzDhg0r1d+lLNNUhS5XB/Pa7S1zHN9uNmaEfTR2U6O/9hvjLV/C+UWMMglVeJpumEzYbiHLULDiYKr1fdqrf5Ni+nNv1vMcNmsC8M6g1sx5qKNMMBXlRoUKPnTD5PVV+8nvY8R5bMLSXaX+QdO0aVMOHDgAwOjRo+nRowfh4eH07NmTV155hfnz5wPg4+ND5cqVURSFkJAQQkJCCAjIzmD4wAMP0KdPH6666io6derE+++/z7Jly1y9JiJbVEQo0+9uS0iQzXVsnXENT9lHYJgK91lW8rTl2xz3kUmowhN0w2TO7/9yTyMdDQfvWqdyo/YHGaaVB7OeYY9Zz1U2JMiXzo2CZYKpKDcqVPCx5UASx1KzCjxvAgnJGaX+QWOapiv3xZo1a+jVqxd16tQhMDCQe++9l1OnTpGWlnbRa/zxxx/079+f+vXrExgYSPfu3QE4ePBgSVe/3ImKCGXD8zfwzcOduDcy+1viUqMzLzuGAfC4ZRGjtIU57rPhnxPS+yFKTXRsAl2nrGZy9F9sPmbylvVj+mqbyTQtPGJ/ii1m0xzlc+/qLERZV6GCj+OpmYUsV7ov5N27dxMeHs6///7LzTffTEREBAsWLGDbtm18+OGHwMWzgaalpdG7d28CAgKYNWsWW7ZsYeHC7A/PrKyCg62KzJmIrI/b2PgsvRcT7UMAeMb6LY9qS13npq7ZJxNQRalwpk5PSM7ANA2GnZvBAG0DdlNjlP2/rDeuyXOfvDlthCjbKlTwUTPQdulClO4LefXq1ezcuZPbb7+drVu34nA4eOutt+jUqRONGzfm6NGjOcr7+Pjk2eJ+z549nDx5ktdee43rrruOpk2bymTTQnLmAnH6VL+F1+3/AeAF6zc8oC1znZMJqKKk6YbJ2MVxrllHY7SvGWxZg24qjLaPZKXRLs99ZFmtKI8qVPDRvkE1agX6UNACNIWSfSFnZmaSmJjIkSNH2L59O5MmTaJ///7ccsst3HvvvTRs2BCHw8EHH3zA/v37+frrr5k+fXqOazRo0ICzZ8/y888/c/LkSdLT06lXrx4+Pj6u+y1ZsoRXXnmlRH4Hb+OeC8TpI30A7zluA+Bl69fcra0EZAKqKDnOpf/vrPyLE6mZgMlLllncb1kBZO/X8qPRKc/9FGRZrSifKlTwoakKz954FUCeAMR5uyRfyNHR0YSGhtKgQQOioqJYs2YN77//PosXL0bTNFq3bs3bb7/NlClTiIiIYPbs2UyePDnHNTp37szw4cMZNGgQNWrU4PXXX6dGjRrMnDmTb7/9lubNm/Paa6/x5ptvlsjv4I3ym4T6juN2PnLcCsCr1hnco61wnTuRmsk7K/+WRGSiWDjndwz+dBNT1+wDTMZbvuQhS3av2wv2B/neuD7P/UIr+zLt7rayrFaUS4ppmmXq3TMlJYXKlSuTnJxMUFBQjnMZGRnEx8cTHh6Or2/Rh0YMwyAlJYXfDqbzyo+7cyy3Da3sy7h+zeWFfJ6zrYKCglDV4o1Rr/TvWFJ0w+SdlX8zdc0/54+YvGiZwyOWHwF4xT6Uz/W+Oe4TWtmXl/s2ISt+GzfffDNWqxVRMLvdzk8//SRtdZ5zfofzTVjB4FXLDIZafsYwFZ53PMR8vUee+43t24xhXcKlx8ONPLcKr6Ta6mKf37lVyHVZUREh3BQRyub4JI6nZlAzMHuoRV7IFZumKnRpFOwWfChMcgwhCwujLIsZa52NDTsfuW3elZicwZPzYpjSwSNVFuVYzvkd2YHHZMtn3GVZi2Eq/J/9URYY150/a+Lsnw2t7CuBhyj3KmTwARdWOwjhzjkB9UKvmMKbjkFkmlaetn7Hs9b52BQ77zjuABRMwDCzs1Dqhol83xKFtTk+6fz8DtDQed36Cbdrv6CbCk/ZR7DY6ArAf67S+Xa/6gpSZI6H8AYVas6HEJeS3wRUgA/025hsHwzAE5aFjLHMxpmaLstQmLDdwpzf/5U5IKJQdMPk++2HAbCRxTTru9yu/YLDVHnCPsoVeGiYdKllYlUhJMjGdJnjIbyEBB9C5JLfBFSAj/V+jLffC8DDlp940/oxFhyASS0/k8nRf0kuEHFJzgmm3247TADpfOkzhd7aNjJNKyPso/nBiHSVtWjZ/84Y1p4Nz98ggYfwGhV22EWIi4mKCKVX8xA2xyex4Z8T51chwEw9irP48ZrlU+7Q1lOVVB63P86xc75oiunKBSKrEISTbpiu+WU7DiXz+YZ4AKqTzEyfKbRUD5Bq+vFQ1jP8bjbLcd+QIF8gTeakCa9T7D0f48ePd+246vwJCQkp7ocRosQ55wU92atJjkRk3+ndeNT+JBmmlRu0P/jKZzKVOYtuSi4QkZP7Mton5sa4Ao+6ygnm+/yPluoBTppB3JU1Nk/gAfB8n6Z5jgnhDUpk2KVFixauXVcTEhLYuXNnSTyMEKUiv3kgPxvXMjTrRZJNf9qpe5nv8z9COeU6L7lAhHuadHctlf0s9HmZhmoCh81g7swaR5zZIEcZ5/yOG5vVKsUaC1F6SiT4sFgsrl1XQ0JCqFGjRkk8jBClJr95INvMJtyZNY5EsypN1MMstr1MS2W/6/zUNf8w+NNNMg+kAtINkwlLd+XZQftGdRvzfF6hhpLMbqMed2SOI97MOTw3qkdDmd8hvF6JzPnYu3cvtWvXxmaz0bFjRyZNmsRVV12Vb9nMzEwyMy9s+JaSkgJkJ0HJvZma3W7HNE0Mw8AwjCLXy5lPzXkNUbCSbCvDMDBNE7vdjqZpxXrtknRDk2C6P309c37/lynL/8Iw4W8jjNsyJ/CFzxs0VQ8x3+d/POUYyQqjHbqpoCkmp8+eY/Q323hnUOsK/03W+Zq+2EaJ3mBzfBKnUs9hPf/1zm7AMG05L1u+RlVM1urXMMo+irNUwpnDw6qaqApEhlfF0B0YesVpr+IgbVV4JdVWRblesWc4XbZsGenp6TRu3Jhjx47x6quvsmfPHuLi4qhePW9ejfHjxzNhwoQ8x+fMmYO/v3+OY84elbCwMHx8fIqz2oLsNn/hhRf4999/S/RxsrKyOHToEImJiTgcjhJ9rJJ06Cy8udOCgkklzjHV+j7dtR0YKLytD2aqvS8BFhjeLHsjwEpWqFa4vQ2FF3hiY/Z3Ow39/D4tywGY4+jJWMf96Gj85yqdjcdUzmTBMy11qsjzQ5Rj6enpDBkypFAZTks8vXpaWhoNGzbk2Wef5amnnspzPr+ej7CwME6ePJlvevVDhw7RoEGDy0rLbZomqampBAYGoiiemTl+6NAhJkyYQHR0NCdPniQ0NJT+/fszduxYV3DWs2dPrrnmGt555518r7FmzRpeffVV/vzzTzIyMqhTpw6RkZF89tlnWCwX78xau3YtN9xwQ57jL774Ii+++CKpqanUrFkT0zQZM2YM0dHRbN++/cp/cTcZGRkcOHCAsLCwMpVevShW7T7G6LkxZBrZY5d3XGWwYD+Mt3zJPZZVAHzj6ME4x31kcSFQtmkmIUG+PN+naYXsBbHb7axcuZJevXp5XQrsVbuP8dqyPRxNzsBhgKpAgHmWD6xTuV7Lnvc22T6Yj/VbANDIXkprmtm5SxWFPL1j3txexU3aqvBKqq1SUlIIDg4uG+nVK1WqRMuWLdm7d2++5202GzZb3nDfarXmaRRd11EUBVVVL2u/EefwgfMaOw6fYfJPe3jh5qa0qlulyNcrqv379xMZGUnjxo355ptvCA8PJy4ujv/7v/8jOjqaTZs2Ua1atRx1zC0uLo6+ffvy3//+lw8++AA/Pz/27t3Ld999B3DJdnGe/+uvv3I8OQICAqhUqRKVKlUCyDHUUtx7u6iqiqIo+f6Ny4s+reqiqBrjl8SRdDaDLrVMFsRrjHXcT7wZykuWWQy2rKGJeojhWU9ynKoAZOoKB09n8ticPyv0ctzy/Ld351xGu3JXIl9sOHD+aPYXm4Yc4lOft6ivHifNtPG0fQTRxoU8/DqgZ3eKXXJvKW9pr9IgbVV4xd1WRblWiQcfmZmZ7N69m+uuu+7ShUvZ99uPsHH/Kb7ffqRUgo+RI0fi4+PDihUr8PPzA6BevXq0adOGhg0bMmbMGKZNm3bRa6xcuZLQ0FBef/1117GGDRsSFRVVpLrUrFmTKlWq5Dg2c+ZMRo8ezZkzZ5g5cyZTpkwBcPUSzZgxg2HDhhXpcbyZMxfIpn+Oc3L3Jiwq6LrCF3of9puhvGedSlv1H36wjWFE1hNsM5sAF3bpmLB0F72ah0j+hnIqOjaBCUt35VnNAnCTupm3rdOopGRy0KjBw/an+cusl6PMkzdeTYPgSrK3lKiQin21yzPPPMO6deuIj4/n999/54477iAlJYX77ruvuB/qshxNzmDnkWRijySz9M+jACz98yixR5LZeTiZw6fTS+Rxk5KSWL58OY899pgr8HAKCQlh6NChzJs3j0uNgoWEhJCQkMD69etLpJ5OgwYNYtSoUTmWTQ8aNKhEH7M80lSFDuHZvVXZCaGyrTVac2vWq+wxwqipnGGuz6vcra3EmZLdBBKSM2Q5bjlV0DJaDZ1nLPP42OddKimZ/Kq34NasV3MEHqoCkwdG8MSNjenfug6RDatL4CEqnGLv+Th8+DCDBw/m5MmT1KhRg06dOrFp0ybq169f3A91WW6ets31f+fLPSkti1s++NV1/MBrfSlue/fuxTRNmjXLm0gIoFmzZpw+fZoTJ05c9Dp33nkny5cvp1u3boSEhNCpUyduuOEG7r333kuOsbmrW7dujtu5J5n6+flRqVIl1yRfcWnP92nKiDl/um7/a4ZwW9YEXrd+wi3aJl61zuBa9W9esj9AGtkB6NQ1/zB1zT+X7HYXZUdBy2hrcpoPfD6go7oHgM8dfZjkGIJO9oqueyPrE9UihNb1quDvI8mlRcVW7D0fc+fO5ejRo2RlZXHkyBEWLFhA8+Z5N+rylIn9rsZy/luG883D+a9FVXh3UGtPVMvV43GpibCapjFjxgwOHz7M66+/Tu3atZk4caKrh6KwfvnlF2JiYlw/VatWvaL6C7ixWa08uUDS8WWU/XEm2QfjMFUGahtY6jOGFsqBHPd1pmWXfCBl3+b4pDw9HtepO/jJ9gId1T2kmn6MzPovrzjucQUeAH0iQuncKFgCDyGogBvL9W1Rk+9HROZ7btHILgxoU6dEHrdRo0YoisKuXbvyPb9nzx6qVq1KcHBwoa5Xp04d7rnnHj788EN27dpFRkYG06dPL3R9wsPDadSokeunuCeVVlRREaFseP4GZj/YkSp+zslXCp/o/RiUNZYjZnWuUhP53udl7tWW4z4MYwLPLdjBwj+OyFBMGaQbJhv3nWKZW4BowcHTlvl8aZ1CsJJCnFGfflmv8qPRKcd9Qyv7uobnhBAVMPhw5+xkKI1Vt9WrV6dXr1589NFHnDt3Lse5xMREZs+ezaBBgy5rCXDVqlUJDQ0lLS2tuKoLZM9c1p3T8UWhaapCl6uDee32ljmObzObcHPmZFbq12JTHPzP+iWfWN+mGimuMsnnHDw5L0Yyo5Yx7nu0fLUxe4jyKuUo3/mM53HLIlTF5GvHjdyWNYEDZt6hs3H9msu8DiHcVMjgo3qADzUCbLSsU5mJAyNoWacyNQJsVA8o2cRlU6dOJTMzk5tuuon169dz6NAhoqOj6dWrF3Xq1GHixImusidOnMgxLBITE0NiYiIff/wxI0aMYMWKFezbt4+4uDiee+454uLi6NevX7HWt169esTHxxMTE8PJkydz5GMRlxYVEcrkgRG4f+YkE8DD9qcYb7+XTNNCb20bK2zP0lvdkuf+MhRTNuSdXGpyj7aCH31epLW6n2TTn5FZ/2Ws4wEyyfke4tyjRebyCJFThRx8DK3sx6/P98BHy843MaRDPbJ0A5ulZFN9X3311WzdupXx48czaNAgTp06RUhICAMGDGDcuHGuHB+QnW10zpw5Oe4/btw4+vfvz6+//srw4cM5evQoAQEBtGjRgkWLFtGtW7dire+tt95KdHQ0PXr04MyZM7LU9jIM7lif/m3qEHPwDL/tO8nUNfsAhZl6FJuNprxlnUYz9RCf+LzDAr0rE+z3kUJ2rhXnoMuz3+4g0NdKp6tkVURp0g2TTftO8ey3O1x/i1ok8Yb1Y1fSsPV6S561P0IiObM3P9ilATc2D5EltEIUoEIGH0COQENRlBIPPJzq16/PjBkzLlpm7dq1Fz3/9ddfX/bjd+/evcDlvMOGDcsRXNhsNr799luZD3KF/H0sdG4UTMerqrNg+xHXN+hdZgP6Z73KaMsCHtWWcrv2K53VXbxof5A1RhvX/VMyHQz97HdZEVOKcufwUDAYqv3Ms5a5BCnnyDCtTHYM4Su9F2auDuQX+zTlkW4NPVFtIcqNCht8CFHaNFVhXL/mDJ91IV19FlZed9zFKr0tb1mnEa4eY4bPG/yod2CC/T5XZlTIHoYZPmu7JKcqAc5MpcdTMzhwMp13V/3t6u1oqBzhNeuntFf/BiDGaMjT9uHsMy9MTu90VTWOp2SSkuGgX+vaHvgNhChfJPjwMn369OGXX37J95xz/xbhOVERoUy/uy3jl8SRmHJhDs12szE3Z01mtGUBD2rL6Ktt5np1J687BjFbvxED1fVh+M6qC1sVSG/IlSsoU6mNLIZrS3nMshib4iDNtPGGYxBf6b0xcvV2PHFDYzpdVa1Uhm+F8AYSfHiZzz77LM9qGif3OSXCc5xp2TfHJ5GYfI5XftxNUloW5/BlsmMoi/UuTLJ+Tmt1H69YZ3K7tp7/2e9lu9k4z7Wck1Ir8j4xV8I5mTTnQKTJTeoWXrLMJkzNTvr3s96Gsfb7OUrepfDOZbSlOXwrRHknwYeXqVOnZPKUiOKlqQqRDbMnKfr5aDmGYnaZDbgtawJDzs8xaK3u53vbeJbokUyx38URarjKOj80n1+wUyalFkF+k0kBmigHednyNV20OAASzGpMtA/lB6MTF3Ii5yTLaIUoOgk+hPCw/IZiDFRm6b1Yrrfnact8/qOt41ZtI73VrXym38w0x62uFO0AZ87ZZVJqIeU3zFKLJP5rWchd2mo0xSTTtPKx3pdpjls5h2++1wkJsjH+1hbS1kJcBgk+hCgDXDvk7jvFyDnbOXPODsAJqvC84xG+0nsz1jKLSG0XoyyLGaL9zMeOfnyl98rx4SjDMHldbDJpVVIYYVnKvdoKfJXsNv9R78Bkx1AOmzXyvV4VPysfDm0rvUxCXAEJPoQoI9wzo7oPw0D2UMxg+xh661t5zjKXhmoCL1i/4SHLj0x33Mos/UYy8XF9qD797Z+8ueJvhnVuwOAO9SrUh2TuYOObzQdJTMk5mbQyZ3nAEs2D2k8EKNnnNhtNeMM+iC1m04te/7XbW9KlUeG2QRBC5E+CDyHKmIJWxIDCCqM9P2e1ZYC6gf9avqe+epyx1lk8avmBmY6bmKXfQAoBpGXq/HP8LC8tiuXDNf9UmKGYglauOIVyigctPzFYW00lJbttY40GvOEYxDqjFQXN6wAZZhGiOEnwIUQZ5L4ixvkN/p1V2XkmdDQWGNezOKszt2u/8LhlIXWVkzxrncdIyyLm6935XO/DYbMmAAnn84PUCLDxxI1XM7hDPQDXtb0lX0j+K1eyNVYO8YjlR/qrG7Aq2fsV7Tbq8YFjAMuMDnkShTlJThUhSoYEH2WIoigsXLiQAQMG5Hv+wIEDhIeH88cff9C6detSrVt+GjRowOjRoxk9erSnq+KV3FfEANQM9GHMolicm906sDBP78H3+nXcom7kEcuPNFMPcr9lOfdqK1hlXMs3eg/WG9dgoHLibCYvLYrlzRV/AXAm3e66dkiQjcEd6pW7D1rdMPlm80FmbjjA8dSMHIGHD3b6qL8z1PIzHdS/XMc36s2Zrve7aE+H9HIIUbIk+Cglw4YN48svvwRA0zRq165N3759mTRpElWrZmexTEhIcP2/rOjZsyfr1q3Lc9xut7NlyxYqVarkOnap4ElcGec+MdsPnGbUN3+4JqXasbDQuI6FWV3pqsbysPYj3bQd3KRt5SZtK0fM6sx3dGe+3p0EqucIOpwSUzLLfPIy97kcNQN9OZ2WxSs/5h1iuVo5zO3aeu7U1lFdSQXAYaqsMNrxseMW/jQbFfgYMplUiNIhwUcpioqKYsaMGTgcDnbt2sUDDzzAmTNn+OabbwAICQnxcA3z9/DDD/O///0vxzGLxUKNGvmvBhAlx9/HQtfGNfKdlAoKvxot+dVoSSPHYQZra7hdW08d5RRPWhfwhOV7tphN+EHvxDK9IyepXODjFDaVe+6A4FJlgv0L/5ZTmImjTnWVE/RTN3KrtoFm6iHX8aNmNb5x9GSe3iNHqvqCyGRSIUqHBB+lyGazuQKMunXrMmjQIGbOnOk6n7vnYPPmzTz66KPs3r2biIgIxowZk+eaS5Ys4emnn+bw4cN06tTJtTnc6dOnqVKlCgC//fYbzz//PFu2bCE4OJiBAwcyefLkHL0WF+Pv759vYOQ+7NKgQQMABg4cCGRvoHfgwIHCNYwosoInpWb7x6zLK457eN0xiJvULQzW1hCp7aKjsoeO6h7GW77kd6MZy432rDWu4V8z59+3MKnc85vcGVrZl7F9m1G1ki3foMGmmbzeAVbEJVIt0D9H0BJ3NJnJP+3hhZubcvTMuYtOHAWTFsq/9FD/4AbtD9qo/7jOZJka64zWzNO7s8Zojc6ls47KMIsQpav8Bx+mCfb0wpU1jOyyWRoUx06tVn9QLq9rdv/+/URHR2O1WvM9n5aWxi233ELPnj2ZNWsW8fHxPPHEEznKHDhwgDvuuIMnnniChx56iD/++INnnnkmR5mdO3dy00038corr/D5559z4sQJRo0axahRoy65u25RbNmyhZo1azJjxgyioqLQNEkzXdIKStPuLhMflhhdWGJ0obb9JDdrv3OLtonW6j46a7vorO0C4IBRi3VGK9YbrdhqNCGZgDyP5+wNGdwhjDV/nSAxn8AgITmDx+b8UWCdnfNVnvr2TzL1C6+dkCAbIZV9iTmUzMuL4/jz0Jk8E0dDOUV7dQ+R6i56aDGEKKfdrquw0WjOEqMzy/T2pORT/9xkMqkQnlP+gw97Okwq3C6SKlClOB/7xaPgU7jeA4AffviBgIAAdF0nIyP7jfvtt9/Ot+zs2bPRdZ0vvvgCf39/WrRoweHDhxkxYoSrzPTp02nSpAlvvPEGAE2aNCE2NpaJEye6yrzxxhsMGTLENSn06quv5v3336dbt25MmzYNX9/8sze6++ijj/jss89ctx999FHeeuutHGWcQzBVqlQps8NH3ih3mvaCVnsAHCWYz/S+fKb3pa5ynJvV3+mu/kk79S8aqMdooK7kPlYCsM8IZbtxNdvNq9lhNOQfszaZ+ADwzeZDBTzCpRkFVC4xJdPVgxNz6AxBpNFYOUQz9SBt1b20V/+irnIyx33STRsbjAhWG635WW9bqGEVkF4OIcqC8h98lCM9evRg2rRppKen89lnn/H333/z+OOP51t29+7dXHPNNfj7+7uORUZG5ijz119/0b59+xzHOnTokOP2tm3b+Oeff5g9e7brmGmaGIZBfHw8zZo1u2S9hw4dmmPIxzmcI8qWqIhQpt3d9hLDFdkOmzX5RO/HJ3o/KnGOSHUX3dQ/6azG0VBNcP3cyXogu2fhoFmTvWZd/jFrc8QMJsGsRuL5nyQCC1yu6k43YcsJUA0HwaQTopyijnKSuspJ6ignqa8co4l6KE+gAdmTRuPMBmwxmrDeaMXvRjNXQHQxqgIPdGlAy7pVpJdDiDKi/AcfVv/sHohCMAyDlNRUggIDUYtr2KUIKlWqRKNG2TPt33//fXr06MGECRN45ZVX8pQ1zYK+v+Yso+Qa9sl9P8MwePTRR/nvf/+b5/716tUrVL0rV67sqrco29yHYlbuSuSLDQcueZ80/FhlXMsq41oAqpBKG/Uf2qp7aavspZn6L9WUszRQjtGAY/RiW55rGKZCGr7ZP2b2vwYKKiYKJiomVhwEKelUPpTGq7asPNfI7bAZzF9GGDvNcDYbTYkxGpFewD4ruVX1t/ByvwhqBdpoXa8K/j7l/61OCG9S/l+RilL4oQ/DAKueXb44go8rNG7cOPr06cOIESOoXTvn0FHz5s35+uuvOXfuHH5+2RuIbdq0KUeZpk2b8tNPP+U4tnXr1hy327ZtS1xcXKkED1arFV3XS/xxxMU5h2IiG1anQ3i1PD0hVfytmKZJ8jlHvvc/QyBrjDasMdq4jlUnmavVI1ytHKahcpRQJYkQJYlQJYlgklEVk0DOEci5iyUJzeO4WYUjZjBHzGAOm8EcNmvwlxHGX2YYKRR+SNOdAky+rZUMqwhRhpX/4KMc6969Oy1atGDSpElMnTo1x7khQ4YwZswYHnzwQV566SUOHDjAm2++maPMo48+yttvv81zzz3Hgw8+SExMjGv1jLNH5LnnnqNTp06MHDmShx9+mEqVKrF7925WrlzJBx98UKy/T4MGDfj555/p0qULNputzOUsqYhyZ0p1DjvAhQynOw8n88WG+ALnYwCcojKnjMpsonmecxYcVOUslZRzVCKDADKopJwDwEABFAwU7FhIMSuRjD8pZiXO4odRiKGaoiiL+UmEEHlJ8OFhTz31FPfffz/PPfdcjuMBAQEsXbqU4cOH06ZNG5o3b86UKVO4/fbbXWXCw8P57rvvePrpp3nvvfeIjIxkzJgxjBgxApvNBkCrVq1Yt24dY8aM4brrrsM0TRo2bMigQYOK/Xd56623eOqpp/j000+pU6eOLLUtI3JnSnVyHuvfug5P9W5MzMEznDibSfyJNN79eW+e8gVxYOEEVThhVrlwsMBAxqRIXSMXERJk4+VbmruW9cp8DiHKDwk+Sol7Pg93Q4YMYciQIUDe+RqdOnUiJiYmx7HcZW699VZuvfVW1+2JEydSt27dHKtY2rdvz4oVKy6r3qtXry5wfkzu4KJfv37069fvsh5HeJa/j4XObsm1moYGFphD5EoogL/FJM3hvFV4tQJ9GNKxviyPFcILSPBRzn300Ue0b9+e6tWrs2HDBt544w1GjRrl6WqJci6/je3eXfV3gR0aD3ZpwA3NanEiNZPJy3bnCFpqBlgZ0rE+4TUCqO6ncWzX77z3lz8HzxQ86VQhu49EcnEI4Z0k+Cjn9u7dy6uvvkpSUhL16tXj6aef5oUXXijUfX/55Rf69OlT4PnDhw8XVzVFOZR7uKZJSEC+GU1zz7G45ZraBaZct9vt/LQHVjzZjT8OpxaYOj1E5m4I4dUk+Cjn3nnnHd55553Lum+7du3yDOs4GYZxBbUS3qigyau5eyMKmmNysTKjeja65HWFEN5Dgo8KzM/Pr8AluIZhkJKSUso1EmVdYQKLsnRdIUTZ5PlkF0IIIYSoUMpl8CFDAuWb/P2EEKJiK1fDLj4+PqiqytGjR6lRowY+Pj550otfjGEYZGVlkZGRUTzp1b1YSbSVaZpkZWVx4sQJVFXFx+fS+3IIIYTwPuUq+FBVlfDwcBISEjh6tHD7ubgzTdOVrrwoQUtFVJJt5e/vT7169SQAFEKICqpcBR+Q3ftRr149HA5HkfcRsdvtrF+/nuuvvx6r1VpCNfQOJdVWmqZhsVgk+BNCiAqs3AUfkL1vidVqLfKHoqZpOBwOfH19Jfi4BGkrIYQQJUX6vYUQQghRqiT4EEIIIUSpkuBDCCGEEKWqzM35cO7aWhLZNe12O+np6aSkpMg8hkuQtio8aavCk7YqGmmvwpO2KrySaivn53bu3dfzU+aCj9TUVADCwsI8XBMhhBBCFFVqaiqVK1e+aBnFLEyIUooMw+Do0aMEBgYW+3LMlJQUwsLCOHToEEFBQcV6bW8jbVV40laFJ21VNNJehSdtVXgl1VamaZKamkrt2rUvmcepzPV8qKpK3bp1S/QxgoKC5MlZSNJWhSdtVXjSVkUj7VV40laFVxJtdakeDyeZcCqEEEKIUiXBhxBCCCFKVYUKPmw2G+PGjcNms3m6KmWetFXhSVsVnrRV0Uh7FZ60VeGVhbYqcxNOhRBCCOHdKlTPhxBCCCE8T4IPIYQQQpQqCT6EEEIIUaok+BBCCCFEqaqwwcett95KvXr18PX1JTQ0lHvuuYejR496ulplzoEDB3jwwQcJDw/Hz8+Phg0bMm7cOLKysjxdtTJr4sSJdO7cGX9/f6pUqeLp6pQpH330EeHh4fj6+nLttdfyyy+/eLpKZdL69evp168ftWvXRlEUFi1a5OkqlUmTJ0+mffv2BAYGUrNmTQYMGMBff/3l6WqVWdOmTaNVq1au5GKRkZEsW7bMI3WpsMFHjx49mD9/Pn/99RcLFixg37593HHHHZ6uVpmzZ88eDMPg448/Ji4ujnfeeYfp06fz4osverpqZVZWVhZ33nknI0aM8HRVypR58+YxevRoxowZwx9//MF1111Hnz59OHjwoKerVuakpaVxzTXXMHXqVE9XpUxbt24dI0eOZNOmTaxcuRKHw0Hv3r1JS0vzdNXKpLp16/Laa6+xdetWtm7dSs+ePenfvz9xcXGlXxlTmKZpmosXLzYVRTGzsrI8XZUy7/XXXzfDw8M9XY0yb8aMGWblypU9XY0yo0OHDubw4cNzHGvatKn5/PPPe6hG5QNgLly40NPVKBeOHz9uAua6des8XZVyo2rVquZnn31W6o9bYXs+3CUlJTF79mw6d+4sWzEXQnJyMtWqVfN0NUQ5kpWVxbZt2+jdu3eO47179+a3337zUK2Et0lOTgaQ96dC0HWduXPnkpaWRmRkZKk/foUOPp577jkqVapE9erVOXjwIIsXL/Z0lcq8ffv28cEHHzB8+HBPV0WUIydPnkTXdWrVqpXjeK1atUhMTPRQrYQ3MU2Tp556iq5duxIREeHp6pRZO3fuJCAgAJvNxvDhw1m4cCHNmzcv9Xp4VfAxfvx4FEW56M/WrVtd5f/v//6PP/74gxUrVqBpGvfeey9mBUn4WtS2Ajh69ChRUVHceeedPPTQQx6quWdcTnuJvBRFyXHbNM08x4S4HKNGjWLHjh188803nq5KmdakSRNiYmLYtGkTI0aM4L777mPXrl2lXg9LqT9iCRo1ahR33XXXRcs0aNDA9f/g4GCCg4Np3LgxzZo1IywsjE2bNnmkC6q0FbWtjh49So8ePYiMjOSTTz4p4dqVPUVtL5FTcHAwmqbl6eU4fvx4nt4QIYrq8ccfZ8mSJaxfv566det6ujplmo+PD40aNQKgXbt2bNmyhffee4+PP/64VOvhVcGHM5i4HM4ej8zMzOKsUplVlLY6cuQIPXr04Nprr2XGjBmoqld1mBXKlTy3RPYb3rXXXsvKlSsZOHCg6/jKlSvp37+/B2smyjPTNHn88cdZuHAha9euJTw83NNVKndM0/TI555XBR+FtXnzZjZv3kzXrl2pWrUq+/fv5+WXX6Zhw4YVotejKI4ePUr37t2pV68eb775JidOnHCdCwkJ8WDNyq6DBw+SlJTEwYMH0XWdmJgYABo1akRAQIBnK+dBTz31FPfccw/t2rVz9aAdPHhQ5g/l4+zZs/zzzz+u2/Hx8cTExFCtWjXq1avnwZqVLSNHjmTOnDksXryYwMBAV89a5cqV8fPz83Dtyp4XX3yRPn36EBYWRmpqKnPnzmXt2rVER0eXfmVKfX1NGbBjxw6zR48eZrVq1UybzWY2aNDAHD58uHn48GFPV63MmTFjhgnk+yPyd9999+XbXmvWrPF01Tzuww8/NOvXr2/6+PiYbdu2lSWRBVizZk2+z6H77rvP01UrUwp6b5oxY4anq1YmPfDAA67XX40aNcwbbrjBXLFihUfqophmBZlhKYQQQogyoeIN3gshhBDCoyT4EEIIIUSpkuBDCCGEEKVKgg8hhBBClCoJPoQQQghRqiT4EEIIIUSpkuBDCCGEEKVKgg8hhBBClCoJPoQQQghRqiT4EEIIIUSpkuBDCCGEEKVKgg8hhBBClKr/B2kH3YcG+45DAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "%matplotlib inline\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from sklearn.linear_model import LinearRegression\n", + "from sklearn.preprocessing import PolynomialFeatures\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn.preprocessing import StandardScaler\n", + "\n", + "def MSE(y_data,y_model):\n", + " n = np.size(y_model)\n", + " return np.sum((y_data-y_model)**2)/n\n", + "\n", + "def OLS_fit_beta(X, y):\n", + " return np.linalg.pinv(X.T @ X) @ X.T @ y\n", + "\n", + "def Ridge_fit_beta(X, y,L,d):\n", + " I = np.eye(d,d)\n", + " return np.linalg.pinv(X.T @ X + L*I) @ X.T @ y\n", + "\n", + "\n", + "np.random.seed(2018)\n", + "n = 100\n", + "d = 3\n", + "Lambda = 0.01\n", + "true_beta = [2, 0.5, 3.7]\n", + "\n", + "# Make data set.\n", + "x = np.linspace(-3, 3, n)\n", + "y_real = 2 + 0.5*x + 3.7*x**2\n", + "\n", + "y = np.sum(\n", + " np.asarray([x ** p * b for p, b in enumerate(true_beta)]), \n", + " axis=0) + 0.1 * np.random.normal(size=len(x))\n", + "\n", + "\n", + "#Design matrix X includes the intercept and scaling is made\n", + "X = np.zeros((len(x), d))\n", + "for p in range(d): \n", + " X[:, p] = x ** (p) \n", + "\n", + "\n", + "#Split data\n", + "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n", + "\n", + "\n", + "#Calculate beta, own code\n", + "beta_OLS = OLS_fit_beta(X_train, y_train)\n", + "beta_Ridge = Ridge_fit_beta(X_train, y_train,Lambda,d)\n", + "print(beta_OLS)\n", + "print(beta_Ridge)\n", + "\n", + "#predict value\n", + "ytilde_test_OLS = X_test @ beta_OLS\n", + "ytilde_test_Ridge = X_test @ beta_Ridge\n", + "\n", + "\n", + "#Calculate MSE\n", + "\n", + "print(\" \")\n", + "print(\"test MSE of OLS:\")\n", + "print(MSE(y_test,ytilde_test_OLS))\n", + "print(\" \")\n", + "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,y_real,label='no noise')\n", + "plt.plot(x, X @ beta_OLS,'*', label=\"OLS_Fit\")\n", + "plt.plot(x, X @ beta_Ridge, label=\"Ridge_Fit\")\n", + "plt.grid()\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "f1352336", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[4.97815077e-01 3.70043176e+00 8.59040228e-04]\n", + "[4.97813734e-01 3.70043111e+00 8.59257518e-04]\n", + "2.003056414623952\n", + "2.0030583797561565\n", + " \n", + "test MSE of OLS:\n", + "0.008737853210440363\n", + " \n", + "test MSE of Ridge\n", + "0.008737821720695461\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAACGGklEQVR4nO3dd3gUVdvA4d/M7GZDSEKHhB66EEIRkCJdkCqKHRvqiw1QRF8VUYEPsHcR7IiF8ipFpQRCBwHp0nvooQYSSEh2d2a+P0I2mWyABJJNe+7r4jJz5szkeAjJk1Oeo5imaSKEEEII4SNqXjdACCGEEEWLBB9CCCGE8CkJPoQQQgjhUxJ8CCGEEMKnJPgQQgghhE9J8CGEEEIIn5LgQwghhBA+JcGHEEIIIXzKltcNyMgwDI4fP05QUBCKouR1c4QQQgiRBaZpcuHCBSpWrIiqXn1sI98FH8ePH6dKlSp53QwhhBBCXIcjR45QuXLlq9bJd8FHUFAQkNL44ODgHH23y+ViwYIFdO3aFbvdnqPvLmykr7JO+irrpK+yR/or66Svsi63+io+Pp4qVap4fo5fTb4LPlKnWoKDg3Ml+AgICCA4OFi+OK9B+irrpK+yTvoqe6S/sk76Kutyu6+ysmRCFpwKIYQQwqck+BBCCCGET0nwIYQQQgifkuBDCCGEED4lwYcQQgghfEqCDyGEEEL4lAQfQgghhPApCT6EEEII4VMSfAghhBBFhG6YrI2OBWBtdCy6YeZJOyT4EEIIIYqAyG0x3PreYp6YtA6AJyat49b3FhO5LcbnbZHgQwghhCjkIrfF8OwvG4mJS7KUn4hL4tlfNvo8AJHgQwghhCjEdMPktRlbSZ1gMUwYt13FMPGUjfprh0+nYCT4EEIIIQqxcYv3cj7R5bl2GbA3PiX4gJQAJCYuybMWxBck+BBCCCEKKd0wmfj3QUvZY9oCGioH0E3rSMepC9Ypmdxk89lnEkIIIYRPrY2O5fyltFGPipxhlH0ShqnQKvlzTlLGc698kL/P2iUjH0IIIUQhlXE0o5e2GoC1Zj1L4BHgp9EirLTP2iXBhxBCCFFIlQ10WK77aKsA+FNvbSl/qm0NNFXxWbsk+BBCCCEKqTUHzng+rqkco4F6CJepMVdvYal3c9VSPm2XBB9CCCFEIaQbJl8vO+C5vuPylMtyI4LzBFnq/nPwrE/bJsGHEEIIUQitOXAWp566o8XkDvVvAP7UW2VS23dTLiDBhxBCCFEoPT9lo+fjhko0YepJLpl+RBnNvOq2qlnGqyw3SfAhhBBCFDKjZ2/jbELaFts+Wsqox0KjKYlYt9T621Va1pDgQwghhBDXyek2+H7lIc+1ikEvbQ2QcZdLypTMM+1q+nSnS0qbhBBCCFFo/Lz6oOW6hbqLEOUccWYAy4xG6e4o2FSFwZ1r+7R9IMGHEEIIUajM2HTUcp260HSe3gIndk+5hsmdjSv6fNQDJPgQQgghCo25W46z/fgFz7UdNz20tQD8aVgTi9k0eLtvhE/b5/ncefJZhRBCCJGjdMPkjT+2WcraqlsoqSRwyizJGqO+p7xHZZ1K1cLws+XNGIQEH0IIIUQhsDY6lth0O1wgLZ36bL0lRrrJjlvKm/TrVs+n7UsvWyHPhAkTiIiIIDg4mODgYFq1asW8efM89/v374+iKJY/LVu2zPFGCyGEEMIq4yFyxUiii7oB8D7LpaT1yBefy9bIR+XKlXn33XepVasWAJMmTaJPnz5s2rSJBg0aANCtWzcmTpzoecbPzy8HmyuEEEKIzJQPsubv6KquJ0BJ5qBRgc1mzTxqVeayFXz07t3bcj127FgmTJjAmjVrPMGHw+EgJCQk51oohBBCiGs6m2Hk487LicX+MNqQPn26TTHJa9e95kPXdX777TcSEhJo1SotT/zSpUspX748JUuWpH379owdO5by5ctf8T3JyckkJyd7ruPj4wFwuVy4XK4rPXZdUt+X0+8tjKSvsk76Kuukr7JH+ivrinpf6YbJmDnb0BQT3YSyxNNW3QrALL1NupomgfaU4CO3fsZmhWKaZrZCoK1bt9KqVSuSkpIIDAxk8uTJ9OjRA4Bp06YRGBhItWrViI6O5s0338TtdrNhwwYcjswnmEaOHMmoUaO8yidPnkxAQEB2miaEEEIUWbHJMGpjyphCfy2Skfaf2GzU5E7naEu9z1q5c+XzJyYm0q9fP+Li4ggODr5q3WwHH06nk8OHD3P+/HmmT5/Od999x7Jly6hfv75X3ZiYGKpVq8bUqVPp27dvpu/LbOSjSpUqnDlz5pqNzy6Xy0VUVBRdunTBbrdf+4EiTPoq66Svsk76Knukv7KuqPfV+/N28e2qw57rWX5v0Fg9wAjXY0zSb/eUa5j8+FgTYvesz/G+io+Pp2zZslkKPrI97eLn5+dZcNqsWTPWrVvHZ599xtdff+1VNzQ0lGrVqrF3794rvs/hcGQ6KmK323PtCyg3313YSF9lnfRV1klfZY/0V9YVxb7SDZPJ69OymtZQjtNYPYDbVJmtW3eclijuoHmNcszfk/N9lZ133XB2EdM0LSMX6Z09e5YjR44QGhp6o59GCCGEEJkYt3gvCU7Dc516gu1yI4KzlLDUfbRltTxJp55RtkY+Xn/9dbp3706VKlW4cOECU6dOZenSpURGRnLx4kVGjhzJ3XffTWhoKAcPHuT111+nbNmy3HXXXbnVfiGEEKLI0g2TiX8fTFdicufls1xm6bd61Q8rV9w3DbuGbAUfJ0+e5JFHHiEmJoYSJUoQERFBZGQkXbp04dKlS2zdupWffvqJ8+fPExoaSseOHZk2bRpBQUG51X4hhBCiyFobHcv5S2m7TJoqe6mmnuKi6c8C42av+hlzgeSVbAUf33///RXvFStWjPnz599wg4QQQgiRNSfiM8/tMd9oRhLW9ZSli/vRIqw0hp47u12yo8icaut0G0xaFQ3AsOlbWLHnNLqR94lWhBBCiOv1997Tno/tuOmlrQYyn3K5s3HFfLHeA4rIwXLvzN3BN8ujsakm1QNVDl6M4ffNJyjup/HRfY3oFi4LYoUQQhQsumHy15bjnut26r+UVi5yyizJKqOBV/0u9fNP9vFCP/LxztwdfL08GhPQDdgbr6JfXhSc4NR55peNRG6LydM2CiGEENm15sBZkt1pI/h3XZ5y+VNvhY5mqVvm8pRLflGogw+n2+Dr5dGe69rKYe5Wl2NkqDfqrx0yBSOEEKJAWbX/jOfjQBK57fIJtjMzmXIZ3Sc830y5QCEPPn5efdDzcU3lGPMdw3jb/h0lSLDUi4lLYm10rI9bJ4QQQly/lenWe3TT1uGvuNhrVGK7Wd1Sr1a54vSIyF/LCwp18HEoNtHz8X6zEjuNqjgUN720NV51T2U4DVAIIYTIr3TD5N+j8Z7rO9WVAMzUrSfYAnRtUMGXTcuSQh18VCttPZhuut4WgL7aCq+6B88kepUJIYQQ+dEXi/Z4Pq5ALK3VHQD8abTxqtumZjmftSurCnXw8Uir6pb47w+9NbqpcLO6l+qKdZHpxFXRsu5DCCFEvqcbJuOX7vNc36GtQlVM1hp1OWpaAw1/m0rLmmV83cRrKtTBh59N5T9tq3uuT1OKFUYEAHdpKy11zye6WLP/rC+bJ4QQQmTbmgNncepp16m7XDLL7dGxXrl8tdA0VaEOPgA61bPua06derlbW4GSYd/LL/8c9FWzhBBCiOvy7fIDno/rKEeorx7CaWrM0W/xqvvwLdV92LKsK/TBR8aFpAuMZsSbxaisnKGFsttyb/HOUzL1IoQQIt/SDZOle9J2uaSuYVxiNCGOQEvd/DrlAkUg+Mh4iE4yfsy9HB1mXHiarJsy9SKEECLfSv8zSkP3LCFIHdVP75n2NfPllAsUgeCjRVhpQoKth+vMuPyX1EP7B3+SLfc+XLDLZ20TQgghsuPVGf96Pr5V3UYF5TyxZiBLjCaWenZVYXDn2r5uXpYV+uBDUxVG3pGS414lZUplnVmXw0Y5gpRLdL2cES7VpiNxON0Zc6AKIYQQeWvsnO0cPZe2lOBubTkAf+htcGU4qu22+hXy7agHFIHgA6BbeChDOtfCdvn/1kRlppG68HS5V/3XZ2zxZfOEEEKIq3K6Db5dcdBzHUQiXdX1QOZTLg+3rOarpl2XIhF8AAzuXAeHlhYFzri8JelWdSvlOWep+8e/x2XhqRBCiHwj/XEhAD21NfgrLnYbldlmhlnu+dtVWtbInwtNUxWZ4ENTFXo0DOXu6imbow+ZIawz6qApJn0u75FO5ZKFp0IIIfKR9MeFQNqofcqoh3V6pWPd/JnbI70iE3wAvNW7Ae1C00Y0ZqTL+QHWkQ7J+SGEECK/2Hsy7RyXasoJmqt70E0l0xNs82tuj/SKVPDhd3nRh3I50Jij30KyaaeeeoT6yiFL3YU7JOeHEEKIvOd0G6w+kLY8IDVNxAojgtOUstQtWcyeb3N7pFekgo9U5uUhqngCiTKaAqmjH2lchskXi/b6vG1CCCFEeuk3QSgYnp9XMzJZaPp4m7B8P+UCRTT4sClpIxrT9XYA3KH9jQ23pd53Kw/I6IcQQog8oxsmMzce81zfou6isnKGeLMY841mXvWrlw3wKsuPimTwoaX7v15hNOS0GUw5JZ626lZLvYvJOmujY33cOiGEECLFF4v2kO4MOe5WUxaaztFbkoyfV/2MWb3zqyIZfKTnxsafehvAe+oFvM+GEUIIIXxBN0y+XRntuS5GEt21tUDmuT1CS/jTIqy0z9p3I4pk8FEhKPN0613UDQRz0XIv+nSCz9olhBBCpFobHUtCctq4Rzd1HYFKEgeNCqw363rVH9G7foFY7wFFNPgY1uMmy/V2sxo7jSo4FBe9tTWWe9+skHUfQgghfG/hjhOW69TcHjMyye0xpHNtuoWH+qppN6xIBh+33VSBJ9tUT1ei8Pvlhaf3asssdROdOuMW7/Nd44QQQhR5umEyeW1aCohQztJa3QHADMM65VLcT8vXh8hlpkgGHwC31Q+xXM/Sb8VlajRW91NbOWq5N27JXhn9EEII4TNrDpzlkivt585d2gpUxWSNcRNHzXKWugPa1igw0y2pimzw0SKsNEH+muf6LCVYfPlI4oyjHy7dZNXeMz5tnxBCiKJrteWID9OzISKzhaZh5Yr7qFU5p8gGH5qqcE/Typay3/T2QEqEmTHnx/RN1tEQIYQQIrccOH3B83ETZR811RgSTQdz9Vu86haU7bXpFdngA6BrA+vinKVGI06bJSinxNNR3Wy5l+jUEUIIIXKbbphEbjvpuU5daBppNCeBYpa6pYv7FZjttekV6eCjRVhpAh1pUy9ubMy4fEhPxqmX0/GS70MIIUTue37KBozLH/uTzB3aKiBtdD69R1tWK3DrPaCIBx+aqvDuXQ0tZal/uZ3UTZQlzlO+6Wgcc7fE+LR9Qgghihan22DO1rRRj27qOoKVSxwxyrHGuMmrfkFc7wFFPPgA6NW4Ek2rlvBc7zMrs8mohU0xuFNbaan74rRNsutFCCFErvl59UHL9X3aUiDlF2Mzkx/ZBXG9B0jwAUDrmmUt1/+7PPqR8peeFmwk63LSrRBCiNwzdd1hz8dVlJO01nZgmGm5qNIrqOs9QIKPy6zzZbP1Vlwy/aijHqORst9y75sV+2X0QwghRI6bu+U4e0+lHelxz+WFpiuNcI5T1qv+mD7hBXK9B0jwAUCrmmUs1xcIINJoDsB9XhlPDTnpVgghRI7SDZPhs7Z5rlUMT/DxP72DV/0nbw2jR0TBSaeekQQfQMsaZSjuZ+2K1L/s3toq/Em23IvKkG9fCCGEuBFro2M5l+jyXLdRt1FJOct5szhRxs2WutVLF+PNXvV93cQcla3gY8KECURERBAcHExwcDCtWrVi3rx5nvumaTJy5EgqVqxIsWLF6NChA9u3b8/xRuc0TVX44J5GlrI1xk0cMcoRrFzidnWd5d7UdUdk6kUIIUSOyXiIXOpC01l6G5Lxs9zr1aiij1qVe7IVfFSuXJl3332X9evXs379ejp16kSfPn08Acb777/Pxx9/zLhx41i3bh0hISF06dKFCxcuXOPNea9HREV6NqzguTZR5bA5IYQQuU43TH7bkJZFuyQX6KquB+C3TKZcWtXwXv9R0GQr+Ojduzc9evSgTp061KlTh7FjxxIYGMiaNWswTZNPP/2U4cOH07dvX8LDw5k0aRKJiYlMnjw5t9qfozJmPJ1upAQfrdUdVFZOW+59v+KAjH4IIYS4YWujY4lPSjvSo4+2CofiZrtRje1mdUvdQIeNlhnWKRZEtut9UNd1fvvtNxISEmjVqhXR0dGcOHGCrl27euo4HA7at2/PqlWrePrppzN9T3JyMsnJaWsq4uPjAXC5XLhcrkyfuV6p77vSe8sG2HBoJm4DdFPhqFmOv/UGtNG2c7e6nM/0u9Pamexmzb5TBXab07Vcq69EGumrrJO+yh7pr6wryH114nyC5Tp1yiXjQlM/1eSBmyti6G6MGzjxI7f6KjvvU0zTzNav71u3bqVVq1YkJSURGBjI5MmT6dGjB6tWraJNmzYcO3aMihXT5qOeeuopDh06xPz58zN938iRIxk1apRX+eTJkwkICMhO03LUC6tT4rI+6ko+8xvPEaMc7ZyfWJK8fNbKfaXHhRBCiCxbGqMw66BKfeUQcxyvk2zaaJE8njgCAbilnE6/Wvl7tD0xMZF+/foRFxdHcHDwVetme+Sjbt26bN68mfPnzzN9+nQee+wxli1LWxOhKNY9x6ZpepWlN2zYMIYOHeq5jo+Pp0qVKnTt2vWajc8ul8tFVFQUXbp0wW63Z1pnwtJ9fLk0LbfHfKM58WYxqqinaanuZLXRwHPv9fU2NrzRpcDus76arPSVSCF9lXXSV9kj/ZV1BbWvnG6DJqOjcBopP0fu1ZYCsMBo5gk8AP45rbHvosrfr3W+4Z85udVXqTMXWZHt4MPPz49atWoB0KxZM9atW8dnn33Gq6++CsCJEycIDU1bO3Hq1CkqVKiQ6bsgZWrG4XB4ldvt9lz7Arrauwd2rscXS6NJzWyahIO/9NY8ZFvE/doSS/CR4IKvVxzihdtq50o784Pc/HsobKSvsk76Knukv7KuoPXVT2sOeAIPB07u1P4GMs/t0aJmOfwdfl7l1yun+yo777rhPB+maZKcnExYWBghISFERUV57jmdTpYtW0br1q1v9NP4jKYq9MmwjWmq3hGA7uo6SnDRcm/i39Gy8FQIIcR1SZ9OvYu6gZJKAsfMMvxthHvVffiW6j5sWe7KVvDx+uuvs2LFCg4ePMjWrVsZPnw4S5cu5aGHHkJRFIYMGcLbb7/NzJkz2bZtG/379ycgIIB+/frlVvtzxdt9IyzXW80wthvVcCgu+morLPfOX3JJxlMhhBDZljGdeupC09/1dhgZfjwXs6uFYpdLqmwFHydPnuSRRx6hbt26dO7cmX/++YfIyEi6dOkCwCuvvMKQIUN47rnnaNasGceOHWPBggUEBQXlSuNzi59NpWPd9PuoFabonQB4QFtC+sPmAL5Zbj3/RQghhLga3TB5Yeomz3VFznCrmpJePbND5N6/p1GhWl+YrTUf33///VXvK4rCyJEjGTly5I20KV94ql0tluw+47n+U2/NcNuv1FWP0lTZy0azjufekt2nmbslpkDn2RdCCOE7ny/ag8tIu77PthRVMVml1+eIaV0nGVYmgN6FIKtpenK2yxW0CCtNiWJpsVk8xZlr3AKkjn5YvfHHNln7IYQQ4pp0w+S7FdGeaxXDM+Uy9fIoe3o9IwpX4AESfFyRpio80SbMUjbFnbLwtJe2hkASLfdiE5yy9kMIIcQ1rY2OJcGZliWsvfovFZVYYs1A5hvNvOpnPHm9MJDg4yoGdaqNw5Y2x7berMteoxIBSjJ9tFVe9U/EXfJl84QQQhRAGQ+R66ctBmCG3tbrELlAh42WNST4KFI0VeHZ9jXTlShMvbz3+oHLXyzp/b3vjFeZEEIIkSrjIXIViKWjmrLwdEomUy7v3x1RqBaappLg4xoGd66Dn5b2Fz9Db0uyaaOhepAGSrSl7uytMbLuQwghxBVlPETuXm0ZNsVgrVGX/WYlS91eDUML7UYGCT6uQVMVBnas5bk+RzALLs/JPZhh9CPJZTBu8T6ftk8IIUTBsWB7jOdjBYP7taUATHF7j3p0aXDl7OAFnQQfWTCoU2380639SB0a66OtohhJlrqS8VQIIURmdMPk59WHPNe3qtuoop4mzgzw7KZMr3yQvy+b51MSfGSBpiq0CEtb8LPaqM8hozxByiV6av9Y6krGUyGEEJl5fsoG3Ol+N00dPZ+p3+q10LRkMTstwkr7snk+JcFHFlUrU9zzsYnKtMvnvWSW80N2vQghhEhv7pbjzNl60nNdlji6qBuAzBea3nZT+UK50DSVBB9ZVL1MgOX6N70dblOlmbqH2spRy70F263bqIQQQhRdumHyxh/bLGX3aMuwKzqbjFrsNqt6PdOmVlmvssJEgo8seqRVdcv1aUqxyGgKeI9+zNt+UtZ9CCGEAFJ2uMQmuDzXKQtNU35uTM5k1AMgpEQxn7Qtr0jwkUV+NpXeESGWsimXp176aitw4LTce2HKJoQQQohTF6wbE1qqOwlTT3LBLMZsvaVX/dAS/oV6vQdI8JEtnz7QFHu6ObjlRiOOmWUopVzkdnW9pe7srTE43UbGVwghhChiDp5JsFynLjT9Q2/NJbx3tIzoXb9Qr/cACT6yRVMVutQv77k2UPlNbw945/wAeH3GFp+1TQghRP6TcojcAc91KeK5XV0HZL7Q9Ik21ekWXjgTi6UnwUc2PXRLdcv1NHdHdFOhlbaDmsoxyz3JeCqEEEXbmgNnuZCcdohcX20FDsXNVqM6280wr/pd6od4lRVGEnxkU8uaZfC3p3VbDGVYfHnh6UPaIkvdJJchOT+EEKIIW7U//ZlfJg9eXmg6Re/sVbdMcb9Cv9YjlQQf2aSpCh/eHWEp+0W/DYC7teX4k2y5d/xcos/aJoQQIn9Zl+4X0FuUXdRSj5NgOvhTb+VVt0/jioV+rUcqCT6uQ6/GlahcKm2R0HKjIYeNcpRQEumtrbbUnbBsv6+bJ4QQIh/QDZNNh895rh+2RQHwh96GiwR41S8qUy4gwcd161g37cAfE5XJl4fQMk697DudwNwtMQghhChahkzdiOvypseyxHkWmqaOlqdXughNuYAEH9fNO+Npe5ymRmN1Pw2UaMu9YTO3ysJTIYQoQuZuOc5fW9KyXd+nLcFP0dlo1GKHWd2r/pg+4UVmygUk+LhuGTOenqUE8y6fSviwttByL04OmxNCiCIjYzp1FYN+tpR0DL+4vUc9WoWVpkdE4d9em54EH9fJz6bStEoJS9mv7pSplz7aKoKwLjRdsF2mXoQQoijImE69g7qZysoZzpmBzDG8M5re17yKL5uXL0jwcQNeur2e5XqtWY/dRmUClGTu0lZY7k3+55BMvQghRBFwIt6aTj11NPw3vT3J+HnVL+znuGRGgo8b0LJGGYrZ03ehwq+WhadpwUayDl8s2uvbBgohhPC5n1elZTStrJyig/ovkPkhckVtoWkqCT5ugKYqPNO+pqVspt6WRNNBXfUozZXdlnufLtpL5DaZfhFCiMJq7pbjbDwS77nupy1GVUyW6w05aHqv67izCOX2SE+Cjxs0qFNtHLa0L5wLBPCH3hqAh2wLveqP+muHTL8IIUQhpBsmQ6amnWjuh4v7tKUA/JrJ9looWrk90pPg4wZpqsKzGUY/UqdeuqtrKU285V5MXJLsfBFCiELoi0V7cKY7zLybuo6ySjwxZmkWXj6GI73QEv5FcsoFJPjIEYM71yHd4AfbzBpsNmrgUNzcqy3zqi8p14UQonDRDZNvV1pzPKWOfk91d0RH83pmRO/6RXLKBST4yBGaqnBnk0qWstQhtn7aIhQMy73NR8/7qmlCCCF8YG10LAnpTq+toxzhFnUXblNlqt7Rq/6QzrXpFl60cnukJ8FHDnm7r/Wwub/0VsSZAVRTT9FO3Wq5t3q/TLsIIURhcuqCdXvtQ5e310YZN3MS69RKcT+NwZ1r+6xt+ZEEHznEz6bSqkYpz3USDqbr7QB4WIuy1N13+iJOt3U0RAghRMEVffqi5+MAkuirrQQyP8dlQNsaRXa6JZUEHzlo0hPWzHU/610A6Kxuoopy0nLvga9W+axdQgghco9umIxfss9z3Uf7myDlEgeMEFYZDbzqh5Ur7svm5UsSfOQgP5tKg4pBnutoM5RlegSqYnqd97LxaJycdiuEEIXA81M2pNvlYvLI5e/3v+qdMTP5MVs+yN93jcunJPjIYX2bVLZcT9K7AnC/thR/ki33Xpm+RXJ+CCFEATZ3y3HmbE0b2W6m7Ka+eohLph+/6+296pcpohlNM5LgI4dlPO12qdGYw0Y5SioJ9NGsUy0Xk92s2X/Wh60TQgiRU3TDZOj//rWU9bctAGCW3oY4Ar2eGd0nvMiv9wAJPnJcxtNuDVR+ujz68Zi2gPTnvQD88s9BH7ZOCCFEThkydSNJ6TYPVCCWbupaACbpt3vV79UwlB4RRXd7bXoSfOSCjKfd/qa355LpR331EM0ynPeyeOcpmXoRQogCxuk2mL3lhKXsIdtCbIrBP0Y9dplVLffsKnz2YBNfNjFfy1bw8c4779C8eXOCgoIoX748d955J7t3W3+Y9u/fH0VRLH9atmx5hTcWTi1rlKG4X1rXxhHILL0NkDYklypZNxm3eB9CCCEKjp9XH7SMY/vh4kFtMQCT3F296j/csppMt6STreBj2bJlDBw4kDVr1hAVFYXb7aZr164kJCRY6nXr1o2YmBjPn7lz5+Zoo/M7TVX44J5GlrLUqZfb1XVUwJpkbMKyfTL6IYQQBciKvact1z3Ufyh3+RyXBUYzr/pdG8h0S3q27FSOjIy0XE+cOJHy5cuzYcMG2rVr5yl3OByEhBTNk/pS9YioSKs1B1l94BwAO81q/GPU4xZ1F/1si/nEfY+nbpLLYM3+s7SpXTavmiuEECKLdMNkTYYDQvvb5gPwi/s23Bl+tJYMsMsOlwyyFXxkFBcXB0Dp0tZOXbp0KeXLl6dkyZK0b9+esWPHUr58+UzfkZycTHJy2hbU+PiUU2BdLhcul+tGmucl9X05/d4rue/mymw4GIvTSBlq+8ndlVv8dtFPW8Q495240nX/5DUHaFG9xJVe5XO+7quCTPoq66Svskf6K+t82VcTlu7D0HUg5Xt7I2UfjdX9JJs2r3NcHJrJE62qYuhuDD2Tl+WB3Oqr7LxPMU3zusb7TdOkT58+nDt3jhUrVnjKp02bRmBgINWqVSM6Opo333wTt9vNhg0bcDgcXu8ZOXIko0aN8iqfPHkyAQEB19O0fOeF1SlBhg03Kx0vEKKc43nnQP402njqjGjqprR39wghhMiHUr+vA3xsH09fbSXT9Vt5yfWcp7xvdZ32oUVnSj0xMZF+/foRFxdHcHDwVeted/AxcOBA5syZw8qVK6lcufIV68XExFCtWjWmTp1K3759ve5nNvJRpUoVzpw5c83GZ5fL5SIqKoouXbpgt9tz9N2Z0Q2T1u8sIi5ZRzdTIuTB2gxesv/OBqM2dzutQdevTzTLN0Nzvu6rgkz6Kuukr7JH+ivrfNVXa6NjeWLSOpL1lO/pZYljlWMQforOHcmj2WLW9NR1aCYDO9Ti2Q41r/S6PJFbfRUfH0/ZsmWzFHxc17TL4MGD+fPPP1m+fPlVAw+A0NBQqlWrxt69ezO973A4Mh0RsdvtufYFlJvvtnwe4NE2Nfh0Udpulql6JwbbZnKzupdw5QDbzBqee7M2x9CmToVcb1d2+KqvCgPpq6yTvsoe6a+sy+2+Wrz7jCfwAHhAW4yforPJqGUJPACSdYVq5YLy7d9dTvdVdt6Vrd0upmkyaNAgZsyYweLFiwkLC7vmM2fPnuXIkSOEhhbNlb6DO9fBnm571WlKMte4BUhNOpbmt43HZNeLEELkU7phMmPTMc+1DTcP21LOcfkxk+21IOe4XEm2go+BAwfyyy+/MHnyZIKCgjhx4gQnTpzg0qVLAFy8eJGXX36Z1atXc/DgQZYuXUrv3r0pW7Ysd911V678D+R3mqrwcEtrspmfLn+R3qGtphTxlnsdPlzss7YJIYTIurXRsZxLTFtUebu6nhDlHKfNYOZd/qUyPdnlcmXZCj4mTJhAXFwcHTp0IDQ01PNn2rRpAGiaxtatW+nTpw916tThscceo06dOqxevZqgoKBrvL3wyri/e6NZmy1GGA4lLSlNqiOxSYyevcOXzRNCCJEFC3dYM5o+dnl77WS9M068pxwebx0micWuIFtrPq61NrVYsWLMnz//hhpUGLUIK02Qv8aFpNR9VgoT3d34xG8Cj9qi+FbvZdl2+/3KaF7tVg8/m2S/F0KI/EA3TCavO+y5rq8cpIW6G5ep8av7Nq/6gQ4bgzrV8mUTCxT56eYDmqpwT1PrwtzZRitOmSUJUc7RQ13j9czPqw/6qHVCCCGuZc2Bs1xyph0il7pmL9JozilKedV//+4IGfW4Cgk+fCTj1IsLGz+5uwDwhC2SjKfdLt97xldNE0IIcQ3pfyEsQxx3an8DMNHdzavu0+3C5PTaa5Dgw0dahJWmQpCfpWyy3pkk004j9QA3K3ss9/7ee1p2vgghRD6gGyYLtp/0XD+sLcShuNhs1GSjWdtS167CK91u8nUTCxwJPnxEUxVG9Qm3lMUSzEz9VgCetM2z3HObcN/Xq3zWPiGEEJl7fsoGUidcHDh52BYFwPfu7qSmWE/lMlJ2xYirk+DDh7qFh/Jkm+qWsol6ypDd7eo6KivWUxI3HDrPX/8e91XzhBBCZDB3y3HmbE0b9bhDW0U5JZ7jZmnmGS0yfebUhSRfNa/AkuDDx26rbz3td49ZheV6QzTF5NEMSccAhs/cKtMvQgiRB3TDZPisrelKTJ7QUkapf3J39Tq9NpUkFrs2CT58rEVYaUoFWL9gf7g8+vGAtpjiXLLci09yyxCeEELkgZSkYm7PdWt1OzepR0g0HUzWO2X6TGgJf0kslgUSfPiYpir0bWLddrvMaMR+I5Rg5RL3aMu9npEhPCGE8L0F22Ms16mjHr/r7Ygn0Ku+AozoXV+22GaBBB95IOPUi4nqWfvxuBaJgmG5/5Pk/BBCCJ/SDZP/bTjiuQ5TYrhN2wSkrdVLL9BhY8LDTekWLltss0KCjzzQIqw0oSWsc4LT9bbEmQFUV0/SWd1kuScLT4UQwrfWRseSkJz2i+DjWiQAC/UmRJvWAMOuKmx8s4sEHtkgwUce0FSFEb3rW8ou4c8UvTOQNrSXniw8FUII35mRbtSjBBc9U+Lf6z286g7qVFuOw8gm6a08ktm220nurrhNldbaDm5SDlnuycJTIYTwDd0w+W3jMc/1g9piApRkdhpVWW1Yf3EsZlflDJfrIMFHHsq49iOGMp5945mNfpyIu+RVJoQQImd1/ijttHEbbh6zpaRB+F73Tir2QPMqssD0OkjwkYcyS7n+g7s7cDmRDect91buk/NehBAiN42evY2DZ9N2GPZQ1xKqxHLaLMGfemuv+hnP7RJZI8FHHsos5fomszbrjTo4FDeP2eZb7s3aeEzWfQghRC5xug2+X5l+ytvkCdtcAH52d8GJ3VLfYVMlp8d1kuAjj3ULD6V7uHX65Vt3TyDl8KIA0iJwHTnvRQghcsvPGdIa3KzsobF6gGTTzq+XNwSk1ysiVKZcrpMEH/nAwy2rWa6jjJs5YIRQUkngfm2J5Z5suxVCiNxxKDbRcv2UbQ4As/Q2nKWEV/13+kb4pF2FkQQf+UDLGmXwt6f9VRionu1cT9rmoaFb6r84bZNMvwghRA6btzUto2mYEkMXdQMA3+g9veo+eWuYbK+9AdJz+YCmKjzYvIql7He9HWfMYCorZ+ih/mO55zbghSnWRGRCCCGu3+jZ2zh90em5HqDNQVVMovSm7DcrWeo2rBTMm73qZ3yFyAYJPvKJjCumk/HjJ3dXAJ6yzQasIx1ztsbgdFvTsAshhMi+jAtNyxLH3doKAL5x97LULRtg56/BbX3avsJIgo98IrNttz/rt3HJ9KOhepBW6g7LPRPvxVFCCCGyL+P30sds83EoLjYZtVhn1rXc69Goog9bVnhJ8JFPZLbt9hzB/Ka3B+ApbbbXM9PTZeATQghxfZbvOe35OIAkHtGiAPja3YuMScWqlQ7wZdMKLQk+8pFu4aGM79fEUvad3gPDVOio/Usd5Yjl3o6YeJl6EUKIG6AbJn/vT0vgeJ+2lJJKAtFGBRYYzSx1FeCRVtV92r7CSoKPfKZHREU+v6+x5/qwWYFIozmQtu0rvWEztviqaUIIUeiMW7yX1N/hNHSevHy0xfd6D4wMPyLb1yknO1xyiPRiPnRH00qUCkjLpJe64OkO9W8qYD1c7o/Nx2XbrRBCXAfdMBm3eK/nuof6D1XU05w1gzxT3uk93b6mL5tXqEnwkU91qFPO8/Fmsxb/GPXwU3Qez5By3W2YjFu8z9fNE0KIAu/5KRtweWauzcs7C+End1eSsW4ACPa3SSr1HCTBRz7Vt2lly/U3l1Ou99MWEog1C9/EVdEy+iGEENkwd8tx5mw96blupe6goXqQS6YfP+ldvOrf3bSSpFLPQRJ85FOta5XFP93c4mKjCfuMigQrl7xSrp9PdLE2OjbjK4QQQmRCN0ze+GObpezpyzsKf9Pbc45gr2fk9NqcJcFHPqWpCh/ek3ZugInKt5dT/P7HNg87bkv9zxft8Wn7hBCioFobHUtsgstzXU85TAftX3RT4bvLR1ukV7q4n0y55DAJPvKxMkH+luuZ+q2cNEsSqsRyp7bScm/1gVjZdiuEEFlwIj7Jcj3g8lqPSKM5h80KXvXH9AmXKZccJsFHPnbqgvUfiBM737lTovJntL9QsQYbD3y1ymdtE0KIgmrBtrQD5Corp+mjpnzv/Nrd26tux3rl6BEhUy45TYKPfKx8hpEPgMl6Z86bxampxtBNXWu5t/FonIx+CCHEVeiGybztaQtNn9JmY1MMVujhbDG9t9I+1Va21+YGCT7ysRZhpQkJdljKEijGJD3lwLnnbH+S8cC5SasO+qh1QghR8AyevMHzcTnOc7+2FIDxeh+vumVkrUeukeAjH9NUhZF3NPAq/9F9O4mmg3D1IO1Ua4bTH1cd8FXzhBCiQJm75Thzt6WNejxhm+c5QG61Ud+r/mhZ65FrJPjI57qFh/LVw00tRxudI5gpeicgdfQjzbHzyTL1IoQQGeiGycu//+u5DiaBh7WFAHzp7oPXAXJlAmStRy6S4KMA6BYeyo/9m1vKvnX3wGlqtFR30lSxbrO9/ZNlvmyeEELke18s2kOiM+0Xs0e0KIKUS+wyqrDIaOJV/+FbqvqyeUWOBB8FxK11yqGlC8xPUIaZelsAns0w+hF9NpG//j3uy+YJIUS+Fbkthk8XpR1D4U8yT9hSDpCb4O6NmcmPwsdah/msfUVRtoKPd955h+bNmxMUFET58uW588472b17t6WOaZqMHDmSihUrUqxYMTp06MD27dtztNFFkaYq3NWkkqXsa70XhqnQRdtIXeWw5d7QaZsl5boQosjTDZMRGbKZPqAtoYxygcNGOWYbrbye6dUwVE6vzWXZ6t1ly5YxcOBA1qxZQ1RUFG63m65du5KQkOCp8/777/Pxxx8zbtw41q1bR0hICF26dOHChQs53vii5u2+EZbrA2ZF5hotAO/RD5dh8sWivQghRFG2NjqWkxecnms7bs8Bcl/rvdHRLPXtqsJnD3pPw4icla3gIzIykv79+9OgQQMaNWrExIkTOXz4MBs2pGxdMk2TTz/9lOHDh9O3b1/Cw8OZNGkSiYmJTJ48OVf+B4oSP5tKqxqlLGUT3Cnbw3qrq6minLTeW7ZPRj+EEEVaxmSNd2orqajEcsosye96O6/6nz3QRHa4+IDtRh6Oi4sDoHTplH3Q0dHRnDhxgq5du3rqOBwO2rdvz6pVq3j66ae93pGcnExycrLnOj4+HgCXy4XL5fKqfyNS35fT7/Wl7x65maZjonDpYKCw3azOMj2C9toWntZm84b7SU/dZLfJ37tP0KpW2Wx/nsLQV74ifZV10lfZI/2VdVfqq9mb0qakVQye0f4C4Dt3d5Lx89zTMGldoxRdbipb6Ps7t76usvM+xTTN6/rV2DRN+vTpw7lz51ixYgUAq1atok2bNhw7doyKFSt66j711FMcOnSI+fPne71n5MiRjBo1yqt88uTJBAQEXE/TioTzyTBiY0rseIuyk2mO0SSbdtomf8op0kZHRjR1U9pxpbcIIUTh9sLqtN+xu6v/MMHvM+LMAFonf0ECxQCwKyZvNNEpKd8rb0hiYiL9+vUjLi6O4GDvk4HTu+6Rj0GDBrFlyxZWrlzpdU9RrENWpml6laUaNmwYQ4cO9VzHx8dTpUoVunbtes3GZ5fL5SIqKoouXbpgt9tz9N2+pBsmTUYvQMXEQOEfsx7rjDo0V/fwlG02Y9yPeOqO2mjDoZl8en9jbrvJ+8CkKyksfeUL0ldZJ32VPdJfWZexr3TDpOXbUelqmAy0/QHAj3o3T+CR6u/ECnx5180+bHHeya2vq9SZi6y4ruBj8ODB/PnnnyxfvpzKlSt7ykNCQgA4ceIEoaFpyVlOnTpFhQqZ/+BzOBw4HN7hpt1uz7V/bLn5bl+wA0+1r81nnq1jCp+7+/Kz37s8pC3iK/cdnKGEp36yrvB/c3bTNbxStucyC3pf+ZL0VdZJX2WP9FfWpfbV+IV7OJ82o09HdTPh6kESTQc/urtannGZCo+3rV3k+jinv66y865sLTg1TZNBgwYxY8YMFi9eTFiYdR90WFgYISEhREWlRZtOp5Nly5bRunXr7HwqcQ3Pd65jyce3wmjIZqMmxRQn/7HN8aofE5fE2uhY3zVQCCHyiG6Y/PB3dLoSkxdsMwD4Se/COayj6g6bKme4+Fi2go+BAwfyyy+/MHnyZIKCgjhx4gQnTpzg0qVLQMp0y5AhQ3j77beZOXMm27Zto3///gQEBNCvX79c+R8oqjRVoW3t9AtJFT5z9wXgUS2K0ngPf2Vc9S2EEIXR2uhY4i65Pdft1C00VvdzyfTjO3dPr/q9IkJlh4uPZSv4mDBhAnFxcXTo0IHQ0FDPn2nTpnnqvPLKKwwZMoTnnnuOZs2acezYMRYsWEBQUFCON76o+/qRZpbrJUZjthrVCVCSedI216t+9OkErzIhhChsFmyPSXeVNurxq97ZMiWd6p0MOZRE7sv2tEtmf/r37++poygKI0eOJCYmhqSkJJYtW0Z4eHhOt1sAxfw0Otezjn58fnn04zFtASW4aKk/Yank/RBCFG4Ld55k4qpDnus26jZuVveSZNr52t3Lq36vCMlmmhekxwu47/vfQqmAtHXDUcbN7DCqEagk8YQt0lI3WTd5YeomXzdRCCF85q1ZW9NdpY16TNE7cRprkkY/m8pnD0g207wgwUchMP6h9NMvCp+77wLgcS2SYKxTLbO3xOB0GwghRGF0/pLu+biVuoMW6m6STRtfuXt71R3YoZas9cgjEnwUAi3CSlPckfZXOd9oxm6jMsFKIo9p3ondfl590IetE0II33GZacHE89pMAKbpHTmJdTeLn01lUKdaPm2bSCPBRyGgqQr33VzFc22i8sXl0Y8nbfMIJNFSf/rGYz5tnxBC5LaFO08Smy6vRwtlJ620HThNjQnuO7zq3yE7XPKUBB+FRNcGoZbrucYt7DMqUlJJ4FEtynJvR0y8TL0IIQoN3TB5e84ORm1MW/822JYy6vGb3oEYyng9k/GUcOFbEnwUEi3CSlOyWNo/PAOVL9x3AvAf2xwCsOb4eOz7f3zZPCGEyDVro2M5ecHpuW6q7KGttg2XqTFB9x716Fi3nOxwyWPS+4WEpirc3bSypWy20YoDRgillYs8qi2w3FsdHSujH0KIQuFEfBJOI20KJXWHy3S9LUfNcl71n2pX02dtE5mT4KMQua1+iOVaR2Pc5dGPp2yzKc4ly/0Hvlrlq6YJIUSu+XnVAc/HjZV9tNe24DZVvtT7eNUN9rdJKvV8QIKPQqRFWGkqBPlZymYZt7LfCKW0cpHHNWvej41H45i7JQYhhCio5m45zsYjacdJDLX9BsBM/VaOmN4Hmo69q6EsNM0HJPgoRDRVYVQfazZZA5VP3XcD8JRtDsEZsp6+PnOrZD0VQhRIumEyJF3ixObKLtppW3GZGp/pfb3q31ytJL0bVfRlE8UVSPBRyHQLD2XcA40tZbONluwyqhCsJPKfDGe+nL/kktNuhRAF0tjZ23F6lq6ZvGz/HwD/0ztw1CxvqWtT4H9Py+nq+YUEH4VQr8aVaBmWlkbYROWTy6MfT2iRlMpw4u2IP7f5tH1CCHGjdMPkhwxnuNyi7iLZtHnWuqXXpUGITLfkIxJ8FFL3t6hmuZ5vNGebUZ1AJYmnbXMs9/acvMglp44QQhQUXyzak+7K5KXLaz0m650zzevxcMtqXmUi70jwUUiFBPtnKFH42H0PAI9qCyhLnOXu3RP+9lHLhBDixuiGydfL0na4dFA301TdxyXTj/Fu7x0ugQ4bLWt4ByQi70jwUUilJB2zW8oWG03YZNQiQEnmWduflns7Yi7IzhchRIHwxaK9XPLkKUob9Zikd+U0Jb3qv393hEy55DMSfBRSmqrweJvqGUrTRj8e1hZSAetC0xenbZKdL0KIfO2duTv4dNFez/Xt6noaqge5aPrztbuXV/3GVUrQIyLUq1zkLQk+CrFBnWoT5LD+Fa8wGrLWqItDcTHQ9oflXrJu8kK6bWtCCJGfzN1ynK+XR3uuFQxetP0OwA96N84R7PXMf7vW81n7RNZJ8FGIaarCB/c2zlCq8Mnl0Y8HtMVU4rTl7uwtMZJ2XQiR7+iGyRt/WHfm9VLXUE89QrwZwHfuHl7PBPhptKwpaz3yIwk+Crlu4aEM6VzLUrbaaMAqvT5+iu45+TG9Xp+v8FXzhBAiS9ZGxxKb4PJca+gMsU0H4Bt3T+IJ9Hrm6XY1Za1HPiXBRxEwuHMdgv01S9mH7vsAuFdbRk3lmOXenlMX+XD+Lp+1TwghruW7Ffst13eqf1NTjSHWDGSi3i3dnZR1a8XsKoM6WX/xEvmHBB9FgKYqPHlrDUvZRrMOUXpTNCVtpXh6P64+5FUmhBB5wek2WLQrbYrYDxcv2lPWenzt7k0CxTz3Un/Nev+eRjLqkY9J8FFEVC9b3KvsA/f9GKZCD20tEYr1t4pkXf7RCiHyh0e/X2O5fkhbSGXlDCfMUkzSu3rKbZhoqpzhUhBI8FFElA/KmHQM9phVmGncCsCrtqmkDlcKIUR+4XQbrIk+57kOJJFBtlkAfOq+myQcnnvP1dfx0xQ5w6UAkOCjiGgRVprQEt4ByCfue0g2bbTRttNW3Wq5F5uM5P0QQuSp12dssVz/xzaXMsoF9huh/Ka3t9yrWQKekkWmBYIEH0WEpiqM6F3fq/yoWY5f9dsAeMU2FYW0bbajNtp4dfoWr2eEEMIXdMPkj81pC+LLEMd/tJSTuT9034dO2kJ67fLI7VPtrOvbRP4kwUcR0i08lK8ebkpxh3Xnyzj3nVwwi9FQPUhP9R/LvXnbTkjadSFEnhgydSOudGmHBtlmEagksdmowTyjhaWu7fK3NRn1KBgk+ChiuoWHsunNrti1tH+gsQTzrbsnAC/Z/ocNt+ee01Ak7boQwufmbjnOX1tOeK6rKCd5SFsIwHvuBwFrkPGonFpboEjwUQT52VQevqWqpex7vTtnzGDC1JPcry213EvWTT5buAchhPAF3TAZPsuazXSo7Xf8FJ3lekNWGw0s96qXLsYr3SSNekEiwUcR1bWB9aClBIrxhfsuAF6wzcCfZMv9zxfvk+kXIYRPrI2O5VxiWjbTm5RD9FFXAfCe+wGv+i/cVsdnbRM5Q4KPIqpFWGlKBdgtZVP0ThwxylFeOc8TWqTXM89N3kjkNglAhBC560TcJcv1K7apqIrJX3pLtpthXvVDShTzKhP5mwQfRZSmKoy9M9xS5sTOR+57AXjG9hclueD13Ki/dsj6DyFErpq0Ku3k2luUnXTU/sVlap5jIdIrXdyPFmGlfdk8kQMk+CjCekRUZEDb6payP4zW7DCqEawk8nwmh87FxCWxNjrWRy0UQhQ1T/64ls1H4wFQMBhm/xWAqXpHDpkhXvXH9AmXHS4FkAQfRdzwng3o3qC859pEZay7HwCPaFFUV7ynWU5dSPJZ+4QQRcd/Jq21nOHSW11NY/UAF01/PnPf7VW/VVhpekSEepWL/E+CD0HN8kGW67+NhizRG2FX9Mtp162idpz0VdOEEEXE7M3HWLgzLfBw4OQV+zQAJrjv4AwlvJ65r3kVn7VP5CwJPgQZ98sDvO1+CN1U6K6to5myy3Jv9pYYnG7D6xkhhLgeumHy4v/+tZT11+ZTWTnDcbM03+vdM31OFpoWXBJ8CFrVLONVtteszDS9IwBv2H8l46Fzbd5d7IumCSGKgM8X7cGVbiF7aeIZePnwuA9d91kOj/PUkYWmBVq2g4/ly5fTu3dvKlasiKIozJo1y3K/f//+KIpi+dOyZcucaq/IBS1rlPFKuQ4ph84lmA4aq/vpra623Dt9MZnRs3f4qolCiEJKN0zGL9lnKXveNoNg5RLbjOqek7czkoWmBVu2g4+EhAQaNWrEuHHjrlinW7duxMTEeP7MnTv3hhopcpemKnxwd4RX+WlK8pW7NwCv2qfiwGm5//3KaJl+EULckC8W7bGc31JDOc5D2iIAxrofwszkx9SAtmGy0LSAy3bw0b17d8aMGUPfvn2vWMfhcBASEuL5U7q0DI3ld97bblOGQL/Ve3LCLEVl5QyPagu8nvt59UGftE8IUfjohsmEZfstZa/ZpmBXdKL0pl5p1AEGtK3O8J7eJ3SLgiVX1nwsXbqU8uXLU6dOHQYMGMCpU6dy49OIHDa8ZwMGtE3JHmi7XJaEw5PYZ7BtFqWItzwzf8cJhBDienwWtYdkd9paj1uUnXTVNuA2Vd51P+hVf9wDTRje0zsgEQWP7dpVsqd79+7ce++9VKtWjejoaN588006derEhg0bcDi8Fw0lJyeTnJx2jkh8fMoPN5fLhcvl8qp/I1Lfl9PvLUxe6Vqb8AqBDJ/1LxfdKfOpM/S2PK5F0kA9xPO2mYxyP+apvzb6HPO2HOW2myrkVZPznHxdZZ30VfYU5v76eMEuflh1iNTddgoGw+2/ADBZ78x+s5KnroLJPU1Cub1BuSv2RWHuq5yWW32Vnfcppmled65sRVGYOXMmd9555xXrxMTEUK1aNaZOnZrpVM3IkSMZNWqUV/nkyZMJCAi43qaJHLDkmMKswykLUVur25js9zYuU6Ob813LN4YRTd2U9o4rhRDiqg7Ew2fbU34HvktdwSd+E7hgFqND8secTZfXY1RTNyXle0y+l5iYSL9+/YiLiyM4OPiqdXN85COj0NBQqlWrxt69ezO9P2zYMIYOHeq5jo+Pp0qVKnTt2vWajc8ul8tFVFQUXbp0wW63X/uBImzNvlN0ZD1/HTbQUVllhBOl30wXbQNv2X7mMderpP7GMmqjjXrli/HX4LZ52+g8Il9XWSd9lT2Fsb90w6TDh0uITXDhNFK+hwSQxGv2KQCMd/exBB5g8va/Gj881vyqW2sLY1/lltzqq9SZi6zI9eDj7NmzHDlyhNDQzFcmOxyOTKdj7HZ7rn0B5ea7C4vmNcoxf0/K6FZqio8x7odor26mvbaFTvomFhtNPfV3nbrEu/P38mavorsQTL6usk76KnsKU3+t33+WExfcpE9uONA2iwrKeQ4aFTJJKKZQOtCflrXKZ2lrbWHqq9yW032VnXdle8HpxYsX2bx5M5s3bwYgOjqazZs3c/jwYS5evMjLL7/M6tWrOXjwIEuXLqV3796ULVuWu+66K7ufSuSh1H/k6f+tHzJD+EHvAcCbtp+x47Y8I1tvhRDX8t0K6+6WKspJ/qOlpGMY634IJ94/wEb0ri85PQqZbAcf69evp0mTJjRp0gSAoUOH0qRJE9566y00TWPr1q306dOHOnXq8Nhjj1GnTh1Wr15NUFDQNd4s8iMlw7/3L9x3csosSZh6kv5apFf9YdP/9SoTQggAp9uwHBwH8IbtVxyKm+V6Q6KMm72e+eLBJnQLl5wehU22p106dOjA1daozp8//4YaJPKX4n4ayZfSRjMSKMb77vv50P41z9tmMku/ldOU9Nyfvuk4nW8KkQRAQggvw2ZssVy3Ubdyu7Yet6nyf+5HyHjOVMe6ZendqKIPWyh8Rc52EVc1uk+4V9l0vS2bjRoEKZd42fY/r/vPTd5I5LYYXzRPCFFAzN1ynOkbj3muNXTesv0MwM96F/aZlb2eeapdLZ+1T/iWBB/iqro2COHpdmGWMhOV/3M9CsC92jIaKge8nnttxlZ047p3cQshCpHIbTE8N3mTpewhbSF11aPEmoF84r7b65kycnBcoSbBh7imYT3qM75fU0vZRrMOM/RbURWTkfZJZDz19nyii3GLrYdFCSGKHt0wGfnndktZSS4w1PY7AB+57yOeQK/nRsvBcYWaBB8iS3pEhPJ693qWsvdcD5BgOrhZ3Usf9W+vZ75cuk9GP4Qo4tZGx3IiPtlS9pLtN0oqCew0qjJF7+T1TK+IUFk3VshJ8CGyrH8b6/TLSUrzpbsPAK/bJ1OcS5b7TrfBC1OtQ61CiKJlYYbzn25SDtHv8qm1o9yPYmT4MVTMrvLZA0181j6RNyT4EFnmZ1N58tZqlrLv9R4cNCpQQTnPENt0r2dmb4mR3B9CFFG6YVpOvlYwGG2fiKaYzNZvYY3hnZTwmfa1ZLqlCJDgQ2TLm73CqVE27cydZPwY4e4PwONaJHWVw17PTFoV7avmCSHykeenbMCZ7nePe7TlNFP3kGA6GON62Kt+oMPGoE6yw6UokOBDZFvU0A7Y033lLDMaMU9vjk1J+a0m4+LTn9cc8m0DhRB5bu6W48zZetJzXYKLvGZLOb/lU/fdnKCM1zPv3x0hox5FhAQfIts0VWFQp9qWsv9zPUqi6aCFupu71RWWe4djLzF6tnW1uxCi8NINkyHTNlvK/mubRhnlAruNykzUu3k9I4tMixYJPsR1GdSpNgHphj9iKMPn7pTze4bZJxPMRUv971ceZOycHT5toxAib3y+aA9OPW0ENELZTz9tMQBvuR7HnSG5tkNTZJFpESPBh7gumqrw4b2NLGXf6z3Ya1SirBLPy7bfvJ75dkU0c7dI5lMhCjPdMPlmeVriQRWDMfYfUBWTGfqt/GPe5PXMsx1kkWlRI8GHuG49IirStEoJz7ULG29dXnz6sLaQ8Ewyn775xzbJ/SFEITZu8V4uudJWmfbTFhGhRhNvFuMdVz+v+v42lcGda3uVi8JNgg9xQ1663Zp4bLXRgD/01qiKyRj7RBSs22zPJjhZGx3ryyYKIXwkclsMnyzc67kuQxz/tU0DUjKZpj+EMtXH9zWWUY8iSIIPcUNa1ihDMbv1y2iM6yEumMVorO7nQW2J1zPvRe6U0Q8hChndMHlhykZL2Wu2KZRQEtlmVOcX/TavZwa0DZNFpkWUBB/ihmiqwjPta1rKTlOKT9z3APCqbQrlOG+5v/lIHG3eXSwn3wpRiHwWtZtkPe26hbKTe23LAXjT9Tg6mqV+hzplGd7TO8mYKBok+BA3bFCn2hR3WL+xTNK7ssUIo4SSyAj7T17PnIhP4tlfNkoAIkQhELkths+X7PdcO3Dyjv07ACa7O7HJ9F7T8XR7SSZWlEnwIW6Ypip8cHeEpUxHY5hrAG5TpZe2hk7qRq/nTGDUXztkCkaIAkw3TF75/V9L2XO2P6ipxnDKLMm77ge9ngn2t9EirLSvmijyIQk+RI7oEVGRAW2rW8q2m9X5Xu8OwGj7RAJI8nouJi5JFqAKUYCNW7yX+KS0+ZbaylGe1f4EYITrMeIp7vXM2LsayiLTIk6CD5FjhvdsQPcG5S1ln7rv5ohRjkrKWV7KJPcHQFSGUy+FEAWDbpiMX7rPc61g8I79O/wUnSj9ZuYZLbye6VyvPL0bVfRlM0U+JMGHyFHjHmpG+l9oLuHPG+4nAOivRRKh7Pd65tc1h2TqRYgC6ItFe0h2p/3bfUhbRDN1DxdNf95y9QesoxshQQ6+79/ct40U+ZIEHyJHaarC4AynUi4zGjFLb42mmLxr/w4bbsv9ZN3kvq9X+bKZQogbFLkthk8XpY16VCCWV2xTAXjffT8xmRwct/zVTj5rn8jfJPgQOe75znXwt1l/4xnteoRzZiD11UM8qc3zembDofP89e9xXzVRCHEDMltkOso+iWDlEpuMWvyid/F6pldEKH42+ZEjUshXgshxmqrwaYZDos5SgrHuhwB40fY7VZWTXs+98vu/Mv0iRAGQcZHp7eo6umnrcJkar7n+g5HhR4u/TZWD44SFBB8iV3QLD+XJNtUtZb/r7Vil18dfcfGe7Vuv1OuXXAbjFu9DCJF/ZVxkGkwCo+w/AvC13ovdZlWvZySFushIgg+Ra26rH5KhROE19wAumX600nbwkLbI65lxS/bK6IcQ+diqfWcsi0zfsP1CiHKOA0YIX7jv8qpfq1xxSaEuvEjwIXJNi7DSlCxmt5QdNivwnvsBAIbZJlNZOW2579JNnp/snZBMCJH3IrfF8MgPaz3XHdTN3GdbhmEq/Nf1NMn4eT0zoncDXzZRFBASfIhco6kKj2eYeoGU1Ov/GPUoriTznu0br+mXOdtOMHbODh+1UgiRFZHbYnjml7RfDIJI9KRQ/0Hvxgazrtcz/naV1rXK+qyNouCQ4EPkqkGdalMywDr6YaLyiuspLpl+tNG2009b7PXctyuimbtFzn0RIj/I7MTa4bZfCFViiTYq8KH7vkyf+/heWeshMifBh8hVmqrwbt+GXuWHzBDed98PZD79AvBf2f0iRL7wxaI9lhNr26n/8oBtqWe6JQmH1zMD2obJWg9xRRJ8iFzXLTyUF2/zPtXyR/121hp1CVSSeNf2DSlHzaVJcOq8MHWTj1ophMiMbph8ni6ZWBCJvGv/FoCJejfWm/W8nmlcOZjhPev7rI2i4JHgQ/jEoE61CQn2t5Sln365VdvOg5lMv8zeEiPTL0LkocGTN1hWZQ2z/UpFJZaDRgU+uMJ0y39vv8k3jRMFlgQfwic0VWHkHd6/CR00Q/ng8vTLcNuvVMJ7+uWV6Vtk+kWIPDB2znbmbktLCNhO/Zd+tiVXnW4pGWCnZU3v1OpCpCfBh/CZbuGhjO/nneXwR/121hl1CFSS+ND+tdful4vJbkk+JoSPzd1ynG9XHPRcl+QCH9i/BlJ2rK3LZLoF4N2+DWWRqbgmCT6ET/WIqMi4DGmWDVRedj1DgumglbYj07NfvlgsyceE8BXdMHltxtZ0JSZj7D9QQTnPPqOiJ1dPRuP7NaVbuCwyFdcmwYfwuV6NK9I7wpr99JAZwmj3IwD81zaNesphy323YfLCFFl8KoQvpJzdknb6dB/1b3pp/+AyNV50PZfpdMu4B5rI7haRZRJ8iDzx6QNNCfDTLGVT9Y5E6U1xKG4+sX+JA6fl/uytMTjd1ikZIUTOitwWwycL93quK3KG0ZfPbvnM3ZetZg2vZ3o1DKVX44q+aqIoBCT4EHlCUxWebpfxm5jCMNcAzpjB3KQeYajtN6/nHv3+H980UIgiSDdMXvn9X8+1gsGH9q8IVhLZaNRign6H1zN2TeGzB+XEWpE92Q4+li9fTu/evalYsSKKojBr1izLfdM0GTlyJBUrVqRYsWJ06NCB7du351R7RSEyqFNtShSzWcrOUILXXAMAGKDNpaVqTbO+JjqWj+fvlvUfQuSCIVM3Ep+Ulk3sCW0erbUdJJgOXnQ9h47m9cygjrVlganItmwHHwkJCTRq1Ihx48Zlev/999/n448/Zty4caxbt46QkBC6dOnChQsXbrixonDRVIX37o7wKl9o3MwUd0dUxeQj+wSCSLTc/3zJPm4eHUXkNsn/IUROmbvlOH9tOeG5rqMc4RXbNABGux/hkJnxlGoIdNgY1KmWz9ooCo9sBx/du3dnzJgx9O3b1+ueaZp8+umnDB8+nL59+xIeHs6kSZNITExk8uTJOdJgUbh0Cw/lq4ebYs/wlTja/QiHjPJUUs7yf/aJXs+dv+TimV82SgAiRA5wug2em5y2oNuBk8/sX+JQ3CzUmzBV75jpc+/fHSGjHuK65Oiaj+joaE6cOEHXrl09ZQ6Hg/bt27Nq1aqc/FSiEOkWHsr2/+uOku57WCL+KcO8psJd2t/0VZdn+uywGVtlCkaIGxC5LYa6b1i3t79u+5Wb1MOcNoN5zfUU4B1g9AgPkd0t4rrZrl0l606cSBmyq1ChgqW8QoUKHDp0KNNnkpOTSU5O9lzHx8cD4HK5cLlcOdk8z/ty+r2Fka/7SgGGdAxj3JIDuMyUko1mHT5x38PL9t8YbZ/IJmdtos20b3Z2xSQx2cmqvSdpWSPvMirK11XWSV9lT27318KdJxkybTNmuuCiq7qOx2xRALzkepYzlEj3hImfCjYVPr63Yb76e5SvrazLrb7KzvsU0zSv+9dGRVGYOXMmd955JwCrVq2iTZs2HD9+nNDQtB8SAwYM4MiRI0RGRnq9Y+TIkYwaNcqrfPLkyQQEBFxv00QBtuc8fLlTAxRUDH61v00rbQfbjOr0dY7Cid1T94UGbmoE51lThSjwYpNh1MaU30NDOcs8x2uUVBL42t2Td9wPpatp8nx9nZolMn+PEImJifTr14+4uDiCg6/+jTlHRz5CQlIWJJ04ccISfJw6dcprNCTVsGHDGDp0qOc6Pj6eKlWq0LVr12s2PrtcLhdRUVF06dIFu91+7QeKsLzqqwlL9/Hd2v1ogE5K9tMXXAOZp75GuHqQ12xT+D/3o576n2234dBM+reqxsu3Z57uObfJ11XWSV9lT27219Bpm5i/4xROI2XUQ0PnU78vKakksNmowYeXz1xKpQFf79YoWczOsv92zHdrPeRrK+tyq69SZy6yIkeDj7CwMEJCQoiKiqJJk5R9306nk2XLlvHee+9l+ozD4cDh8M6WZ7fbc+0LKDffXdj4uq8Gdq7Hd38fIllP2+53ilK87HqGiX4f8IQtkpVGOIuNpp77ybrC1ysPg6oxrEfeHeMtX1dZJ32VPTndX2PnbOevbadJv5ZjsG0mt6i7uGAW43nXYFwZfjzoKOg6jLozAn+HX461JafJ11bW5XRfZedd2V5wevHiRTZv3szmzZuBlEWmmzdv5vDhwyiKwpAhQ3j77beZOXMm27Zto3///gQEBNCvX7/sfipRBGmqwvv3NPIqX2I04Xt3dwA+tH9FBWK96ny9PFoyoApxDRkPjAO4RdnJYG0mAMNdT3DY9B6pLu6n8dXDcnaLyBnZDj7Wr19PkyZNPCMbQ4cOpUmTJrz11lsAvPLKKwwZMoTnnnuOZs2acezYMRYsWEBQUFDOtlwUWqnbbzMO677nfoCtRnVKKxf5zO9LVLwDjcckA6oQV6QbJi9O22wpK0U8n/h9iaaY/K6340+jjddzIcF+bBl5uwQeIsdkO/jo0KEDpml6/fnxxx+BlEWoI0eOJCYmhqSkJJYtW0Z4eHhOt1sUct3CQ/npiRaWMid2BrsGc9H0p6W6kyG2372eWx0dK6MfQlzBF4v2kKyn7TFQMfjUPp6KSiz7jVDecvX3ekYF/n7ttny3xkMUbHK2i8i3WtYoQ+ni1jnEg2Yor7v+A8Dztll0VL1Puh02/V+vMiGKOt0wmbBsv6VskDaL9toWLpl+POd6gUT8vZ4bn8kopBA3SoIPkW9pqsKYPt6jZn8arZnk7gLAJ/bxVFZOW+5P33ScsXPkPCEh0hs8eSPJ7rRRj7bqFobYpgPwhusJdptVvZ558bY6MtUicoUEHyJf6xFRkQFtq3uVj3U/zCajFiWVBMbbP8WB03L/2xUHGTtnh9dzQhRFT0xcy9xtaee2hHKWz+zjUBWTye6OTDfaeT0TEuyQc1tErpHgQ+R7w3s2YEDbMEuZEzsDnc8TawYSoUbzlu1nr+e+XRHN3C1y9oso2np/vpzFu9NGB+24+dLvM0orF9lqVGeU+7FMnxt5RwOZbhG5RoIPUSAM71mfx1tXt5QdpywvuAZhmAoP2RZlev7L4Ckb5ewXUWQ9+eM/bD1uPVH8dduvNFX3EWcG8KxrCMl45+yQ6RaR2yT4EAVG1wbeR3qvMCL4zJ1ywvJY+w/UUw5b7usm3PeVHGooip7Zm4+xaNcZS1kvdTWP2+YDMNT1LEfN8l7PlQ6wy3SLyHUSfIgCo0VYaUoFeGfQ+1y/i6V6I4opTr6xf0RJrL/pbTh8nktO3es5IQor3TB5efoWS9lNyiHet38DwHj3HSwybs702TF3NpTpFpHrJPgQBYamKoy903v3i4nKENdzHDbKUVU9zTj752hYg4327y/yVTOFyHPPT9lIkist300p4vnW7yMClGSW6w350H1fps8NaBtGjwiZbhG5T4IPUaD0iKjI0+3CvMrPE8QA10skmA5u1bbzum2y5f6piy5ufXehrP8Qhd7o2duZszVtZ4sNN+Ptn1NZOUO0UYFBrsEYmXzr79kwhOE98+5sJFG0SPAhCpxhPeozvl9THJp1aHi3WZWhrmcBeNI2j3u0ZZb7R88nU/+tSCK3yQ4YUTiNnr2d71cetJS9YfuFVtoOLpr+DHC9RDyBXs+VKGbj8webepULkVsk+BAFUo+IUHaM7k73+tYFc/ONFnyaugDV9j2NlX2W+8lug2d+2SgBiCh0xs7xDjzu05bQ37YAgBddz7HPrJzps+/dHSHrPIRPSfAhCixNVZjwaHPKB1m3Cn7m7st8vRkOxc3Xfh9TnnNez770279yBowoNDI7qbapsocxth8A+Mh1D1FGs0yflW21Ii9I8CEKvE/ub2K5NlEZ6nqW3UZlKijn+cbvY68MqAnJOk1HR8kIiCjwMjuptiJn+MrvU/wUnbl6C8bpd2b6rGQxFXnFltcNuF66ruNyubL1jMvlwmazkZSUhK7L1surya2+stvtaJqWY++DlAPogv014pPS2plAMQa4XuIPvzdprO7nY/t4Brmex0wXb19MdvPMLxv56uGm8pufKLDunbDKclJtIIl87/cB5ZXz7DSq8LLrGcvXfXqSxVTklQIXfJimyYkTJzh//vx1PRsSEsKRI0dQFPkHdzW52VclS5YkJCQkx96rqQrv39OIZ37ZaCk/bFbgGeeL/Oz3Nj21tRw2p/Ge+0Gv51+bsZUu9UPkm7AocJ6Y+A8bj5z3XGvojLN/wU3qEU6ZJXnS+d9MT6r1t6t8en9jCbpFnilwwUdq4FG+fHkCAgKy9QPMMAwuXrxIYGAgqiozTleTG31lmiaJiYmcOnUKgNDQnPvG1y08lK8ebsqzv2wk/Wbaf8ybeNX1FJ/4TeBZ218cMiswVe9kefZ8ootxi/fxwm21c6w9QuS2J39cy+Ld6TOYmoy0TaKD9i+XTD+edL7Mccp6PedvU9ky4nb8bPI9UOSdAhV86LruCTzKlCmT7ecNw8DpdOLv7y/BxzXkVl8VK1YMgFOnTlG+fPkcnYLpFh7KT0+04JEf1lrKZxptqeY+yRDbDMbYfuCYWZYVRoSlzicL9zCoUy0Z/RAFwujZ21i067Sl7EltHo/YFmKYCkNcA9lq1sj02Y/vayyBh8hzBeorMHWNR0BAQB63RNyI1L+/7K7ZyYrWtcpS3OEd0HzqvpsZ+q3YFIPx9s+ooxzxqnPP+L9zvD1C5LS5W47z/cpDlrIu6nqG234F4B33g8w3mmf67NPtJIOpyB8KVPCRStZrFGy5+fenqQof3B2RyR2F11wD+MeoR5ByiYl+71MuwxbcTUfjGD17R661TYgb5XQbvDpjq6WsibKXz+xfoiomk92d+Fbvmemz4x5owrAeksFU5A8FMvgQ4mp6RFRkQNvqXuVO7DzlHMp+I5RKylkm+b1PEImWOt+vjOaDubskDbvIdyK3xXDL2wu5kOT2lNVUjvGD3wcEKMks0yN4y90f8A7ux/drQq/GFX3XWCGuQYIPUSgN79mAAW29z4CJI5DHXa9w2ixBffUQ3/p95JUD5Mvl+yUNu8hXFu48yTO/bORcYtpUZQVimeT3HqWUi2w2avCsawjuTJbxje/XlB4REniI/EWCDx/p378/iqKgKAp2u50KFSrQpUsXfvjhBwwj65k2f/zxR0qWLJl7DS1EhvdMOQMm4wqQw2YF+jtfJd4sRkt1J5/bx3mdgitp2EV+MjzDVEswF5nk9x6VlTMcMEJ4wvlKpltqxz3QRNZ4iHypyAYfumGyev9Z/th8jNX7z/pkmL1bt27ExMRw8OBB5s2bR8eOHXnhhRfo1asXbrf72i8Q2dYjIpQ9b/egcaVgS/l2szpPuV4i2bRzu7b+chpq76+BV3/fIlMwIs9dSE4Ljh04+dbvY+qpRzhpluRR1zBiCfZ6plO9cjLVIvKtIhl8RG47wa3vLebBb9fwwtTNPPjtGm59b3Gu/5brcDgICQmhUqVKNG3alNdff50//viDefPm8eOPPwLw8ccf07BhQ4oXL06VKlV47rnnuHjxIgBLly7l8ccfJy4uzjOKMnLkSAB++eUXmjVrRlBQECEhIfTr18+TT6Oo01SFWYPbUrmk9TfDNUZ9nncNRDcVHrQt4SXbb17PxiW5eejb1RKAiDyxYPsJYpPBZaas49DQ+dw+jlvUXcSbxXjM+RpHzXKZPjugbU1fNlWIbClywcei3WcZOHkTMXFJlvITcUk8mwfD7J06daJRo0bMmDEDAFVV+fzzz9m2bRuTJk1i8eLFvPLKKwC0bt2aTz/9lODgYGJiYoiJieHll18GwOl0Mnr0aP79919mzZpFdHQ0/fv39+n/S373wb2NvcrmGy0Y7n4SgMG2WTypzfGqsyb6HBEj58sUjPApp9tg6G//MmpjyjoOFYOP7BO4XVtPsmlngPNldplVM302tIQ/LcJK+7K5QmRLgUoydqN0w+T9hQcyGVxPGXBXgFF/7fB5qu169eqxZcsWAIYMGeIpDwsLY/To0Tz77LOMHz8ePz8/SpQogaIohISEWN7xxBNPeD6uUaMGn3/+OS1atPBkKRXQIqw0IcEOTsQnW8qn6p0oTTyv2P/Hm/ZfScaPX/QuljoJTl3OgRE+E7kt5vJxAanfh0zG2r7nTm0VLlPjOdfz/GPedMXnR/SuLwnzRL5WpEY+1h2M5eQF5xXvm0BMXBJro2N91yhS0o6n5r5YsmQJXbp0oVKlSgQFBfHoo49y9uxZEhISrvqOTZs20adPH6pVq0ZQUBAdOnQA4PDhw7nd/AJDUxVG3tEg03vj9T586b4DgDH2idyrLc203pCpm7jklEMJRc5LXYf2f39tz3BOkckI2088aFuCfjl76SLj5iu+Z3w/CZBF/lekgo9TF5KvXQk4dSHp2pVy0M6dOwkLC+PQoUP06NGD8PBwpk+fzoYNG/jyyy+Bq2cDTUhIoGvXrgQGBvLLL7+wbt06Zs6cCaRMx4g03cJDGd+vSSaZEBQ+cN/P9+7uALxn+5Y7VO+Mp0luk5veiuSduZKMTOScyG0xtHk3ZR3aD38ftNz7r20aj9vmA/CK62nmGC2v+B7Z3SIKiiI17VI+yJHFet5b1nLL4sWL2bp1Ky+++CLr16/H7Xbz0Ucfec5T+d///mep7+fn53XE/a5duzhz5gzvvvsuVapUAWD9+vW++R8ogHpEVORLFJ6bvDHDHYXR7odx4ORh2yI+tk/A6bITabTwesfXy6MxTJPhPTMfSREiq9KmWLwN0mYy0PYnAG+4Hme60e6K73m6XZjsbhEFRpEa+WhevTQVgvwy+a03hULuLtRKTk7mxIkTHDt2jI0bN/L222/Tp08fevXqxaOPPkrNmjVxu9188cUXHDhwgJ9//pmvvvrK8o7q1atz8eJFFi1axJkzZ0hMTKRq1ar4+fl5nvvzzz8ZPXp0rvw/FBY9IlJOwS3hnzH+VnjT/Ti/udthUwy+sH9BV3Vdpu/4dsVB/vr3eO43VhRaumHy0v/+zfTe89oMXran7MAa43rIax1SKodNYXw/SZ0uCpYiFXxoqsIrt6Wc9JgxAEm9zs2FWpGRkYSGhlK9enW6devGkiVL+Pzzz/njjz/QNI3GjRvz8ccf89577xEeHs6vv/7KO++8Y3lH69ateeaZZ7j//vspV64c77//PuXKlePHH3/kt99+o379+rz77rt8+OGHufL/UJh0Cw9l41td6dnQunjXROVV91P8obfGruh8af+cnuqaTN8xeMomFmw/4YvmikLo+SkbSPBaQ2Tyou03htp/B+A91wN8d4XzWrqFV2DH/3WXDKaiwFFM08xXCQzi4+MpUaIEcXFxBAdbE+ckJSURHR1NWFgY/v7ZnxoxDIP4+HhWHU5k9Jydlu22oSX8GdG7vizUuiy1r4KDgz1TQDnlRv8ec8Oov7Yx8W/rSaEaOh/Yv6avthLdVHjJ9SyzjFu9nnVoJu+30Lm9W3f8HX6+anKB5HK5mDt3Lj169MBut+d1c/LU3C3HeW7ypgylJq/YpvHc5amWMa6Hrhh4dKxblomP35LLrSw45Gsr63Krr6728zujIrXmI1W38BBuDw9lbXQspy4kUT4oZapFtqYVXSN6h3P4bCKLdp32lOlovOx6BrepcZ9tGR/bJ2B3u/lN72B5NjV8b/f+Yv7vrkYSwIpr0g2ToV7TLSav2ybzlC0l18wo1yNM1Ltn+nzVUv4SeIgCrUhNu6SnqQqtapahT+NKtKpZRgIPwff9W9C5njVbpIHKq+4B/OLujKqYfGD/hoe0hZY6TkMhNjklG6qcByOy4uMFu0lyp53ppGAwwvaTJ/B409X/ioFHhUA/lr/a2SftFCK3FNngQ4jMfN+/BY+3qW4pM1F5w/0EP7i7ATDW/gNPaX9Z6ozaaPOMgLw6Xc6DEZnTDZO+41fy5dL9njIbbj6yf+XZTjvM9SQ/610zPJn29bTiNQk8RMEnwYcQGYzo3YABbcMylCr8n/sRJrh7A/C6fQqv235FIe23V6eRMnoWd8nNF4v2+qq5ogDQDZNPo/ZQZ/hcNh6O85T7k8zX9k/oq63Ebaq86HyWKbo1uHAoJvbLA7MD2obhZ5Nv26Lgk69iITIxvGd9xj3QJEOpwnvuB3nb9SAAT9nm8KH9K2x4n0j86aK9Mv0igJQ8Hg1HzufTRXvR0w2IBXORn/3eobO2iUumHwNcLzHTaJvhaZNXGumoKnSpX57hPWU7rSgccjz4GDlypOfE1dQ/Gc8hEaIg6NW4IuP7NfUq/0bvzUvOZ3CbKndrK/nG/jHF8M6KO3jyRpzp5vVF0ZOaQCwxw3ba8pxjmt9omqt7iDMDeNg5jCVGxmAX7AqULQYf3h3Bt48291Wzhch1uTLy0aBBA8+pqzExMWzdujU3Po0Qua5HRChDOtfyKp9utOMp11AumX500jbzq9/blCLeUsdlQN035/HZwj2yBqQIcroNXpiacSst1FaOMsMxgpvUI5w0S3K/8y02mHUzfUeZ4inbILs1lB1UonDJleDDZrMREhLi+VOuXLlrPyREPjW4cx1KFPPelb7YaMpDztc5bxanqbqPmX4jqKkcs9QxTfhk4V7qvzWPT6N2SxBSRERuiyF8xDyS3da/71vVrUz3G0Fl5QzRRgXucY5gl1k103cowIgrHIQoREGXK3k+9u7dS8WKFXE4HNxyyy28/fbb1KhRI9O6ycnJJCenHfgWH5/y26PL5fI6TM3lcmGaJoZhYBjZH85OzaeW+g5xZbnZV4ZhYJomLpcLTdNy9N255d27GjBk2mYMA1yenycKG8063OMcwUT7B1RXTzLDbwTPuF5ktZHhh4ZpMGHpXn76+wBj+zbktpsq+Pp/IV9I/Td9tYMSC7oF208y9LfNlzeopG3hf0BbzBjbD9gUg3+MejztfJHzBKV70kQFbCoEOWyMuSuc9rVKExVduPsrpxSFr62cklt9lZ335XiG03nz5pGYmEidOnU4efIkY8aMYdeuXWzfvp0yZcp41R85ciSjRo3yKp88eTIBAQGWstQRlSpVquDnJ5kkc9rkyZMZNmwYhw4dunblG+B0Ojly5AgnTpzA7fZerJnfuQ345yT872Ba7F6GOL7x+5ib1b24TI033E8wTe/ouR9sN7i/hkF47hwbJPKZ1Sdh6oGUrw8Fg1dtU3nGNhuAGfqtvOYagJP0mSVNHqup07R8HjRWiBySmJhIv379spThNNfTqyckJFCzZk1eeeUVhg4d6nU/s5GPKlWqcObMmUzTqx85coTq1atfV1pu0zS5cOECQUFBKEreJBU7cuQIo0aNIjIykjNnzhAaGkqfPn148803PcFZp06daNSoEZ988kmm71iyZAljxozh33//JSkpiUqVKtGqVSu+++47bLarD2YtXbqUzp298wS8/vrrvP7661y4cIHy5ctjmibDhw8nMjKSjRszP3HzeiUlJXHw4EGqVKmSb9KrZ1Xk1hO8PP1fTBOcngGhlK8lB07et39DH20VAF+7e/Ke+0GMdLObdsVEVaFCkIMFL7YvcsntXC4XUVFRdOnSpdCkwF648yQj/9zO+Uspv/W5ddAvf00U5xIf2ydwu5ZyyvQnrrv5TO+L9XSplBEP++VBwE/vb+wZGSuM/ZVbpK+yLrf6Kj4+nrJly+aP9OrFixenYcOG7N2bed4Dh8OBw+F91L3dbvfqFF3XURQFVVWv67yR1OmD1HdsOXqed+buYliPekRULpnt92XXgQMHaNWqFXXq1GHKlCmEhYWxfft2/vvf/xIZGcmaNWsoXbq0pY0Zbd++nZ49e/L888/zxRdfUKxYMfbu3cvvv6ccQnWtfkm9v3v3bssXR2BgIMWLF6d48eIAlqmWnD7bRVVVFEXJ9O84v+vdtApbYuL5dsXBdKUmoJCMHy+4BnLACOVF+3Sets2hnnKEF1wDPcPrLlMBHQ6fd9Jn/Gre6Fmf1rXKFrkgpCD+3WcmclsMz05OTZNu/TusoRzna/sn1FaPkWzaeMX1FH9kcjYQl7PF2DSNj+7LPD1/YekvX5C+yrqc7qvsvCvX83wkJyezc+dOQkPz32rtGRuPsfrAWWZsPHbtyjlg4MCB+Pn5sWDBAtq3b0/VqlXp3r07Cxcu5NixYwwfPvya74iKiiI0NJT333+f8PBwatasSbdu3fjuu++yNRVVvnx5y6LgwMBAfvzxR0qWLAnAjz/+yHvvvce///7r2TL9448/Xuf/eeEyvKc1CVnKL6ypA4gKn+l3M9g5iEumH+21Lfzl9wYNlGiv9+w+eZFHflhL/TfnMXfLcV80XeQg3TB5dfqWTO91Vjcwy+9NaqvHOGGW4n7nW1cIPFI0rVqSLSNvl3OBRJGR48HHyy+/zLJly4iOjuaff/7hnnvuIT4+nsceeyynP9V1OR6XxNZjcWw7Fsdf/6Z8w//r3+NsOxbH1qNxHD2XmCufNzY2lvnz5/Pcc89RrFgxy72QkBAeeughpk2bxrVmwUJCQoiJiWH58uW50s5U999/P4MGDbJsm77//vtz9XMWJMN71md8v6b421VsGgy8SSd9Cuy/jNbc5fw/DhoVqKKeZrrfSO7RlmX6rmTd5LnJm3hn7g4ftV7khC8W7SXuknXNkoLBC9p0vvf7iGDlEmuNuvROHstm03u7dqonb63GjOfaFLnRL1G05fi0y9GjR3nwwQc5c+YM5cqVo2XLlqxZs4Zq1arl9Ke6Lj0mbPB8nPpPPTbBSa8vVnrKD76b+RHWN2Lv3r2YpslNN92U6f2bbrqJc+fOcfr06Uzvp7r33nuZP38+7du3JyQkhJYtW9K5c2ceffTRa86xpVe5cmXLdcZFpsWKFaN48eKeRb7CW4+IUG6rX4F27y2kTslLaED6VFK7zKrc4RzNJ/YJdNY28aH9axor+xjtfoRkvEepvl4eTTG7jcGda8sPonzM6TZ4bfq/zNhkHa0qQxyf2MfTTkvJazTJ3YUx7kdwXeXb7LgHmtCrccVcba8Q+VGOj3xMnTqV48eP43Q6OXbsGNOnT6d+/fyTEnhs79rYLn9jT/09NfW/NlXh0/sb50WzPCMe11oIq2kaEydO5OjRo7z//vtUrFiRsWPHekYosmrFihVs3rzZ86dUqVI31P6iys+mMqJ3yte3LZNdw/EE8h/XS3zsugfDVHjYtijTfCCpPl20l0Yj5zN7s2+mAkX2jJ2zgzpvzPMKPFqp25nnGEY7bSuXTD9edj3NCPfjEngIcQVF7myXng3KM+PZVpnemzWwDXc2qZQrn7dWrVooisKOHZkPre/atYtSpUpRtmzZLL2vUqVKPPLII3z55Zfs2LGDpKQkvvrqqyy3JywsjFq1ann+5PSi0qIkdWdCyWKZL7YyUflc78sTrv9y1gyivnqIv/ze4D5tCemnalJddOoMmrqZvl+ukKRk+YTTbdD142V8u8K6dkfF4EXb7/xqf5vyynn2GJW4wzmG3/X2V33f0+3CJPAQRVqR/omTOsjgi123ZcqUoUuXLowfP55Lly5Z7p04cYJff/2V+++//7q2AJcqVYrQ0FASEhJyqrlAysplXdevXVEAsOy/HXnxtjrYtcz/DpcajemW/C4r9QYEKMm8b/+WcfYvCCbzv7eNR+KpM3yuLEbNI7phsnr/WR6fuJY6b8xjz6mLlvtVlJNM8RvDC7YZqIrJFHdH7nCOYa9Z+QpvhECHxvh+TRjWI/+MBguRF3J9q21+VCbQj3KBDkJL+nN/8ypMW3eEmPNJlAnM3cRl48aNo3Xr1tx+++2MGTPGstW2UqVKjB071lP39OnTbN682fJ8SEgIf/zxB5s3b+auu+6iZs2aJCUl8dNPP7F9+3a++OKLHG1v1apViY6OZvPmzVSuXJmgoKBMt0WLFJqq8MJttRnUqRYvTNnE7K3e02CnKcUjrmE8bczmJdtv9NLW0FTdw6uup1hhRHjV1014bvImemw5zu3hoZQP8qdFWGlZE5KLnG6DYTO2MHtLDMmZHgxo8qC2mDdsv1BcSeai6c/rrif502hzxXf6aQoDO9ZiUCdZzyMEFNHgI7REMVa+1hE/LSXfRL8WVXHqBo7MJu1zUO3atVm/fj0jR47k/vvv5+zZs4SEhHDnnXcyYsQIT44PSMk2OnnyZMvzI0aMoE+fPqxcuZJnnnmG48ePExgYSIMGDZg1axbt2199qDe77rjjDiIjI+nYsSPnz59n4sSJ9O/fP0c/R2GkqQrjHmpKt83HGZTJwWImKl/pd7DaqM9n9i+prp7kZ793meLuyFj3Q1wkwOuZudtOMnfbSQBCS/gzond92ZaZC96Zu4Ovl3tvi05VnnO8Z/+GjlpKbo81xk287Hqao+aVU5M+37EmL3SpK0GHEOkUyeADsAQaiqLkeuCRqlq1akycOPGqdZYuXXrV+z///PN1f/4OHTpccTtv//79LcGFw+Hgt99+k/Ug16lX44qoqsJzkzPPEPuvWYvuznd4xTaNx23zedC2hHbaFl51PcVKo+EV3xsTl8Szv2xkwsNNJQDJIbphMujXDczbfvIKNUzu1ZYx3PYrJZUEkk0777vv5we92+X8pJn74sEm9G4kazuEyEh+qgiRi3pEhPLVw00JCc58uuoS/oxyP8YDzjc4bJSjknKWX/ze4SP7eMoQd8X3msDIP7fLgtQcMHdLDHVen3vFwKOmcoxpfqP5wP4NJZUE/jVq0MP5Nt/rPa4aeAxoW10CDyGuQIKPQqZ79+4EBgZm+uftt9/O6+YVSd3CQ/n7tc68eFudK9ZZY9Snm/M9Jrm7YJgKd2srWex4iYe0hZeTb3s7EZ9Mhw8W8/68Xfy974wEItngdBt8vWwft4yN4rnJG8lsWbUDJy/afmee32vcou4i0XQwxvUQfZ2j2G9efVfc0+3CGN6zwVXrCFGUFdlpl8Lqu+++89pNkyr9mhLhW6mLUeuGBPLajK2cT/Q+ejoRf0a4H2em3pax9u9poB5irP0H7tWW8YbrcbaZNbyeOXIuifHL9jN+2X6C/TXevyfzs0GKMt0wWRsdy6kLSZQP8mfhzpN8v/LK6zrApIu6geG2X6mupoyGLNYb85b7cY6a5a76uW6pXoqf/9MSP5v8XifE1UjwUchUqpQ7eUpEzugWHkqX+iGs2X+W12dt4dBZ70Bxs1mLO5xjeESL4iXbbzRW9zPb8QbT9Vv50HU/MZTJ9N3xSTrP/LKRdrXL0r5OOR5pVb3I/xCM3BbDqL92EBOXlKX69ZTDvGn7mTbadgBOmiUZ5XqUucYtZDw4Lj27pvDZ/Y3pESHTLEJkRdH+ziREHtBUhTa1y7Lsv5344sEmFPfz/meoo/Gj3o3OyR8yU0/Zwnm3tpIljqG8bJtGcTIf3QJYvvcMo+fspN6b84r0eTGR22J49peNWQo8ynGOt23fMcdvGG207SSbdsa5+9Ap+SPmGi25WuDhsKlsH9VNAg8hskGCDyHyUO9GFdkyshv3NM08MdUpSvGiayB3JI/mH6Me/oqLQbY/WOp4kSe1OfiTfMV3G2bKeTFP/riW1fvPFqk1IU63weszt2WSP9aqNPG8bvuVFY4h9LMtRlNMZust6ez8kA/d95NAsWu8AT57oHGRH2ESIrtk2kWIPKapCu/dE8HCXSczXQsCsMWsyf3ON+mibuA12xRqqjG8af+VZ2x/8ZW7N7/qt5FE5jtqFu06zaJdpylRzM4TbaoXmkRXGddy3FytFOuiY/lpzUEW7zyJK/N1ugCU4CIDbHN4XIukuJISwK036vCe6wHWmfWy9Pkl34oQ10+CDyHyAU1VeLdvQ579ZeNVfltXiDKascTZmLu0lQzWZlJVPe0JQr5z92SK3ol4imf6dNwlF58s3MtXy/bTs2EobWqXIyS4YGZMjdwWw8g/t3Mi/sojP5mpxGmesEXygLbYE3RsMcL4yH0fy4wIrja9AqAp8GiranRtEFog+02I/EKCDyHyiW7hoUx4uOk1f6i6sfGb3oGZ+q2WIGSYfQrP22bwm96eiXo3DpkhmT5/yWXw+8Zj/L4x5eTcksXsPN6mOs92qMWGQ+c4dSGJssUdoMCZi8n5LqX73C0xV0zcdiX1lYM8ZZtNL3UNNiVlSGSnUZWP3fcQZdzMtYIOu6bwbLsakqlUiBwiwUc+oigKM2fO5M4778z0/sGDBwkLC2PTpk00btzYp23LTPXq1RkyZAhDhgzJ66YUGqm7YcYt3scnC/dctW76IORO7W+e1OZyk3qE/rYFPKpFsdBoymS9E8uNRhhXWd51/vKIyKcL915x1KV0cTtj+oT7bFFl+imVssUdGKbJ6gNnWHvgLOsPXzn5WnoOnPRU19DPtphmalpfrtQb8I3ei+VZGOkIsKs83b5moZmqEiK/kODDR/r378+kSZMA0DSNihUr0rNnT95++21KlSoFQExMjOfj/KJTp04sW7bMq9zlcrFu3TqKF08b4r9W8CSyJn1OkKxsE3Vj43e9Pb/r7Witbuc/2lw6aZvpqm2gq7aBGLM0v+vt+J/eniNmhSu+52qLM2MTXDw3eRM9tx7ntptCiE1wUjrQQflA6wjJzdVKeUZPsjNikhpsnIi7xN/7zhC18xRxlzJf/3J1JvWUI9yvLaGvtoISSiIAblNljtGSb9w92W6GZelNQzrXZnBnCTqEyA0SfPhQt27dmDhxIm63mx07dvDEE09w/vx5pkyZAqScWpsfDRgwgP/7v/+zlNlsNsqVu3rCJXFjUkdB1kbHcvxcIhuPnOPE+SR2xsRxPN6ZyRMKq4xwVhnh1HQf4yFtEXdpKwlVYhlsm8Vg2yz+MeoxV7+FSL05J8l+0rk5W08yZ+uVzj8BVUnZZZPK367SvnZZmlUvQ+nifpy9kEAFYPCvG3ArGuWDHByOTWTdwViS3de/GydMiaGXupo7tNXUVo95yo8Y5Ziid+Q3vT2nyXpgP75fE9k6K0QukuDDhxwOhyfAqFy5Mvfffz8//vij537GkYO1a9fy9NNPs3PnTsLDwxk+fLjXO//8809eeukljh49SsuWLT2Hw507d46SJUsCsGrVKl577TXWrVtH2bJlueuuu3jnnXcsoxZXExAQkGlglH7apXr16gDcddddQMoBegcPHsxax4gr0lSFVjXLAGW4u1kVT/m11j3sNyvxf+5Hedf9IF3UDdynLaWtupVb1F3cou5ilH0S64w6ROotWGo0Yr9ZkWtNQWRFxt28SS6D+TtOMX/HKQAcmsn7LWDJ3jMk69f/+RQMGirRdNI2cZu6kXD1oOdesmljsdGEKXonVhgNr3r+Skayg0UI3yj4wYdpgisxa3UNI6WuU4OcOKnVHgDK9X0DPXDgAJGRkdjt9kzvJyQk0KtXLzp16sQvv/xCdHQ0L7zwgqXOwYMHueeee3jhhRf4z3/+w6ZNm3j55ZctdbZu3crtt9/O6NGj+f777zl9+jSDBg1i0KBB1zxdNzvWrVtH+fLlmThxIt26dUPTfHNKcFHVIyKUr9RrL051YmeO0ZI5RktCOEsPbS09tH9opu6h+eU/b/ILR82yrNAbssJoyBqjPrEE50q7k3V4cbV2hdNqrsakqnKKFuouWqk7aK/+S1kl3nPXZWqsNML5S29FlNGMCwRk6+1PtqnObfVD8tXCWiEKs4IffLgS4e2sDY+qQMmc/NyvHwe/rI0eAMyePZvAwEB0XScpKWUe/+OPP8607q+//oqu6/zwww8EBATQoEEDjh49yrPPPuup89VXX1G3bl0++OADAOrWrcu2bdsYO3asp84HH3xAv379PItCa9euzeeff0779u2ZMGEC/v7+12z3+PHj+e677zzXTz/9NB999JGlTuoUTMmSJfPt9FFhk35a5utl+1m65/RV65+gDD/o3flB704FYumureU2dQPN1d1UVs7woG0JD7IEgINGBTaZtdhk1GKbEcZes3K2f6Bn7krH5FmVJp766iHqKweJUKNpru6ignLeUifeLMYKoyFLjcYs1Jty7joCplIBdt7p21BGOoTwsYIffBQgHTt2ZMKECSQmJvLdd9+xZ88eBg8enGndnTt30qhRIwIC0r7ht2rVylJn9+7dNG/e3FLWokULy/WGDRvYt28fv/76q6fMNE0MwyA6Opqbbrrpmu1+6KGHLFM+qdM5Iu+lTsu0qlmGd+bu4OvlVzswLc1JSvOj3o0f9W74k8wt6i7aqltoq26lrnqU6upJqnOSu7S/Pc8cN0uz16jMATOUY2ZZjptlOGaW5ZRZijiKk4iDrE3dmARyiTJKPGWIJ0SJpZpyiqrKSaopJ6mhxhCinPN6ymlq/GvWZJ1Rj+VGBOuNOriv41uYv12lw+Wzb1rWKCMjHULkgYIffNgDUkYgssAwDOIvXCA4KAg1p6ZdsqF48eLUqlULgM8//5yOHTsyatQoRo8e7VXXNK+9+M40TZQM0z4ZnzMMg6effprnn3/e6/mqVatmqd0lSpTwtFvkX8N61OelrvWY+PcBpq49QvTZrE1HJuFgmdGIZUYjAIJJoJG6nybKPpqoe7lJPUyIco6KSiwVtVjasyXT9zhNjTiKk2AWQ0fFjYaOhoKJAyf+ihN/nBQnCYfivma7Dhgh7DCrscOoxgajLpvNmiTjl/UOSeemkEAGtK1JaMliMrUiRD5Q8IMPRcn61IdhgF1PqZ8TwccNGjFiBN27d+fZZ5+lYkXr1FH9+vX5+eefuXTpEsWKpZwvsWbNGkudevXqMXfuXEvZ+vXrLddNmzZl+/btPgke7HY7uq7n+ucRV+ZnU3m6fS2ebl+LyG0xDJ22mcSr5RnPRDzFWWFEsIIIuPzXGcxFaivHqKsepapyiorKGSoqZ6mknKEscdgVHT9FpxzxlEu3FuNqLpr+nDWDOU1JDpvlOWyW55BRgUNmBXabVbJ0rsrVODSFl7rWpX+bMDl7RYh8puAHHwVYhw4daNCgAW+//Tbjxo2z3OvXrx/Dhw/nySef5I033uDgwYN8+OGHljpPP/00H3/8Ma+++ipPPvkkmzdv9uyeSR0RefXVV2nZsiUDBw5kwIABFC9enJ07dxIVFcUXX3yRo/8/1atXZ9GiRbRp0waHw5HvcpYUNalrQr5YtJdvVx4gIfn6A8N4Atlg1mWDXjeTuyYBJFOCBEooCQSQhIaBTdGxXY5ekkw/kkj5k2g6OEvwFc+iyQkK8NmDTWQthxD5lPw6kMeGDh3Kt99+y5EjRyzlgYGB/PXXX+zYsYMmTZowfPhw3nvvPUudsLAwfv/9d2bMmEFERAQTJkzwrM1wOFK+sUdERLBs2TL27t1L27ZtadKkCW+++SahoTn/Tfmjjz4iKiqKKlWq0KRJkxx/v8g+TVUY0qUOW0bczpQBLXm8dTWC/HP6dw6FRPyJoQy7zKpsNOuwzqzHaqNBygiKEcE6sx5bzRrsNStzjHK5GniElvBnwsNNJfAQIh+TkQ8fSZ/PI71+/frRr18/wHu9RsuWLdm8ebOlLGOdO+64gzvuuMNzPXbsWCpXrmzZxdK8eXMWLFhwXe1evHjxFdfHZMzj0bt3b3r37n1dn0fkrvQLU9/o1cCTujz6dAJfLdtPkjv7m1+vT/qv3xtfd1HcT+M/bcNoEVYmX55DI4TInAQfBdz48eNp3rw5ZcqU4e+//+aDDz5g0KBBed0skY+lJS5LMbhzbVbtPcPvG49wJDaRZLeBw6ZSzM/GxWQ3/x7N2lkqWTWsgc47268vD0yAXSW8Ugmah5Wmdc2ysltFiAJKgo8Cbu/evYwZM4bY2FiqVq3KSy+9xLBhw7L07IoVK+jevfsV7x89ejSnminyMU1VaFu3HG3rZp4u3+k2mLQqmnUHzxFgV6lfsQRlAh3EJiSz/tA5Fu86hUvPWmp0hwYhwSn/zbgEpWQxG4+2qoZTN9h6NJ4Ah0azaqWpHxpMbKJTRjWEKEQk+CjgPvnkEz755JPrerZZs2Ze0zqpDMNXw/Aiv/OzqQxoV5MB7bzvDSDlULg1+8/y9/7THD+fRGhJf0oHOCgb5KBsgB+7Tl7gyLlEqpUO4IFmlVi4IJIfHmvOqYsuzwF1IcESWAhRlEjwUYQVK1bsiltwDcMgPj5rWyZF0aapCm1ql6VN7bKZ3k8/ouJypZxU2yKs9BWPFhBCFH6y20UIIYQQPlUggw+ZEijY5O9PCCGKtgI17eLn54eqqhw/fpxy5crh5+fnlV78agzDwOl0kpSUlDPp1Qux3Ogr0zRxOp2cPn0aVVXx87u+VNlCCCEKtgIVfKiqSlhYGDExMRw/nrXzXNIzTdOTrjw7QUtRlJt9FRAQQNWqVSUAFEKIIqpABR+QMvpRtWpV3G53ts8RcblcLF++nHbt2slit2vIrb7SNA2bzSbBnxBCFGEFLviAlHNL7HZ7tn8oapqG2+3G399fgo9rkL4SQgiRW2TcWwghhBA+JcGHEEIIIXxKgg8hhBBC+FS+W/ORemprbmTXdLlcJCYmEh8fL+sYrkH6Kuukr7JO+ip7pL+yTvoq63Krr1J/bmc8fT0z+S74uHDhAgBVqlTJ45YIIYQQIrsuXLhAiRIlrlpHMbMSoviQYRgcP36coKCgHN+OGR8fT5UqVThy5AjBwcE5+u7CRvoq66Svsk76Knukv7JO+irrcquvTNPkwoULVKxY8Zp5nPLdyIeqqlSuXDlXP0dwcLB8cWaR9FXWSV9lnfRV9kh/ZZ30VdblRl9da8QjlSw4FUIIIYRPSfAhhBBCCJ8qUsGHw+FgxIgROByOvG5Kvid9lXXSV1knfZU90l9ZJ32Vdfmhr/LdglMhhBBCFG5FauRDCCGEEHlPgg8hhBBC+JQEH0IIIYTwKQk+hBBCCOFTRTb4uOOOO6hatSr+/v6EhobyyCOPcPz48bxuVr5z8OBBnnzyScLCwihWrBg1a9ZkxIgROJ3OvG5avjV27Fhat25NQEAAJUuWzOvm5Cvjx48nLCwMf39/br75ZlasWJHXTcqXli9fTu/evalYsSKKojBr1qy8blK+9M4779C8eXOCgoIoX748d955J7t3787rZuVbEyZMICIiwpNcrFWrVsybNy9P2lJkg4+OHTvyv//9j927dzN9+nT279/PPffck9fNynd27dqFYRh8/fXXbN++nU8++YSvvvqK119/Pa+blm85nU7uvfdenn322bxuSr4ybdo0hgwZwvDhw9m0aRNt27ale/fuHD58OK+blu8kJCTQqFEjxo0bl9dNydeWLVvGwIEDWbNmDVFRUbjdbrp27UpCQkJeNy1fqly5Mu+++y7r169n/fr1dOrUiT59+rB9+3bfN8YUpmma5h9//GEqimI6nc68bkq+9/7775thYWF53Yx8b+LEiWaJEiXyuhn5RosWLcxnnnnGUlavXj3ztddey6MWFQyAOXPmzLxuRoFw6tQpEzCXLVuW100pMEqVKmV+9913Pv+8RXbkI73Y2Fh+/fVXWrduLUcxZ0FcXBylS5fO62aIAsTpdLJhwwa6du1qKe/atSurVq3Ko1aJwiYuLg5Avj9lga7rTJ06lYSEBFq1auXzz1+kg49XX32V4sWLU6ZMGQ4fPswff/yR103K9/bv388XX3zBM888k9dNEQXImTNn0HWdChUqWMorVKjAiRMn8qhVojAxTZOhQ4dy6623Eh4entfNybe2bt1KYGAgDoeDZ555hpkzZ1K/fn2ft6NQBR8jR45EUZSr/lm/fr2n/n//+182bdrEggUL0DSNRx99FLOIJHzNbl8BHD9+nG7dunHvvffyn//8J49anjeup7+EN0VRLNemaXqVCXE9Bg0axJYtW5gyZUpeNyVfq1u3Lps3b2bNmjU8++yzPPbYY+zYscPn7bD5/DPmokGDBvHAAw9ctU716tU9H5ctW5ayZctSp04dbrrpJqpUqcKaNWvyZAjK17LbV8ePH6djx460atWKb775Jpdbl/9kt7+EVdmyZdE0zWuU49SpU16jIUJk1+DBg/nzzz9Zvnw5lStXzuvm5Gt+fn7UqlULgGbNmrFu3To+++wzvv76a5+2o1AFH6nBxPVIHfFITk7OySblW9npq2PHjtGxY0duvvlmJk6ciKoWqgGzLLmRry2R8g3v5ptvJioqirvuustTHhUVRZ8+ffKwZaIgM02TwYMHM3PmTJYuXUpYWFheN6nAMU0zT37uFargI6vWrl3L2rVrufXWWylVqhQHDhzgrbfeombNmkVi1CM7jh8/TocOHahatSoffvghp0+f9twLCQnJw5blX4cPHyY2NpbDhw+j6zqbN28GoFatWgQGBuZt4/LQ0KFDeeSRR2jWrJlnBO3w4cOyfigTFy9eZN++fZ7r6OhoNm/eTOnSpalatWoetix/GThwIJMnT+aPP/4gKCjIM7JWokQJihUrlsety39ef/11unfvTpUqVbhw4QJTp05l6dKlREZG+r4xPt9fkw9s2bLF7Nixo1m6dGnT4XCY1atXN5955hnz6NGjed20fGfixIkmkOkfkbnHHnss0/5asmRJXjctz3355ZdmtWrVTD8/P7Np06ayJfIKlixZkunX0GOPPZbXTctXrvS9aeLEiXndtHzpiSee8Pz7K1eunNm5c2dzwYIFedIWxTSLyApLIYQQQuQLRW/yXgghhBB5SoIPIYQQQviUBB9CCCGE8CkJPoQQQgjhUxJ8CCGEEMKnJPgQQgghhE9J8CGEEEIIn5LgQwghhBA+JcGHEEIIIXxKgg8hhBBC+JQEH0IIIYTwKQk+hBBCCOFT/w+wva6mXOOYSAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "\n", + "np.random.seed(2018)\n", + "n = 1000\n", + "d = 3\n", + "Lambda = 0.001\n", + "true_beta = [2, 0.5, 3.7]\n", + "\n", + "# Make data set.\n", + "x = np.linspace(-3, 3, n)\n", + "y_real = 2 + 0.5*x + 3.7*x**2\n", + "\n", + "y = np.sum(\n", + " np.asarray([x ** p * b for p, b in enumerate(true_beta)]), \n", + " axis=0) + 0.1 * np.random.normal(size=len(x))\n", + "\n", + "\n", + "#Design matrix X does now include the intercept. \n", + "X = np.zeros((len(x), d))\n", + "for p in range(d): # (d-1)\n", + " X[:, p] = x ** (p+1) # (p+1 if not intercept included)\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", + "#Center by removing mean from each feature\n", + "X_train_scaled = X_train - X_train_mean\n", + "X_test_scaled = X_test - X_train_mean\n", + "#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered, note)\n", + "y_scaler = np.mean(y_train)\n", + "y_train_scaled = y_train - y_scaler\n", + "\n", + "\n", + "#Calculate beta\n", + "beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)\n", + "beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d)\n", + "print(beta_OLS)\n", + "print(beta_Ridge)\n", + "\n", + "interceptOLS = y_scaler - X_train_mean @ beta_OLS\n", + "interceptRidge = y_scaler - X_train_mean @ beta_Ridge\n", + "print(interceptOLS)\n", + "print(interceptRidge)\n", + "#predict value\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", + "print(\"test MSE of OLS:\")\n", + "print(MSE(y_test,ytilde_test_OLS))\n", + "print(\" \")\n", + "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,y_real,label='no noise')\n", + "plt.plot(x, X @ beta_OLS+interceptOLS,'*', label=\"OLS_Fit\")\n", + "plt.plot(x, X @ beta_Ridge+interceptRidge, label=\"Ridge_Fit\")\n", + "plt.grid()\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "14449052", + "metadata": {}, + "source": [ + "\n", + "\n", + "" + ] + } + ], + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/doc/src/week37/programs/LinearReg.py b/doc/src/week37/programs/LinearReg.py index be84e77f5..90d86a253 100644 --- a/doc/src/week37/programs/LinearReg.py +++ b/doc/src/week37/programs/LinearReg.py @@ -4,7 +4,7 @@ from sklearn.model_selection import train_test_split class linregOwn: """ - A class of linear regressions. Perform ordinarly least squares (OLS) and Ridge regression manually. Lasso + A class of linear regressions. Perform ordinary least squares (OLS) and Ridge regression manually. Lasso is performed using scikit-learn functionality. """ def __init__(self, method = 'ols'): @@ -267,4 +267,4 @@ class linregSKL: if self.yHat is None : self._sklPredict() self._R2 = metrics.r2_score(self.y_test, self.yHat) - return self._R2 \ No newline at end of file + return self._R2 diff --git a/doc/src/week37/programs/codeexamplesscaling.do.txt b/doc/src/week37/programs/codeexamplesscaling.do.txt new file mode 100644 index 000000000..e392d9812 --- /dev/null +++ b/doc/src/week37/programs/codeexamplesscaling.do.txt @@ -0,0 +1,370 @@ +TITLE: Scaling examples with own code and the library Scikit-Learn +AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo & Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University +DATE: today + + + + + +===== This note contains code examples with a simple scaling ===== + +The programs here use both ordinrary least squares and Ridge regression with one value only for +the hyperparameter $\lambda$. The first example has no scaling and includes the intercept as well and we are trying to fit a second-order +polynomial. + +!bc pycod +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 + +# Same random numbers for each test. +np.random.seed(2018) +n = 100 +d = 3 +Lambda = 0.01 +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 includes the intercept and scaling is made +X = np.zeros((len(x), d)) +for p in range(d): + X[:, p] = x ** (p) + + +#Split data, no scaling is used and we include the intercept +X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) + + +#Calculate beta, own code +beta_OLS = OLS_fit_beta(X_train, y_train) +beta_Ridge = Ridge_fit_beta(X_train, y_train,Lambda,d) +print(beta_OLS) +print(beta_Ridge) +#predict value +ytilde_test_OLS = X_test @ beta_OLS +ytilde_test_Ridge = X_test @ beta_Ridge + +#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,'*', label="OLS_Fit") +plt.plot(x, X @ beta_Ridge, label="Ridge_Fit") +plt.grid() +plt.legend() +plt.show() + +!ec + + +In this example we do not include the intercept and we scale the data by subtracting the mean values. This follows the discussion in the "lecture material":"https://compphysics.github.io/MachineLearning/doc/LectureNotes/_build/html/chapter3.html#more-on-rescaling-data". +see also the weekly slides "for week 36":"https://compphysics.github.io/MachineLearning/doc/pub/week36/html/._week36-bs029.html". + +Before we discuss the code, we repeat some of the basic math from the slides of week 36. + +Let us try to understand what this may imply mathematically when we +subtract the mean values, also known as *zero centering*. For +simplicity, we will focus on ordinary regression, as done in the above example. + +The cost/loss function for regression is +!bt +\[ +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,. +\] +!et + +Recall also that we use the squared value. This expression can lead to an +increased penalty for higher differences between predicted and +output/target values. + +What we have done is to single out the $\beta_0$ term in the +definition of the mean squared error (MSE). The design matrix $X$ +does in this case not contain any intercept column. When we take the +derivative with respect to $\beta_0$, we want the derivative to obey + +!bt +\[ +\frac{\partial C}{\partial \beta_j} = 0, +\] +!et + +for all $j$. For $\beta_0$ we have + +!bt +\[ +\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). +\] +!et +Multiplying away the constant $2/n$, we obtain +!bt +\[ +\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. +\] +!et + +Let us specialize first to the case where we have only two parameters $\beta_0$ and $\beta_1$. +Our result for $\beta_0$ simplifies then to +!bt +\[ +n\beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} X_{i1} \beta_1. +\] +!et +We obtain then +!bt +\[ +\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \beta_1\frac{1}{n}\sum_{i=0}^{n-1} X_{i1}. +\] +!et +If we define +!bt +\[ +\mu_{\bm{x}_1}=\frac{1}{n}\sum_{i=0}^{n-1} X_{i1}, +\] +!et +and the mean value of the outputs as +!bt +\[ +\mu_y=\frac{1}{n}\sum_{i=0}^{n-1}y_i, +\] +!et +we have +!bt +\[ +\beta_0 = \mu_y - \beta_1\mu_{\bm{x}_1}. +\] +!et +In the general case with more parameters than $\beta_0$ and $\beta_1$, we have +!bt +\[ +\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. +\] +!et + +We can rewrite the latter equation as +!bt +\[ +\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \sum_{j=1}^{p-1} \mu_{\bm{x}_j}\beta_j, +\] +!et +where we have defined +!bt +\[ +\mu_{\bm{x}_j}=\frac{1}{n}\sum_{i=0}^{n-1} X_{ij}, +\] +!et +the mean value for all elements of the column vector $\bm{x}_j$. + + + +Replacing $y_i$ with $y_i - y_i - \overline{\bm{y}}$ and centering also our design matrix results in a cost function (in vector-matrix disguise) +!bt +\[ +C(\boldsymbol{\beta}) = (\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta})^T(\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta}). +\] +!et + + + +If we minimize with respect to $\bm{\beta}$ we have then + +!bt +\[ +\hat{\bm{\beta}} = (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^T\boldsymbol{\tilde{y}}, +\] +!et + +where $\boldsymbol{\tilde{y}} = \boldsymbol{y} - \overline{\bm{y}}$ +and $\tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj}$. + +For Ridge regression we need to add $\lambda \boldsymbol{\beta}^T\boldsymbol{\beta}$ to the cost function and get then +!bt +\[ +\hat{\bm{\beta}} = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}}. +\] +!et + +Now we try to implement this. + +!bc pycod + +np.random.seed(2018) +n = 100 +d = 3 +Lambda = 0.01 +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 does not include the intercept. +X = np.zeros((len(x), d)) +for p in range(d-1): + 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) + +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() + +!ec + +We see that we get the same values for the parameters! As it should be. The MSE may however change (not the case here). + +Finally, instead of using our own function we repeat the same example using the _standardscaler_ functionality of the library _Scikit-Learn_. + +#!bc pycod +np.random.seed(2018) +n = 100 +d = 3 +Lambda = 0.01 +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 does not include the intercept. +X = np.zeros((len(x), d)) +for p in range(d-1): + 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) + +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() + + +#!ec + + + + + + + + + diff --git a/doc/src/week37/programs/codeexamplesscaling.ipynb b/doc/src/week37/programs/codeexamplesscaling.ipynb new file mode 100644 index 000000000..436fd939b --- /dev/null +++ b/doc/src/week37/programs/codeexamplesscaling.ipynb @@ -0,0 +1,649 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "718b0cee", + "metadata": { + "editable": true + }, + "source": [ + "\n", + "" + ] + }, + { + "cell_type": "markdown", + "id": "cdce7555", + "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", + "\n", + "Date: **Sep 11, 2023**\n", + "\n", + "Copyright 1999-2023, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license" + ] + }, + { + "cell_type": "markdown", + "id": "67634fc9", + "metadata": { + "editable": true + }, + "source": [ + "## This note contains code examples with a simple scaling\n", + "\n", + "The programs here use both ordinrary least squares and Ridge regression with one value only for\n", + "the hyperparameter $\\lambda$. The first example has no scaling and includes the intercept as well and we are trying to fit a second-order\n", + "polynomial." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "e75a3307", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from sklearn.linear_model import LinearRegression\n", + "from sklearn.preprocessing import PolynomialFeatures\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn.preprocessing import StandardScaler\n", + "\n", + "def MSE(y_data,y_model):\n", + " n = np.size(y_model)\n", + " return np.sum((y_data-y_model)**2)/n\n", + "\n", + "def OLS_fit_beta(X, y):\n", + " return np.linalg.pinv(X.T @ X) @ X.T @ y\n", + "\n", + "def Ridge_fit_beta(X, y,L,d):\n", + " I = np.eye(d,d)\n", + " return np.linalg.pinv(X.T @ X + L*I) @ X.T @ y\n", + "\n", + "# Same random numbers for each test.\n", + "np.random.seed(2018)\n", + "n = 100\n", + "d = 3\n", + "Lambda = 0.01\n", + "true_beta = [2, 0.5, 3.7]\n", + "\n", + "# Make data set.\n", + "x = np.linspace(-3, 3, n)\n", + "y_real = 2 + 0.5*x + 3.7*x**2\n", + "\n", + "y = np.sum(\n", + " np.asarray([x ** p * b for p, b in enumerate(true_beta)]), \n", + " axis=0) + 0.1 * np.random.normal(size=len(x))\n", + "\n", + "\n", + "#Design matrix X includes the intercept and scaling is made\n", + "X = np.zeros((len(x), d))\n", + "for p in range(d): \n", + " X[:, p] = x ** (p) \n", + "\n", + "\n", + "#Split data, no scaling is used and we include the intercept\n", + "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n", + "\n", + "\n", + "#Calculate beta, own code\n", + "beta_OLS = OLS_fit_beta(X_train, y_train)\n", + "beta_Ridge = Ridge_fit_beta(X_train, y_train,Lambda,d)\n", + "print(beta_OLS)\n", + "print(beta_Ridge)\n", + "#predict value\n", + "ytilde_test_OLS = X_test @ beta_OLS\n", + "ytilde_test_Ridge = X_test @ beta_Ridge\n", + "\n", + "#Calculate MSE\n", + "print(\" \")\n", + "print(\"test MSE of OLS:\")\n", + "print(MSE(y_test,ytilde_test_OLS))\n", + "print(\" \")\n", + "print(\"test MSE of Ridge\")\n", + "print(MSE(y_test,ytilde_test_Ridge))\n", + "\n", + "plt.scatter(x,y,label='Data')\n", + "plt.plot(x, X @ beta_OLS,'*', label=\"OLS_Fit\")\n", + "plt.plot(x, X @ beta_Ridge, label=\"Ridge_Fit\")\n", + "plt.grid()\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "25d58cba", + "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", + "\n", + "Before we discuss the code, we repeat some of the basic math from the slides of week 36.\n", + "\n", + "Let us try to understand what this may imply mathematically when we\n", + "subtract the mean values, also known as *zero centering*. For\n", + "simplicity, we will focus on ordinary regression, as done in the above example.\n", + "\n", + "The cost/loss function for regression is" + ] + }, + { + "cell_type": "markdown", + "id": "5b820b6d", + "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", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "9d3946fb", + "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", + "output/target values.\n", + "\n", + "What we have done is to single out the $\\beta_0$ term in the\n", + "definition of the mean squared error (MSE). The design matrix $X$\n", + "does in this case not contain any intercept column. When we take the\n", + "derivative with respect to $\\beta_0$, we want the derivative to obey" + ] + }, + { + "cell_type": "markdown", + "id": "11228e0d", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\frac{\\partial C}{\\partial \\beta_j} = 0,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "6dd5cd47", + "metadata": { + "editable": true + }, + "source": [ + "for all $j$. For $\\beta_0$ we have" + ] + }, + { + "cell_type": "markdown", + "id": "4bf8184e", + "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", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "c8c7dc9d", + "metadata": { + "editable": true + }, + "source": [ + "Multiplying away the constant $2/n$, we obtain" + ] + }, + { + "cell_type": "markdown", + "id": "3c6a746e", + "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", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "00b74b33", + "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" + ] + }, + { + "cell_type": "markdown", + "id": "20d80ba8", + "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", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "6cf3aa5d", + "metadata": { + "editable": true + }, + "source": [ + "We obtain then" + ] + }, + { + "cell_type": "markdown", + "id": "3492889a", + "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", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "6f6526ea", + "metadata": { + "editable": true + }, + "source": [ + "If we define" + ] + }, + { + "cell_type": "markdown", + "id": "a34541ec", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\mu_{\\boldsymbol{x}_1}=\\frac{1}{n}\\sum_{i=0}^{n-1} X_{i1},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "8d08709c", + "metadata": { + "editable": true + }, + "source": [ + "and the mean value of the outputs as" + ] + }, + { + "cell_type": "markdown", + "id": "0efce920", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\mu_y=\\frac{1}{n}\\sum_{i=0}^{n-1}y_i,\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "5a0488c6", + "metadata": { + "editable": true + }, + "source": [ + "we have" + ] + }, + { + "cell_type": "markdown", + "id": "9727879c", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\beta_0 = \\mu_y - \\beta_1\\mu_{\\boldsymbol{x}_1}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "d0297a4d", + "metadata": { + "editable": true + }, + "source": [ + "In the general case with more parameters than $\\beta_0$ and $\\beta_1$, we have" + ] + }, + { + "cell_type": "markdown", + "id": "5b4f7606", + "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", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "cd6fa027", + "metadata": { + "editable": true + }, + "source": [ + "We can rewrite the latter equation as" + ] + }, + { + "cell_type": "markdown", + "id": "73827116", + "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", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "83069a4e", + "metadata": { + "editable": true + }, + "source": [ + "where we have defined" + ] + }, + { + "cell_type": "markdown", + "id": "1fd7a5df", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\mu_{\\boldsymbol{x}_j}=\\frac{1}{n}\\sum_{i=0}^{n-1} X_{ij},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "e6e597a6", + "metadata": { + "editable": true + }, + "source": [ + "the mean value for all elements of the column vector $\\boldsymbol{x}_j$.\n", + "\n", + "Replacing $y_i$ with $y_i - y_i - \\overline{\\boldsymbol{y}}$ and centering also our design matrix results in a cost function (in vector-matrix disguise)" + ] + }, + { + "cell_type": "markdown", + "id": "6ef1dd83", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "C(\\boldsymbol{\\beta}) = (\\boldsymbol{\\tilde{y}} - \\tilde{X}\\boldsymbol{\\beta})^T(\\boldsymbol{\\tilde{y}} - \\tilde{X}\\boldsymbol{\\beta}).\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "190376f5", + "metadata": { + "editable": true + }, + "source": [ + "If we minimize with respect to $\\boldsymbol{\\beta}$ we have then" + ] + }, + { + "cell_type": "markdown", + "id": "7bfde500", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\hat{\\boldsymbol{\\beta}} = (\\tilde{X}^T\\tilde{X})^{-1}\\tilde{X}^T\\boldsymbol{\\tilde{y}},\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "9fdcab64", + "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", + "\n", + "For Ridge regression we need to add $\\lambda \\boldsymbol{\\beta}^T\\boldsymbol{\\beta}$ to the cost function and get then" + ] + }, + { + "cell_type": "markdown", + "id": "1ce3f20e", + "metadata": { + "editable": true + }, + "source": [ + "$$\n", + "\\hat{\\boldsymbol{\\beta}} = (\\tilde{X}^T\\tilde{X} + \\lambda I)^{-1}\\tilde{X}^T\\boldsymbol{\\tilde{y}}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "0431630c", + "metadata": { + "editable": true + }, + "source": [ + "Now we try to implement this." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "9da6e338", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "\n", + "np.random.seed(2018)\n", + "n = 100\n", + "d = 3\n", + "Lambda = 0.01\n", + "true_beta = [2, 0.5, 3.7]\n", + "\n", + "# Make data set.\n", + "x = np.linspace(-3, 3, n)\n", + "y_real = 2 + 0.5*x + 3.7*x**2\n", + "\n", + "y = np.sum(\n", + " np.asarray([x ** p * b for p, b in enumerate(true_beta)]), \n", + " axis=0) + 0.1 * np.random.normal(size=len(x))\n", + "\n", + "\n", + "#Design matrix X does not include the intercept. \n", + "X = np.zeros((len(x), d))\n", + "for p in range(d-1): \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", + "#Center by removing mean from each feature\n", + "X_train_scaled = X_train - X_train_mean\n", + "X_test_scaled = X_test - X_train_mean\n", + "#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered, note)\n", + "y_scaler = np.mean(y_train)\n", + "y_train_scaled = y_train - y_scaler\n", + "\n", + "\n", + "#Calculate beta\n", + "beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)\n", + "beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d)\n", + "print(beta_OLS)\n", + "print(beta_Ridge)\n", + "\n", + "interceptOLS = y_scaler - X_train_mean @ beta_OLS\n", + "interceptRidge = y_scaler - X_train_mean @ beta_Ridge\n", + "print(interceptOLS)\n", + "print(interceptRidge)\n", + "#predict value\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", + "print(\"test MSE of OLS:\")\n", + "print(MSE(y_test,ytilde_test_OLS))\n", + "print(\" \")\n", + "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,y_real,label='no noise')\n", + "plt.plot(x, X @ beta_OLS+interceptOLS,'*', label=\"OLS_Fit\")\n", + "plt.plot(x, X @ beta_Ridge+interceptRidge, label=\"Ridge_Fit\")\n", + "plt.grid()\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "c2271c6b", + "metadata": { + "editable": true + }, + "source": [ + "We see that we get the same values for the parameters! As it should be. The MSE may however change (not the case here).\n", + "\n", + "Finally, instead of using our own function we repeat the same example using the **standardscaler** functionality of the library **Scikit-Learn**.\n", + "\n", + "\n", + "np.random.seed(2018)\n", + "n = 100\n", + "d = 3\n", + "Lambda = 0.01\n", + "true_beta = [2, 0.5, 3.7]\n", + "\n", + "\n", + "x = np.linspace(-3, 3, n)\n", + "y_real = 2 + 0.5*x + 3.7*x**2\n", + "\n", + "y = np.sum(\n", + " np.asarray([x ** p * b for p, b in enumerate(true_beta)]), \n", + " axis=0) + 0.1 * np.random.normal(size=len(x))\n", + "\n", + "\n", + "X = np.zeros((len(x), d))\n", + "for p in range(d-1): \n", + " X[:, p] = x ** (p+1)\n", + "\n", + "\n", + "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n", + "\n", + "\n", + "\n", + "X_train_mean = np.mean(X_train,axis=0)\n", + "\n", + "X_train_scaled = X_train - X_train_mean\n", + "X_test_scaled = X_test - X_train_mean\n", + "\n", + "y_scaler = np.mean(y_train)\n", + "y_train_scaled = y_train - y_scaler\n", + "\n", + "\n", + "beta_OLS = OLS_fit_beta(X_train_scaled, y_train_scaled)\n", + "beta_Ridge = Ridge_fit_beta(X_train_scaled, y_train_scaled,Lambda,d)\n", + "print(beta_OLS)\n", + "print(beta_Ridge)\n", + "\n", + "interceptOLS = y_scaler - X_train_mean @ beta_OLS\n", + "interceptRidge = y_scaler - X_train_mean @ beta_Ridge\n", + "print(interceptOLS)\n", + "print(interceptRidge)\n", + "\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", + "\n", + "print(\" \")\n", + "print(\"test MSE of OLS:\")\n", + "print(MSE(y_test,ytilde_test_OLS))\n", + "print(\" \")\n", + "print(\"test MSE of Ridge\")\n", + "print(MSE(y_test,ytilde_test_Ridge))\n", + "\n", + "plt.scatter(x,y,label='Data')\n", + "\n", + "plt.plot(x, X @ beta_OLS+interceptOLS,'*', label=\"OLS_Fit\")\n", + "plt.plot(x, X @ beta_Ridge+interceptRidge, label=\"Ridge_Fit\")\n", + "plt.grid()\n", + "plt.legend()\n", + "plt.show()\n", + "\n", + "" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/doc/src/week37/programs/noscaling.py b/doc/src/week37/programs/noscaling.py new file mode 100644 index 000000000..7eb9769e9 --- /dev/null +++ b/doc/src/week37/programs/noscaling.py @@ -0,0 +1,72 @@ +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 = 100 +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) # (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) + + +#Calculate beta, own code +beta_OLS = OLS_fit_beta(X_train, y_train) +beta_Ridge = Ridge_fit_beta(X_train, y_train,L,d) +print(beta_OLS) +print(beta_Ridge) + +#predict value +ytilde_test_OLS = X_test @ beta_OLS +ytilde_test_Ridge = X_test @ beta_Ridge + + +#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,'*', label="OLS_Fit") +plt.plot(x, X @ beta_Ridge, label="Ridge_Fit") +plt.grid() +plt.legend() +plt.show() diff --git a/doc/src/week37/programs/pfaffian.py b/doc/src/week37/programs/pfaffian.py new file mode 100644 index 000000000..b563cbe66 --- /dev/null +++ b/doc/src/week37/programs/pfaffian.py @@ -0,0 +1,10 @@ +from pfapack import pfaffian as pf +import numpy.matlib + +A = numpy.matlib.rand(100, 100) +A = A - A.T +pfa1 = pf.pfaffian(A) +pfa2 = pf.pfaffian(A, method="H") +pfa3 = pf.pfaffian_schur(A) + +print(pfa1, pfa2, pfa3) diff --git a/doc/src/week37/programs/scaling.py b/doc/src/week37/programs/scaling.py new file mode 100644 index 000000000..c1f0c364f --- /dev/null +++ b/doc/src/week37/programs/scaling.py @@ -0,0 +1,86 @@ +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() diff --git a/doc/src/week37/programs/sklearnscaling.py b/doc/src/week37/programs/sklearnscaling.py new file mode 100644 index 000000000..44c51e37b --- /dev/null +++ b/doc/src/week37/programs/sklearnscaling.py @@ -0,0 +1,88 @@ +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 does include 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) + +scaler = StandardScaler() +yscaler = StandardScaler() +scaler.fit(X_train) +yscaler.fit(y_train) +X_train_scaled = scaler.transform(X_train) +X_test_scaled = scaler.transform(X_test) +y_train_scaled = yscaler.transform(y_train) +y_test_scaled = yscaler.transform(y_test) + + + +#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 +ytilde_test_Ridge = X_test_scaled @ beta_Ridge + + +#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,'*', label="OLS_Fit") +plt.plot(x, X @ beta_Ridge, label="Ridge_Fit") +plt.grid() +plt.legend() +plt.show()