From ab0ec583d35f6aee6752b0531edea1631eafb81a Mon Sep 17 00:00:00 2001 From: Kees van Kempen Date: Fri, 18 Feb 2022 18:03:05 +0100 Subject: [PATCH] 07: Compare algos (task 8), explain conditioning, plot that (task 9) --- Week 3/7 Linear Equation Systems.ipynb | 128 ++++++++++++++++++++++--- 1 file changed, 117 insertions(+), 11 deletions(-) diff --git a/Week 3/7 Linear Equation Systems.ipynb b/Week 3/7 Linear Equation Systems.ipynb index 0ef5a4f..51eac6f 100644 --- a/Week 3/7 Linear Equation Systems.ipynb +++ b/Week 3/7 Linear Equation Systems.ipynb @@ -269,7 +269,7 @@ }, "outputs": [], "source": [ - "def GS(A, b, eps, k_max = 10000):\n", + "def GS(A, b, eps, k_max = 1000000):\n", " \"\"\"\n", " Return the Gauss-Seidel algorithm estimate solution x to the problem\n", " Ax = b and the number of iterations k it took to decrease maximum\n", @@ -563,7 +563,7 @@ }, "outputs": [], "source": [ - "def SD(A, b, eps, k_max = 10000):\n", + "def SD(A, b, eps, k_max = 1000000):\n", " \"\"\"\n", " Return the Steepest Descent algorithm estimate solution x to the problem\n", " Ax = b and the number of iterations k it took to decrease maximum\n", @@ -723,7 +723,7 @@ }, "outputs": [], "source": [ - "def CG(A, b, eps, k_max = 10000):\n", + "def CG(A, b, eps, k_max = 1000000):\n", " \"\"\"\n", " Return the Conjugate Gradient algorithm estimate solution x to the problem\n", " Ax = b and the number of iterations k it took to decrease maximum\n", @@ -872,7 +872,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "metadata": { "deletable": false, "nbgrader": { @@ -887,7 +887,20 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEMCAYAAAAvaXplAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA93UlEQVR4nO3dd3hUVf7H8fdJJyEJhCQk9NBJA0LoAhFEUBEUUURUsGFZbPuzgLq2tcuuLpZVdkVEEVRAhQUVFTGggNIJNYAgoaQBCenJzPn9ccOQCglk5s4k39fzzJOZOyWfGcL5zj3n3nOU1hohhBCiLDezAwghhHA+UhyEEEJUIsVBCCFEJVIchBBCVCLFQQghRCVSHIQQQlQixUEIIUQlUhyEEEJU4mF2gIqUUgnA34EdwAKt9arzPSc4OFi3a9fOrrmEEKK+2bhxY4bWOqSq+xxSHJRSs4FRQJrWOrrM9pHAvwB34L9a61cADeQAPkBKTV6/Xbt2bNiwoc5zCyFEfaaUOlTdfY7qVpoDjCy7QSnlDrwDXAFEAhOUUpHAaq31FcDjwHMOyieEEKIMhxQHrXUicKLC5j7APq31Aa11EbAAGKO1tpbefxLwdkQ+IYQQ5Zk55tASOFzmdgrQVyk1FhgBNAHeru7JSqkpwBSANm3a2C+lEEI0QE43IK21XgwsrsHjZimljgFXe3l59bJ/MiGEaDjMPJT1CNC6zO1WpdtqTGu9VGs9JTAwsE6DCSFEQ2dmcfgd6KSUilBKeQE3AktMzCOEEKKUQ4qDUmo+sBboopRKUUrdobUuAaYC3wG7gM+11jtq+bpXK6VmZWVl1X1oIYRwRkV5cPIQpGyA3cvh1OHzP+cCqPqwElx8fLyW8xyEEC5JayjIgtx045KTVvl62W1FOeWfP/otiLv1gn61Umqj1jq+qvucbkC6NpRSVwNXd+zY0ewoQghxltUCeZnVN/YVr1uKqngRBb7NwC8EGodAy15nr/uFnr0e1N4ub0H2HIQQoiZKCmve2Odlgu2UrTLcPKtu4Ku67tsM3O37/b3e7jkIIcRFKzwNWUcgO8X4efpYFQ1/OhRWM7bp6Xe2UW8aAa37VN/w+zQBpRz69i6USxcH6VYSQpxTcT5kH4WsFMg+Ur4InLldVaPfqKnRqDcOhbCYc3zLDwEvP8e/LweQbiUhhGuyFBsNf5WNfmkxyMus/DzfYAhsCQGtSn+2hMBWpT9bQuMw8PBy/PsxgXQrCSFci9UCOanVN/pZR4z7qfDl1jvQaOgDWxoDuBWLQEBL8PQx5S25GpcuDtKtJIQL0tr4Rl9VV8+ZbaePgbWk/PM8fc9+u+94WZlv/WUKgLe/Oe+pHpJuJSFE3SvKhYxkyNhr/Mw6fLbhzz4KJQXlH+/uBQEtKnT1VPjW36ipywzmugrpVhJC1D2tjSN6MvZWuJQWgzOUG/iHGw18eHfociUEti7f3+8bDG6yarEzkeIghDg3SwmcOgTpe842/hml1wvKHOnj6QfBnaBNfwieBCGdIbizcZKWhyzN4mpcujjImIMQdagwBzKTjca/bCE4sb/8GbyNmxuNfvQ44+eZIuDfQr791yMy5iBEQ1KuK2hPmUKQbAwKn6HcISjCaPSDO0Fwl7PXGzUxLb6oWzLmIERDYymBkwerGA+opiuo3cDSItDZKARBEdIV1MBJcRDClZ3pCkqvUAAy94O1+OzjGocZjX/0OAjpcrYQBLSUI4BElaQ4COEqigsg5Tc48DMc2WgUgewyiyfauoK6QOeRpeMBXaBZR+kKErXm0sVBBqRFvWa1wLEtRjH442f4c51xfoByh7BoaHdJaTdQ2aOCGsa0D8L+ZEBaCGehtTE4/MfPRkE4uObspHChURAxGNoPgbYDwSfA3KyiXpABaSGc1anDZ4vBH4mQc9zY3qQtRI2BiCFGUWgcam5O0eBIcRDCkXIzjCJwpiCc/MPY7hdiFIGIIcbeQdN2psYUQoqDEPZUmAOHfj1bDFK3G9u9/I0xg753GwUhtJscNSScihQHIepSSSGk/H52EPnIRmN2UXdvY4WwoU9BRAK06Gn3JSCFuBjy1ynExbBa4Pi2s8Xg0FooyTcmm2vREwY8YHQTte4Lno3MTitEjbl0cZBDWYXDaW1MNfHHz3BglXFEUcEp476QrhB369kjiuTcAuHC5FBWIc4n60j5I4pOHzW2B7aB9oONbqKIweDf3MyUQtSaHMoqRG3knYCDq892FWXuM7b7NqtwRFGEDCKLekuKgxBgHFG0Z7lREI5vBzR4NTa6h+JvLz2iKFKmpBYNhhQHIX75F3z/tLFUZas+cOkTRjFoGQfunmanE8IUUhxEw6U1rHwBVs+AqLEw5m3w8jM7lRBOQYqDaJisVvh2Gvz2vnGE0ag3wc3d7FRCOA0pDqLhsZTAkvth66fQfypc/oIMLAtRgRQH0bCUFMKiO2DXUrj0SRj8qBQGIarglIdeKKX8lFIblFKjzM4i6pGiXJh/o1EYRr4CQx6TwiBENRxSHJRSs5VSaUqppArbRyql9iil9imlppW563Hgc0dkEw1E/in4eKxxVvOYd6DfvWYnEsKpOWrPYQ4wsuwGpZQ78A5wBRAJTFBKRSqlhgM7gTQHZRP1XU46fDTKmARv3IfQ82azEwnh9Bwy5qC1TlRKtauwuQ+wT2t9AEAptQAYAzQG/DAKRr5SarnW2lrxNZVSU4ApAG3atLFjeuHSso7A3DGQlQITFkCny8xOJIRLMHNAuiVwuMztFKCv1noqgFJqMpBRVWEA0FrPAmaBMbeSfaMKl5S5H+ZeY0yMd8tiaDvA7ERCuAynPVpJaz3nfI+RWVlFtVJ3wsfXGGspTFoKLXqYnUgIl2Lm0UpHgNZlbrcq3VZjWuulWuspgYGBdRpMuLiUjTDnSmNNhdu+kcIgxAUwszj8DnRSSkUopbyAG4EltXkBpdTVSqlZWVlZdgkoXNAfq2HuaPAJhNu/hZAuZicSwiU56lDW+cBaoItSKkUpdYfWugSYCnwH7AI+11rvqM3ryp6DKGfPt/DJdRDYCm77Fpq2MzuREC7LUUcrTahm+3JguSMyiHpu+0L48m4Ii4GJi8CvmdmJhHBpTnmGdE1Jt5IAYMOHsOhOY53mW5dIYRCiDrh0cZBuJcEvM+F/D0Gn4XDzIvAJMDuREPWCSxcH2XNowM6sxfD93yDqWhg/DzwbmZ1KiHrDpYuD7Dk0UFYrfPM4JL5urMVw3Qfg4WV2KiHqFac9CU6IKllKYOkDsGWerMUghB1JcRCuo6TQGHjetQQSnpApt4WwI5cuDjJ9RgNSlAuf3Qz7VxprMciU20LYlYw5COdXkHV2LYbRb0thEMIBXHrPQTQAuRnwyVhjIr1xs40jk4QQdufSew5yKGs9l30UPrwC0vcaazFIYRDCYVy6OEi3Uj124gDMHgHZx4y1GGSRHiEcSrqVhPNJ3QkfXwuWIpi0BFrGmZ1IiAbHpfccRD10pHQtBjDWYpDCIIQppDgI5/HHavhoNHgHGGsxhHY1O5EQDZYUB+Ec9n4H88YZazHc/i0ERZidSIgGzaWLgxytVE8kLYIFN0FIV5i8HAJamJ1IiAbPpYuDHK1UD2ycAwvvMNZimLRU1mIQwkm4dHEQLu7Xt2Dpg9DxMpi4UNZiEMKJyKGswvG0hp9egsTXIPIaGPsfmXJbCCcjxUE4ltUK302H9e9Bz1vg6n+Bm7vZqYQQFUhxEI5Tdi2Gfn+BES/KlNtCOCkpDsIxSopg8Z2w82tImA5DHpfCIIQTc+niIOs5uIiiPPj8Ftj3A4x4Cfr/xexEQojzcOmjleRQVhdQkAWfXAf7foTRb0lhEMJFuPSeg3ByuZnwybWQusNYiyF6rNmJhBA1JMVB2Ef2MZg7Bk4dghvnQ+fLzU4khKgFKQ7CPpY/AlkpcPMiaHeJ2WmEELXk0mMOwkkd2QS7/weXPCSFQQgXJcVB1L2fXoRGQdD3HrOTCCEukBQHUbcOrTUOWb3kIZkrSQgX5nTFQSnVTSn1nlJqoVLqXrPziFrQGla+AH6h0Psus9MIIS6CQ4qDUmq2UipNKZVUYftIpdQepdQ+pdQ0AK31Lq31PcANwEBH5BN15I+f4dAaGPwIePmanUYIcREctecwBxhZdoNSyh14B7gCiAQmKKUiS+8bDSwDljson7hYZ/YaAlpBr8lmpxFCXCSHFAetdSJwosLmPsA+rfUBrXURsAAYU/r4JVrrK4CJjsgn6kDyCkj53dhr8PA2O40Q4iKZeZ5DS+BwmdspQF+lVAIwFvDmHHsOSqkpwBSANm3a2C2kqIEzew1N20HPm81OI4SoA053EpzWehWwqgaPmwXMAoiPj9f2TSXOaddSOL4NrnkP3D3NTiOEqANmHq10BGhd5nar0m01ppS6Wik1Kysrq06DiVqwWoxV3YI7Q+wNZqcRQtQRM4vD70AnpVSEUsoLuBFYUpsXkFlZnUDSYkjfZazRICu6CVFvOOpQ1vnAWqCLUipFKXWH1roEmAp8B+wCPtda76jl68qeg5ksJbDqJWgebawFLYSoNxwy5qC1nlDN9uVcxOGqWuulwNL4+Hg548oMW+fDiQPGrKtuTnc+pRDiIrj0/2jZczBRSSH8/Cq0iIMuV5idRghRx1y6OMiYg4k2zYWswzD0KVkLWoh6yKWLgzBJcT4kzoA2A6DDULPTCCHswKWLg3QrmeT3DyDnOAx9UvYahKinXLo4SLeSCQpzYM0b0D5BFvIRoh5z6eIgTPDb+5CXAZc+ZXYSIYQduXRxkG4lB8s/Bb/8CzqPhNa9zU4jhLAjly4O0q3kYOvehYIsuPQJs5MIIezMpYuDcKDcTFj7LkSOgfDuZqcRQtiZFAdRM7/+C4pyIEH2GoRoCFy6OMiYg4OcToX1s4xZV0O7mp1GCOEALl0cZMzBQdb8EyxFMORxs5MIIRzEpYuDcICsFNgwG3rcBM06mJ1GCOEgUhzEuSXOMJYBHfKY2UmEEA4kxUFU78QfsPlj6DUZmsg63UI0JC5dHGRA2s5+fg3cPGDQ/5mdRAjhYDUuDkqpF6rYZuq6kDIgbUfpe2HbAuh9JwSEm51GCOFgtdlzaKmUsq3oppQKBX6o+0jCKax6GTwawSUPm51ECGGC2iwTejfwnVJqP6CBDwE5trE+Op4EOxbDoEfAL9jsNEIIE5y3OCil5gKbgM3AX4BPgRLgGq31PvvGE6b46SXwDoQBU81OIoQwSU26leYACrgN+ARoB5wEblZKjbNbMmGOIxthzzIYcD80amp2GiGESc6756C1XgmsPHNbKeUBdAO6A32BhXZLJxxv5YvQKAj63WN2EiGEiWoz5gCA1roE2F56+aTOE9WCUupq4OqOHTuaGaP+OLQW9v8Iw58Hb3+z0wghTOTS5znIoax1SGtY+QL4hULvu8xOI4QwmUsXB1GH/vgZDq2BwY+Al6/ZaYQQJpPiIM7uNQS0MqbKEEI0eFIcBCSvgJTfYcij4OFtdhohhBOQ4tDQWa3GXkPTdtBjotlphBBOotZHK4l6ZvdSOL4Nrn0f3D3NTiOEcBKy59CQWS3G2dDBnSHmerPTCCGciOw5NGRJiyB9N1w/B9xMnWBXCOFknLI4KKWuAa4CAoAPtNYrzE1UD1lKjJlXm8dAtzFmpxFCOBmHdSsppWYrpdKUUkkVto9USu1RSu1TSk0D0Fp/pbW+C7gHGO+ojA3K1vlw4gBc+gS4Se+iEKI8R7YKc4CRZTeULhb0DnAFEAlMUEpFlnnIU6X3i7pUUgg/vwot4qDLFWanEUI4IYcVB611InCiwuY+wD6t9QGtdRGwABijDK8C32itNzkqY4OxaS5kHYahT4FSZqcRQjghs/sTWgKHy9xOKd12P3AZME4pVeX0oEqpKUqpDUqpDenp6fZPWl8U50PiDGgzADoMNTuNEMJJOeWAtNZ6JjDzPI+ZBcwCiI+P147IVS/8/gHkHIdxs2WvQQhRLbP3HI4ArcvcblW6rUaUUlcrpWZlZWXVebB6qTAH1vwT2l8K7QaanUYI4cTMLg6/A52UUhFKKS/gRmBJTZ8sU3bX0vr3IC/TGGsQQohzcOShrPOBtUAXpVSKUuqO0oWDpgLfAbuAz7XWO2rxmrLnUFP5p+DXmdD5CmgVb3YaIYSTU1q7fnd9fHy83rBhg9kxnNtPLxmHr969GsJjzU4jhHACSqmNWusqvy2a3a10UWTPoYZyM2HtuxB5jRQGIUSNuHRxkDGHGvr1X1CUAwnTzU4ihHARTnkoq6hDp1Nh/SyIvQFCu5qdRghRBatVk1tUQm6hhZzCEnIKS8gt89O4brFtK3v/lMHtGdQppM4zuXRxUEpdDVzdsWNHs6M4rzX/BEsRDHnc7CRC1CuFJRZyCs426LlF5Rvz06X3VdxuXC9fBPKKLDX6nW4KGnt70NjbA7/SS1GJ1S7vz6WLg9Z6KbA0Pj7+LrOzOKWsFNgwG3pOhGYdzE4jhKlKLFajUS4621DbGvbCEvKKyn87zy1t8I1GvnwRyC0sodhSs4N5vD3c8Pcpbcy9jIY9uLEXbZv5Gtu9jPvONvju5QpA2Z8+nm4oB5286tLFQZxH4uvGz8GPmZtDiAtQbLHaGuq8Iku5b9+Vv5FbbI15TqGFvDMNfJmumpp+w3ZTGI24z9lv54293Qlu7G002j5lGm0v9/INuE+Z615GQ+/h7ppDuy5dHKRb6RxOHIDNn0D87dCk9fkfL0Qd0FqTW2QhK7+YrLxiTuUXkZ1fvo88r6j8N/bcIku5LpczhaBWjbn3mW/g7rbGuZWvL429yzfevl7u5b6J+1Zo3P283Wnk6e6wb+cXIy0vje3p24kKjiLML6zOX9+li4N0K53Dz6+BmwcM+j+zkwgXo7Umv9ho4E/lFRsNfeklu/RnVdtPlf4ssZ67u+VMY16239zPy50gP9/Sbe62rpYz39p9vSp3u5zZ5siuFrMUlBSw68QutqVvMy4Z2zieexyAp/o+xfiudb/sjUsXB1GN9L2w7TPo/xfwr/tvFMI1FJQ28FU15ll5ReVvl2ncs/KLz9mfrhQE+HjSxNeTwEbGpWXTRgQ28qRJo7PbzlwCGnna+twbe3vg7VH/G/OLobXmUPYhtmVssxWD5JPJlOgSAFo2bknPkJ7ERMYQGxJL1yD7HIUoxaE+WvUyePrCwIfMTiLqQEGxhRO5ReUuJ0sb91N5Zxv0so18Vn7xebtlAnw8CCzTwIcF+hDYyKtS4162CAQ08sTf2wM3N2nc60pWYRbbM7bb9gi2p28nuygbAD9PP6KDo7kt+jZigmOICYkhuFGwQ3K5dHGQMYcqHE+CHYth0CPg55g/IlFzVqsmu6C4UmOfmVvEyTO388rfd67DHP29PQgo05B3DG1criGv2Lifufj7eOIuDbzDFVuL2XtyL9vSjSKwLWMbh7IPAeCm3OjQpAPD2w4nNiSW2OBYIgIjcHdzNyWrzK1U38y/CQ6ugYe2QqOmZqep9wpLLJzMLdPY5xVxIqewykbe+MZfjKWaPvlGnu4E+XnV6NLU14sAHw+XPRKmIdBaczz3uK17aHvGdnZm7qTQUghAcKNgYoNjiQmJITY4lqjgKPw8/Rya8VxzK7n0noOo4MhG2LMMLn1KCsMF0FpzurCEEzlnGvmqG/iyl5zCkipfSylo0sjT1phHBPvRq23T0tveBPl5Gj99vQhq7EWQrxeNvMz5hijqRl5xHjsyd9jGCbZnbCc931il0tvdm25B3RjfZTwxITF0D+5OmF+YU4+9SHGoT1a+CI2CoF+VK6s2WBarJjO3kLTsQtJOF5CabVxPPV1AWnYBaacLSc0u4ERuUbUDsV4ebjQr8829bTNfmvp60czPi6Z+XuXuC/Lzoomvl3Tb1GNWbeXAqQPGWEHpnsG+U/uwamOcp21AW/qG97V1D3Vu2hlPd0+TU9eOFIf64tCvsP9HGP538PY3O41DWK2aE3lFpGYXGI19mYY+NbuQ9NJCkJ5TWGVXTpCfF6H+3jQP8KFLc3+aNfautrH39XKNY9+FfWTmZ5YbNN6RsYOc4hwA/L38iQ2OZVibYcagcXAMTXyamBu4Drh0cZAB6VJaw8oXoHFz6H2n2WkumtWqOZlXRGrZb/el11OzC0k7XUhadgHppwurPKa+qa8nzQN8CA3woVNzf5oHGAUg1N+H0NLrIY298fKQ/npRWZGliN0ndtsKwbb0bRzJMVYvdlfudG7amavaX0VsSCwxwTG0DWiLm6p/f0suXRzkJLhSB1bBoV/gitfBy9fsNNXSWnMyr7j0m32BrZFPLfOt/0w3T1WNfhNfT5qXNvAdQ4JtjX7zAG9C/M/89MbbQ/ruzZaRn8G6Y+tYd3QdG1M3UmApMDtSjWUVZlFsLQYgzC+MmOAYbuxyI7EhsXRr1o1GHo1MTugYLl0cBGf3GgJaQa9JpsUoLLFwPKuAY1kFHMvK51hWAcezCsr07Rv9/VX16Qc28rQ19O1DmhkNfml3T2iAN6H+PoT4e+PjKY2+s8orzmND6gajIBxbR/LJZAACvALoE9bHpbpZArwCbEcRhfqGmh3HNFIcXF3yCjiyAa6eCR7edvkVhSUWUrMKbY1+2QJwLCufY6cKyMwtqvS8AB+P0m/2PvSN8CO09Fu+0cVT2r0jjb5LKrGWsCNzB+uOrmPtsbVsTd9KibUELzcvejbvyYNxD9K/RX+6Nu1q2nH64uJIcXBlVqux19A0AnrcdEEvUVRiJTW7QoN/qmwRKCAjp7DS8wJ8PAgPbER4Ex9iWjYhPNCn9GJsCwvwwc9b/rzqC601B7MP2rqKfj/+O6eLTwPQLagbt0TeQr/wfsSFxuHj4WNyWlEX5H+vK9u9FI5vg2vfhyoOkyu2nG34j57Kr9Ttc6bhr3gepL+Ph62hj2oRYGvwbY1/oDT8DUFGfgbrj623dRWdmeithV8LLm93Of1a9KNvWF+a+sg5NfWR/A93UcXFxfDDC5QEduAHPZBjifs5eqqgtAAYjX96VQ2/twdhgT6EN2lEt7AAwpv40CKwEWGBPrRo4kNYYCMaS8PfIOUV57EpbZOtq2jvyb2A0QffN7wvd8XcRf/w/rTybyWH9TYALt0K1PdDWbXWpJ0uZOfRbHYey2bXsWwOn8zneFY+A3J/5A3PvTxY9ADLF2wDwM/LnfAmxjf7LmH+tm/54U0a0SLQh7BAH/x9XOtEHGE/FqvFGDco3TPYkraFYmsxnm6exIXG8WDcg/QL70e3oG4ybtAAydxKTqLEYuVARi47jxpFYOexbHYezS430Ns6qBHtmvnRwt+DaQcmgacvW69aQngTP8Kb+BAgDb84B601f57+k7VH17Lu2Dp+O/4bp4uMcYOuQV3pF96P/uH96dm8Z4M5XLOhk7mVnMzpgmJ2Hz9tFIHSvYLdx0/bplj2cnejc1hjhnULJTI8gMgWgXQN9z/b+G+aCzsPw7ULSOgi6zWI6p0oOGEbN1h7dC3Hco8BEO4XzvC2w+kX3o++4X0J8gkyOalwNlIc7EhrzbGsgnJFYOexbA5l5tke09TXk8gWAUzq35bIFgFEhgfSPsQPz+pm2ywpNFZ5a9kLOo900DsRriK/JJ9NqZtsXUW7T+wGjCke+ob15Y7oO+jXoh9t/NvIuIE4JykOdaTYYmVfWk6lQnAqr9j2mHbNfIlqEcD1vVoR2SKAbuEBhAX41O4/6aa5kHUYRs80pv4UDZrFamHXiV22rqLNaZspthbj4eZBz9Ce3N/zfvqH9yeyWaSMG4hakeJwAbLyi9ldZlxg57FsklNzKLIY3ULeHm50DfPniugwIsONItA1PODijwIqzofEGdB2ILS/tA7eiSi2FnM85zgpOSm2KRNcwfHc46w9upb1x9fbxg26NO3CTV1vol8L43wDX0/nnUpFOD8pDuegtSblZH65AeKdx7JJOZlve0wzPy8iWwRw2yXtjPGB8AAigv3sswjL7x9AznEYN1v2GmohrziPw6cPk3I6hcOnD5e7HMs9hkVXv9KaM2vu25xhbYbRP7w/fcL7OGz5SNEwSHEoVVRiJTnt9NkuodKjhrILjMVclIKIYD+6t27ChD5tiGwRQFR4ACH+3o7pu83cD6v/AR2GQruB9v99LkRrzYmCE7YGv2IRyCzILPf4QO9AWjduTUxwDFdEXEFr/9a08m/lUkfoBHgF0Nq/tYwbCLtp0MXh573pfL3lCDuPZrMvLcc2E2gjT3e6hPkzqnuL0qOFAuga5o+vl0kf16k/Ye4Yo0KNfNWcDCYrsZZwPPd4tQUgr+TsIL9C0dyvOa39WzOk9RBb49/avzWt/VsT4BVg4jsRwjU4XXFQSrUHngQCtdbj7Pm7klNPszo5g8jwAC7tGmorBO2a+TnPKl7Zx+Cj0VCYDZOWQkhnsxPZTX5JfqVG/8ztozlHKdFnl+T0dPO0NfjxYfG2hr+VfytaNm6Jt7t9JiEUoqFwyElwSqnZwCggTWsdXWb7SOBfgDvwX631K2XuW1jT4nChJ8FZrNp5ikBVcjPgwysh+wjc8hW07m12oouitSarMKtSv/+ZIpCWn1bu8f5e/rZGv+Il1De0Xi6wUheKi4tJSUmhoMB11lAQ9uXj40OrVq3w9Cx/oqwznAQ3B3gbmFsmlDvwDjAcSAF+V0ot0VrvdFAm5y4M+Sdh7jVGl9LNC12qMBRZithzYg97T+6tVADOzOR5RmijUFr5t2JAywGVCkCgd6BJ78C1paSk4O/vT7t27WRMQqC1JjMzk5SUFCIiImr8PIcUB611olKqXYXNfYB9WusDAEqpBcAYoEbFQSk1BZgC0KZNm7oL6wwKsuGT6yBjD0yYD+0uMTtRtSxWC/uz9rMjYwdJGUkkZSax9+ReSqxGF5CHmwctG7eklX8ruod0L9f4t/JvJdM720FBQYEUBmGjlKJZs2akp6fX6nlmjjm0BA6XuZ0C9FVKNQNeBHoqpaZrrV+u6sla61nALDC6lewd1mGK8uDT8XB0C4z/GDpeZnYiG601KadTSMpMMgpBRhK7Tuwiv8Q4tLexZ2OigqOYFDmJ6OBougZ1JdwvXE6+MoEUBlHWhfw9ON2AtNY6E7inJo+td7OyFhfAgpvg8Dq47r/Q9SpT46Tnpdv2BnZk7CApM4mswiwAvN296RrUlbGdxhLVLIro4Oh6u9C6EA2RmcXhCNC6zO1WpdtqTGu9FFgaHx9/V10GM4WlGL6YDAd+gjHvQvR1Dv312UXZ7MjYwY5Mo3toe8Z20vKMAWJ35U7HJh25rM1lRAVHEd0smo5NO+LpJrPAiuqlpqby8MMPs27dOpo2bYqXlxePPfYY1157rcMyrFu3jgcffJDCwkIKCwsZP348zz77bLWP37BhA3PnzmXmzJmV7mvXrh0bNmwgOLj6kw1r8hhXYWZx+B3opJSKwCgKNwK1Wuuy3uw5WC2w+C7Y+w1cOQN6TrTrr8svyWfPiT1sz9hOUkYSOzJ3cCj7kO3+tgFtiW8eT3RwtK17yJVOEBPm01pzzTXXMGnSJD799FMADh06xJIlSxyaY9KkSXz++ed0794di8XCnj17zvn4+Ph44uOrPHinwXFIcVBKzQcSgGClVArwjNb6A6XUVOA7jENZZ2utd9TmdevFnoPVCl9PhR1fwvC/Q5+6fSvF1mL2n9rP9ozttkHjfaf22aaMCPUNJbpZNGM6jCEqOIqoZlFylJC4aCtXrsTLy4t77jnbQ9y2bVvuv/9+Dh48yC233EJubi4Ab7/9NgMGDGDVqlXMmDGD//3vfwBMnTqV+Ph4Jk+ezLRp01iyZAkeHh5cfvnlzJgxgy+++ILnnnsOd3d3AgMDSUxMrJQjLS2N8PBwANzd3YmMjAQgNzeX+++/n6SkJIqLi3n22WcZM2ZMuQyZmZlMmDCBI0eO0L9/f8oe9v/JJ58wc+ZMioqK6Nu3L++++y7u7vVrbM1RRytNqGb7cmC5IzI4Ja1h+SOw9VNIeAIGPnBRL2fVVg5lH7LtDSRlJLH7xG4KLYWAMeVCdHA0g1sNtu0VhPqG1sU7EU7quaU72Hk0u05fM7JFAM9cHXXOx+zYsYO4uLgq7wsNDeX777/Hx8eH5ORkJkyYwLnOU8rMzOTLL79k9+7dKKU4deoUAM8//zzfffcdLVu2tG2r6OGHH6ZLly4kJCQwcuRIJk2ahI+PDy+++CJDhw5l9uzZnDp1ij59+nDZZeUP/njuuee45JJLePrpp1m2bBkffPABALt27eKzzz7jl19+wdPTk/vuu4958+Zx6623nvMzcTVONyBdGy7draQ1rHgKNnwAAx+EIY/V8uma1LxU21FDSZlJ7MzYaTuPoJFHI7oFdeOGLjcQ3SyamOAYWftXmOYvf/kLa9aswcvLix9++IGpU6eyZcsW3N3d2bt37zmfGxgYiI+PD3fccQejRo1i1KhRAAwcOJDJkydzww03MHbs2Cqf+/TTTzNx4kRWrFjBp59+yvz581m1ahUrVqxgyZIlzJgxAzAO//3zzz/LPTcxMZHFixcDcNVVV9G0aVMAfvzxRzZu3Ejv3sa5R/n5+YSG1r8vWS5dHFy6W2nVy7D2begzBS57rkazrG5J28K6Y+tsBeHMhHIeyoNOTTsxMmKkbY+gfWB7PNxc+p9X1IHzfcO3l6ioKBYtWmS7/c4775CRkUF8fDxvvPEGzZs3Z+vWrVitVnx8jHNdPDw8sFqttuecOcPbw8OD3377jR9//JGFCxfy9ttvs3LlSt577z3Wr1/PsmXL6NWrFxs3buSRRx5h8+bNtGjRguXLjU6JDh06cO+993LXXXcREhJCZmYmWmsWLVpEly5dyuVOTU0973vTWjNp0iRefrnKo+zrDZc+7lApdbVSalZWVpbZUWpnzRvw86vQ82ZjIr3zFIaCkgL+vvbv3PLNLby75V3+PP0nA1oMYFqfacy7ch7rJq7j86s/5+n+TzO201g6N+0shUGYaujQoRQUFPDvf//bti0vz5gcMSsri/DwcNzc3Pj444+xWIzxr7Zt27Jz504KCws5deoUP/74IwA5OTlkZWVx5ZVX8sYbb7B161YA9u/fT9++fXn++ecJCQnh8OHDfPjhh2zZssVWGJYtW2YbK0hOTsbd3Z0mTZowYsQI3nrrLdt9mzdvrvQeBg8ebBtM/+abbzh58iQAw4YNY+HChaSlGUfznThxgkOHDlV6vqtz6RbEJfcc1r8PPzwL0ePg6pngdu76vO/kPh5NfJR9p/YxOWoyU2Kn4O/l75isQlwgpRRfffUVDz/8MK+99hohISH4+fnx6quvEhcXx3XXXcfcuXMZOXIkfn5+ALRu3ZobbriB6OhoIiIi6NmzJwCnT59mzJgxFBQUoLXmn//8JwCPPvooycnJaK0ZNmwY3bt3r5Tj448/5uGHH8bX1xcPDw/mzZuHu7s7f/vb33jooYeIjY3FarUSERFhGwg/45lnnmHChAlERUUxYMAA20wMkZGRvPDCC1x++eVYrVY8PT155513aNu2rT0/UodzyMR79nahE+853Ka5sOR+6DoKrp8D7tWfJ6C15ou9X/Da76/h5+nHS5e8xMCWso6DOL9du3bRrVs3s2MIJ1PV34UzTLwntn0BSx4wpsMYN/uchSGrMIvn1j7H94e+p394f14a9JKs8iWEcCiXLg4uc7TSrqXw5d3GBHrjPwGP6tca2JS6icdXP05GXgZ/7fVXJkVNkikphBAO59KtjtZ6qdZ6SmCgE5+0lfw9fHEbtOxlzLDqWfWZxharhfe2vsdt392Gh/Lg4ys/5rbo26QwCCFM4dJ7Dk7vj0T47GZoHgkTvwDvqgeSj+ceZ/rq6WxI3cCVEVfyt35/o7FXYweHFUKIs6Q42Muf6+HTG6FpBNz8JTRqUuXDfvrzJ/72698oshTxwsAXGN1htJyoJoQwnUsXB6cdczi6GeaNA/8wuPUr8GtW6SGFlkL+seEfzN89n25B3Xht8Gu0C2zn8KhCCFEVl+7Qdsoxh9Sd8PG14NMEJi0xCkQFB04d4KZlNzF/93xuibyFT678RAqDqFdefPFFoqKiiI2NpUePHqxfvx6AN99803YynFlOnTrFu+++W+397u7u9OjRg6ioKLp3784//vGPcmdum+Gll15y/C/VWrv8pVevXtoppCdr/VpHrWd00Tpzf6W7rVarXrR3ke79SW89aP4g/fPhn00IKeq7nTt3mvr7f/31V92vXz9dUFCgtdY6PT1dHzlyRGutddu2bXV6erqZ8fQff/yho6Kiqr3fz8/Pdj01NVUPGzZMP/30046IVq2ymS5UVX8XwAZdTbvq0nsOTuXkIZg7GrQVbl0CQe3L3Z1dlM1jiY/xzK/PEBscy8LRCxncarBJYYWwn2PHjhEcHIy3t3HIdnBwMC1atGDmzJkcPXqUSy+9lEsvvRSAFStW0L9/f+Li4rj++uvJyckBYOPGjQwZMoRevXoxYsQIjh07BkBCQgIPPvggPXr0IDo6mt9++w0wpuC+/fbb6dOnDz179uTrr78GjNlh+/TpQ48ePYiNjSU5OZlp06axf/9+evTowaOPPnrO9xIaGsqsWbN4++230VpjsVh49NFH6d27N7Gxsbz//vu29zx48GBbrtWrVwPw7bffEhcXR/fu3Rk2bNg5s86ZM4exY8cycuRIOnXqxGOPGZNxTps2jfz8fHr06MHEifZd66UsOUO6LmQfhdkjoSALJv8PwmLK3b0lbQvTVk/jeO5xpvacym1Rt8m6ysJuyp0J+800OL69bn9BWAxc8Uq1d+fk5HDJJZeQl5fHZZddxvjx4xkyZAhQfqW0jIwMxo4dyzfffGObWqOwsJDp06czZMgQvv76a0JCQvjss8/47rvvmD17NgkJCXTq1In//Oc/JCYmct9995GUlMQTTzxBZGQkN998s20K7s2bNzNt2jT69evHxIkTKSoqwmKxkJqayqhRo0hKSqoyf+PGjW1F6owmTZqwZ88evv76a9LS0njqqacoLCxk4MCBfPHFFyxevJiCggKefPJJLBYLeXl5FBQUEBcXR2JiIhEREZw4cYKgoKBqs37xxRc8//zzbN68GW9vb7p06cKaNWto3bp1lZlqq0GdIe0UA9I5afDRaMg7AZO+LlcYLFYLH+74kLc3v02YXxgfXfER3UMqz/8iRH3SuHFjNm7cyOrVq/npp58YP348r7zyCpMnTy73uHXr1rFz504GDjSmhSkqKqJ///7s2bOHpKQkhg8fDoDFYrEt2AMwYYKxPMzgwYPJzs7m1KlT1U7B3b9/f1588UVSUlIYO3YsnTp1uqj3tmLFCrZt28bChQsBYxLB5ORkevfuze23305xcTHXXHMNPXr0YNWqVQwePJiIiAgAgoKCbK9R3XThw4YN48wYamRkJIcOHaJ169YVYziESxcHbfbEe3knYO41kH0Ebl5snOhWKi0vjSdWP8H64+sZ2W4kT/d/WibME453jm/49uTu7k5CQgIJCQnExMTw0UcfVSoOWmuGDx/O/Pnzy23fvn07UVFRrF27tsrXrniot1Kq2im4u3XrRt++fVm2bBlXXnkl77//Pu3bl+/yPZ8DBw7g7u5OaGgoWmveeustRowYUelxiYmJLFu2jMmTJ/PXv/7Vtv5DRdVlXb9+va0rDozPsKSkpFZZ65KMOVyogiz4ZCxk7oMbP4W2/W13JaYkMm7JOLZlbOP5Ac/z2uDXpDCIBmPPnj0kJyfbbm/ZssU2Y6m/vz+nTxsLUvXr149ffvmFffv2AUZf/N69e+nSpQvp6em24lBcXMyOHWdXEP7ss88AWLNmDYGBgQQGBlY7BfeBAwdo3749DzzwAGPGjGHbtm3lMpxPeno699xzD1OnTkUpxYgRI/j3v/9NcXExAHv37iU3N5dDhw7RvHlz7rrrLu688042bdpEv379SExM5I8//gCMqb2BGk0XXpGnp6ftdzqKS+85mKYoF+bdYPTljp8HHYzBtSJLEW9sfINPdn1C56adeX3w67RvUrtvKUK4upycHO6//35OnTqFh4cHHTt2ZNasWQBMmTKFkSNH0qJFC3766SfmzJnDhAkTKCw0lrJ94YUX6Ny5MwsXLuSBBx4gKyuLkpISHnroIaKijIWLfHx86NmzJ8XFxcyePRug2im4P//8cz7++GM8PT0JCwvjiSeeICgoiIEDBxIdHc0VV1zB66+/Xi7/mcHf4uJiPDw8uOWWW/jrX/8KwJ133snBgweJi4tDa01ISAhfffUVq1at4vXXX8fT05PGjRszd+5cQkJCmDVrFmPHjsVqtdqWR63JdOEVTZkyhdjYWOLi4pg3b16d/ntVRwaka6u4AD69AQ6uNmZXjboWgINZB3ks8TF2ndjFTV1v4q/xf8XbvfoJ9oSwl/o8ZXdCQgIzZswgPr7KMVRxDg1qQNrhSorg81uNOZOufQ+irkVrzdf7v+al9S/h7e7NW0PfIqF1gtlJhRDiokhxqClLCSy6A5K/g1FvQPcbySnK4e/r/s7yP5bTO6w3L1/yMs39mpudVIh6a9WqVWZHaDBcujg47FBWqxW+vg92LYERL0P87WxP385jiY9xLPcYU3tM5c6YO+XcBSFEveHSRytpR8ytpDUsexi2fQZDn8La7x5mJ83m1m9uxaItfDjyQ+7ufrcUBiFEveLSew52pzV8Ox02zoFB/0dG79t44vt7WHtsLcPbDueZ/s8Q6O1Ek/4JIUQdkeJwLiv/Duv/DX3v5ZcuQ3liyXXkFufydP+nGddpnKy7IISot1y6W8muEl+H1f+gOO5W/hEczD0/3kuQTxALrlrA9Z2vl8IgxDkcP36cG2+8kQ4dOtCrVy+uvPJK9u7de0Gvdeedd7Jz5846zXe+aburk5OTw7333kuHDh2Ii4ujV69e/Oc//7moLHPmzGHq1KkAvPfee8ydO/eCXufgwYN8+umnF5WlLCkOVVn7Lqx8gT+jR3OLWzpzdn7E+C7jmX/VfDo2dbKFhYRwMlprrr32WhISEti/fz8bN27k5ZdfJjU19YJe77///S+RkZF1mvFCi8Odd95J06ZNSU5OZtOmTXz77be2M5/LutBpL+655x5uvfXWC3quFAd72zAbvpvO0s6XcH3hXg6fPsybCW/yVL+n8PHwMTudEE7vp59+wtPTk3vuuce2rXv37gwaNAitNY8++ijR0dHExMTYpsJYtWoVCQkJjBs3jq5duzJx4kTb9BIJCQmcOcm1ceOza6svXLjQNl/T/v376devHzExMTz11FO2x+Xk5DBs2DDi4uKIiYmxTY9d1bTdr7/+um0q7meeeabS+9q/fz+//fYbL7zwAm5uRtMZEhLC448/bnsPgwYNYvTo0bZids0119CrVy+ioqJsZ4kDfPjhh3Tu3Jk+ffrwyy+/2LY/++yztgn59u/fz8iRI+nVqxeDBg1i9+7dAEyePJkHHniAAQMG0L59e9skgNOmTWP16tX06NGDN954o5b/apXJmENZWxeQu+z/eKl9DEuK/yQuNI5XBr1CeOPw8z9XCCf06m+vsvvE7jp9za5BXXm8z+PV3p+UlESvXr2qvG/x4sVs2bKFrVu3kpGRQe/evRk82FjXZPPmzezYsYMWLVowcOBAfvnlFy655JIaZXrwwQd58MEHmTBhAu+9955tu4+PD19++SUBAQFkZGTQr18/Ro8ezSuvvEJSUhJbtmwBjJlSk5OT+e2339BaM3r0aBITE23ZwFgbonv37rbCUJVNmzaRlJRkm4l19uzZBAUFkZ+fT+/evbnuuusoKirimWeeYePGjQQGBnLppZfSs2fPSq81ZcoU3nvvPTp16sT69eu57777WLlyJWCsH7FmzRp2797N6NGjGTduHK+88gozZsw471QcNSXF4YwdX7Jj+QM81q4dKZzmvu73cVfsXXi4yUckRF1Zs2YNEyZMwN3dnebNmzNkyBB+//13AgIC6NOnD61atQKgR48eHDx4sMbFYe3atXz11VcA3HTTTTzyyCOA0cX1xBNPkJiYiJubG0eOHKmye2vFihWsWLHC1kjn5OSQnJxcrjhU9OKLL/LFF1+QlpbG0aNHAejTp4+tMADMnDmTL7/8EoDDhw+TnJzM8ePHSUhIICQkBIDx48dXGo/Jycnh119/5frrr7dtOzP/FBh7JG5ubkRGRl5wd935OF3Lp5TyA94FioBVWmu7zzJl3b2cj1c8yJvhzWnWKIgPBr9KfJjM3SJc37m+4dtLVFSUraujNmoyXXXZA0EKCgrO+5rz5s0jPT2djRs34unpSbt27ap8ntaa6dOnc/fdd1f7WpGRkWzduhWr1YqbmxtPPvkkTz75ZLmuLj8/P9v1VatW8cMPP7B27Vp8fX1JSEioUWYAq9VKkyZNbHs2FZX9rOw1P55DxhyUUrOVUmlKqaQK20cqpfYopfYppaaVbh4LLNRa3wWMtne2zF1L+Muqh5gRFMjgVoNYNGaxFAYhLsLQoUMpLCws18e+bds2Vq9ezaBBg/jss8+wWCykp6eTmJhInz59avzazZs3Z9euXVitVts3cjCm/160aBEACxYssG3PysoiNDQUT09PfvrpJw4dOgRQadruESNGMHv2bNtqa0eOHCEtLa3c7+7YsSPx8fE89dRTWCwWwChQ1TXOWVlZNG3aFF9fX3bv3s26desA6Nu3Lz///DOZmZkUFxfzxRdfVHpuQEAAERERtvu01mzduvWcn01tpiKvCUcNSM8BRpbdoJRyB94BrgAigQlKqUigFXC49GEWe4Zat2U249ZO5zcfb57s+RBvDntHTmoT4iIppfjyyy/54Ycf6NChA1FRUUyfPp2wsDCuvfZaYmNj6d69O0OHDuW1114jLCysRq8J8MorrzBq1CgGDBhQbnW4N998k3/+85/Exsayb98+22pqEydOZMOGDcTExDB37ly6du0KQLNmzWzTdj/66KNcfvnl3HTTTfTv35+YmBjGjRtXZUP73//+l8zMTFuhGD58OK+99lqVmUeOHElJSQndunWzLVcKEB4ezrPPPkv//v0ZOHBgtTPozps3jw8++IDu3bsTFRVlG0yvTmxsLO7u7nTv3r1OBqQdNmW3Uqod8D+tdXTp7f7As1rrEaW3p5c+NAU4qbX+n1Jqgdb6xmpebwowBaBNmza9znwjqI1FKx7m45QfeW3ELDq37Ffr5wvhjOrblN0xMTEsWbKkXF9+RXl5eTRq1AilFAsWLGD+/PnnbUwbGleasrslZ/cQwCgKfYGZwNtKqauApdU9WWs9C5gFxnoOFxJg7PB/MqowG28f2VsQwhkNHz6cmJiYcxYGgI0bNzJ16lS01jRp0sS2CJC4cE43IK21zgVuq8ljL3ZWVqWUFAYhnNj3339fo8cNGjTovH3yonbMPAnuCNC6zO1WpdtqzCGzsgrhgurDCo+i7lzI34OZxeF3oJNSKkIp5QXcCCypzQsopa5WSs3KysqyS0AhXJGPjw+ZmZlSIARgFIbMzEx8fGo3w4NDBqSVUvOBBCAYSAWe0Vp/oJS6EngTcAdma61fvJDXd+ga0kI4ueLiYlJSUmp8TL2o/3x8fGjVqhWenp7ltp9rQNphRyvZQ5kxh7uSk5PNjiOEEC7lXMXBpSfekzEHIYSwD5cuDkIIIezDpYuDDEgLIYR9uPSYwxlKqXTgEBAI1LZSBAMZZZ5b9jVqcp1abqvN/bXJe65857vt6KwVX08+W3M+2+p+r3y2Deezbau1DqnyVbXW9eYCzLqA52wo+9yyr1GT67XdZq+858p3vtuOziqfrXN8ttX9XvlsG95nW9XFpbuVqlDtdBu1eO7SKrad63ptt9Xm/po893y5a3K7Jlnks636dV35s63u98pnW/65DeGzraRedCtdDKXUBl3NoVzOyJXyulJWcK28rpQVXCuvK2UF++Wtb3sOF2LW+R/iVFwprytlBdfK60pZwbXyulJWsFPeBr/nIIQQojLZcxBCCFGJFAchhBCVSHEQQghRiRSHc1BKtVdKfaCUWmh2lppQSl2jlPqPUuozpdTlZuc5F6VUN6XUe0qphUqpe83Ocz5KKT+l1Aal1Cizs5yPUipBKbW69PNNMDvPuSil3JRSLyql3lJKTTI7z/kopQaVfq7/VUr9anaec1FKtVFKfaWUmq2Umlbb59fb4lD6gaQppZIqbB+plNqjlNp3vg9Ma31Aa32HfZPactVF3q+01ncB9wDjnTzrLq31PcANwEBnzlrqceBz+6Qsl6su8mogB/DBWH7XmbOOwVjoq9ieWUtz1cXf7erSv9v/AR85c1YgBliotb4d6FnrELU9G89VLsBgIA5IKrPNHdgPtAe8gK1AZOmH+L8Kl9Ayz1voYnn/AcQ5e1ZgNPANcJMzZwWGYyxGNRkY5ex/B4Bb6fOaA/OcPOs04G5H/D+r4/9jnwP+zpwVaAb8BKwEbqt1Bnv+Y5h9AdpV+HD7A9+VuT0dmF6D17F7caiLvIACXgUuc/asFV5rmTNnBV7EWJRqBfD1mcbXWfOWeZyXAxrci/1sbwZuKL3+mT2z1tVnC7QB/uPsWYFHgMGl12v9d+BBw9ISOFzmdgrQt7oHK6WaYTQMPZVS07XWL9s5X0W1ygvcD1wGBCqlOmqt37NnuApq+9kmAGMBb2C5PYNVoVZZtdZPAiilJgMZWmurXdNVVtvPdiwwAmgCvG3XZJXV9m92MfCWUmoQkGjPYNWobV6AO4AP7ZaoerXN+i3wrFLqJuBgbX9ZQysOtaK1zsTov3cJWuuZwEyzc9SE1noVsMrkGLWitZ5jdoaa0Fovxmh0nZ7WOg+jsXUZWutnzM5QE1rrJGDchT6/3g5IV+MI0LrM7Val25yVK+WVrPbjSnldKSu4Vl6HZm1oxeF3oJNSKkIp5YUxyLjE5Ezn4kp5Jav9uFJeV8oKrpXXsVntPahi1gWYDxzj7CFyd5RuvxLYizHq/6TZOV0xr2SVvK6W1dXyOkNWmXhPCCFEJQ2tW0kIIUQNSHEQQghRiRQHIYQQlUhxEEIIUYkUByGEEJVIcRBCCFGJFAchhBCVSHEQQghRiRQHIeykdGGWLaWX9Uop+f8mXIacIS2EnSilkjHm0z9mdhYhaku+yQhhP8uBbUqpN80OIkRtyXoOQtiBUmoAxsp84VrrErPzCFFbsucghH1cD+zVWpcoQ4DZgYSoDRlzEMIOlFJ9gA8ADeQD92mtN5qbSoiak+IghBCiEulWEkIIUYkUByGEEJVIcRBCCFGJFAchhBCVSHEQQghRiRQHIYQQlUhxEEIIUYkUByGEEJX8PwHirA85xCEuAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "A = np.array([[ .2, .1, 1.0, 1.0, 0.0],\n", " [ .1, 4.0, - 1.0, 1.0, - 1.0],\n", @@ -895,8 +908,52 @@ " [ 1.0, 1.0, .0, 8.0, 4.0],\n", " [ .0, -1.0, - 2.0, 4.0, 700.0]])\n", "b = np.array( [ 1 , 2 , 3 , 4 , 5 ] )\n", - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "x_exact = linalg.solve(A, b)\n", + "\n", + "eps_list = np.logspace(-8, -1, 8)\n", + "\n", + "fig, ax = plt.subplots()\n", + "\n", + "for alg, alg_name in [(GS, \"Gauss-Seidel\"), (SD, \"Steepest Descent\"), (CG, \"Conjugate Gradient\")]:\n", + " #diff_list = []\n", + " k_list = []\n", + " for eps in eps_list:\n", + " x, k = alg(A, b, eps)\n", + " #diff_list.append(diff(x_exact, x))\n", + " k_list.append(k)\n", + " #ax.plot(eps_list, diff_list, label=alg_name)\n", + " ax.plot(eps_list, k_list, label=alg_name)\n", + " \n", + "ax.set_xscale(\"log\")\n", + "ax.invert_xaxis()\n", + "ax.set_yscale(\"log\")\n", + "ax.set_xlabel(\"$\\epsilon$\")\n", + "#ax.set_ylabel(\"$||\\\\vec{x}^* - \\\\vec{\\\\tilde{x}}||_\\infty$\")\n", + "ax.set_ylabel(\"$k$\")\n", + "ax.legend()\n", + "\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The condition number for A is K(A) = 12269.877667702964 >> 1,\n", + "so A is ill-conditioned, so the Conjugate Gradient method is highly susceptible to rounding errors.\n", + "This explains the great difference in order of required iterations k as observed above.\n" + ] + } + ], + "source": [ + "print(\"The condition number for A is K(A) =\", linalg.cond(A), \">> 1,\")\n", + "print(\"so A is ill-conditioned, so the Conjugate Gradient method is highly susceptible to rounding errors.\")\n", + "print(\"This explains the great difference in order of required iterations k as observed above.\")" ] }, { @@ -935,7 +992,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": { "deletable": false, "nbgrader": { @@ -950,10 +1007,59 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The number of iterations is brought down by a lot due to the conditioning.\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEMCAYAAAAvaXplAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABDWElEQVR4nO3dd3xUVfr48c/JpEx6IQlJCL0mIQkldIGIougiKLIi2HBVVl1s+7Miu7q7Yv+u3VVcUVkFe8EVV9eCERdQUIRAgACChJ6EJKQnM+f3x50M6SSQzJ1JnvfrNa/M3HvnzpMhnOeecs9RWmuEEEKI2rzMDkAIIYT7keQghBCiAUkOQgghGpDkIIQQogFJDkIIIRqQ5CCEEKIBSQ5CCCEakOQghBCiAW+zA6hPKZUO/A3YAryptV51svdERkbqXr16tWtcQgjR0WzYsCFXax3V2D6XJAel1BJgKnBEaz241vYpwFOABfin1vphQAPFgBXIacn5e/Xqxfr169s8biGE6MiUUnub2ueqZqVXgSm1NyilLMBzwHlAIjBbKZUIfKu1Pg+4C/iLi+ITQghRi0uSg9Y6A8ivt3kksFNrvVtrXQm8CUzXWtsd+48Bfq6ITwghRF1m9jl0A/bVep0DjFJKzQDOBcKAZ5t6s1JqHjAPoEePHu0XpRBCdEJu1yGttX4feL8Fxy1WSh0ELvD19R1ef39VVRU5OTmUl5e3R5jCA1mtVuLj4/Hx8TE7FCHcnpnJYT/QvdbreMe2FtNafwx8nJaWdl39fTk5OQQHB9OrVy+UUqcXqfB4Wmvy8vLIycmhd+/eZocjhNsz8z6HH4D+SqneSilf4FJgRVudvLy8nC5dukhiEAAopejSpYvUJIVoIZckB6XUcmANMFAplaOUukZrXQ3MBz4DsoC3tdZbWnneC5RSiwsLC5vaf5qRi45E/h5Eh1BZCsf2Qs562LYSCvad/D2nwCXNSlrr2U1sXwmsPI3zNtms5C4OHz7Mbbfdxtq1awkPD8fX15c777yTiy66yGUxrF27lltuuYWKigoqKiqYNWsW999/f5PHr1+/nqVLl/L000832FdzT0lkZGST72/JMUIIB62hvBBKjhqP4iMNn9feVllc9/3TnoFhV7Z5WG7XId0aSqkLgAv69etndiiN0lpz4YUXctVVV7Fs2TIA9u7dy4oVbdZ61iJXXXUVb7/9NqmpqdhsNrZv397s8WlpaaSlpbkoOiE6ILsNSvOaLuzrP7dVNnISBQFdIDAKgqKg2/ATzwOjTzyP6NMuv4JHJwd3rzl89dVX+Pr6cv311zu39ezZk5tuuok9e/ZwxRVXUFJSAsCzzz7L2LFjWbVqFY8//jj//ve/AZg/fz5paWnMnTuXu+++mxUrVuDt7c0555zD448/zjvvvMNf/vIXLBYLoaGhZGRkNIjjyJEjxMbGAmCxWEhMTASgpKSEm266iczMTKqqqrj//vuZPn16nRjy8vKYPXs2+/fvZ8yYMdRec/z111/n6aefprKyklGjRvH8889jsVja7fsUwlTVFS0v7EvzwHnLVi1ePnUL+OjEhoV9zfOALmAxr4j26OTQUn/5eAtbDxS16TkT40K474KkZo/ZsmULw4YNa3RfdHQ0//3vf7FarWRnZzN79uxmpwDJy8vjgw8+YNu2bSilKCgoAOCvf/0rn332Gd26dXNuq++2225j4MCBpKenM2XKFK666iqsViuLFi1i0qRJLFmyhIKCAkaOHMnZZ59d571/+ctfOOOMM/jzn//MJ598wssvvwxAVlYWb731Ft999x0+Pj7ceOONvPHGG1x5ZdtXb4VoVxXHoXA/FOUYP48fbKTgPwoVjfdt4hN4olAP7w3dRzZe2AdFgTUMPKTvy6OTg7s3K9X3hz/8gdWrV+Pr68sXX3zB/Pnz2bhxIxaLhR07djT73tDQUKxWK9dccw1Tp05l6tSpAIwbN465c+dyySWXMGPGjEbf++c//5nLLruMzz//nGXLlrF8+XJWrVrF559/zooVK3j88ccBY4TXr7/+Wue9GRkZvP++cdvJb37zG8LDwwH48ssv2bBhAyNGjACgrKyM6OjoU/9yhGgPVWVQdAAKc6Bof90kUPO6sULfP9wo1IOiISa56cI+MAp8A13/e7mARyeHljYrnewKv70kJSXx3nvvOV8/99xz5ObmkpaWxhNPPEHXrl35+eefsdvtWK1WALy9vbHbT1RHa4Zeent78/333/Pll1/y7rvv8uyzz/LVV1/xwgsvsG7dOj755BOGDx/Ohg0buP322/npp5+Ii4tj5Uqjv79v377ccMMNXHfddURFRZGXl4fWmvfee4+BAwfWifvw4cMn/d201lx11VU89NBDp/09CXFKbFVGwd9ooe9IBqV5Dd8XEAmh3Yyr/F5nQEg3CI13/OwGQTHg7ev638fNeHRycHeTJk1iwYIF/OMf/+CGG24AoLS0FIDCwkLi4+Px8vLitddew2azAUafxNatW6moqKCsrIwvv/ySM844g+LiYkpLSzn//PMZN24cffoYnVC7du1i1KhRjBo1ik8//ZR9+/bxyiuv1Injk08+4fzzz0cpRXZ2NhaLhbCwMM4991yeeeYZnnnmGZRS/PTTTwwdOrTOeydMmMCyZctYuHAhn376KceOHQPgrLPOYvr06dx2221ER0eTn5/P8ePH6dmzZ7t+p6KTsNug+HDThX7hfmM/uu77/EKNgj60m9GBG9oNQhyvQxwPH6spv5Kn8ejk4O7NSkopPvzwQ2677TYeffRRoqKiCAwM5JFHHmHYsGFcfPHFLF26lClTphAYaFRNu3fvziWXXMLgwYPp3bu3s7A+fvw406dPp7y8HK01f//73wG44447yM7ORmvNWWedRWpqaoM4/vWvf3HbbbcREBCAt7c3b7zxBhaLhT/96U/ceuutpKSkYLfb6d27t7MjvMZ9993H7NmzSUpKYuzYsc55rBITE3nggQc455xzsNvt+Pj48Nxzz0lyECentXFF31hTT8224wfBXl33fT4BJ67u+519osCvnQD8gs35nTogVXv0iadKS0vT9Ttzs7KySEhIMCki4a7k78JFKksgNxtydxg/C/edKPiLDkB1vTvVLb4QElf3Kr/+Vb9/uMd05noKpdQGrXWj49Y9uuYghDCR1saIntwd9R6OZFBDeUFwrFHAx6bCwPMhtHutJBBv9AN4yarF7kSSgxCiebZqKNgLR7efKPxzHc/La4308QmEyP7QYwxEXgVRAyBygHGTlrcszeJpPDo5uHufgxAepaIY8rKNwr92IsjfVfcO3qCuRqE/eKbxsyYJBMfJ1X8H4tHJwd3vkBbC7dRpCtpeKxFkG53CNZQFInobhf6AcyByoPE8sj/4h5kWvnAdj04OQogm2Krh2J5G+gOaaArqNc74GTnASAQRvaUpqJOT5CCEJ6tpCjpaLwHk7QJ71YnjgmKMwn/wTIgaeCIRhHSTEUCiUZIc2tGiRYtYtmwZFosFLy8vXnzxRUaNGsWTTz7JvHnzCAgIMC22goICli1bxo033tjofovFQnJyMlVVVXh7e3PllVdy22234WVim/KDDz7IggULTPt801WVQ873sPsb2L/BSAJFtRZPdDYFDYQBUxz9AQOhSz9pChKt5tH3OdTqkL4uOzu7zj6zx7OvWbOGP/7xj6xatQo/Pz9yc3OprKwkLi7OLdY72LNnD1OnTiUzM7PR/UFBQRQXG/PGHzlyhDlz5jBu3Dj+8pe/uDLMJmM6VWb/XbSK3QYHNxrJ4Jdv4Ne1xv0BygIxgyFqkKMZqPaoIJn2QbRcc/c5oLX2+Mfw4cN1fVu3bm2wzZXee+89PXXq1Abbn3rqKe3j46MHDx6s09PTtdZaf/bZZ3r06NF66NCheubMmfr48eNaa63Xr1+vJ0yYoIcNG6bPOeccfeDAAa211hMnTtQ333yzTk1N1UlJSXrdunVaa62Li4v11VdfrUeMGKGHDBmiP/zwQ6211pmZmXrEiBE6NTVVJycn6x07duhZs2Zpq9WqU1NT9e23394gzsDAwDqvd+3apSMiIrTdbtfV1dX69ttv12lpaTo5OVm/8MILWmutDxw4oMePH++MKyMjQ2ut9aeffqqHDh2qU1JS9KRJk5qN9ZVXXtEXXXSRPvfcc3W/fv30HXfcobXW+q677tJeXl46NTVVz5kz55T/Xcz+u2iW3a714Syt176g9bLZWj/YXev7QozHc2O0XnmX1ttWal1WaHakooMA1usmylWPrjnUOOkd0p/eDYc2t+2HxiTDeQ83ubu4uJgzzjiD0tJSzj77bGbNmsXEiROBuiul5ebmMmPGDD799FPn1BoVFRXcc889TJw4kY8++oioqCjeeustPvvsM5YsWUJ6ejr9+/fnpZdeIiMjgxtvvJHMzEwWLFhAYmIil19+uXMK7p9++om7776b0aNHc9lll1FZWYnNZuPw4cMtrjnUCAsLY/v27Xz00UccOXKEhQsXUlFRwbhx43jnnXd4//33KS8v595778Vms1FaWkp5eTnDhg0jIyOD3r17k5+fT0RERJOxvvPOO/z1r3/lp59+ws/Pj4EDB7J69Wq6d+/eMWsOBfuMWsHub+CXDCg+ZGwP6wl9JkLvidB7gjE7qBBtTO6QNkFQUBAbNmzg22+/5euvv2bWrFk8/PDDzJ07t85xa9euZevWrYwbNw6AyspKxowZw/bt28nMzGTy5MkA2Gw254I9ALNnGyuvTpgwgaKiIgoKCpqcgnvMmDEsWrSInJwcZsyYQf/+/U/rd/v888/ZtGkT7777LmBMIpidnc2IESP43e9+R1VVFRdeeCFDhgxh1apVTJgwgd69ewMQERHhPEdT04WfddZZhIaGAsYcTnv37qV79+6nFbPbKMk1kkBNQjj2i7E9MMpIAr0nGkkhvJepYQrROZJDM1f47clisZCenk56ejrJycm89tprDZKD1prJkyezfPnyOts3b95MUlISa9asafTcqt4IE6VUk1NwJyQkMGrUKOfsrC+++KJzVteW2r17NxaLhejoaLTWPPPMM5x77rkNjsvIyOCTTz5h7ty5/PGPf3Su/1BfU7GuW7cOP78TQygtFgvV1dX13+45Koph7/9OJIPDjhqsb7AxXfSo3xsJITpBRg0JtyK3M7aT7du3U7uTfOPGjc4ZS4ODgzl+/DgAo0eP5rvvvmPnzp2AsXTnjh07GDhwIEePHnUmh6qqKrZs2eI831tvvQXA6tWrCQ0NJTQ01DkFd01T4U8//QQYBXufPn24+eabmT59Ops2baoTw8kcPXqU66+/nvnz56OU4txzz+Uf//gHVVXGUMkdO3ZQUlLC3r176dq1K9dddx3XXnstP/74I6NHjyYjI4NffjGukPPz8wGajLU5Pj4+zs90W9UVsGc1fLUIXj4HHukJy34L379kjBiatBCu+QLu2gNz3oTRN0DXREkMwu10jpqDCYqLi7npppsoKCjA29ubfv36sXjxYgDmzZvHlClTiIuL4+uvv+bVV19l9uzZVFRUAPDAAw8wYMAA3n33XW6++WYKCwuprq7m1ltvJSnJWLjIarUydOhQqqqqWLJkCUCTU3C//fbb/Otf/8LHx4eYmBgWLFhAREQE48aNY/DgwZx33nk89thjdeIvKytjyJAhzqGsV1xxBX/84x8BuPbaa9mzZw/Dhg1Da01UVBQffvghq1at4rHHHsPHx4egoCCWLl1KVFQUixcvZsaMGdjtdufyqC2ZLry+efPmkZKSwrBhw3jjjTfa9N/rlNltcGjTiRFFe9dAdZkx2VzcUBh7s9FM1H0U+PibHa0QLebRHdLuPJS1PaWnp/P444+Tltb4CDTRtNP+u9DamGril29g9yqjllBeYOyLGnSiz6DnOLm3QLi9DtshrWVuJeEKhfvrjig6fsDYHtoDEqZC73SjMzm4q5lRCtGmPDo5dFarVq0yO4SOrTQf9nx7oqkoz+gPIqBLvRFFvaWvQHRYkhyEAGNE0faVRkI4tBnQ4BtkNA+l/c4xoihRpqQWnYYkByG+ewr++2djqcr4kXDmAiMZdBsGFh+zoxPCFJIcROelNXz1AHz7OCTNgOnPgm+g2VEJ4RYkOYjOyW6H/9wN378Iw66EqU+Cl8XsqIRwG9KA2o4OHTrEpZdeSt++fRk+fDjnn38+O3bsOKVzXXvttWzdurVN4ysoKOD5559v9fuKi4u54YYb6Nu3L8OGDWP48OG89NJLpxXLq6++yvz58wF44YUXWLp06SmdZ8+ePSxbtqz5g2zV8NEfjMQwZj5c8LQkBiHqkeTQTrTWXHTRRaSnp7Nr1y42bNjAQw89xOHDh0/pfP/85z9JTExs0xhPNTlce+21hIeHk52dzY8//sh//vMf553PtZ3qtBfXX389V1555Sm996TJQWt4dy78vAzOvBfOeUBGHAnRCLdMDkqpQKXUeqXUVLNjOVVff/01Pj4+XH/99c5tqampjB8/Hq01d9xxB4MHDyY5Odk5FcaqVatIT09n5syZDBo0iMsuu8w5vUR6ejo1M88GBQU5z/nuu+8652vatWsXo0ePJjk5mYULFzqPKy4u5qyzzmLYsGEkJyfz0UcfAXD33Xeza9cuhgwZwh133AHAY489xogRI0hJSeG+++5r8Hvt2rWL77//ngceeMC58E9UVBR33XWX83cYP34806ZNcyazCy+8kOHDh5OUlOS8SxzglVdeYcCAAYwcOZLvvvvOuf3+++93Tsi3a9cupkyZwvDhwxk/fjzbtm0DYO7cudx8882MHTuWPn36OCcBvPvuu/n2228ZMmQITzzxRN3g7TYoOQpZH8OUh2HinZIYhGiCS/oclFJLgKnAEa314FrbpwBPARbgn1rrmhny7gLebqvPf+T7R9iWv62tTgfAoIhB3DXyrib3Z2ZmMnz48Eb3vf/++2zcuJGff/6Z3NxcRowYwYQJEwBjjqEtW7YQFxfHuHHj+O677zjjjDNaFNMtt9zCLbfcwuzZs3nhhRec261WKx988AEhISHk5uYyevRopk2bxsMPP0xmZiYbN24EjJlSs7Oz+f7779FaM23aNDIyMpyxAWzZsoXU1NRmV4T78ccfyczMdM7EumTJEiIiIigrK2PEiBFcfPHFVFZWct9997FhwwZCQ0M588wzGTp0aINzzZs3jxdeeIH+/fuzbt06brzxRr766isADh48yOrVq9m2bRvTpk1j5syZPPzwwzz++OMNp+KwV0PebmPuo+nPwdDLW/SdCtFZuapD+lXgWcDZkKyUsgDPAZOBHOAHpdQKoBuwFbC6KDaXW716NbNnz8ZisdC1a1cmTpzIDz/8QEhICCNHjiQ+Ph6AIUOGsGfPnhYnhzVr1vDhhx8CMGfOHG6//XbAaOJasGABGRkZeHl5sX///kabtz7//HM+//xzZyFdXFxMdnZ2neRQ36JFi3jnnXc4cuQIBw4Ydw6PHDnSmRgAnn76aT744AMA9u3bR3Z2NocOHSI9PZ2oqCgAZs2a1aA/pri4mP/973/89re/dW6rmX8KjBqJl5cXiYmJzTfX2aqMNZWryyGwC6Sc0/SxQgjARclBa52hlOpVb/NIYKfWejeAUupNYDoQBAQCiUCZUmql1tpe/5xKqXnAPIAePXo0+/nNXeG3l6SkJGdTR2u0ZLrq2tN1l5eXn/Scb7zxBkePHmXDhg34+PjQq1evRt+nteaee+7h97//fZPnSkxM5Oeff8Zut+Pl5cW9997LvffeW6epKzDwxHDQVatW8cUXX7BmzRoCAgJIT09vUcwAdrudsLAwZ82mvtrfVZNzhFVXGnc426qMZTQL9zd+nBCiDjP7HLoB+2q9zgG6aa3v1VrfCiwDXmosMQBorRdrrdO01mk1V5/uZNKkSVRUVNRpY9+0aRPffvst48eP56233sJms3H06FEyMjIYOXJki8/dtWtXsrKysNvtzityMKb/fu+99wB48803ndsLCwuJjo7Gx8eHr7/+mr179wI0mLb73HPPZcmSJc7V1vbv38+RI0fqfHa/fv1IS0tj4cKF2Gw2wEhQTRXOhYWFhIeHExAQwLZt21i7di0Ao0aN4ptvviEvL4+qqireeeedBu8NCQmhd+/ezn1aa37++edmv5s6v1N1OeRlG01KXfqCNaTZ9wohTnDLDmkArfWrWutm53BWSl2glFpcWFjoqrBaTCnFBx98wBdffEHfvn1JSkrinnvuISYmhosuuoiUlBRSU1OZNGkSjz76KDExMS06J8DDDz/M1KlTGTt2bJ3V4Z588kn+/ve/k5KSws6dO52rqV122WWsX7+e5ORkli5dyqBBgwDo0qWLc9ruO+64g3POOYc5c+YwZswYkpOTmTlzZqNrPvzzn/8kLy/PmSgmT57Mo48+2mjMU6ZMobq6moSEBOdypQCxsbHcf//9jBkzhnHjxjU5U+obb7zByy+/TGpqKklJSc7O9KakpKRgsVhITUnhiUV/Bm2HLv3AL6jZ9wkh6nLZlN2OZqV/13RIK6XGAPdrrc91vL4HQGv9UGvPfdI1pDuA5ORkVqxYUactv77S0lL8/f1RSvHmm2+yfPnykxamHVJlidHHoLyMxOBzovuqo/1dCHE63HXK7h+A/kqp3sB+4FJgTmtOUGs9h3YIz31MnjyZ5OTkZhMDwIYNG5g/fz5aa8LCwpyLAHUqFcchfzd4eRuJwdvv5O8RQjTgqqGsy4F0IFIplQPcp7V+WSk1H/gMYyjrEq31lmZO00BnWc/hv//9b4uOGz9+/Enb5Du08kLI/8VICF36GhPpCSFOiatGK81uYvtKYKUrYhAdXGk+FPxqLMUZ0RcsMm2YEKfDbTukW8KdO6SFC5XkQsFeY0bVLv0kMQjRBjw6OWitP9Zaz6sZlSM6oeLDULgP/EKMGoNMoCdEm/Do5CA1h05Mayg6YDysYRDRW1ZpE6INefT/JnevOciU3S3Xqim7tYai/UatIaALhPcyhq3Swim7hRAn5dHJwZ3JlN3tNGW31kbHc8lRCIyG0O51ZlaV5CBE25Dk0E5kyu52mrJ73RdQls/cOx7m5oUPM3bcuJZP2S2EaDGPHtbR0pvgDj34IBVZbTtlt1/CIGIWLGhyv0zZ3cZTdvftw7ovPuDG/7eArz77BHz8OXjoUMun7BZCtIpH1xzcvc+hKU1N2Q04p+z28vJyTtndUmvWrHFObz1nzombzWum7E5JSeHss89u0ZTdw4YNY9u2bWRnZzf7mYsWLWLIkCHExcU5tzU2ZXdqaiqjR492Ttm9bt0655Tdvr6+zJo1q8G5T0zZPZMhKYP5/W33cjC3EIKigVZM2S2EaDWPrjm0VHNX+O1Fpuw2tMmU3f99C6rKIbwn+Ic797doym4hxCnx6JqDOw9llSm7T3z2KU/ZHWild3wM77z/MUT0QVvDWjdltxDilHl0cnDnZiWZsttwylN226shN5s3nl3Ey+99Tuqo8a2bsjs1VTqkhTgNLpuyuz3JlN2GDjNld1WZMeW2thsT6PkGnvw9LdTR/i6EOB3uOmW3aKFONWV37bUYIvsbE+kJIVxOkoMH6DRTdstaDEK4DUkOwj3IWgxCuBWP7pA+2WiljtCf0imUHTMSg48VuvRvt8Qgfw9CtJxHJ4fmRitZrVby8vKkQHB3JblwbI9jLYb+7bYWg9aavLw8rFbryQ8WQnTcZqX4+HhycnI4evSo2aGIplQUQVkBePtDoC8cPbUZa1vKarUSHx/frp8hREfRYZODj4/PSUf3CJNoDV8/CBmPQuKFMOMl8JY+BiHcSYdNDsJN2e3w2T2w7gUYegVc8JSs3iaEG5LkIFzHVg0f3wwb34DRf4BzF9VZi0EI4T4kOQjXqK6E96+FrR9B+j0w8S5JDEK4MY9ODi1dz0GYrLIU3r4Cdn4B5z4IY/5gdkRCiJPosENZhZsoL4TXL4adX8K0ZyQxCOEhPLrmINxcSR68fhEc3gIzl8DgGWZHJIRoIUkOon0UHYSl06FgL1y6HAacY3ZEQohWkOQg2sfK26EwBy5/D3q1bA1sIYT78Og+B+Gm9v8I2/4NZ9wqiUEIDyXJQbS9rxeBfwSMut7sSIQQp0iSg2hbe9cYQ1bPuBWsIWZHI4Q4RW6XHJRSCUqpF5RS7yqlbjA7HtEKWsNXD0BgNIy4zuxohBCnwSXJQSm1RCl1RCmVWW/7FKXUdqXUTqXU3QBa6yyt9fXAJcA4V8Qn2sgv38De1TDhdvANMDsaIcRpcFXN4VVgSu0NSikL8BxwHpAIzFZKJTr2TQM+AVa6KD5xumpqDSHxMHyu2dEIIU6TS5KD1joDyK+3eSSwU2u9W2tdCbwJTHccv0JrfR5wmSviE20g+3PI+cGoNcjaz0J4PDPvc+gG7Kv1OgcYpZRKB2YAfjRTc1BKzQPmAfTo0aPdghQtUFNrCO8FQy83OxohRBtwu5vgtNargFUtOG4xsBggLS1N1gI1U9bHcGgTXPgCWHzMjkYI0QbMHK20H+he63W8Y1uLKaUuUEotLiwsbNPARCvYbcaqbpEDIOUSs6MRQrQRM5PDD0B/pVRvpZQvcCmwojUnkFlZ3UDm+3A0y1ijQVZ0E6LDcNVQ1uXAGmCgUipHKXWN1roamA98BmQBb2utt7TyvFJzMJOtGlY9CF0HG2tBCyE6DJf0OWitZzexfSWnMVxVa/0x8HFaWprccWWGn5dD/m5j1lUvt7ufUghxGjz6f7TUHExUXQHfPAJxw2DgeWZHI4RoYx6dHKTPwUQ/LoXCfTBpoawFLUQH5NHJQZikqgwyHoceY6HvJLOjEUK0A49ODtKsZJIfXobiQzDpXqk1CNFBeXRykGYlE1QUw+onoE+6LOQjRAfm0clBmOD7F6E0F85caHYkQoh25NHJQZqVXKysAL57CgZMge4jzI5GCNGOPDo5SLOSi619HsoL4cwFZkcihGhnHp0chAuV5MGa5yFxOsSmmh2NEKKdSXIQLfO/p6CyGNKl1iBEZ+DRyUH6HFzk+GFYt9iYdTV6kNnRCCFcwKOTg/Q5uMjqv4OtEibeZXYkQggX8ejkIFygMAfWL4Ehc6BLX7OjEUK4iCQH0byMx41lQCfeaXYkQggXkuQgmpb/C/z0Lxg+F8JknW4hOhOPTg7SId3OvnkUvLxh/P8zOxIhhIu1ODkopR5oZJup60JKh3Q7OroDNr0JI66FkFizoxFCuFhrag7dlFLOFd2UUtHAF20fknALqx4Cb3844zazIxFCmKA1y4T+HvhMKbUL0MArgIxt7IgOZcKW92H87RAYaXY0QggTnDQ5KKWWAj8CPwF/AJYB1cCFWuud7RueMMXXD4JfKIydb3YkQgiTtKRZ6VVAAVcDrwO9gGPA5Uqpme0WmTDH/g2w/RMYexP4h5sdjRDCJCetOWitvwK+qnmtlPIGEoBUYBTwbrtFJ1zvq0XgHwGjrzc7EiGEiVrT5wCA1roa2Ox4vN7mEbWCUuoC4IJ+/fqZGUbHsXcN7PoSJv8V/ILNjkYIYSKPvs9BhrK2Ia3hqwcgMBpGXGd2NEIIk3l0chBt6JdvYO9qmHA7+AaYHY0QwmSSHMSJWkNIvDFVhhCi05PkICD7c8j5ASbeAd5+ZkcjhHADkhw6O7vdqDWE94Ihl5kdjRDCTbR6tJLoYLZ9DIc2wUUvgsXH7GiEEG5Cag6dmd1m3A0dOQCSf2t2NEIINyI1h84s8z04ug1++yp4mTrBrhDCzbhlclBKXQj8BggBXtZaf25uRB2QrdqYebVrMiRMNzsaIYSbcVmzklJqiVLqiFIqs972KUqp7UqpnUqpuwG01h9qra8DrgdmuSrGTuXn5ZC/G85cAF7SuiiEqMuVpcKrwJTaGxyLBT0HnAckArOVUom1Dlno2C/aUnUFfPMIxA2DgeeZHY0Qwg25LDlorTOA/HqbRwI7tda7tdaVwJvAdGV4BPhUa/2jq2LsNH5cCoX7YNJCUMrsaIQQbsjs9oRuwL5ar3Mc224CzgZmKqUanR5UKTVPKbVeKbX+6NGj7R9pR1FVBhmPQ4+x0HeS2dEIIdyUW3ZIa62fBp4+yTGLgcUAaWlp2hVxdQg/vAzFh2DmEqk1CCGaZHbNYT/QvdbreMe2FlFKXaCUWlxYWNjmgXVIFcWw+u/Q50zoNc7saIQQbszs5PAD0F8p1Vsp5QtcCqxo6Ztlyu5WWvcClOYZfQ1CCNEMVw5lXQ6sAQYqpXKUUtc4Fg6aD3wGZAFva623tOKcUnNoqbIC+N/TMOA8iE8zOxohhJtTWnt+c31aWppev3692WG4t68fNIav/v5biE0xOxohhBtQSm3QWjd6tWh2s9JpkZpDC5XkwZrnIfFCSQxCiBbx6OQgfQ4t9L+noLIY0u8xOxIhhIdwy6Gsog0dPwzrFkPKJRA9yOxohBCNsNs1JZXVlFTYKK6opriimpJaP43nNue22vvnTejD+P5RbR6TRycHpdQFwAX9+vUzOxT3tfrvYKuEiXeZHYkQHUpFtY3i8hMFekll3cL8uGNf/e3G87pJoLTS1qLP9FIQ5OdNkJ83gY5HZbW9XX4/j04OWuuPgY/T0tKuMzsWt1SYA+uXwNDLoEtfs6MRwlTVNrtRKFeeKKidBXtFNaWVda/OSxwFvlHI100CJRXVVNlaNpjHz9uLYKujMPc1CvbIIF96dgkwtvsa+04U+JY6CaD2T6uPF8pFN696dHIQJ5HxmPFzwp3mxiHEKaiy2Z0FdWmlrc7Vd8MrcpuzMC+usFFaU8DXaqpp6RW2l8IoxK0nrs6D/CxEBvkZhba1VqHta6lbgFtrPfc1Cnpvi2d27Xp0cpBmpWbk74afXoe030FY95MfL0Qb0FpTUmmjsKyKwtIqCsoqKSqr20ZeWln3ir2k0lanyaUmEbSqMPeruQK3OAvn+IAAgvzqFt4BvpY6V+IB9Qr3QD8L/j4Wl12dn44jpUfYfHQzSZFJxATGtPn5PTo5SLNSM755FLy8Yfz/MzsS4WG01pRVGQV8QWmVUdA7HkWOn41tL3D8rLY339xSU5jXbjcP9LUQERjg2GZxNrXUXLUH+DZsdqnZ5sqmFrOUV5eTlZ/FpqObjEfuJg6VHAJg4aiFzBrU9sveeHRyEE04ugM2vQVj/gDBbX9FITxDuaOAb6wwLyytrPu6VuFeWFbVbHu6UhBi9SEswIdQf+PRLdyfUH8fwvxPbKt5hPj7ONvcg/y88fPu+IX56dBas7doL5tyNzmTQfaxbKp1NQDdgroxNGooyYnJpESlMCiifUYhSnLoiFY9BD4BMO5WsyMRbaC8ykZ+SWWdxzFH4V5QeqJAr13IF5ZVnbRZJsTqTWitAj4m1Eqov2+Dwr12Egjx9yHYzxsvLync20phRSGbczc7awSbj26mqLIIgECfQAZHDubqwVeTHJlMclQykf6RLonLo5OD9Dk04lAmbHkfxt8Oga75IxItZ7drisqrGhT2eSWVHKt5XVp3X3PDHIP9vAmpVZD3iw6qU5DXL9xrHsFWHyxSwLtclb2KHcd2sOmokQQ25W5ib9FeALyUF33D+jK552RSolJIiUyhd2hvLF4WU2KVuZU6muVzYM9quPVn8A83O5oOr6LaxrGSWoV9aSX5xRWNFvLGFX8Vtiba5P19LEQE+rboER7gS4jV22NHwnQGWmsOlRxyNg9tzt3M1rytVNgqAIj0jyQlMoXkqGRSIlNIikwi0CfQpTE2N7eSR9ccRD37N8D2T+DMhZIYToHWmuMV1eQX1xTyjRfwtR/FFdWNnkspCPP3cRbmvSMDGd4z3PHaj4hAH+NngC8RQb5EBPji72vOFaJoG6VVpWzJ2+LsJ9icu5mjZcYqlX4WPxIiEpg1cBbJUcmkRqYSExjj1n0vkhw6kq8WgX8EjG50ZdVOy2bX5JVUcKSogiPHyzlcZDw/fLycI0XlHDleweGicvJLKpvsiPX19qJLrSv3nl0CCA/wpUugL+GBvnX2RQT6EhbgK802HZhd29ldsNvoK3DUDHYW7MSujX6eniE9GRU7ytk8NCB8AD4WH5Ojbh1JDh3F3v/Bri9h8t/AL9jsaFzCbtfkl1ZyuKjcKOxrFfSHiyo46kgER4srGm3KiQj0JTrYj64hVgZ2DaZLkF+ThX2Ar2eMfRftI68sr06n8ZbcLRRXFQMQ7BtMSmQKZ/U4y+g0jkwmzBpmbsBtwKOTg3RIO2gNXz0AQV1hxLVmR3Pa7HbNsdJKDte+unc8P1xUwZHjFRwpKufo8YpGx9SHB/jQNcRKdIiV/l2D6RpiJIDoYCvRjudRQX74ekt7vWio0lbJtvxtzkSw6egm9hcbqxdblIUB4QP4TZ/fkBKVQnJkMj1DeuKlOt7fkkcnB7kJzmH3Ktj7HZz3GPgGmB1Nk7TWHCutclzZlzsL+cO1rvprmnkaK/TDAnzo6ijg+0VFOgv9riF+RAXX/PTDz1va7s2WW5bL2oNrWXtgLRsOb6DcVm52SC1WWFFIlb0KgJjAGJIjk7l04KWkRKWQ0CUBf29/kyN0DY9ODoITtYaQeBh+lWlhVFTbOFRYzsHCcg4WlnGwsJxDheW12vaN9v7G2vRD/X2cBX2fqC5Gge9o7okO8SM62EpUsB9WHyn03VVpVSnrD683EsLBtWQfywYgxDeEkTEjPaqZJcQ3xDmKKDog2uxwTCPJwdNlfw7718MFT4O3X7t8REW1jcOFFc5Cv3YCOFhYxsGCcvJKKhu8L8Tq7biytzKqdyDRjqt8o4nH0bwjhb5HqrZXsyVvC2sPrGXNwTX8fPRnqu3V+Hr5MrTrUG4Zdgtj4sYwKHyQaeP0xemR5ODJ7Haj1hDeG4bMOaVTVFbbOVxUr8AvqJ0EysktrmjwvhCrN7Gh/sSGWUnuFkZsqNXxMLbFhFgJ9JM/r45Ca82eoj3OpqIfDv3A8arjACREJHBF4hWMjh3NsOhhWL2tJkcr2oL87/Vk2z6GQ5vgohehkWFyVbYTBf+BgrIGzT41BX/9+yCDrd7Ogj4pLsRZ4DsL/1Ap+DuD3LJc1h1c52wqqpnoLS4wjnN6ncPouNGMihlFuFXuqemI5H+4h6qqqoIvHqA6tC9f6HEczNjFgYJyRwIwCv+jjRX8ft7EhFqJDfMnISaE2DArcaH+xIRaiQuzEhPqT5AU/J1SaVUpPx750dlUtOPYDsBogx8VO4rrkq9jTOwY4oPjZVhvJ+DRpUBHH8qqtebI8Qq2Hihi68Eisg4Wse9YGYcKyxhb8iVP+OzglsqbWfnmJgACfS3EhhlX9gNjgp1X+bFh/sSFWokJtRJs9awbcUT7sdltRr+Bo2aw8chGquxV+Hj5MCx6GLcMu4XRsaNJiEiQfoNOSOZWchPVNju7c0vYesBIAlsPFrH1QFGdjt7uEf706hJIXLA3d+++CnwC+Pk3K4gNCyQ2zEqIFPyiGVprfj3+K2sOrGHtwbV8f+h7jlca/QaDIgYxOnY0Y2LHMLTr0E4zXLOzk7mV3Mzx8iq2HTpuJAFHrWDboePOKZZ9LV4MiAnirIRoEmNDSIwLZVBs8InC/8elsHUfXPQm6QNlvQbRtPzyfGe/wZoDazhYchCA2MBYJveczOjY0YyKHUWENcLkSIW7keTQjrTWHCwsr5MEth4sYm9eqfOY8AAfEuNCuGpMTxLjQkiMDaVPVCA+Tc22WV1hrPLWbTgMmOKi30R4irLqMn48/KOzqWhb/jbAmOJhVMworhl8DaPjRtMjuIf0G4hmSXJoI1U2OzuPFDdIBAWlVc5jenUJICkuhN8OjycxLoSE2BBiQqyt+0/641Io3AfTnjam/hSdms1uIys/y9lU9NORn6iyV+Ht5c3Q6KHcNPQmxsSOIbFLovQbiFaR5HAKCsuq2FarX2DrwSKyDxdTaTOahfy8vRgUE8x5g2NIjDWSwKDYkNMfBVRVBhmPQ89x0OfMNvhNRJW9ikPFh8gpznFOmeAJDpUcYs2BNaw7tM7ZbzAwfCBzBs1hdJxxv0GAj/tOpSKap7WG6mp0ZSX2ykp0RYXxqKzEXlGJrjSe64oK/AYMwCc2ts1jkOTQDK01OcfK6nQQbz1YRM6xMucxXQJ9SYwL4eozehn9A7Eh9I4MbJ9FWH54GYoPwcwlUmtohdKqUvYd30fO8Rz2Hd9X53Gw5CA23fRKa+6sa0BXzupxFmNixzAydqTLlo/s6LTWUFV1olB2/KxfKNsrKtCVVca22q8rKpzH2StqFexVtV7XnKOy0vm8/mvszS/zWiN20SLCLp7R5t+DJAeHymo72UeOn2gScowaKio3FnNRCnpHBpLaPYzZI3uQGBdCUmwIUcF+rmm7zdsF3/4f9J0Evca1/+d5EK01+eX5zgK/fhLIK8+rc3yoXyjdg7qTHJnMeb3Po3twd+KD4z1qhE6Ibwjdg7tLv0ETtN2OvaQEW2EhtsJC7EVFjudFjtcnntuKirAVFWIvMJ7bS0pocINQaymF8vND+fqi/Hzx8q157ud87RUUhMX52hfl6+fcr3x98fLzM7bVnKPR1774dO/eNl9aPZ06OXyz4ygfbdzP1gNF7DxS7JwJ1N/HwsCYYKamxjlGC4UwKCaYAF+Tvq6CX2HpdCNDTXnEnBhMVm2v5lDJoSYTQGn1iU5+haJrYFe6B3dnYveJzsK/e3B3ugd3J8Q3xMTfRLSU1hpdVmYU3oVF2AoL6hbyRY6Cv7DIcUzhiUL++PFmr7yVry9eoSFYQkOxhITi0zUGS/8BxragIJSftUEhbBTejbx2FvrGdi9fX/Dx8fjE7XbJQSnVB7gXCNVaz2zPz8o+fJxvs3NJjA3hzEHRzkTQq0ug+6ziVXQQXpsGFUVw1ccQNcDsiNpNWXVZg0K/5vWB4gNU6xNLcvp4+TgL/LSYNGfBHx8cT7egbvhZ2mcSQtF6urLyROHtKNTthXWv3O1FhdgcV+7ObYWF6Kpm+oG8vLCEGAW8V2golrAwfHv0wBIaYrwOCTUK/9AQLCGObY6Hl1XmfzoZl9wEp5RaAkwFjmitB9faPgV4CrAA/9RaP1xr37stTQ6nehOcza7dJwk0piQXXjkfivbDFR9C9xFmR3RatNYUVhQ2aPevSQJHyo7UOT7YN9hZ6Nd/RAdEd8gFVjyBvbSU6vxj2I7lY8vPdzw3Xlfn52PLP2ZsP2b8tBcXN3s+r6AgRwF/4kreEhKCJSwUr5Ba2xz7vUJCjX2BgR5/dW42d7gJ7lXgWWBpraAswHPAZCAH+EEptUJrvdVFMbl3Yig7BksvNJqULn/XoxJDpa2S7fnb2XFsR4MEUDOTZ41o/2jig+MZ221sgwQQ6hdq0m/QeWitsR8/XquQr1/AO54fO+Z8rsubWLjHxwfv8HAsERF4R4TjHx+PJTwcS0Q4lrCwxq/kg4NR3m7XgCFwUXLQWmcopXrV2zwS2Km13g2glHoTmA60KDkopeYB8wB69OjRdsG6g/IieP1iyN0Os5dDrzPMjqhJNruNXYW72JK7hczcTDLzMtlxbAfVdqMJyNvLm25B3YgPjic1KrVO4R8fHC/TO7cxbbNhKyho+or+mGN7TcF/rACqqxs9l/L3dxb2li4R+PXrZzyPCMc7IsIo+MMdzyMi8AoKkiv5DsTMlN0N2FfrdQ4wSinVBVgEDFVK3aO1fqixN2utFwOLwWhWau9gXaayFJbNggMbYda/oN/ZZkfkpLUm53gOmXmZRiLIzSQrP4uyamNob5BPEEmRSVyVeBWDIwczKGIQsYGxcvNVK+iqKuwlJdhLS41H7ec1r0tKsR0vwnasoO7VfX4+tsLCJkfaeAUHGwV7eAQ+8fFYU5LxDo9wXulbIiKwhEfgHR5mFPb+njN6S7Q9t6vPaa3zgOtbcmyHm5W1qhzenAP71sLF/4RBvzE1nKOlR521gS25W8jMy6SwohAAP4sfgyIGMaP/DJK6JDE4cnCHXWi9Kbq6umEBXlKKvbTE8bN+AV9S95jSUnRpKbaSErTj+GY7YGvz8jKaasLD8Q4Pd1zV11zR17q6d1zhe4eFoXx92/cLER2KmclhP1B7gG68Y1uLaa0/Bj5OS0u7ri0DM4WtCt6ZC7u/hunPw+CLXfrxRZVFbMndwpY8o3loc+5mjpQaHcQWZaFfWD/O7nE2SZFJDO4ymH7h/fDx8vxZYHVVFRW7dlGxfTvVx47VLbBrF/iNXMHryoZLozZFWa14BQSceAQGYgkKxiu6q/O1V2Dd/V4BAah6r70CHMf5+6MsUiMT7cfM5PAD0F8p1RsjKVwKtGqtyw5Tc7Db4P3rYMencP7jMPSydv24suoytudvZ3PuZjJzM9mSt4W9RXud+3uG9CStaxqDIwc7m4c86QaxpthLSijfvp3yrCzKs7Ko2JpFRXZ2g6t15efXoCD3CgzEOzqq8YK7fuFdv7D395dOV+FxXDWUdTmQDkQCh4H7tNYvK6XOB57EGMq6RGu96FTO79HrOdjt8NEf4OdlMPlvMO7mNj19lb2KXQW72Jy72dlpvLNgp3PKiOiAaAZ3MZJAUmQSSV2SOsQooer8fMq3ZlGetZWKrCzKt2ZRuXevsz3eEh6ONSEBa2ICfgkJWBMS8I6MNApyH8+vEQnREqYPZdVaz25i+0pgpSticEtaw8rbjcSQvuC0E4Nd29lbtNdZG8jMzWRb/jYqbBWAMeXC4MjBTIif4KwVRAdEt8VvYhqtNVX791O+dSsV27Y5EkIW1YcPO4/xiYvDLzGBkAumYk1IxJqYgHfXrjKyRohmeHRd16OblbSGzxfC+pdh3C0w8c5Wvl1zuPSwc9RQZl4mW3O3Ou8j8Pf2JyEigUsGXsLgLoNJjkz2+LV/dXU1Fbt3O2sC5VlZlG/bhr2oyDjAywu/vn0IGDXSSAIJCVgHDcQSFmZq3EJ4Ilkm1CxfPwjfPAIj58F5j7ZoltWNRzay9uBaZ0KomVDOW3nTP7y/szYwOHIwfUL74O3lubnfXlZGxfbtlNeqDVTs2GHMVonRL+A3cKCRAGqahwYMkGkRhGgF05uV2ovH1hxWP2EkhqGXGxPpnSQxlFeX89gPj/H2jrdRKHqF9mJs3FiSIpNIjkxmYMRAj55LyFZQYNQCamoDWVlU/vKLc+I0r5AQrAkJhM+ejTXRSAa+vXtLJ68Q7UhqDq627kX49E4YPBNmLIaT3CC289hO7si4g50FO5mbNJd5KfMI9g12UbBtS2tN9cGDjgSwzfFzK9UHDjqP8Y6JqVsbGJSAT7c4j24OE8Jdddiaw+myl5RgLys7+YFtZfM78NkC6HcujF8E+ceaPFRrzce7P+bZn54jwCeAF0c+wsjYkVBUQTUVrov5NNgKC2vVBrZSkbUNW0GBsVMpfHv1ImDIUKxz5pwYMRQhC90L4Q46dc3h6PPPk/v0M+0QkahN+fjg178/fo4mIWtCItaBA/AKDDQ7NCE6tQ5bczjdPoeg8RPwDg9v26Aac2gzbHgNIvrAyOvA0vQ4+l+L9vHhzg8orizmzB5nMip2NF4e2qTiFRCA36BB+PXpI/cOCOFhOnXNwSWy/wvLZ0PcULjiffBrvL/AZrfx0uaX+MfP/yAuMI7HJj7G4MjBjR4rhBBtocPWHNzeLxnw1uXQNREue6fJxHCo5BD3fHsP6w+v5/ze5/On0X8iyDfIxcEKIcQJkhzay6/rYNmlEN4bLv8A/MMaPezrX7/mT//7E5W2Sh4Y9wDT+k6TkTlCCNN5dHJw2/scDvwEb8yE4Bi48kMI7NLgkApbBf+3/v9Yvm05CREJPDrhUXqF9nJ5qEII0RiPnnxfa/2x1npeaKgbTRR3eCv86yKwhsFVK4wEUc/ugt3M+WQOy7ct54rEK3j9/NclMQgh3IpH1xzcTu5OWDodvK1w1UcQGl9nt9aaD3Z+wMPfP4zVYuW5s55jQvwEk4IVQoimSXJoK8f2wtJpoO1w5SfGsNVaiiqL+Nuav/GfPf9hVMwoHhz/oMfPiCqE6LgkObSFogPw2gVQWQJz/w1RA+rs3nhkI3d/ezeHSg5xy7BbuDrpallXWQjh1jw6ObhFh3TxEXhtGpTmG01JMcnOXTa7jVe2vMKzPz1LTGAMr533GqlRqebFKoQQLSQd0qejNB+WXghF+437GLoNd+46UnqE3//39zz141NM7jmZdy54RxKDEMJjeHTNwVTlhfD6DMjbCXPegp5jnLsycjJYuHoh5bZy/jr2r1zY70K5d0EI4VEkOZyKyhJ44xJjzqRZb0DfM43Ntkqe2PAEr2e9zoDwATw24TH6hPU5ycmEEML9SHJorapyY66knO9h5hIYOAWAPYV7uDPjTrLys5gzaA5/TPujRy/AI4To3CQ5tEZ1Jbx9pTFn0kUvQNJFaK35aNdHPLjuQfwsfjwz6RnSu6ebHakQQpwWSQ4tZauG966B7M9g6hOQeinFlcX8be3fWPnLSkbEjOChMx6ia2BXsyMVQojT5tHJwWVDWe12+OhGyFoB5z4Eab9j89HN3JlxJwdLDjJ/yHyuTb5W7l0QQnQYMpT15B8Cn9wGm96CSQuxj76eJZlLuPLTK7FpG69MeYXfp/5eEoMQokPx6JpDu9Ma/nMPbHgVxv8/ckdczYL/Xs+ag2uY3HMy9425j1A/N5r0Twgh2ogkh+Z89TdY9w8YdQPfDZzEghUXU1JVwp/H/JmZ/WfKvQtCiA5LkkNTMh6Db/+PqmFX8nRkJK9+eQP9wvrx8jkv0y/czdaPEEKINibJoTFrnoevHuDXwdO40+soW7auYtbAWdyedjtWb6vZ0QkhRLuT5FDf+iXw2T18POAMHqjYgXeVN0+mP8lZPc8yOzIhhHAZSQ61/fwmJZ/8Px7sk8yKql8ZFj2Mh8c/TGxQrNmRCSGES0lyqLHlA7asvJk7e/Uih+PcmHoj16Vch7eXfEVCiM7H7Uo+pVQg8DxQCazSWr/R3p9p37aSf31+C0/GdqWLfwQvT3iEtJi09v5YIYRwWy65CU4ptUQpdUQplVlv+xSl1Hal1E6l1N2OzTOAd7XW1wHT2ju2vKwV/GHVrTweEcqE+PG8N/19SQxCiE7PVXdIvwpMqb1BKWUBngPOAxKB2UqpRCAe2Oc4zNaeQa3duISZa+7he6sf9w69lSfPek5uahNCCFzUrKS1zlBK9aq3eSSwU2u9G0Ap9SYwHcjBSBAbaSZ5KaXmAfMAevTocUpx7T+ymVAsvDh5MQO6jT6lcwghREdkZp9DN07UEMBICqOAp4FnlVK/AT5u6s1a68XAYoC0tDR9KgHMmPx3plYU4WeV2oIQQtTmdh3SWusS4OqWHHu6s7IqpSQxCCFEI8yclXU/0L3W63jHthZzyaysQgjRCZmZHH4A+iuleiulfIFLgRWtOYFS6gKl1OLCwsJ2CVAIITorVw1lXQ6sAQYqpXKUUtdorauB+cBnQBbwttZ6S2vOKzUHIYRoH64arTS7ie0rgZWnel6XrQQnhBCdjKwEJ4QQogGPTg5CCCHah0cnB+mQFkKI9qG0PqX7x9yKUuoosBcIBVqbKSKB3FrvrX2Oljynldtas7818TYX38leuzrW+ueT79ac77apz5XvtvN8tz211lGNnlVr3WEewOJTeM/62u+tfY6WPG/ttvaKt7n4Tvba1bHKd+se321Tnyvfbef7bht7eHSzUiOanG6jFe/9uJFtzT1v7bbW7G/Je08Wd0tetyQW+W4bP68nf7dNfa58t3Xf2xm+2wY6RLPS6VBKrddae8wc3Z4UryfFCp4VryfFCp4VryfFCu0Xb0erOZyKxWYH0EqeFK8nxQqeFa8nxQqeFa8nxQrtFG+nrzkIIYRoSGoOQgghGpDkIIQQogFJDkIIIRqQ5NAMpVQfpdTLSql3zY6lJZRSFyqlXlJKvaWUOsfseJqjlEpQSr2glHpXKXWD2fGcjFIqUCm1Xik11exYTkYpla6U+tbx/aabHU9zlFJeSqlFSqlnlFJXmR3PySilxju+138qpf5ndjzNUUr1UEp9qJRaopS6u7Xv77DJwfGFHFFKZdbbPkUptV0ptfNkX5jWerfW+pr2jdQZV1vE+6HW+jrgemCWm8eapbW+HrgEGOfOsTrcBbzdPlHWiast4tVAMWDFWH7XnWOdjrHQV1V7xuqIqy3+br91/N3+G3jNnWMFkoF3tda/A4a2OojW3o3nKQ9gAjAMyKy1zQLsAvoAvsDPQKLjS/x3vUd0rfe962Hx/h8wzN1jBaYBnwJz3DlWYDLGYlRzganu/ncAeDne1xV4w81jvRv4vSv+n7Xx/7G3gWB3jhXoAnwNfAVc3eoY2vMfw+wH0KvelzsG+KzW63uAe1pwnnZPDm0RL6CAR4Cz3T3Weuf6xJ1jBRYBTwKfAx/VFL7uGm+t43xdUOCe7nd7OXCJ4/lb7RlrW323QA/gJXePFbgdmOB43uq/A5cs9uNGugH7ar3OAUY1dbBSqgtGwTBUKXWP1vqhdo6vvlbFC9wEnA2EKqX6aa1faM/g6mntd5sOzAD8OI0Fn05Rq2LVWt8LoJSaC+Rqre3tGl1Drf1uZwDnAmHAs+0aWUOt/Zt9H3hGKTUeyGjPwJrQ2ngBrgFeabeImtbaWP8D3K+UmgPsae2Hdbbk0Cpa6zyM9nuPoLV+Gnja7DhaQmu9ClhlchitorV+1ewYWkJr/T5Goev2tNalGIWtx9Ba32d2DC2htc4EZp7q+ztsh3QT9gPda72Od2xzV54Ur8TafjwpXk+KFTwrXpfG2tmSww9Af6VUb6WUL0Yn4wqTY2qOJ8UrsbYfT4rXk2IFz4rXtbG2d6eKWQ9gOXCQE0PkrnFsPx/YgdHrf6/ZcXpivBKrxOtpsXpavO4Qq0y8J4QQooHO1qwkhBCiBSQ5CCGEaECSgxBCiAYkOQghhGhAkoMQQogGJDkIIYRoQJKDEEKIBiQ5CCGEaECSgxDtxLEwy0bHY51SSv6/CY8hd0gL0U6UUtkY8+kfNDsWIVpLrmSEaD8rgU1KqSfNDkSI1pL1HIRoB0qpsRgr88VqravNjkeI1pKagxDt47fADq11tTKEmB2QEK0hfQ5CtAOl1EjgZUADZcCNWusN5kYlRMtJchBCCNGANCsJIYRoQJKDEEKIBiQ5CCGEaECSgxBCiAYkOQghhGhAkoMQQogGJDkIIYRoQJKDEEKIBv4/V1MtDKK/vF8AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "def CG_cond(A, b, eps, k_max = 1000000):\n", + " D = np.diag(np.diag(A))\n", + " C = np.sqrt(linalg.inv(D))\n", + " A_tilde = C @ A @ C\n", + " b_tilde = C @ b\n", + " \n", + " x_tilde, k = CG(A_tilde, b_tilde, eps, k_max)\n", + " x = C @ x_tilde\n", + " return x, k\n", + "\n", + "# Sorry for copying.\n", + "\n", + "fig, ax = plt.subplots()\n", + "\n", + "for alg, alg_name in [(GS, \"Gauss-Seidel\"), (SD, \"Steepest Descent\"), (CG, \"Conjugate Gradient\"), (CG_cond, \"Conjugate Gradient\")]:\n", + " k_list = []\n", + " for eps in eps_list:\n", + " x, k = alg(A, b, eps)\n", + " k_list.append(k)\n", + " ax.plot(eps_list, k_list, label=alg_name)\n", + " \n", + "ax.set_xscale(\"log\")\n", + "ax.invert_xaxis()\n", + "ax.set_yscale(\"log\")\n", + "ax.set_xlabel(\"$\\epsilon$\")\n", + "ax.set_ylabel(\"$k$\")\n", + "ax.legend()\n", + "\n", + "fig.show()\n", + "\n", + "print(\"The number of iterations is brought down by a lot due to the conditioning.\")" ] } ],