{ "cells": [ { "cell_type": "markdown", "id": "81b1d2f0", "metadata": {}, "source": [ "---\n", "title: \"P2\"\n", "subtitle: \"Problem Class 2: Iterative Methods in Matrix Analysis\"\n", "abstract: \"Exercises based on content from Lectures 8-11\"\n", "date: 2026-03-03\n", "format:\n", " html:\n", " other-links:\n", " - text: This notebook\n", " href: P2.ipynb\n", "---" ] }, { "cell_type": "code", "execution_count": 19, "id": "db6fe346", "metadata": {}, "outputs": [], "source": [ "#| echo: false\n", "\n", "using Plots, LaTeXStrings, LinearAlgebra, Polynomials\n", "\n", "function splitting(A,b,x; P=0., Niter=100, tol=1e-10) \n", " \n", " X = [x]\n", "\n", " if (P==0.)\n", " P = D\n", " end\n", "\n", " if ( isapprox(det(P), 0; atol=1e-12) )\n", " @warn \"P is singular\"\n", " return X\n", " end\n", "\n", " iP = inv(P)\n", " r = b - A*x\n", " for i ∈ 1:Niter\n", " push!( X, X[end] + iP*r ) \n", " r = b - A*X[end]\n", " if (norm(r) < tol)\n", " return X\n", " end \n", " end\n", "\n", " @warn \"max iterations exceeded\"\n", " return X\n", "end;\n", "\n", "function SD(A,b,x0; α=0., Niter=50, tol=1e-3, xlims=[-25,25], ylims=[-25,25])\n", "\n", " g(x,y) = dot([x;y], A*[x;y])/2 - dot([x;y], b)\n", "\n", " x_range = range(xlims[1], xlims[2], length=100)\n", " y_range = range(ylims[1], ylims[2], length=100)\n", "\n", " P = contour(x_range, y_range, g, \n", " xlims=(xlims[1],xlims[2]),\n", " ylims=(ylims[1],ylims[2]),\n", " fill=true, \n", " levels=50, \n", " xlabel=\"x\", ylabel=\"y\", legend=false,\n", " colorbar=true)\n", "\n", " if (α==0.)\n", " adaptive_step = true\n", " title!( P, \"Gradient descent (with line search)\" )\n", " else\n", " adaptive_step = false\n", " title!( P, \"Gradient descent (with fixed α=$α)\" )\n", " end \n", " \n", " r = b - A*x0\n", " anim = @animate for i ∈ 0:Niter\n", "\n", " if (i==0)\n", " scatter!(P, [x0[1]], [x0[2]], markersize=5, color=:white)\n", " else \n", " \n", " if (norm(r) 0$ such that $(**)$ converges, \n", "* You can write this method as a splitting method - what is $P$ and $N$ in this case? \n", "* Explain how this method can be written as a gradient-based method.\n", "\n", "Here, we plot the iterates of this method with two different choices of $\\alpha$ for a simple $2\\times2$ system:" ] }, { "cell_type": "code", "execution_count": null, "id": "a5a0ea06", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[33m\u001b[1m┌ \u001b[22m\u001b[39m\u001b[33m\u001b[1mWarning: \u001b[22m\u001b[39mmax iterations exceeded\n", "\u001b[33m\u001b[1m└ \u001b[22m\u001b[39m\u001b[90m@ Main In[1]:28\u001b[39m\n" ] } ], "source": [ "#| include: false \n", "b = [6, 25, -11, 15]; x0 = [0.,0.,0.,0.];\n", "A = [10 -1 2 0; -1 11 -1 3; 2 -1 10 -1; 0 3 -1 8];\n", "L = -tril(A, -1); D = Diagonal(diag(A));\n", "ξ = A \\ b\n", "\n", "ω=1.15\n", "α=.0005\n", "\n", "J = splitting(A,b,x0) # Jacobi\n", "GS = splitting(A,b,x0; P=D-L) # Gauss-Seidel\n", "SOR = splitting(A,b,x0; P=D-ω*L) # SOR\n", "alpha = splitting(A,b,x0; P=Matrix(I,4,4)/α);" ] }, { "cell_type": "code", "execution_count": null, "id": "3572dbb9", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd0AT9/sH8OcyIWFvwgp7gwxxgCLgApW69957rypq3RMnddS9R90DF1oXKKIIOEFl7xkISci83x+xkX5/tlUM0tbn9Vee3N3nHgL6zudydyFIkgSEEELoe0Vp6gYQQgihpoRBiBBC6LuGQYgQQui7hkGIEELou4ZBiBBC6LuGQYgQQui7hkGIEELou4ZBiBBC6LuGQYgQQui7hkGIEELou/afCsLMzMyYmJg/W0qSJN5P7luSy+VN3cJ3RKFQNHUL3xd8wb+lxv7P5D8VhBkZGVeuXPmzpVKpVCKRfMt+vnNCobCpW/iOiEQifJ/3zZAkKRKJmrqL70hj/2fynwpChBBC6EthECKEEPquYRAihBD6rmEQIoQQ+q5hECKEEPqufadBOH78eI4px5JjuWPHDgBIS0vbumXr4UOH+Xw+AKSnpx85ciQuLk55zm5eXl5sbOzr16+V25aXlycmJlZUVChLgUCQmZkpk8mUpVwuVw6CEELoX4HW1A00AUc7R6aEPaXl/DqpaMXCVft+OaAhZbW1ap8hzl23Yn3nLp0Tbj5qbhZYIrq5aN6S4JC2v12542nqm8V7Z2Srb+9gf/3iTQcj54zSV5F9u4nrxLHnYzn6VrlV2VE/LUx9lnr5/GVtDV05Vbp5x6ZL5y5fOn+ZSlAtbDhbdmw+uPfQpfOXCSBatA5YsXb5z1u2X7lwhSCI3v17TZo2acvGrdcuXWOz2WMmj+7arWvM5pjfbt4xMjGaNmequ7v77l17Eh8kWnGtJk2baGxsfPTw0ZSnqS7uzsNGDKPT6efOnnub/tbLx6tLly4KheL69etFhUV+/n7NmjWTy+UJCQlVVVXNmzc3NzdXKBRpaWkikcjLy4vNZpMkmZ2drVAobG1tKRQKAFRWVtLpdG1tbeVrJZFIGAxGk/66EEKocRH/pWuPrl27tmXLlqtXr35yqUQiIUmSz+e72bmfHHGNTqUDQHbF+wWXph0ZdpFCUADg7ru448n7d/Q9QgABAAeTdr4sTl3XbYdyhJW3FigUiqgOqwkgSCCnnR/JNbCb0TaKAEIsE487M6CFdeD4VjMJIMpqS2dcGBXmHD7MbzyFoLwqTlt9J6qDY9fBvqMpBOVG+uVTrw6F2Yb38RpCkoqDT3c9q3jcltOhm2svgaR2e+KGKrK8tUlIO9uOZbWlu55spGvRmxsENrdonV+de+zFXh09bW+d5u7Gzd5Vpt8vuanB0HBhe9vrOr6sTK3SLOVV8RxYbmYaFs/KEl1aOiYlPrHRtNdnGiYVJgwZM/DE0VOmDA6brvWiJGXRioWb1m3RJnUpFEqFpHTt5jWL5y0hRYREJjbnmi1esWjutHnV5TV1UlH78PZjJ46ZOWlWUUExQSXGTRob2iFk9tS5eTl5Orrai1ZE2djYLJq3+G3GWwdHx6Wrl9BotDXL1r5/l+kX4Ltg8Y88Hm9r9LaSopLWbVtNmDwhPz9/3y/7eZW89hFh3bt3z8nJOXn8lLhO3PWHLj4+Pvn5+bGxsRQKJSIigsPhFBcX379/n8VihYSEsFisioqK5ORkAwMDX19fgiD4fH56erq5ubmFhQUASKXSvLw8DoejoaGh/JXV1NTo6Oh8g7+9JicQCFgsFkEQTd3Id4EkSaFQyGazm7qR7wWfz1e9O28M/5oZoVQqTUlJMTMzs7Ky+ppx9u/fb2Nkp0xBAOCJqppZ+itTEABkcmkgt50yBQGAQtI6OHZVbUsnGe1dwpVLCSCYoNnDfYCyZNKYdDljoM8oZWmsZULKiKG+45Qju5l5yerkQ/3GKscJdex88NGuAc1GKMvenoMfnXvQL3woALAY7H6eI/YmbOvbfigAmGibdXHolZz/uL/3cACwN3LKKHmlAHKwzxgA8LFs/vraC2cd134+wwAgGDpMOzeqs0u3cNfuAPAD9BlxvPfIFhPb2IUCwADvEUNWR84J+8nXMgAAqkVVY8b3XxK+3t3MCwAyy98O7jN0eeeNzqZuAHAvM65HeK/lnTc6tHAmgTz4ZGeXsG5LO26w93WSyMSrdy/auHZTVMhqBw/nSkH53ImzxDLRzMDFUzq7p5e86hnRi0JQJwXM7e07MiX3SdtWwYScMspnSkv9Dg+v3Gt3rB2/SjDQYxRHw/Hw6pOH9x1+9+r9D079mTStKb9O9wtpdvf6/Q623UhSsW1tl95Dep05fDbIOlQkE86f8ePoCaP3bd/vy2lRVVdRpigaNnLYtugYFxP3wpp8a2fLsI5hm9ZtttK3ya/K7Tuoj4WlRfSaaF1Nfb64OmpZFJ/P3xa9DRSEli57Y0x0WurznVt3ymRye0f7TT9Hx924tXv7bqlE2iakzdKVP104e+HQ3sNyubxH3+4Tp0w8deLUmRNn6XT6kNGDu3TpcvrU6diLV3X1dEeNH+np6Xnx4sV7t+6bmBkPHzXcxMTk5s2bTxKfWHOte/fpzWQy4+PjX7967eDoEBwcTBDEs2fPsrOz3d3dnZycACAjI6OoqMjDw8PQ0BAACgsLeTyeg4ODchZeW1tbV1dnZGT04Y9TJgMAGu1f828WoX+Lf8c/KrFY3KFDB29v75SUlGHDho0ePbrBQw0cODB65SYFqVBGlBZTO6v8nWopjULPr879uDZJlgtK65VQK66tP5pEJqm3kIQ/TK9J1dtzEkhKvY9jRRIhm6GlKisFFSZaZqqSJ6y00bdVlWKZ2NXEU1XKFXJPjs/HUi7zsQhQlTKprIV1UP2yNTf49x+NRiGpyhQEADZTm0nVUKYgABiyjfSZhsoUBABLHS5X197B2BkACCDcjZtV1FTYGzkBAIPG9Db2s9N1Ui41YBs563o4GDu7mnoAgIupuznDKsK9u7eFLwC0tQs79fTwxDYz3cy8AMBanxt75NzKrlss9awBwMvCt+/+zjv6HjZkGwNAa9vg/gfC9w86o6OhCwBt7MLGRg84POQ8i8EGAP/C1qtWLDww6CyTxgSASy9Ob1i+cVefo3QqAwC2xa+LWbMjJvIgncpQkIoFp6bKSXlM10MMGlMkFU6bP9JUh7Op8z4mjVnKLx7Ya5Cjseu69rs0aBovi9LC2nRwM/FcHrRFg64ZlxEbFNDWUddlut9iOpV+4uyBdsfamVGs+3mMlimkW6O2b92wTVdkGO7Yo6akeljPEZYOlho8dqBFaPnb0pA9YU4ujpRyho9JwP24pI1rN9twreUl4KzreZV/ay17naampqhIwtV12FS2tVmQV052Tm2RyFzb4kVJytjJo29cu1meW2nANsyqfLcqeuX+Xw7kvsvTZLCl1LotOzdvWBWd8eotADi42K+OXvXTgqUvUl+QAJ0iOk6fMy1q7qLkJ8kMBmP0xNE9e/dYPH9J0uMkQ0PDeYvmenl7rfxpVcrTFBs77oIl8w0NDTet3/zmxRtXT9eZc2cAwK7tu7Le5fi38B0xekRtbe3Rw0eLC0rahAZ16tSpurr6wvkL1dXVIaEhHh4efD4/Li5OKpUGBwebmpoKBILExEQKhdKyZUsNDY26urqXL19qaWk5OTkRBCGXy3NycvT19fX19ZV/VFVVVbq6usoj8ACgUChUjxFqWv+OIDxz5oyLi8u2bdtqamq8vLyGDRtGp9MbNpSFhQVNkzL/8pShfmOEUmHM/fUkXX4weVd7u4gqUeWZN0cUdHlcxpVW3OCi6oJHxfckCrEnx8fDvFluVfY7/pvsl++cTd1Mtc2LqgvKJMUHn+1YFLZWk84qrikUE6LDKXsmtpxFISjFNYVAg6tvzndx7QkAJTVFckKWWpjszfEFgOo6XlVdRXFNoZkOBwAUpPxtxRuRVKhJZymbTC5MUkW1VC5JKXzS3buvcpFULntVkhZo1+7Dz0NAZsVbJxNXZUWjUAuq8wzYH+YQNCqtSlhhpGWiLBWgkClkNAoNAKgUqkwhU70sFIKqUHy8m5+ClFMpVFUplotZv/cGAGK5RF/TQFVK5GIjtnG9UmKuY6EqRRKRjYGdqpTJZcoUVKITDMPft5XIxfqahsoUBIBaMd/ByEmZggAgU0h9LVooUxAANGisEPtOyhQEAG2aTjeX3sqSQlD06IbtncMZNCYAaNJZulT9Eb4TlNuaaJsxFJqTWs3RoGkAgLu5l0womxG4ULlyR+eue+N/nh6xUPnjj/CfOOBAl9Ujdih/F5NbzJ1+dsyhIdHKnTKpzF8ebt3We/+Hl0Va9+T9w2URG5Ul716ltFA6J3ihslx0dYaZgc3M4EnKcuKZIe0cOvQNHQoAcoV8xKpe/XyGdunYEwBq6qrHjRwwqtXkmRFLASCr4l3vrn3Ht54xrdtiAEjIutuxbeeJrWZN7PYjCeSx5H3BLduND5g5uuusOlndpj0rNq3bNNp36sAOE0r5xfMmLhDKa4d6jpvv1zer4l2P8F5Agf6uI/qajnrz8kXbFsEkQHeH/v76bVNOP2m3P6S6qibctrsJm7Pzt717d+17lfYq1Dpcm6E7btfENhGtr5y/2soimElhrly4esTEYXu37/cxD1CQ8mml02cvnLVuxQYXI/daMZ9PrZ49f+bi+Us4OlaVteX27nbDRg2dP+tHNlW7SljZd2CfVkEtF86OkoplQIX5UXP1jfR/+nGpUCDSN9Bbs3F1dU31qiWrq3nVLm6uazetzsjI2LRmc2VlZbvQdguW/Pgk6cmubb8IBcLwyPBxE8c+uP/g6P5jcoW8z6De4eHh8Q/iz5++wGAwBg4b4O7unpSUdPPqTT1DvT59+xgbG6empibEJ5iamXbp0oXJZKanp6emplpZWbVq1QoAcnNzMzIyHB0dbWxsAKC8vDw/P9/e3l55RE4oFFZUVFhYWCjDW6FQiEQiPDT6n/HvCMLk5GTlH6uOjo6xsXF+fr6tre3fbvVn8ovyO3XqtOTmbKo+VUMnNO3Z/qOH9568ssfAwGD3iV0ODg6rl61ZcX+uubn5nhO/GBoaLpq7eOfVaA6Hs//EXqlUGjV3UVUlz9DI4Mivh95lvJuzepxCRuroax8/e+zyhSvjTg7Q1tAmmHD09OFtG2Omx17W1tCtllbuPbpnadTyi291GVRmfm321l1blkXN8TDxkZOy9MqXUcsWzN46LtA6pFbCTy1L6jW4x8Jr04I4oeV1JY9LHzg6O21OWOVv2iqfn/NakKpFsI+n7vcwbva24k0pFJ1/d4JOozsYOqcVP6vTEO1Iih7vP9NKn5uY94ChTVv/YOmM1gsN2ca3313T1tfalrB2YsvZDBrjRvplTW3Ns8+P9fAcQADxIOt2HSFKyn3Y3LoVALwoSc3kZeTzcpWhlVP9/lHe/aGSscqJbKmo6F3umwj37sp4qJXV3Mm60dymlfLllYLkUe4DrqG9sqTQiJSCJ4G27ZQlQSFyKjOV0UgCKVGIa8V8LaY2AGjQNXmiKrniQwbraeqX8ItUvzUWnV0h/Dg7p1CImrpqValQkFL5x9m5XPGHW/TKFDJVZAKAXCGtn+skCYzf81VBKuhUhupNQLWoyljLRHXkvKy21M7QQbVhdR3Pw7xZ/Z02tw6sP2wrm+CPTciJtnZh9Xoi2jtEKB9SKVRSBp2cI5WljoYuoaC0d/yw1FrflqqgBdu3V5bu5t6aBLu1bTAAEEAEcUNS85624rYFAA2aRii3M1VGU5ZmOpxWnGChpDbILhQAPDk+LrqedoaOoY6dAMBCz+rm69jOrt3aOXYEADczr/unbo9tPV15zKCVbZvBhyOXdF7naOwCAB1duvbfE7Gpx27l30MX157Dfuq5o+9hE20zAMitypo5ddwv/Y4p34E9yPxt6tjpO/oc1tXUB4DDT3dPGjVla/d9upr6ClKx8tKCU4d+Xd9lh66mvkgqnLN4AkHAsg6b9FkGJfyiYf2G62rqR4WsNmQbvyhK6dC2oxHbZGbrRUaexrfTr7dp0VaXqj/Wf7qWqfal02fCTrQnayiDvUZTCdrmH3/eu3NfRSYv0qmvVC4Z1nOkdwvPzJSc9jYR2eKSdhtD24UFv3z0ppV58NO6FysWr2wb3Db5Xoq3sX++4Fchi+/q5pp4J8nFyD2j/HWzQC8qlZpw+6G1ge27svRxk8e8fv3m3q37JtqmJbVFazauvnbl2u0bv2kyWEwtRswv2w7sPXjr+m0KQXF2dVq/dd3Wjdvirt8iCKJj5w7zouZuWB0ddz2OzqAPGz1s8NBB61dvuH3ztoGBwZRZk1sHtt4cveXhvYccS87M+TOsra13bd/1LCnVwcV+0tSJLBbryOEjr5+/8fTxGDR4kEKhOHP6TF5Onl+AX/v27SUSybVr1yorK1u0aOHq6iqVSuPj42tra1u0aGFsbCyTyVJTU+VyuZeXl4aGhlwuz8rKotFoNjY2ysNU5eXlLBaLxfrwr6Curk71sfr36d8RhNXV1ao3X1paWjwe7ysHvH79OgCsPD6h+JXV/v20qVPHjZ84TrV0TfTq+isfPX2kfnnn4W+qxy1btRw8bHD9cvnqZSKRSPkXFtQmiM/n19bWmpubA0DHTh0zMzMlEomjoyONRuvarWtycjKVSvX19WUwGH369bl//762tvbu4Bgmk/li2Iv79+77mbhu7bJeU1Pz+vXrD+8/8rVzW9d/GY1GO7j/YOrTx04tnR6NSRCJRFuit156c9w7yDvx9MPc3NwNq6LPvTraMrDFw+MJCQkJm9ZsrKqqahfaLuFo/JGDR+fvmSCTydqEtLm/7+7q5WsmXxkKAK2CWt7Yem3OtHkHL++QK2QBrQNOX/p16vhpNBmDL6r2CfDZOHXDrLljufp2xdWF7r5uQ/oMnLpruJupV05VpqO3PVODGRU3zdnQI6PypXtLlzf5qVserrZm26ZVPG3Rzv/4s305vPcmLPNHRfeCwgJX31/Y03mQroberezYtqFBi+NmDPIco0HXPJ9+3C/Ad929JQO9RipIxZHUPaY2poeSf+nh1l8gqT3+Yl81pfLOuxtt7dvX1FXfzrpaxC/IKI1wMnGtFfNfVqY8KqxuyQ0y0jIRiGsLhTmnXhz05PgwaUyhRFAj4515eXRa0I8EECKpUE6RXU0/38NjAAAIJQKSIk/Of6wMAKFEIFaICnh5FnpWAEChUIv5hQJxLZupBQCaDM2MsteqyTqdSn9d8lz125eT8syytwC9lCUJinxeNsDHI9XltWXKXAEAKoVSXcdTzd0JAiQyMY1B+33bj8fYSSDrnwEjkYmVc1mlWnGtjoaeqhRIak21zeutLLHU56pKqUzCrTc7l8okjsauqlIkFrmbeX/cr4JUdUuhUBgURv2pvIGmoTIFAYAkwc7QUfWzsBlavhYByhQEAFOWWXv7CGVJISgWLOvm5q2VpSadxWFZtXeO0GcZAICptrk+1Wh8wHTlQQIP82YUCW1u56XGWqYA0Mm524GHOzcM/EV5kGCk/8Q++zoeGXpJOdGfb7R8wIEuJ0bEKg942Bk6zjo/7vDQC8rflJm25Z4rW3f1O6b8CJ/1TCvxesL6yA8nwW26s+L1vbfrO+1UlnMuTeDoWW7usg8A5Ar5uOiBLaxbb488AgC1Yv6EMYM6OHfZ1f0EAORWZXeP6NnNrdeOyKMEEI+yH7QP6hDp2vfnLocB4Eza0TYBbbs69opuv0csq9u5f9P2Lds7c7svClhfIShfPG2piBC0t+o6nDulsDq/d0RfuiatHadTiGnXzKS3bVsGa2hotjIOdjRwffQ8ZVfMbqlE4mfY2lLLZseFPbu2/fL+faa3gb8Bw2jPhv1tI4KuxV530fNg07QWFiwaP33s7h17uNoOdAr9ddnzJasWr166xphpJlPIhNTa5WuWzp/5I5NkCer4ju4OM+ZOnzl5lrhWKpQIevXr2bNPj5mTZ1WWV1FplBlzp3t4e8ybPr+0pMzAUH/52mUsNuunH5fm5eW5ubktW7uUz+evW7G+qKCoResW8xfNy83N3bbxZ15FVUinkLHjx7x58+bAnoNikTiie3h4ePirV6/O/XpOQZLde/3g6emZkZFx88ZNBpMRGRlpamqalZWVkJCgq6vbvn17DQ2NwsLCp0+fmpub+/n5EQRRWVn59OlTV1dXS0tLABAIBLm5uTY2Nqog/3r/6CBUnSlkYmJSXl6ufLKsrMzU1FQt4/uYN7spil8TBWPGgKamWoYEgiDq/3q0tbVVJzsRBGFvb69apKmpGRj4cQJhYmLSq1cvVenh4eHh4aEqO3Xq1KlTJ1U5ZtwY1WM2m7189TJV6ebmtu/IXlXZuXPnzp07q8rJ0yZNnjZJVW6O2VS/+fOxZ+VyOYVCUb5tfPTsYWVlpZaWlvLcjYiuEZmZmSYmJspPfYaOHPrq1SsbGxvl6Uvp6enp6ekznCa4uLiQJHnnzp03b94MatHT19dXJBKdPXO2qLB4cdCC1q1bFxUVHT54uKIib86UGSEhIc+ePTu894hAIJy2alLnzp3PnT138vBxgiAmLx/XuXPnbZu2rTsXxWKxJi4eGxIasmj+4unXR2hpa01ZMtnH12fm5Fk5SblMDcbspbNMzUznz/xRUCOkMagLV/8oFIqmrhzGZmiL5IKVm5Y/fZw84dwgc12Lwuq81RtXnTx26sebk020TNPLXq+KXhmzeftvudd0GLrJhYmLVy5auXl+c04gnaDH5/02fd60+QcmhdlESEnp7ezYnv27L46b0YHbtUZcc+X9aUd3x58frQ+yCiutLbqZc0lbV/vM82P+Fi1zqrLSqpJf1aSa6Vi4m3mnl74qkOQcSt1pxDayM3JKLXzKJ6p3PI5e0G6VnqZ+WmEyoQE7EjdObT2PTmWkFSbTNKin0g739RoCAC+LUuUU+eOc+ACbQADIrcouFRblVeVY6dsAQKWw7GVJqiqqa8S8+OzfhgSMUQaAUMZ/kH27g/OHyaWUlDwteNTM0l9ZkhR4WZyqjHwAYNAZ78rS3c29fv9ThkphhQHLEAAIIKRyqVQuVZ5lpkln8SUfL5ZlM9m1dTWqkkahiWQfv5lBQZJy8uMEXTXIh1ImqT87l8glqgRVLtVnGapKClDqHSqXadJYqkPllcJyC10rZQoCQKWo3MXUQzWVr5MKA6wCVSfBEUC0s+/4cViS2vGP58R1c+mtfEylUGlyei/PD292tZjaIKOoTnOz1ucqJIr+zYYrywCb1lvvrunpOUBZdnKKvPr84g/u/QCAxWD3cOsfc299pHsfZdnNsc/9zFs/uPUFAFNt89aF7QRSQU+PgQDgYOycnPvY1tChn/dQAPC1apF64Wknp+7hrj8AQIhTx7EnBwwLGKc8yvKDZ98h+yLntV/mxfEFgH6ew4ZGdV/RZbPy45ISftHkMcPXRf5sa+gAACkFTwb3HRoduUv5x3Ph5cl+P/Rf13WHhZ6VglRsvb2mz6F+qzpvtdCzEkmFi9fOrBJV/hS23tLfpoRfNG7IeCCI+W2Xc93tXxamhodFaFA1p7VcwPGxTMx4EBjQhq6gj/WbbmBs+NuZG2HH2teWCwZ6jGIx2DELdx3ad/ht2vseTv2pBGX0yXEtQvwTbiV2tO0mkYvbr+3Yve8PV89eD7QM4Ut5C+dE9RvQ78zRc36cFsWCwlo6r3OXzsf2H3c18cjhZbn6Oju7OB3ad4RraJ9d8X74mGGz588CdWjiIDx27NjBgwdTU1P79eu3ZcsW1fPnzp0bP368XC7ncDgnT57s2LHjihUrJkyYkJ6eTqVSORyOWvbewjM8Jje2ZQvZnj20KVPUMuS/HpVKrV8aGHz8LJBOpzs7O9dfFBT0cbrj7OysWkoQREhIiL+/v/JNgKam5qDBg1Rrmpubz50/V1X6+Pj4xHw8/adHzx49evZQlTPnzpw5d6aqjNm1rX57F69dqF8+TE6oXw4fOay2tlZLSwsAevftvWTF4tLSUktLSwqF0m9gv/z8/PLychcXFw0NjUFDBiUnJ/P5/Bj/jdra2kOGDrl3755UKl3d9id9ff0hQ4dcv36dQWdEhc8yNjZ+8uTJ9SvXrQwMb/WNMzExOXf23G837pi7m8Xtuamvr78jZsel+ONcB+6t7TdpNNrqZWtuPb/s5OJ4635cVVXVsqgV2alZnl6et+7HpaWlrVm2gM+vdXVzvf3g1rnT52bsGU3KFfZODjfv3di8bvP4CwNoVLol1/LijfM/zlp44tV+IAh9U92jZ47MmjzbiGkikooYurSNP2+YvXCsm4lXVV2ljCUeOWXEjP2j/TktiwT5tYxqxwDHpb/N8TT0za3NlOqL8ilZOxKjHXRdM3gvtS01L2WfqpbwbHRtU0qSzOxNtj1ePcB9pImWWXz+Ha6j9fLf5o72naarqXft7QVbR+76+z+N8pvMpGn8+vywmYXJgac7B3iPIEnF6edHpEzJb++uhzh0ksoltzOvvuO9zih97WTiKpVLU0uS3pS9ivTobcA2kilkeYKsVxWpbexDGTSmXCGvklZcTj/jzvEmgFCQCjGI4t7FDvQZCQAKUgFUeJR9P8guBABkCpkMZFkV75T/pysUcqFUUC2qUgYnm6FVUJ2n+ghci6mTWfFW9ZfAoDJyqjJVJUmSRTX5H/9QCKJKVFm/FEgE9Sqi/lF3kiQBPs7QiXqPJXIJk8pUldV1PD3Wx1Dniarqf3YukNTaGTipSplC5mLi9rGUy+rPzqUyqY9Fc1UplogDrFsrH1MICigIZQoCAIPGZFCYqpMGWHS2voaB7e/H8/U1Dbl69soUBABrXbtm5s0/HPwgKM767gZ0I2WpSWc56rpZcW0s9WwAwFTb3JLFbe8coRzKg9NMh9Qf12Ka8vS69o4Rxx7vW9xlnXLbIT5j+h+I2NnvqJ6mPgC4m3v32dfp0JBzyg9W/K1bDjnc/eiwS8oDG804/j/umnZo8Dnl2yO7t04Hd+7Z3e+E8k3MwaRdJ3b/urX7AWW5PO7Ht0+ytkceoRAUBalYdnxuQOvmbSBy8RUAACAASURBVNu2ha/WxEHIYrFGjRp15cqV+l83VV1dPWzYsAsXLoSEhCxfvnzMmDEPHjy4ceOGcv60Z88ede3d0NBaX0HtP/rO4IFvT536zdhYb8GCsf7+/uoaHzU5ZQoqaWhoWFt/PLhnaWmpPNICABQKpf7vXUtLKyIiov6ao0aNUpX+/v71V+7Zq2fPXj1V5YzZM2D2xwY2//xxzm1lZXX64ilVaWNj061bN1U5fdb06bOmq8qYX2IAQC6XK9+aXL11RSAQkCSp/ImePE/Kzs7W0NBQHnXv0q3L48ePzc3Nle9FBg4d8PTpUw6H4+fnBwDPnj1LTU2N5HYIDg4mSTI2NvbNq/QhHv3Cw8OFQuGRw0ey3+VEDujcs2fPoqKinTG7Hhe+CRsVNGjIjuTk5B1bdlWX8ToP6bx31PYzp8/s3721rq6uR98eO8ds2bJx64JfJ1Ep1H6D+24eum7hnKhJV45QqZRBwweu7P7TtAnTCxILgSAHDx80J2jGrEmzQUrUSmpHjB5uZm42ecUwjo5FAS9/1LiRvCrelLPDuIb278syRk0e+TD+4eLbMy11rJ8XPxs/fdyZU8dSyh7rM40eF8SPnzp23dHFbSzDWDTtO7nXho0btvDctC72vagU6pV3Z37o3W1x3Mzuzv0kMsmvrw+5+7nGPFzXySGSJ6w68nyPhg7j7IvjbblhRTUFd3Kvyyly92xvf6uWWRXvX1alvKxKcTJ2szNyyKnMLK7LP5iyY4nBOl1N/XxejgBq9j39eXabxXQqI5+XA3TyZNrBwT6jASCvKoekKuKz7wRy2wFASU2hQCZQRXWdVJTLy1JFtUwuTSl8oopqqVz6KPd+f79hyl+3RC5+Vpik/LwWABSgSC996WH+IQupFEp2Zaby3DoAoNPoRTWF1r8f9CYJsk4q0qBrAgCFoEgVUhJIZUIzaMw6WZ3qj4pKocrrnSInV8gpxMc3vhK5WLPe7FwqE2szP16DK5FJDOudEyeW1pnVy3WJTMLRs/y4I6Dq/T6zl8ml2gwd1Xny1SKetZ6t6vA+X1zjZd6s/kGCNtxQ1VSeSWGGO/2gKtkUdie3rsqSQlAiHHpcvXhNLUH4j7igfs6cOTweb/fu3cpy3759u3btSkxMBAA+n29kZPTq1av6BxX/zIEDB2bPnu3l9eHYDoVCmTt3buvWH946KS+oZzKZ9TeJuTQ392XBL2v85fLBAKX6+rMPHpwfEhL8v0OjL6eajaFvQCgUampq/pMvqBcKhapPDeRyeXFxsampqfKyyOrq6ry8PC6Xq/yDef/+fWFhoZubm6GhoVwuf/jwYVVVlZ+fH4fDqa2tvXHjBr+GHxIaYm1tnZube+XSFblcEd6ls729/ZMnTy6evcRkMnv37+Xk5HT29NmbV2/p6ukOHzOUy+XuiNkZfzfB1Mx04vTxxsbG61asT36abG1tPWvBTAqFsmzhivz8PA6HE7V8YUlRydoV60TCOmMTo+XrlsXfTzi45xCFpOgZ6q7asPLw/iPxdxNYDDZVk1ixbvnaFetqy4QsBosnr4patmB51AoLtjVJQqEwd9KMiZtWb/bhtBDJBO+r03v3733u6PmWFsE8ScXLyhT/AP/slLwA86BiQcGzikRDI0NTmZWXse97XkZS2QMFqQi1CHcycn9Vmnq/OE4skgz1Gmetb/ukICEuL5YqpU5uMc9Y2/Ru5s0bORdNmZypreZrMbVj35y/mnnO3yxwhP8ECkE5/+LE5fQz/T1HdHTuSgJ5KuXQpVenF4SudDPzUl4fHJd+ddMPe4y1TEggtz1Y+7woJabHQeUB52U354ll4hXhm5SZOuviWGcTt7EtP7xXG3dqQHfP/sqjtQAw5Ej3uWFLPM19AEBBKvofiPi5zyFjLRMAkCvk/fZ3PvL7FFAoqR15rO/x4VeUw+ZUZG64vWxbnwPKceLf//Yg6+689j8pyxNPDpIEOcBvuLLccGtZC25QG/tQZXn3XZzYo3Lugo+Hlz6JxWL97YU6/8QgXLhwYUFBwYEDB5SlnZ3dzp07O3bs+Kfb/+7SpUvLly9fuXKl6hkfHx/Vwb1PBuHZ2C0XMm4fmXXp9ycKvb1HJCd/+t406Is09s0gUH14Z5lvSSqVVlVVmZh8uDCpoqJCJBIpDzDI5fI3b94QBOHs7EylUmtra5OSklgslp+fH41Gy8vLS0hIMDQ0DA4OptPpSUlJjx4+srC06Nq1K41Gu3TxUnLSMwdn+779+srl8j2/7E1/ke7h4zFi1PCKioqt0dtysnICWjefMHnCy5cvN6/bWlZaFhzWdvqsabFXYndu/UUoELTv3H7WvJk7YnaePnGGJMkOndrPmDt90fzFSQ+TgCDCOoaOmTh6xsSZJQWlErk4rGNonwF9pk+coQksnogX2rFdy9YtV/60ytbQIb8qt21YkKGR0YVTF51N3N9XZLQMCSgrKSt9X2Gv7/Si5FmL0IDHiUlWdFsLltXjkgSfIM/fbt0JtuyozzS8lx/nEuDw8G5id6cBLDrrauZ5K3dOVlpuP/dhFIJy6tUhHQ6bwWP3chsokooOPNsu15C2MGjbwakrT1gV82htjbx6mMf4ltw2pfzi1XejaiX8hcGr7IwcKwRlC29MoxK0FR0367MMKoUVS2/PPnTugKur61/+ruBzLlf9Jwbh5MmTFQrF9u3blaW3t/eCBQv69ev3t+N8zi3W/icIr984v6l43+3Rv0qlH563smqem5vU8B8G/Q6D8FvCIPyW/u23WJPJZPVvUVRZWamtra28OFsqlWZmZpqbmyvvTVhWVpaRkWFvb29mZgYAr169ev/+vYeHh62trUKhuH37dkFBQUBAgKurK5/PP3f2XGVFZXBIsI+PT25u7vGjx0WCuojI8ICAgEcPH508ckqhUPQZ2DswKPDUiVMXTl9kajCHjxnWomWLzdFbbsTe1NfXmzRjooenx5IFPz1+9NjQ0HBu1ByOBWfOtHl5Obm6erqLVkQBwOK5S4QCEUtLc9napZ8zQfoc/8QgXLp0aUZGxtGjR5WljY3NgQMHQkJC/nacBgRheXn5xNMjsg9MSkrsDAAEkRgauiku7sRX/TwIADAIvy0Mwm/p3x6E/zrf471G3dzcTpz4EEUlJSVFRUUuLi6NtC8jIyOuRJvZctvTpCKSLDY1/XXPnrONtC+EEEL/QE18r7+cnJy4uLicnJyCgoK4uLj3798DQGRkZHV19ZYtW/Ly8ubOndu1a1flqXGNJKJ5hJYXe9MmUVgYd/z4e1wut/H2hRBC6J+miYPw8ePHa9euraqqkkqla9euvXfvHgAwmczY2NgrV66EhoaSJKk6ZNpImntGZNPrxo8bvmzZgGPH8CxHhBD6vjTxodE+ffr06dPn/z/frFmzGzdufJse2FoGHDkz7fWtVq26EQQkJUHz5n+/FUIIof+Gf+JnhN+eXhlzwbo1rRyf9Oo16/BhHQxChBD6fmAQwuAhAy8/uGoUaPE26WhZYjRdcX/DBh8G4+83RAgh9B/wvQdhVlbWhZuXfNYFE1QKAOj5mWbtjLx2LS8ysqk7Qwgh9E18798QffnyZQN3I2UKAoCOs6FUUXX4cNM2hRBC6Nv53oPQzc1NVPLxft9SvphBp1+/XnDs2KUnT578E+42gBBCqFF970EYFhZGrYGiK++ltRJhYe2b6CR/z1YSSb9Ro5LDw7c1b965/tdiIIQQ+u/53oMQAF6npdtWmacvTcw5/WZU/zGpqUKx+G5d3ZLy8oMvXvRduHBtUzeIEEKoEX3vJ8sAgJ6e3vWrcaRYNHNffwMDYz6/C8CH7+gSi3vduvWJyxwRQgj9Z+CM8AOCqekp0ayklGhoFNV7utDU1KTJekIIIdT4MAg/asa2ziaLTEwSqdTLAHKAbF3dSQsXjmvqvhBCCDUiDMKP7MxcZHLJ6Ru/9Ot3zdY20MZmsqPjonbt2jZ1XwghhBoRfkb4EcOc68XTypbmHz0aAwAyGXh4wO3bEBra1J0hhBBqNDgj/IhuxnWvkj0tTlWWNBosWQLz5wNeTIgQQv9hGIQf0Uyt3Ip4T4tTFb9HX//+IJfD+fNQUlIil8ubtj2EEEKNAYPwI4JGN9I206Ox31VlfniGgA4dzvbt28zbe5S5ud+yZRubtkOEEEJqh0H4B3Rzm2ZMc9XR0YyMjL17N8lkD0pKLpeVPdm0Kfn8+ctN2yFCCCH1wiD8A7o511Os8aQ4RVleu3arqmoEgPJr62k83ozDhy81YXsIIYTUDoPwD+jmXLdy0avydLFcAgA0GpUgFPWWy2k0alP1hhBCqDFgEP4B3ZxLK8qz1bV5UfYaACIiOurr7wUoBwAAkb7+qlGjejZthwghhNQLg/APaEacF5k51xee69ymk527/dmL5w4dWm5r29XYOJBCCZw6NbJjx/ZN3SNCCCF1wgvq/6C6pmbU7VTTqc0szdgKqeLnQzuXjlmcmflILBZHRTElkqbuDyGEkLrhjPAPEhISNH1NNc3YAEChU4y7W+89vBcAmEzmmDGwfz9IpU3dIkIIIbXCIPyDuro6gvnxdBgKnVJXV6d87OQEzs5wCU8aRQih/xYMwj9o3bq1MIUnE36Y91XeLu4V+fHsmDFjYPfuJuoMIYRQ48Ag/ANTU9Nty9blLXuUvys9bWlCKwO/mVNnqpb27g3JyZCd3XT9IYQQUjcMwv8V2XdA/OCw69tPtlgUunHLJgrl40vEZMKAAbBvXxN2hxBCSM0wCD9Bw9zGSoPiY+WZVvryfxZFRuZu2DDZ27vT0KHT8vPzm6Q9hBBCaoRB+Al0c1tpUXYzU8+U0uf1ny8pKRkypIdIFJmWdvzIkc4tW/5QXl7eVE0ihBBSCwzCT6CZ29RkZzQz9XhW/Icg3LPnSHHxNICOAAYkGV5cPP7gweNN1SRCCCG1wCD8X48ePgwbM6PnwnW9mofmVxZW1fFUi96+LVAobFWlXG739m1BU/SIEEJIbTAI/6CysnLC0IE/B1qe7ul7pasrPZe36dhm1dKwMH8trauqUls7tn375k3RJkIIIbXBIPyDhISEMAstMy0NAKASRAcFPSHzoWrpoEH9W7R4Z2Awgkr9mUIZ6uZW1KsX3oMbIYT+3fBeo3+gUCgoBKEqHatkcgdNVUmhUOLiTsXHx2dkZLx9OyE+vlW9dRFCCP0r4YzwDwIDA+Py+eVCCQAoSPJ8YhFNl1kuqvyfdUaMGLF8eavycrh+vYkaRQghpCYYhH9gaGi4Zc+BsfGFg66nR5x8zAmKaGHtl/r/riYEACoVli2DefNAofj/CxFCCP1rYBD+rzbBwY9fvjn/ICl2aPDSZUubmXimlDz/5Jq9egGTKRg9eueUKVGnTp1WYCQihNC/EAbhp+kbGmlaO0myXvuYeT77kyCsrKzMzW134IA4JqbN6NFJoaG9SZL8xn0ihBD6ShiEf4pp7yHOfGGnZ1Mj5pcLK/7/CitXbi0rm0GS0wA68flrU1PN4+Livn2fCCGEvgYG4Z9i2HlI3r+QSWXyx4LBY4euWb+mpqam/gpPn76Wy1uqSh6vVWrqq2/eJkIIoa+CQfinmFxXYd67oNCgt3nvShyqD2X86hvoV1n58QxSLy9HCiVFVerqJru7OzVFpwghhBoOg/BPEUzNuxWyKsNa8662Oo4Gxu0smG11N2yJVq0QFTXV3HwFjXYQ4DGVusrR8WXHjh2bsGGEEEINgEH4V94pGDSrjxfUs7g6qS/TVKWJiUlaWtyCBcW9ex9ksUwOHbpMpVKbok2EEEINh0H4V5r5t5K9//i5oDCzxsezWf0VDAwMli6d9+uvP48fP/rAAfo3bxAhhNDXwiD8K5HDxhi8rym5kMN7VV58I0caz589bdYn1xw3Dvbvh7q6b9wgQgihr4VB+Fdo2npnhnab32aYdbGJk7H9s4fJenp6n1zT3h6aNYOzZ79xgwghhL4WBuHfYDl69va2j14fbRxiyWaz/2LNceNg165v1hdCCCH1wCD8G0w7D/H7F7a61gW1RVKF7C/W/OEHeP8eXrz4Zq0hhBBSAwzCv8G085BkvqBT6WZsk9ya/L9Yk0aDQYNqpk7dHRW16u7du9+sQ4QQQl8Dg/BvUA1MgEaXlRXY63HfV2X/xZoFBQWHDrW7c0e4cqVrjx57R4369Gk1CCGE/lEwCP+GXC6/VCyZM2Vi2Zvi9LK3f7HmjBnLS0rWkeQ0gB5VVYcuXkx/8+bNN+sTIYRQw2AQ/hWFQtElrN3rtNQW0kIi/sGJ6yd4PN6frZya+pIkg1SlUNjmBX5giBBC/3gYhH/lypUrXEnlrOY2QVaGI030tc00YjZF/9nKdna2AB+Tj8VKs7e3/yZtIoQQajgMwr/y5uVzb70P94sxlADQKc/fPPuzlaOj55uYjCOISwBpBPGTh4fYx8fnW3WKEEKogTAI/4qTq/uLaqnyMUECu0pi3uxPJ3lubm6JiWfGj0/o2nWLq6t1r14nv1WbCCGEGo7W1A38o3Xp0mXr2lU/pxa2NGGlVwoKQe+Hoa3/Yn0ul7t9+2oASEuDTp1g2DDQ1v5WvSKEEGoQnBH+FRqNdu3uA7v+E6/KDOv0zKaPnFksKf2cDb28IDQUtmxp7AYRQgh9LZwR/g06nT5m3PiR/XqXRk+p5LjefvLgMzdcuRL8/SWurvepVH6rVq1MTU0btU+EEEINgzPCz0LVMwKCsCG0s6tzFaTiczZhMAql0jaDBsUOHvzSy+uHw4dPNXaTCCGEGgCD8HPRrRyphTkGGvr5/KLPWX/48Nm1tdFicbRAsLC09LfZs9fU1tY2dpMIIYS+FAbh52JYOUrz3trrczN52Z+z/uvXbxUK1fX1mnJ5a7y+HiGE/oEwCD8Xw8pBkvfWTo+bycv5nPXZbBbAx2+3p1LzzM3NG607hBBCDYRB+LnoVs6S3Aw7PZvPnBEuXDhRT28oQB6AkEqNcXGh2NjYNHKPCCGEvhieNfq5qDr6BI3OJbT/+jsoVIYM6aelpbl69YTSUn51ddiVK0cbuUGEEEINgUH4BRjWTjoV1VV1PIFUyKaz/nb9Hj0ie/SIBIDmzeHBA+jcufFbRAgh9IXw0OgXoFs6yvLfcXWts3i5X7ThhAmwY0cjNYUQQuirYBB+AYa1ozTvraFU71LcpaKiz7qIQmnAAHj4ELKzG60zhBBCDYVB+AUYVk6zj57Zv3DX7qN7W4UHzZw/8zM31NSEIUNg9+5G7Q4hhFBDYBB+gSu/3btPETvO87Ea4GQzx/1iUuyNmzc+c9uJE2HPHhCLG7VBhBBCXwyD8AvE3rqq1YajKln++rE3Yz9zW3t7cHbOHDx4fVTUqpSUlMZpECGE0BfDIPwCJgbGMoFUVcoEUmNDk8/c9ubNW6mpA0+ftly50qV9+7k7duxvnB4RQgh9GQzCLzBk4BD+9aK6UiEAiIoFgt8qBvUb+JnbTpiwuKbmIsAAgJ4VFRd/+mmzXC5vzGYRQgh9FryO8As4ODic3HloytiBbwkaU0/zzMFfuVzu52xIkiSfLwVQTR81CMIxPz8f7zWDEEJNDmeEX6ZVcOjlwZ3jr10ImB/s7+//mVsRBMFiUerdelSmULy3sLBopCYRQgh9PgzCL/ZGoXX3+AWBWFAqLP/8rdaunWdg0B3gNkACnT5g8uTBNBpOxxFCqOnh/8VfgCTJIX168tNTvA01Re3MJidNOLn+FEEQn7Nt37497O2td+w4IRRK7twZ37ZtWGN3ixBC6HNgEH6BUydP6BSlrwlzBgBLDeoVUemF8+e79+jxmZv7+fnt2eMHAIcOwfz58PAhfF6GIoQQakR4aPQLPLr7WxhHW/nYma+g2bAf3r3dgHEGDwaxGC5eVGtzCCGEGgRnhF/A3MqmMDtR+ZgrJKtZVEMbzl9v8kkUCqxYAdOmFRQVXdPQoEVEhJuYfO71iAghhNQLZ4RfYPDwEfve1SQV8iRyxZOCKlG+0KtT84YNpVBczc3tNWmSeOzYak/PLvfuPVBvqwghhD4TBuEX4HA4p2JvXCAsB1x9c6GC0qtN7zxJYcOGmjBhoVR6VaGYKJVOLS29OHLkXPW2ihBC6DPhodEv4+joeODkaVFavPDxzSyn1gdfnGzAIOXl5TKZKYD+70+Y8/mERCJhMBhqbBUhhNDnwBlhQ9A5dpLCTHdjl7eV7yVyyZdubmBgQJIlAIrfn5BQqSJMQYQQahIYhA1BMzQjhbVMicxa1zK98t2Xbk6hUEaM6KmjMxbgLcBrNnvIjBkjG6NPhBBCfwsPjTYIQdDMudKiLC9jt7TSV57Gbl86wJo1UV5ev+7Zs7i0lMpkDp4zp1tjtIkQQuhvYRA2EINjJy3I9OS6Xc2Ma9gIAwf2GTiwj1gMNjaQkQFOTuptECGE0GfBQ6MNRLewkxZmepm4vSh7oyDJBo/DZMLw4bBrlxpbQwgh9AUwCBuIzrGVFmRpkMyKB0WT5k2+cuUK2dA4nDgRDh4EgUC9DSKEEPosGIQNROfYVua+92nlW1ZWdl/+ZOaOeV16dm1YFlpbQ2AgHD+u9h4RQgj9PQzCBiIYGj+/KqS1Zlt0sTP0NTUbYPeuLvv27YbcehQAJk2CmBj1NogQQuizYBA23IsaEdtRT1VSuMyUtNSGDRUcLCkoiHZw6BgU1Ov8+UtqahAhhNDfw7NGG86Fa/cgr1zDhPWhLpa5dHVu2FC9e4+urnYtLz/5/n3Vy5ezKip4o0YNUVujCCGE/hzOCBtu7viJNWeyKh4VCXJriq/laJdpdOrUqQHjlJeXJyYWSKU/AugD2PF4R1at2qH2bhFCCH0SBmHDWXsHXO7Ttptee3YSxdrA8sGt+zRaQ2bYRUVFBGFT7wm2SCRVV5MIIYT+Gh4abTiqvrEhg1g5b14tndrvwmgKrYHvKpydnQGSAQQAbAAAeGZr25CvOUQIIdQAOCP8CgRBN7eVFmbrMLWtdCxeVqQ3bBgGg7F58yJj4/YaGhsYjCg2e/Thw5vU2ylCCKE/g0H4VegWdpLCTABobu7zpCilweMMGNArJeXMgQNWR4+20tJ6UFNjp74eEUII/RUMwq9C59hJC7MAoLlZs68JQgDgcDj9+vXr3bvLjz9qRkWpqT+EEEJ/B4Pwq9AtbKUFmQDgbuyaW5NfI+F//ZgTJsCbN9DQS/MRQgh9GQzCr0I348pK80i5jE6heRq7JRenff2YDAYsWlQ3dOj6gIBuXbuOSEpK+voxEUII/RkMwq9C0BlUAzNZaT4A+Js3Syp69vVjkiT5yy99ioroSUm/XLkyKTx85t27979+WIQQQp+EQfhV8vPzV917OXDAgHUrl3vquyZ93ceESi9fvnz7VluhmA5gDuBfUXFw/vzorx8WIYTQJ2EQNlxubm5EcJAPXTTFhqa4d2ZMl74KhTyfX/iVw+bn54vF9vWesCkpKf7KMRFCCP0ZvKC+4bZtWDfT06g91wgA7PTZBU/yFDLDpKJnltpfdTm8l5cXk7kKQKF8m0IQt3x9vdTTMUIIof8HZ4QNl/3urZM+W1U6aNE0K2hfeREFAHA4nMmTexgbR1CpOymUxWZmC2Niln3lmAghhP4Mzggbzrt5wMPES9a6msryUYW4t1XA0hNrGfHysLahQUFBDR558eIZffuG378fn5DgKxYvNDNjqqllhBBC/4to2Jeq/zNdu3Zty5YtV69e/eRSiURCkiSTqbZQqa2t7RDUqrUO6aBFiy8V1pk5JGelU3xYDAOm5Lmgo3fo9s0/f+UueDyws4PXr8HUVC0tf1N8Pl9bW7upu/heCAQCFotFEERTN/JdIElSKBSy2ey/XxWpQ2P/Z4KHRhtOS0vr7uOnLSYtKjew7j9oCNPEQCfS1CLc1rgFx2K047WHN9LTG3j3URU9PejTB3bvVku/CCGEPgGD8KswGIy+/frNGD8m0EI39XmqtrOBahHTQev58+dfv4upU2HHDpDi9zIhhFDjwCBUA7oZV1qUY8u1ExZ8vMWavFhsZ6eGe2e7u4ODA5w///UjIYQQ+gQMQjWgc7iyouzVi1eUHcmuSikVFvBLLuXYaFr6+PioZfzJkyEmRvru3Ts+Xw33MkUIIVQfBqEaUFjaBIPpamV258pvbWT+NTdKI33Cr56PVdeZC2Lx2fj4gBYtFtnbdxgxYsZ/6fwmhBBqcnj5hHrQzbnSomw71+Y/b4rZnrxPX0OPRlPPa5uTkzNzZrRc/qCykg0AZ87M9fPbO3nyaLUMjhBCCGeE6kE350qLcpSPbXWts6pz1TXyvXv3a2r6Anw4UZvPn3Ty5DV1DY4QQgiDUD1o5jbS4t+DUM8msypbXSNramrQaKJ6T4g0NTXUNThCCCEMQvWgm3NlRdnKx1xd6zx+gYJUqGXkkJAQHZ3TAO8BAECgo7Nw4sT+ahkZIYQQYBCqC93MRlqSBwoFAGjQmAYa+vn8IrWMbGhoePnybi+v8aamzXV02rm5devevataRkYIIQR4soy6EAwNqra+rKKIZmwBAHZ6Nlm8HGsdC7UM7uvrk5p6EwBqa8HZGZ49AzVdl4EQQghnhOpDN7eRFmUrH9vpcbOqc9S+Cy0tWLgQFixQ+8AIIfT9wiBUG5o5VxWEXD3rLJ7aThytb+xYyMmBa9ck7969q6ura4xdIITQdwWDUG3oZh9PHLXTtc7kqX9GCAA0GrRqtbdbtxaBgVHW1q0XLlzTGHtBCKHvBwah2tQ/cdRKx7JEUCqRS9S+l+Tk5AsXfpXJEktLT5SVPdm+/cWVK5/+2imEEEKf499xsoxCobh9+7bysZGRUbNmzZq2n0+imVrJKopJmZSg0WkUKkfbPLcm30FfDffdru/SpbjKylEADAAAoPB4E48fP9mlS7h694IQQt+Pf8eMUCqV9u/fPy4uLi4uLiUlpanb+TSCSqMZmslK85WlnZ5NZiN8TMhmaxBE/Y8G8b7GLAAAIABJREFURWw2Xl+PEEIN9++YEQKAsbHxmjX/9M/D6Oa20qJsOscWAGx1rbPVd6M1lZ49u65fP7i0tAOAGUCVoeGqUaP+6S8LQgj9k/07ZoQAkJeX17x58/Dw8MTExKbu5U/RzP54ozVettp3YWdnd/z4KienPiYmzen08BkzpgQENFf7XhBC6Pvxz5oRVlZWTps2rf4zurq6MTExdDr99evXVlZW8fHxvXv3Tk9PZ7FYTdXkX6Cbc4WPbygfN9KhUQAIDW2Xnn4fALZtgwcPGmMPCCH0HflHBGF+fj6TyTQ2NtbR0fnxxx/rL1J+mRGFQrGysgKAwMBAKyurzMxMDw+Ppun1L9HNbVTfQWHGNq0R1wikQja9sTJ7xAhYuhSyssDWtpH2gBBC/32fODQqk8ni4uLy8/PVuBsej9elSxczMzOCILKyslTPV1VVBQUFBQYGuru7jxw5kkKhuP2Rk5MTAJSUlNTU1ABAWlpaTk4Ol8tVY29qRDPiyGt5pFgEABSCsNG1aoyPCVW0tGD4cNi5s/H2gBBC/32fCMLS0tIOHTrUj6uvR6VSBw4ceOPGjf95ft26dQYGBtnZ2e/fv3/w4MG5c+c+uXlOTk779u39/PymTp169OhRLS0tNfamTgRRQNW+H3uBz+fD73ccbdQdTpkC+/aBQNCoO0EIof+yTxwaNTIyYrFYtbW1atyNtrb2oEGDxGLx/zx/9OjRXbt2EQShra09ZMiQY8eO9erV6/9vHhAQ8Pjx47/dS35+fkJCgq+vr+qZJUuWhIWFKR9LJBKSJKVS6Vf8HH+lrq5ueP8+0ry3lnGJU2bNnRW1hONrml72vp2ZOl/J/2FoCIGBGrt3y0ePbqyfq8EEAgFBEE3dxfdCKBQqFAp8wb8NkiRFIhFJkk3dyPfia/4zYbFYFMrfnBb6iSBkMBjjxo2Ljo4ODQ1lMpkN2/fnkMlkBQUFDg4OytLR0fH8+fNfMyCHw/Hw8Ni6davqGRcXFzb7w3e7K4Ow8X6i9SuXB1KqhoR7AECdTN5/9fJVF/YlVzxv7Plry5ax8+atWrtWbGSks3Pn8sDA1o26u89HkuQ/d+7+n0MQBIvFwiD8NkiSpFAoqv9bUGNr7P9MPn2yjK6ubkpKiqOjY8eOHY2MjOovUuPFfAKBQKFQaGpqKksWi6X8ILDBKBSKjo6On5+fOrr7Ynfibu5q/uG10qBRQy10yjJK3suyG3WnT58mb9iwTiI5X1xsVFyc27Nnn8ePT9nY2DTqThFC6L/k00F47NgxhUJRW1t79uzZ/1mkxiDU0dHR0NDg8XiWlpYAUFlZaWpqqq7Bvz1dXZ0asUyL8eElrZaRlvoceaH8fWGWPaexTuvcs+fXioooAGUAW5eXTzhz5tLMmZMbaXcIIfTf8+kgTE9P/wb7JgjC29s7MTFReS1EYmLiP/Mmop9p1KRpPy2eHR1kq82kJRZWPa2Slu3c+vDJvVDdRG1gH913xNvLW+07rakRAHw8PqNQaPF4xWrfC0II/Yd9u+sIz58/LxQKASA2NtbExKRnz55UKnXy5MlRUVFOTk7FxcXHjh2Lj4//Zv2o3Q89eohEwpHrVovKil29vD0C3ZLkaZ4/tQIAUbGg16Der5Nf0el09e60T58OsbF7eLyWAASAzNBwf2TkUvXuAiGE/tv+NAglEsmFCxdSUlLy8vLMzc09PDx69er1NfdzuXDhgkAg6NOnz927dwEgMjKSSqUOHjy4trZ20aJFLBbr1KlT/8zL5D9f/4GD+g8cVHksmmnv0WzsFM4cZ+XzmmZsuoXmq1evvL3VPCns3r3b3btPjx5tCeBWVZXWo8cof39/9e4Cof9r7z7Dmsi6OICfNHrvHWkCChYEu4JdUbFgX3Wtay+ra1uxshbW7tp11y72ggUVFbGtutiwgiId6QkEQnreD7ObN0KwBDID4fweP2TuTGb+QZ4cZubOvQhpNuWFMCcnp0ePHgkJCQwGw9zcvKioSCwWL1u27Nq1ax4eHqodaf/+/UrbJ0+ePHnyZNX2WTuxrB3FuRl0OkMmBRrjv1apjMFgfOltqtq0afny5T9nZGS8eOGyaZO+TAbYcxAhhL6d8qcrpkyZkp6efvTo0fLy8tzcXD6fHxUVJRKJRo4cSXK+uohl4yTKSe/fOyQvOgNkAAClqcWyPLG3t7eajmhsbOzj4zNihD6TCWfOqOkgCCGkmZScEfJ4vEuXLu3fv3/EiBFEC4PB6Nu3r76+fpcuXdLS0rB3/pcxrZ3Euemrlu8rnlNy+bcrPODbmdtcPBWlpjNCORoNwsNh+nTo3x+YtWIQWYQQqgOUfF+y2WyxWFy5AycxYkteXh4Wwi9jmttKuByGVLxr604A2BV/QEuLpfIl5e/SrRuYm38MDNzKYn0KCvKfP39a7ZymAyGEag8ll0YtLS11dXVv3LhRof3GjRs0Go2YBQJ9CY3GtLSTT1Xf3M43Ie8NOUdOSkpKTBz84EHvuLg1ERH6rVv3FovF5BwaIYTqKOVDrI0cOXLBggXFxcWDBw+2s7PLz8+/dOnSihUrgoODbWxsyE9Z57CsnUQ56SwHdwDwsfROLPoglkqYdPVeGgWApUu3FBVtBAgEAD5/anr6+9jY2G7duqn7uAghVHcpv5W0efPm/Pz85cuXL1++XN4YGBh44MABcmLVdUxrJ1HuvxMw6bP0bPStktkpnubu6j7ux4/pAJ7yRS7XKyVFvdNfIIRQXae8EOrp6Z07d+7Fixd3795ls9lGRkatW7du1aoVyeHqLpaNEy/+lnzR17JRQv4bEgph69ZNnz2LFYuHE4umpjf9/Oar+6AIIVSnKSmEubm53t7ekZGRPXr0qPEHwOsJxTNCAPCx9L6f+WiwV4i6jxsePi86uuenT2/LyjwZjMs9e1rj8/UIIfRlSjrL6OnpcTgcnGGkOpiWdhJ2nkz87xyBvpbeCflk9JcxNjZ+/TruwIEmERH5Xl4zhgz5g4SDIoRQnabkjNDQ0LBLly5RUVHt27cnP5BmoDGYDDMbcX4Wy7YBANgaWDNojE+lubYGap9eQ0tLa9CgQQDg5AQbNkCI2s9CEUKoblN+j3D27Nnjx4/ncDh9+/a1tbVVnO2Tqtn+6hyWjZM4J50ohADQ2MLzZf4bEgqh3ODBsHgxPHoEeG8XIYS+QHkhHD9+fG5u7t69e/fu3VthlUwmU38qTcCycRLlpuv+t+hr2ehl/tvuLp1IC8BgwPTpsGkTHD9O2jERQqjuUV4IT548KRQKSY6iYZjWTuUvH8gXfa28LyVfJznDhAmwahV8/AiuriQfGSGE6gwlhbCsrOzKlSvDhg2r09PkUo5l7ciN+X/HUXdTl7yy/BIh10jLkLQMhobQtu3pgIBtBga8du0Ctm5dYWFhQdrREUKoTlDSa5TL5UZERAgEAvLTaBKmlaO4IBukEmKRQWN4mru/KUgkM8OJE2fj4o4UFZ1JT3908mT3wMBQHHENIYQqUFIIraysbG1tk5OTyU+jSWgsLYaxhbgwR95C3CYkM0NExN6Skt0A5gA0iaRfXp5HQkICmQEQQqj2U1II6XT6hg0blixZ8uzZM/IDaRKmtZMo5/8jnPlaer/KI7UQlpQUA5jJF0Ui64KCAjIDIIRQ7ae8s8z+/fuLiopatGjh4OBgZWWluCo+Pp6UYJqAeIICfNsCgEAg+DNiz7Go4xd0IlsHtN6xabuZmdlX91BNHTu2Tk8/KxINBQAAvrZ2TIsWc9R9UIQQqluUF0JnZ2cTExOSo2gelrUjP+nfs+oZv8y8w37UfGUHoMGrR8mDfhh8K/qmugNs2bLi2bMBWVmX+XzbsrLYlSvnmpubq/ugCCFUtygvhJUfH0QqYNo4i+9GEa9jYmOcFvkQr81b2aQ+eFtSUmJkZKTWAIaGhk+fxrx48aKoqGjXroWFhaZqPRxCCNVFygshqhEsa0dRbgYxBkGFYQgYLAaPx1N3IQQAGo1GPAbj5ARt2sCUKWCK1RAhhBR81llm1qxZUVFR8sXdu3cr9h29du1ay5YtyYtW99G0dUvo2h9ePAEAT7eGxW8KiXZedqmWkEnyFMfu7tC3L2zZQuYxEUKoDvjsjPD8+fM2NjYhISEAIJVKJ0+eHBkZ6ebmRqwtLCz8559/KMhYN5WVlY0eHJrz9pl59JMsIW3Vxq1LVi/Lvp1UJivX5jIvn7hIfqRlyyAgQNKx4zMmk+fn52dgYEB+BoQQqm3w0qi6LJozqwPkDerXHAA+lfInzpr++NXb3NzcE6/O6VjqN2rUiPxILFa2UDioT5/GurrGWlqzDxxY16NHF/JjIIRQraLkOUJUIx7cuTPQw5J4bWug42um++rVK0dHx76tgh/nPKUk0qhRP3O5EeXle4uK1ufkxIwZMxcHmkEIISyE6sJgMqQKPWT4Eqmuri4ANDRzKxFwc8vyyI/07t1HqbTDf0vmUmmzxERSh3xDCKFaCAuhugwYMnT9kwyJTAYA/2Szs0QMT09PAKABLcC2+aNsCk4KWSwGgOKkInn4WCFCCFW8R7hv376YmBj54m+//bZnzx7idW5uLnm56r75i5euKOf3O3VCxmW7+DY/HnWKTv/3z47Wdi1upd0N8ehJcqSpU0f+9ts0LncTgC6dfsDTU5vknqsIIVQLfVYInZ2ds7Ky0tL+HR7T1dW1vLxcvki0kJquLqPT6StWr12xem3eptkm/X/ScnKSrwqwa77h8Q6RRMRisMiMNG/eND29fVu39iorE3I4QadPHyLz6AghVDt9Vgjv3LlDVQ4NxrSyF+Vlarn8v5uokZahs7Hjq4J3za19yUxCo9GmT584ffpEAAgJgXPnYNIkMo+PEEK1Ed4jVDumlYM4L7NCYys7v0fZTyjJQ5g/H37/HSQSCiMghFCtgIVQ7VjKC2ELagth+/ZgYwPnzlEYASGEagUshGrHtHIUVSqEnmYebD4ntyyfkkgE4qQQIYTqOSyEase0sJMU5YD0s6uQdBqthU2zfz5ROfVxSAjk5Fxu1erH7t1HnjhxhsIkCCFEISyEakdjaTEMTcVFnz1BL5PJih8VzB43s3u/HqfOnKIk2MaNOwoLDz5+PC8mJmzSpEuLF6+hJAZCCFELCyEZKveXWbBk4fnr580HOnACxQu3ha3dEEF+qg0b9vF4RwB8ALyKi/ft23dSKpWSHwMhhKiFhZAMTCsHcV6GfFEikUSePm4/0l3HUk/XWt9+rMeufbtIjiQUCiUSXQCt/xoYNJpNUVERyTEQQohyWAjJwLRyUOwvw+FwtI21gUYjFmkMOrBoQqGwinerhZaWloGBFEA+WhCHycy1sLAgMwNCCNUGWAjJUOEJCnNzc2mpRFwqIhYFbL4eS09LS6uKd6vLvn1rraz6Mpk76PTd2to9t28PJzkAQgjVBjgfIRkq3yPcsXH7T3Mm6wWYisUi4ZPSY/uOkp+qU6fAFy+iLl68zONJV6065ebmSH4GhBCiHBZCMjCMLaSCcimfR9fRI1p69uj56Mbf12JiNsXvuHUr2smWmiJkY2MzceJ4ABAIYPVqOHaMkhQIIUQlvDRKChqNaWEnzv/spNDW1nbM6NFtgtuXMMuoyiU3bRrExsLLl1TnQAgh0mEhJAnLykGcl1W53c20QTInhfw8Fejrw5w5sHo11TkQQoh0WAhJwrSyrzziKAC4mjRIZqeSHkeJkSOLz5+fam3t7+gY8Msvy0UiEdWJEEKIDFgISaJ0DgqoNWeEADBixDihsE1eXnxm5sMdO2izZi2hOhFCCJEBCyFJKjxKKOdu4vKRkyYDGfmRFHG53DdviqTSUQAAwCgvX3r+/A1qIyGEEDmwEJKEZeUozs8CWcWCZ6Clb6hl8Kk0V+m7SFNWVkajGSg00GQyuqxSWoQQ0jxYCElC09al6+hJigsqr3IzaZDMpvjqqI2Njb5+AcBHYpFOv+vh4UD7b+wbhBDSYPgcIXmI24QME8sK7W6mLsmc1A6ObShJJXf27O6QkKHl5T5lZTxt7U/Hjx+nNg9CCJEDCyF5mFYOorws7YbNK7S7mjjHpt2jJJKipk2bJCf//e7du9JS7eBgdzodTwcRQvUCXholT4U5KOTcTV2SOamkx1GCyWT6+Pi0bu0xYgRt2zaq0yCEECmwEJKHVcUTFA6GdoXl7HIxn/xIVZk7F3btguJiqnMghJD6YSEkT1WPEtJpdGdjhxROGvmRquLiAj17wu7dVOdACCH1w3uE5GGa2Ui4bJlISGNVnHHJzaTBB3ZKIwtPSoIptWABdOoUz+XesrExHjo0FKcqRAhpKjwjJIlIJFoW9mvw0XstvBtOnziupKREca2bSW25TSh3+vR6Ljf8t99cZs+W+vp2f/PmDdWJEEJILbAQkiRs3tzSuxejBgdE9fFyy3w6bsRQxbXuprVlxFECm83eseO0UHgOYLBYPCUn58j48QupDoUQQmqBhZAkMVcuzW7hyKDRAGCwp01OchKXy5WvdTVt8JGTSvlAa3Lv3r2TyVoq/Ho0ysj4RGUghBBSGyyEZJDJZDSZTPG5PAMtZmlpqXzRSMtQj6WXW5ZPfjalXFxc6PTXCg3ZZmbGlKVBCCF1wkJIBhqN5uDiGp/z7+MI6cXlRRK6ra2t4jZEfxkq0ilhY2PTrZuXkdHPAAkANw0NB23ZgpNRIIQ0E/YaJcmO/YeG9A22fs/WEvDeS7T/ijxVYQPiNmF7h1aUxKvsyJFthw8fP316C49nnJq6OzDQl+pECCGkFlgISWJvb383/tn7pKSMtZMDd15h6epV2MDVpMGdjL8pyaYUjUYbPXr46NHDAaBlSzh/HgYOpDoTQgipAV4aJQ+dTvf08vLxcIOSinNQSKXSp1f/ORC2d9iY4Q/+fkBJvC9YsgRWrqw8hRRCCGkCLIRkY1rYSgpzKjROmDpx7+X91n2dklyzhk/74XzUeUqyVaVvX9DSgqgoqnMghJAaYCEkG9PcVlzw2aMIHA7n5oNbtsNcdW0NDN1MHKZ4LV4ZRlW8qixaJJ01K3LEiBkLF4ZnZ2dTHQchhGoMFkKyMcxtxJ+fEaanp+vZGcoXmXosXjmP9FxfsWfPqMzM55GRI37/3dPPr8/79++pToQQQjUDCyHZmBa2ksLPzggbNmxYmlosE0uJxfKcMmtLayqiVenFixfx8WKJJAKgjUw2JDd3+5w5q6kOhRBCNQN7jZKNaW5b4YxQR0dn0ZyFv29Zr9vSWMwTCR9zL5++RFU8pT58+MDjKc4n7JeYiGeECCENgWeEZGOa24gLP1Xogjn1pynRRy9NbvyjvoHBrVu3fHx8qIqnlJeXl57eI/kijfbQx8ebwjwIIVSD8IyQbDRtXZqWroTLZhiZKbZ7e3t7e3tnxxZnS/LcwJWqeEo1bty4c2fr69encjj9abQ0M7MdmzbVrn6tCCGkMiyEFGBa2EgKP1UohITGlp6vC951cGxNfqovO3Fi16VLl69ejfv0yYrLvensrCQ8QgjVRVgIKcA0txMX5Gi5NK68ytey0Z8vjpAf6Vv06dO7T5/eAgG4ucHTp+DnR3UghBCqCXiPkAJMCxtxofJZjbzNPd6zU4QSIcmRvp22NsyeDevWUZ0DIYRqCBZCCjAqdRyV02HqOBs5JBV9JDnSd5kyBW7dAnySECGkGbAQUoBpbiOp4owQABpber0qeEtmnu+lrw+TJ8PKlblxcXEfP9bqmo0QQl+FhZAClUdZU+Rj6f06/x2ZeVQgFG46dmxgv34XW7eeERLyo1QqpToRQgipCAshBRjG5tLyUplQoHStj4XXy/xafUb44sWLvXuvS6V3i4vX5+dfjo11+eOPvVSHQgghFWEhpAKNxjCzFhcpv01orW/JpDOzS5WvrQ1u377PZg+R//KUlo64dOkOtZEQQkhlWAipwTSvOOKoIh9Lr1e1+KTQyspMS0txSsV8S0t8rBAhVFdhIaQG0+KLtwktvF/V4tuEPXv2MDc/ARAPAACZZmaLfv75R4ozIYSQqrAQUoNZaTImRbX8jNDU1DQ2NrJz53UODgFmZmP79l0REOBPdSiEEFIRjixDDYa5LT/xWVVr3U1ds0tzykQ8fZYemam+nYeHx82bJwAgORlatwYOB0xMqM6EEEIqwTNCanz5jJBJZ3iaub8pSCQzkmrc3KBnT9i5k+ocCCGkKiyE1GCa20qKcipMxqTI26xhzKObBQUFVW1QeyxdCps3A5dLdQ6EEFIJFkJq0LS06boGkpIipWuvxlwLHxG2Y+P2gO6t+wzqKxAof+KwlvDwgE6dYO3a1GvXrqWmplIdByGEvg8WQsowzG2VDr1dXFw8aeakBnN83Sb5Os1tlGKUHbYijPx430Uk+mXt2p+GDLnbsuWk0aNnUR0HIYS+AxZCyjAtbJQ+QfHkyRP9RiZMAxaxaNbRNjrmKrnRvs/ly1du3SqWSq+XlPyWn38tKkp0/vwFqkMhhNC3wkJIGaa5rURZfxk9PT0Q/v/eoVQk0dHRITHXd4uOvs/hDJYvFhcPvXjxLoV5EELou2AhpAzTXPkZYfPmzaUZQu5HDgBIRdLcU6kTRo8nPd13sLMzp9Hy5Ys0Wq69vQWFeRBC6LtgIaRMVfcItbW1r5y9bPyAlboqISH8wYTuYyZPnEx+vG83atRgK6vNAAkAAPDK0nLduHHDKM6EEELfDB+opwzTQvmlUQBwdXWNjb4JAOOvzO4T0I/cXN/N0dExJuav6dNXpKSklZQ4zpy5t0GDBlSHQgihb4VnhJQp4oujEpJPH4/kcDhVbeNv0/RJzgsyU6nG19c3Lu50evo/J0+ePXq0Gc5OiBCqQ7AQUiP21q0ebVumFZXE71wV5N/sn8ePlW7Wwqbp07pQCOW6dwcDA7iAnUYRQnUHFkJqzJ3y0/4urlNbNJje3GF3YIPZk5R3h2li1Tix6ANfzCc5XnUsWgTh4V8YMwchhGoXLIQUKCwsNNWimetqEYv2hjpiXqlQKKy8pQ5Tu6GZW0L+G3IDVkv//iASwY0bIBaLqc6CEEJfh4WQAiYmJkXl/39UUCyV8SUyLS0tpRu3sGn6NCeBtGzVR6NBixZ/9e3r17BhLze3NjExt6hOhBBCX4KFkAIMBqNnyICVD9PyeYJPpfxF91NGjBlX1cZ+daS/jNz16zeioqIEgvsFBQ8/frwwcuSvaWlpVIdCCKEqYSGkxqr1G1uNnfXrG96SO+/7zA6bv3hJVVs2smiYXZpTLCghM151/PnnOTb7VwBdAACwKioaf/nyNYozIYRQ1bAQUoNOp0+YNDn67sO/BrYaPmggjUaraksGjeFr6f0s9yWZ8apDLBYrPp8qlTJFIrxZiBCqvbAQUopGY1jYifOzvryVX526TThyZG8Tky0AEgAAKDU339+rV1eKMyGEUNWwEFKMZe0ozs388jYtbJrG5zwnJ0/1DRgQMn16Y0vLFpaWfen09nPnzmrYsCHVoRBCqEpYCCnGtHIQ532lELqaOJeL+TlleeREqr7w8PnZ2fH37m1YvfpZQkIo1XEQQuhLsBBSjGnlIMrL+PI2MqmM8UI0bur4tevWstlscoJVE5PJtLW1nT6ddvMmJCZSnQYhhKqGhZBirK+dEcpksh79ej5/+SLbqejQh1N+7fw/fVIyZ0XtpK8P06fD2rVU50AIoaphIaQY08pRXJANVQ9Tfffu3XRRln2om7GXuWVHe8PelstXryAzYTXNnAmXL0NKCtU5EEKoClgIKUZjadENTMTs3Ko2SExMpNn/f9AZfRfjl2/qzKMUAGBkBL6+f/j4+Flb+zdr1u3ly7oUHiFUH+B8hGBnZ0f9xcblh7+8PvGvz3qNfuG5Q9WEhoaePn26ZvdJOHjweHz8Pzze3zyedl5eUq9eI968uWVkZKSOYyGEkArwjBAKCgqEQqGsHouLi8vPz1fTj3fPnpMlJasAtAEAoCGXG3z//n01HQshhFSAhRCpl0Ag+K8KAgBIJDp8fl2aVQohpPGwECL1Gjy4p4HBtv+WOPr659u1a0dlIIQQ+hwWQqRe8+ZNGzCAY2nZ2sJiEJ3eefPmcCsrK6pDIYTQ/2FnGaRedDr90KGtZWVlBQUFK1c6vn+Pf3shhGoX/FZCZNDX13d2dv71V/r27cDlUp0GIYQUYCGsGUKhkOoIdYCbG3TpArt3U50DIYQUYCGslufPnw8dOszIyFRbW1tHR69Xr+Bbt25RHepfAoHg+XMlc1ZkZmbm5OSQn4cQFgYbN0J5OVXHRwihirAQqu7IkSMBAa3OnLnL5Y4F2CwQzIqJ+di1a9eVK1dWZ7fr1q1bs2ZN9eNlZmZ27ty5cntERMSePXuqv3/VNGoEzs7XfHwGNm7caf78cB6PR1UShBAiYGcZFb1582bs2HESSbBMdhRAn2iUSMIB5ixbtszPz69Pnz6q7TkjI0MkEtVc0oqWLFnCYDDUt/8vO3Xq/OvXe7jcPwCskpMPPXgw7O7dCzU+UA5CCH07PCNU0caNG2UyY5nssLwKAgAAE2Azg9F89eqI6h9CKBTOnz+/Y8eO7dq1W7hwYWlpKdHO4/GWLl0aFBTUoUOHjRs3Eo1RUVG9e/cODAwMDw9XrKPnzp0LDAzs1atXbGws0XLp0iX5a/KtWPEHl3sYwA3AUCCYlpSk9eHDB6rCIIQQ4Bmhym7ejJNIegMYVlpDl0iGPHoUJhQKtbS0lLzzm0mlUk9Pz59++olGoy1ZsmTBggXbt28HgMGDB+vp6e3atYtOpxNjWN+6dWvChAnHjh1zdHScMmVKbm7utm3bAKC0tPTEiRMHDhxISEgYMGBAfHy8u7v7s2fPLC0tqxOsOrhcLoCZfFEicf706ZOHhwdVeRBCCM8IK8rMBC0toNG+8i81tQjApop92EmlEm3toi/vgcWCL5+Y6ejojB07tqysLCUlpXPnzteuXQOAd+/e3b17d/87PU0FAAAgAElEQVT+/V5eXg0bNgwNDQWAPXv2zJgxo2vXrp6enlu3bt23b59EIgEAkUi0ceNGFxeXfv36hYaGHjp0qKZ/Wt+tSZNGNFrcf0sCBuNOkyZNqAyEEKr38IywIgcH+JZHIdzcLD5+zK5iZSadzigvN6veCSHk5uZ26dLFxcXFxcVFJBIR42KnpKQ4OTkZGBgobpmRkTFgwADitaenp0AgIPqF6unp2dnZEe0eHh5v376tVqCasHfv2nbt+hUWdigvt5ZIzq9cOcPExITqUAiheg3PCFXUtWsQg3EFoLjSGgmDcbxt23bVvC4KAEeOHGnRosXFixe3bt06YsQIopGYNKpCbxpra+usrCzidWZmJoPBIC5+8ng8DodDtGdlZdnYVHUKSx4bG5t37+6dOtXj8GFnH5+TNjajqU6EEKrvamMhlEqliYmJ586dU5wmMC8v7+DBg8Tlwdpg7ty5dHopjTYcoEShWQgwVSp9FRb2a/UPoaOjk5qayufzORxOeHg40ejr6+vu7v7rr78KBAKxWBwfHw8Aw4YN2759e1pamkAgCAsLCw0NJcownU5fvny5WCx+/fp1ZGTkkCFDqp+q+lgsVrdu3YYMGbJ8uf3KlSCTUR0IIVS/1cZCeO3atZ9//nnhwoVPnz4lWoqLi7t06VJeXn706NGwsDBq4xEaNmx49OhhFiuWwXADmAIQATCXyfSk0fatWbOmR48eKu9ZIBDo6uoCwNixYy0tLV1cXNq1azd8+PBmzZoBAJ1OP3/+fHp6upubm6ur67FjxwBg8ODB06ZN69atm4eHB4PB2LlzJwBoa2sHBgY6Ojq6u7uHhIT8/vvvLVq0AABHR0dbW9ua+SlUT79+IJHA1atU50AI1XNUTwpbpREjRly6dIl4vXPnznnz5slkMqFQ6ODgwOPxlL4lOjq6Z8+eVe1QIBDw+fzK7SwWS+WJed+8efPjjz9aWtoCgLGx2YABA+/fv6/arggikcjPz2/Hjh3V2cn3iouL69ixY43vtqSk5KvbHD8ua926xo9cH5WWlkqlUqpT1BdSqbS0tJTqFPXIt3yZVEdtPCOs7O3bt76+vgDAYrHs7e0zMjKoTvQvb2/vAwcO5OVly2QyDqfw7Nkzbdu2VXlvV69edXZ2ptFoP/zwQw2GrM0GD4bs7OsBAaMCA4fs23dIKpVSnQghVO9Q2Wv03Llzuz8fgLlv377Tpk2rvCWfz2exWMRrbW3t8lo2VKVUKi0uLjY0NGQyq/XzbNmy5d9//+3o6Fh/RlrZv/9wYeGp9PTVAPrPnm159Chh7971VIdCCNUvaiyEZWVlEonEyMhIsZHNZiclJTVo0MDa2rpXr14dOnRQXEvcG6vMzs4uO/vfZxWys7Pt7e3VlPl7PXz4cPWqVTdv3uCV85kMRrt2bX+ZN1/lwdXMzMzMzMy+vp0GCQ/fVlYWC6AHAFzu5gsXWm3ZwtPT06M6F0KoHlHLpdELFy54eHgYGBhUGPT54sWLHh4eCxcubNy48c6dO3V0dCw+p6+vDwACgeDjx4+lpaU5OTnEVdCBAwcePXq0rKzsxo0b9vb2FhYW6oj9vfbt29ehffsX926P97Vb26Xxzy1d8t8l9O3bd+HChdXZ7f79+5s0aaKlpWViYuLv7y8/aX737l23bt2MjY2NjIy8vb1Pnz5dEx+CYny+lKiCBBrNVf4XD0IIkUMthdDLyysyMvLPP/9UbJRIJNOnT9+zZ09sbOzNmzfnzZvHZrOVvj0rK2vhwoXa2trXrl0jxtL09fWdO3du//799+zZc+TIEXVk/l4JCQmTJ03q5WZ5a2Sb+W08fvBxmNnSLXpYq6n+LhEREWfPnlVttxEREWFhYWvWrOFyubm5uVu2bDl37hyxaujQoa1atSooKCguLj5x4kQt6flZTZaWRgAp/y2VAbxxcXGhMhBCqP6hydT2GNeRI0c2b95MPOgGAPfv3+/Xr19eXh6dTgeAgICAGTNmjB5dk89T7927d9asWdbW1vKWdevW9ezZk3hNdA3V1tau8C4zM7OysjL5PchvNHbs2MunT9wb3U6P9dlMDjKAvicf6zk1fPDw4ffmLy4utrOzi4yMDAkJqbBKLBZra2s/f/6c6DRUs+7cubN48eLo6Oia3W1paWmFEXAqS0hIGDhwOpv9g1hswGQe2rFj8tChA2o2Rj3B4/F0dXXrz91laslksvLycryGT5pv+TKpip6eHlF0voC8zjIZGRnOzs7yQA0aNEhPT6/ZQzg6OrZs2fKvv/4iFhkMhoODg3zKoaoKYWVSPg++1n3xTuytTk5mFaogANAAertZrn3wTzm78MvHojEYNO3P7on+888/QqGwV69elTdmMpnt2rUbN27c5MmT27Vr5+Xl9dVP8V3odLrKv2dVkclkX91n27Zt37y5cfnyFQ6nfNWqo61aNajpFPUFjUbT09PDQkgOmUxGp9OJWzmIBN/yZVId5BVCHo+nOOqYrq6uOiZl1dXVdXV1rc4eJJyC3HVTQfaVQliUm2Pho7zPjqWetlQqTVz+o5W+zpd3Yj52ibZHU/kim802MzOTn5uOGjWKGDJ0x44dHh4eUVFRmzZt2rFjx+TJkz08PA4dOuTv7/8dH6y2MjMzGzVqJACUlcHvv0MtGBgcIVS/kFcIbWxsioqK5IsFBQW183ucYWJht+rkVzezPu3xqbRM6apsbjmTwWi8/tz3Xm61sLAoKCiQz9+0YMECgUDQvn17LpcLACYmJitWrFixYgWHwyGuKr958+a79l/LTZ0Kbm7w4QO4u1MdBSFUn5D3QH2zZs1SU1OJUxyRSPTo0aOAgADSjl7junTrfiutiMOvOJW8WCq78D6vY8cO31sFAaBly5aGhoanTp0iFn18fFq0aFH56raJiclPP/2UnJysvvu7lDAygilTYMMGqnMghOoZtRTC1NTUiIiIS5cu5eTkREREEN/sdnZ2Q4YMGT169LVr18aPH+/p6dmmTRt1HJ0cc+fOFcpgcnQCW6EW8sWSeTdevy/khi1ZqsI+9fX1161bN3PmzL/++isjI4PNZt+8eZOYaEIikUyePPnRo0clJSWpqambNm3q2LGj5t0Qmj0bTpx4+8svv69cuS4pKYnqOAihekEtl0bFYjGbzW7QoEGDBg3YbHZpaSnRvmfPnvXr1+/cudPDw2Pz5s3qODRpXF1dT5w8NXzY0PYH7/VwsXA01ssrE8SkFhbyBJs3b+7UqZNqu504caKDg8OmTZvCwsL09PQ8PT1Pnz7dvHlzmUxmYmIyY8aM9PR0U1PT9u3b79ixo2Y/UW1w/foZPn/rhg3TAKR//DFuz54FAwb0pToUQkjDqfHxCfJdvXp1y5YtVT0GUFWvUS0tLRUenyB8/Phx/fr1V69czs3LMzM1DezUec6cOX5+firsikJ37txZsmRJXFzc1zf9Hlwu19DQ8LveYm/vl50dB0C8q9jBoUtGRnzNptJUZWVl2GuUNDKZjMfjYa9R0qjwZfJdcIb6anF1ddXIMzNKFBcXSySm/1VBADAWCvV5PBxxDSGkXnVj9olajs/np6enl5Up70SKvpGxsTGNVgQg+K+Bz2CUYBVECKkbFsJqiY2N7dS5k4GhgbOzs6GhYQv/FidOnKj+bokJFyu3C4VCPp9f/f3XWosWTTE1HQJwH+Curm7o0qUzqU6EENJ8WAhV98cff3Tp2iU+6ZnjQA/vaX4uw72TS9KGDRs2Y8YMlfd5+PBhZ2dnFxcXAwOD5s2bf/jwgWiPj4/v0KGDsbGxubm5n5/fzZs3ifbMzExizgpbW9tWrVrt37+/Bj4YdWbO/Ons2bk//HCyf/8zOjq//vDDWKoTIYQ0H94jVNGTJ09mz55t1c6h4YQmNMa/f0849HJLO5+0bdu2tm3bDh8+/Hv3mZSUNGnSpJs3b7Zp00Yqld67d4+4MJiUlNS5c+clS5bcvHmTTqcfPXo0JCTkxo0bxGZsNjs9PV1PT+/evXtDhgzx9vZu3bp1DX9aEgUFdQwK6ggAI0bArl0wbx7VgRBCmg7PCFW0efNmbRMdj3H/r4IAADRwHtDQ2MN8w0ZVHgt/9+6dpaUl8XglnU7v2LGjnZ0dAERERAQFBc2bN09LS4vJZP7444+jR49evny5/I0mJibm5ub9+vXz9vZ+8uRJNT9aLREWBps2gRqG4UMIoc/gGaESpcIyGXzlqZK4e3eMm1jQWUr+kjD3t356/GkB9yuDbgOAgZY+Df7f371ly5bFxcWhoaEDBgzo2LGjk5MT0f706dPJkycrvrF3796DBw+ucB/x3bt3Hz588Pb2/vJB64pGjaBNG9i3D2bijUKEkDphIawon1cwMfpnsVTy5c1yCnJsfJ2UrmIZa8tksuEnJmobf2nQbTqNvrLDgmbW/59WycbGJj4+fuvWrWvXrh09enS3bt2OHTtmbm6ek5NjaWmp+F4rKys+ny/vONOkSROxWJydnb1mzZoKkyHXaUuWQHBwKpN508xMv2fPniYmJlQnQghpICyEFVnqWZwPPfzVzTyXenIKlD8vwc/nMVnMKz+eUOEhfXd3961btwJAWlpanz59Vq9evWHDBmdn56ysLMXNMjMzjYyMdHX/ncXpzp07DAZj9+7dBw8enDlzpo7OV2a9qCtevDheWLh9xowRLFaBmVnn69cP+fj4UB0KIaRp8B6hinp068F5kS8qEVZol4mlhX9/CgoMUm2oGjlnZ+devXqlpqYCQNu2bY8ePaq49ujRo4rTFpqYmNjZ2S1fvtzMzGyDpoxaLRKJ5s9fKxRek0qnCAQLPn06PmrUHKpDIYQ0EBZCFf3yyy8MYLzdEi8slj8ADhK+OHHXs7Kc0mXLlqmwz9u3by9fvvzZs2f5+fmxsbGRkZHEdc5FixZ9+vRpzJgxr1+/fv/+/S+//HL79u3ffvutwttpNNqKFSvWrVvHZrOr89FqieTkZABfAPkD9Q1zcjThcyGEahu8NKoiJyen8+fODxo8KH5OrGkTCx0bfSGbz3lRIOVLdu/a1b59exX26eLicurUqYkTJ7LZbDs7u/nz50+dOhUALC0t//7775UrVw4aNEgqlbZs2fLhw4fu7u4AoKOj07VrVybz3//Hzp07DxkyJC4urn///jX4YSlhb28vk31UaCjT1cVfV4RQzcNvFtV17979zes3mzdvjr52NetepqWFZciwUbNnz27UqJFqO3R2dt6+fbvSVXZ2drt27arcbmVlFRMTo9iyZ88e1Y5e2xgaGvbo4Xf+/ILS0qkAXH39hcuWYf9RhFDNw0JYLQ4ODuvXr1+/fj3VQTTTwYNbdu3aHxk5t6xMt7h49ujRXalOhBDSQHiPsAZwOJzk5OS8vDyqg2gaOp0+der4u3dPP3162NS065UrVAdCCGkiLITVcuXKldat2pibm7u7u1tbWzfybnTgwAGqQ2mmhQshPJzqEAghTYSFUHURERF9+vRhpxVP7zjvt76bfumyVI9nNHbs2PHjx2vSdMe1xMCBUFoKt25RnQMhpHHwHqGKHj58uGjRohDfQXM6h9Fp//490a/J4MgnB3b8tSEwMHD06NHUJtQwdDrMnw8rV3J1dF6amppqzEhyCCHK4RmhirZu3WppaD0raKG8ChKGtxjTzMF/86bNqu321atXXbp0cXBwmDBhQkFBQb9+/WoirIbQ1b14715Qr14nAgOX+/v34HK5VCdCCGkCPCNUgl8k/NqY2/Dg7t8tndqxGFqVV7V37bT97npOdomO9leGOtM2ZdHo/x90m8/n9+7de/ny5cOHD+/fv/+IESPwvEeuuLh4+vTlEklcSYkBAHA4x6dPX3LwoIp/cCCEkBwWwooEHNGrXalfLYTsQraxrbHSVUa6xjKZ7PEfr8wMLL68E49h9sZu+vLF+Ph4kUg0ZswYGo02duzYYcOGqTZCjUZ68eKFSBQEYEAsikSD7t7dQmkihJCGwEJYkbYJy//Xhl/dzP6I3afiLKWrsosztVhancNbygd8+UZ5eXn29vY0Gg0AfHx8rK2tibkJEQCYmJgwmYpDrJUSsxYjhFA14T1CFfXs1fNh2t3CsvwK7UKxICbxsuKwZ9/O1tY2Ozub6HH64MGD0tLSmsmqERo3bmxpmUynRwMAANfQcOrPP4+lOBNCSCNgIVTR3LlzmVqshVEzckqy5Y0l/OKlV37JLfm0bLkqlzQDAgK0tbU3btz48ePHiIgIe3v7zZs3YzkkMBiMW7dO9u9/1sEhwNy8Z+PGPcaPH0l1KISQJsBLoyqyt7ePunghdGDoiIN9/B3bWBvZsnlF8ekPJDTpgYMHWrZsqcI+mUzm1atXlyxZEhkZOWXKlODg4DVr1mRkZGCXGYK1tfWZM3sBoKQE3Nzg/Xvw8KA6E0Ko7sNCqLrAwMA3b99s3779avS157mPTExNJk6ZOGPGDFdXV5X32bBhwxMnTsgXDx06VBNJNY2REUybBuvWgaYMMI4QohIWwmqxsrJasWLFihUrqA5S78ycCR4esHgxODtTHQUhVMfhPUJUJ5mZwcSJsH69JC0traysjOo4CKE6DAshqqu8vS/v2BEQEDDXxaXz2LE/S6VSqhMhhOokLISoTsrOzp4/P1wqvZ2ffzo//9GZM9qbN++mOhRCqE7CQojqpHv37nG5oQBGxCKXO/PEiWhqIyGE6igshKhO0tHRYTD4Cg18Xd2vjOyKEEJKYa9RAICnT5+qMBCMxkhKSqI6wnfr0KGDoeHK0tIRAG4AfBOTsMmTh1IdCiFUJ9Xfb3+5wMDAadOmUZ2iSmVlZYJitpE2g1gUSWRcGd3K2qaqjQuLi5iGLKCBTCKTlIjs7ewZDMZXjxIYGFiTodXP1NT08uV9P/44KTeXw2LJ5s6dOGxYKNWhEEJ1EhZCiImJoTrClxw8cCDjyIYxPvbE4qdSfngqRN24rXTj4NDehm3tdG3+ndEiLzp9Ub9FQwYPIScqyZo3b5aQcIPqFAihOg/vEdZ27Tt0iM4sE0n+fTbg5PvCbr37VrWxSCSiMf4/wSEwaUKBUN0JEUKoTsNCWNu5ubn9NG9RaHRS2OPMkTHv82y8ps6cXdXGg/sPZt/MISZTFPNEvHhOUFAQaVERQqguwkujdcCP4ycOHDLs3bt39vb2dnZ2X9hy4tgJr9++PvPbWT1LfWGBYOvvmx0cHEjLiRBCdREWwrrB0NAwICDgq5vRaLQt6zavW/W74hy/CCGEvqAeXRr98OFDXXxOQAVaWloODg7UVkGpVHr9+nUKA9Q3Dx8+LCwspDpFfcFmsx8+fEh1inrk+vXrEolEffuvR4Xw1KlTx48fpzpFfZGfnz916lSqU9Qjmzdvvn//PtUp6osHDx6sX7+e6hT1yMyZM3NyctS3/3pUCGUyGdUREFIj/A1HSDX1qBAihBBClWEhRAghVK9pVK/R3Nzc+Pj4bt26KV2bkpIilUofP35Mcqr6SSgUcjicqv4vUI178eJFRkbG9u3bqQ5SLxQWFmZkZOCvN2kKCwuHDx+ura2twnsHDBjw1f4KNE26r8DhcC5cuGBvb1/VWplMZmpqSnKqeislJcXFxYXqFPVFdna2ubm5at8U6HsJBIKCgoKqvmpQjUtJSWnQoIFqPeFdXFzc3Ny+vI1GFUKEEELoe+E9QoQQQvUaFkKEEEL1GhZChBBC9RoWQoQQQvWaRj0+8QUlJSVXr16VyWQ9e/Y0NjamOo4GKi4uvn//fkFBga+vb/PmzeWN//zzj3wbHx8fGxsbigJqlOTk5JSUFPliYGAgi8UiXsfHx798+dLHx+dbRmlH3yIvLy8hIUGxJSAgwNjYOCUlJTk5Wd7YoUMH7LVbHRwOJyEhwdra2tPTU94ok8lu376dmpratm1bxfa8vLyYmBh9ff2ePXvq6OhU99iyeiAnJ8fZ2blfv34DBw50dHTMysqiOpGmycrKMjQ07N69+5gxY2xtbSdNmkS0379/X19fv+t/bty4QW1OjREWFubk5CT/wZaUlBDta9ascXBwmDRpkpOTU3h4OLUhNcadO3fkP2p/f38ajZaSkiKTycLDwx0cHOSr8vPzqU5ah02ePFlLS8vY2HjGjBmK7aNGjWrcuPFPP/1kYWFx8uRJovH169cWFhYjR47s0qVL06ZNS0tLq3n0elEIFy9ePGjQIOL1iBEjFixYQG0ezVNWVpadnU28TktLo9PpiYmJMpns/v37Xl5elEbTTGFhYXPmzKnQyOFw9PX1X758KZPJ3r17p6enV1RUREU6TbZ69epOnToRr8PDw6dMmUJtHo2Rnp7O5/OnTZumWAifP39uYmJC/BqfPXvWzc1NIpHIZLIffvjhl19+kclkEomkXbt2O3furObR68U9wkuXLg0aNIh4HRoaeunSJWrzaB49PT1bW1vitYWFBZPJFAqFxKJAIIiJiXn48CGfz6cuoAbKzc2Njo5+9eqVvCU2NtbR0dHHxwcAPD093dzcbty4QV1ADSSTyfbv3z9u3Dh5S35+fnR0dEJCggwfyK4eR0fHyheWL1261LlzZ2IUlN69e2dlZb19+xYALl68SHyl0+n0gQMHVv8rvV4UwqysLPkYEPb29llZWdTm0WwREREtWrRo1KgRschisf74448JEyZ4eXlVuNGCVEan0z98+LBz586uXbt269atvLwcALKyshwcHOTb4K96jbtz505eXt7AgQOJRTqdnp6evnPnzl69egUFBZWWllIbT/Mo/kpraWlZWlpmZWWVlpaWlJTU7Fd6vSiEEolEPjYPg8EQi8XU5tFgkZGRf/755+HDh+l0OgC0atUqMTExKirq5cuX/fv3nz59OtUBNcTSpUsfPnwYFRWVnJycn5+/ZcsW+Pz3HACYTCb+qtesv/76a8SIEXp6esTiggULHj16FBUV9eHDBz6fv27dOmrjaR6lv9LEDL01+5VeLwqhra1tfn4+8To3N9fOzo7aPJrq7Nmzc+fOvXr1qnxkPwaDQbyg0WjDhg17/vw5dek0ivwHq6+vHxIS8uzZMwCwtbXNy8uTb4O/6jWrpKTk9OnTitdF5f8Lurq6/fv3J/4XUA1S/JWWSqX5+fl2dnbGxsa6uro1+5VeLwphUFDQtWvXiNfXr18PCgqiNI5mio6Onjp16sWLF4l7VJU9ffrU0dGR5FT1wdOnT52cnACgXbt27969y87OBoC8vLyXL1926NCB6nSaIzIy0s3Nzd/fX+la+f8CqkFBQUG3b98mOhwQXdCJey6dOnW6fv06sU2NfKXXi+cIZ82a1apVKxMTEwaDcfDgwb///pvqRJomLS2tf//+/v7+e/bsIVqmTp3atGnTJUuWFBQUuLm5ffz48fDhwwcOHKA0puYIDg5u1qyZiYnJnTt34uPjd+3aBQB2dnY//vhjSEjIqFGjjh079sMPP+BfHjXozz//nDBhgmJL//79vby8zMzMHjx4cO/ePcVHZtH3unz5clRU1P379wFg0qRJ/fr1Cw4ODgwMdHNz69+/f7du3bZt2zZ//nwtLS0AWLBgQUhIiEgk+vTp0+PHj/ft21fNo9eX2Sc+fPgQGRkpk8mGDx/u4eFBdRxNU1RUdPr0acWWHj16ODs7v3v37vr169nZ2VZWVr1791Z8HhZVx82bNx8+fMjj8ZydnYcOHSofI0IqlUZGRhIP1A8fPlx+7Q5Vk1AoPHDgwODBgxXncYuLi3vw4AGXy3Vycho6dChO8VYdT58+jY+Ply8GBAQQ43LweLwDBw5kZma2b98+ODhYvsHz58/PnTunq6s7evTo6l8arS+FECGEEFKqXtwjRAghhKqChRAhhFC9hoUQIYRQvYaFECGEUL2GhRAhhFC9hoUQIYRQvYaFECGNcuXKlbi4OPXt/9GjR48fPyZeP3v2bO/evUo34/F4p06dKikpUV8ShGoKFkKEaoZIJHJ1dV22bBmxmJaWtmfPHjabrb4j5ufn79mzp8LQ+8uXL//jjz/Ud8QePXrIh3m8cuXKlClTlG6pp6e3YcOG1atXqykJQjUICyFCNUMmk6WkpBQUFBCLL168mDRpEjHyp5qkpaVNmjTpzZs3io1TpkwZPny4mo4YHh7u7u7eu3fvb9l48eLFW7ZswamgUO2HhRAhskkkkpycHC6XW9UGIpEoJydHJBLJW9hsNo/H+5adjx07NjQ0tEIjl8vNycmRSqVVvau4uPir+y8tLT148OCYMWO+JQYA9OrVy8TEpKprpwjVHlgIEap5xJjXANC2bVszMzMzMzNiKnmJRLJy5UobGxtbW1tjY+O2bdsqTjHv4eERFha2ePFiMzMzW1vbmJiY8vLyIUOGmJqampmZ6evru7m5yQcuv3XrVpcuXQBg4MCBxCGI+eiDg4MnTpwo32d8fHzr1q2NjIxsbW1tbW03bdokX3XmzBkzM7O7d+8GBQWZmJgYGhq2b98+MzOzqg919uzZkpIS+bS0lT19+tTNzW348OF8Ph8AmExmSEjIoUOHVPwhIkQWLIQI1bygoKD58+cDwMaNG0+ePHny5MkGDRoAwIwZM9auXTt//vwXL17ExsYCQJcuXeRXUzkczt69e+/cuXPs2LG4uLjGjRsLhUIdHZ3Dhw+/fv3677//bt++/bhx44g3Nm3adM2aNQCwZMkS4hDEIMUFBQXyG5MpKSldunThcrnR0dFPnjzp16/fnDlz5LVQKBSy2ezRo0cHBwc/fvx4//79L1++nDt3blUfKjY21snJqaoBjmNiYjp16tShQ4dDhw7p6OgQjW3atElJSUlJSan+jxQhNZIhhGqCQCAAgKlTpxKLFy5cAIBXr17JN0hMTKTRaBs3bpS35OXl6evrr1u3jli0sLCwsLAoKSmp6hBSqdTHx2fChAnEIjHvz/Xr1xW3CQgICA0NJV7PmjWLyWR+/PhRvjYoKMjc3FwoFMpksmPHjgHAqlWr5GvnzxxTQJQAAAQXSURBVJ+vra0tkUiUHt3Pz69r166KLb/99huDwZDJZPv372exWAsWLJBKpYobPHjwAADOnTtX1SdCqDaoF/MRIlQbxMTEyGQyc3Nz4homwd7eXvHqaPfu3Q0NDRXfVV5efuLEicTExKKiIgAoKytLTk7+xiM+f/68ZcuWLi4u8pZhw4bdvn379evXzZo1I1p69eolX9uoUSOBQJCXl2djY1N5b/n5+e7u7pXbly5dumbNml27do0fP77CKnNzcwCQTzKOUO2EhRAhkuTm5gLA7NmzK7Qr9pqxtrZWXJWcnNyxY0eJRNK9e3cLCwstLS1tbe3i4uJvPGJ6erqfn59iC3Fhs7CwUN6iOIseMespMSF4ZUwmUyKRVGiUSqVbtmxp0aIFcU+0AqK/D4vF+sbACFECCyFCJCGmz3327Jmzs3NV29Dpn92237lzp1AoTExMNDMzI1oePnz4he6mFejr68tvQBKIRwDlE/l+F2tr6wp7IwLHxsb26NGjZ8+eFy9erHA6S1TcCtUdodoGO8sgpBYGBgYAUF5eLm8JDAwEgFOnTn37TlJSUlxcXORVMDc3V3EWb+IQRBdNpVq3bv348WPFK5MXL140MDDw8fH59gxyAQEBL1++lFWaytvPzy8uLi4pKSk4OLjCUDIJCQk0Gs3f31+FwyFEGiyECKmFp6cnk8ncsWPHvXv3njx5Ul5e7u/vHxoaunTpUuIxcx6P9/Lly9WrV1+5cqWqnTRr1uzZs2dnzpwRCASvXr0aNGgQjUaTr3V2djY0NNy3b9/du3efPHlS+ZLpzJkzJRLJsGHDEhMTCwsL16xZc+HChZkzZ8p7dX6Xrl27FhUVJSYmVl7VqFGj2NjY1NTUzp07K153ffDgQdOmTa2srFQ4HEKkwUKIkFrY29tv27bt9u3bnTt39vf3T0pKAoAjR45MmDBh4cKFDg4O+vr6TZo0OXLkSIXLiYrmzJnTuXPnQYMG6ejoNG/evHXr1sHBwfK1urq6f/75Z2JiYteuXf39/Ykumop8fX3Pnj2blJTk5eVlYWGxbNmyWbNmrVixQrVPFBwcbG9vf+LECaVrPT097969y2azO3bs+OnTJwAoKyu7ePHihAkTVDscQqShVb7QgRBSjUQiodPpiidtSvF4vMTERKlU6uDg8C33z9LS0vLy8jw8PExMTL66MTF8jOK9RolE8vbtWx6P5+XlZWRk9NU9fMHatWt3796dlJT0Lf1f/vzzz4ULFyYnJ1fzoAipGxZChNC34vP53t7ev/76q+LgNUqJRCIvL685c+ZMmzaNnGwIqQwLIULoOxQUFEgkkq+eyAqFwszMTCcnJyYTu6aj2g4LIUIIoXoNO8sghBCq17AQIoQQqtewECKEEKrXsBAihBCq1/4Hge4HbBn5NB4AAAAASUVORK5CYII=", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# | include: false\n", "\n", "errJ = zeros(length(J))\n", "for (i,x) ∈ enumerate(J)\n", " errJ[i] = norm( x - ξ, Inf ) \n", "end\n", "\n", "errGS = zeros(length(GS))\n", "for (i,y) ∈ enumerate(GS)\n", " errGS[i] = norm( y - ξ, Inf )\n", "end\n", "\n", "errSOR = zeros(length(SOR))\n", "for (i,z) ∈ enumerate(SOR)\n", " errSOR[i] = norm( z - ξ, Inf )\n", "end\n", "\n", "errA = zeros(length(alpha))\n", "for (i,w) ∈ enumerate(alpha)\n", " errA[i] = norm( w - ξ, Inf )\n", "end\n", "\n", "plot( errJ, label=\"Jacobi\", m=3,\n", " xlabel=\"Iteration (k)\", yaxis=(:log10, \"Error\"), \n", " legend=:bottomleft, color=:blue ) \n", "plot!( errGS, label=\"GS\", m=3) \n", "plot!( errSOR, label=\"SOR\", m=3) \n", "plot!( errA, label=\"α\", m=3) " ] }, { "cell_type": "code", "execution_count": 23, "id": "cb83f239", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mSaved animation to c:\\Users\\math5\\Math 5485\\Math5486\\tmp.gif\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "Plots.AnimatedGif(\"c:\\\\Users\\\\math5\\\\Math 5485\\\\Math5486\\\\tmp.gif\")" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mSaved animation to c:\\Users\\math5\\Math 5485\\Math5486\\tmp.gif\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "Plots.AnimatedGif(\"c:\\\\Users\\\\math5\\\\Math 5485\\\\Math5486\\\\tmp.gif\")" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "A = [5 0 ; 0 1]; b = [10; 3]; x0 = [1;4]; ξ = A\\b\n", "display( gif( SD( A, b, x0; α=.005, xlims=[ξ[1]-1,ξ[1]+1], ylims=[ξ[2]-1,ξ[2]+1]) , fps=1 ) )\n", "display( gif( SD( A, b, x0; α=.39, xlims=[ξ[1]-1,ξ[1]+1], ylims=[ξ[2]-1,ξ[2]+1]) , fps=1 ) )" ] }, { "cell_type": "markdown", "id": "8623cda3", "metadata": {}, "source": [ "* For what step size $\\alpha$ does the method converge for the example above? \n", "* Compare this with the method from lectures that used line search: $\\alpha^{(k)} = \\frac{(r^{(k)},r^{(k)})}{(r^{(k)},r^{(k)})_A}$:" ] }, { "cell_type": "code", "execution_count": 21, "id": "8ce8745a", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mSaved animation to c:\\Users\\math5\\Math 5485\\Math5486\\tmp.gif\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "Plots.AnimatedGif(\"c:\\\\Users\\\\math5\\\\Math 5485\\\\Math5486\\\\tmp.gif\")" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gif( SD( A, b, x0; xlims=[ξ[1]-1,ξ[1]+1], ylims=[ξ[2]-1,ξ[2]+1]) , fps=1 )" ] }, { "cell_type": "markdown", "id": "f31b0c56", "metadata": {}, "source": [ ":::\n", "\n", "::: {#exr-4}\n", "\n", "Recall that $(x,y)_A := x^TA y$ and $\\|x\\|_A := \\sqrt{(x,x)_A}$\n", "\n", "* Write down the definition of a vector norm, \n", "* Show that, if $A$ is SPD, then $\\|\\cdot\\|_A$ defines a vector norm.\n", "\n", ":::\n", "\n", "::: {#exr-3}\n", "\n", "Here, we prove Theorem 11.1 from lectures. Suppose $A$ is SPD, $x$ is the solution to $Ax=b$ and let $x^{(k+1)} = x^{(k)} + \\alpha^{(k)} r^{(k)}$ be the iteration given by gradient descent with line search (i.e. $r^{(k)} = b - Ax^{(k)}$ is the residual and $\\alpha^{(k)} = (r^{(k)},r^{(k)})/(r^{(k)},r^{(k)})_A$ is the step size). \n", "\n", "Let $e^{(k)} := x^{(k)} - x$ be the error between the iteration and the exact solution, let $\\|x\\| := \\sqrt{(x,x)}$ be the 2-norm, and define $\\lambda_{\\rm min}, \\lambda_{\\rm max}$ be the minimum and maximum eigenvalues of $A$.\n", "\n", "* Expand $( e^{(k+1)}, e^{(k+1)} )_A = (e^{(k)} + \\alpha^{(k)} r^{(k)}, e^{(k)} + \\alpha^{(k)} r^{(k)})$ and thus show\n", "\n", "\\begin{align}\n", " \\| e^{(k+1)} \\|_A^2 \n", " %\n", " &\\leq \\| e^{(k)} \\|_A^2 \\left( 1 - \n", " \\frac\n", " {\\| r^{(k)} \\|^4}\n", " {\\|r^{(k)}\\|_A^2 \\,\\, \\|r^{(k)}\\|_{A^{-1}}^2} \n", " \\right)\n", "\\end{align}\n", "\n", "Hint: $e^{(k)} = x^{(k)} - x = A^{-1}( Ax^{(k)} - Ax ) = - A^{-1}r^{(k)}$.\n", "\n", "* Use the so-called Kantorovich inequality\n", "\n", "\\begin{align}\n", " \\frac\n", " {\\|x\\|^4}\n", " {\\|x\\|^2_A \\,\\, \\|x\\|^2_{A^{-1}}}\n", " \\geq\n", " \\frac\n", " {4 \\lambda_{\\rm min} \\lambda_{\\rm max}}\n", " {(\\lambda_{\\rm min} + \\lambda_{\\rm max})^2}\n", "\\end{align}\n", "\n", "to prove that \n", "\n", "\\begin{align}\n", " \\| e^{(k+1)} \\|_A \n", " %\n", " &\\leq \\left( \n", " \\frac\n", " {\\frac{\\lambda_{\\rm max}}{\\lambda_{\\rm min}} - 1}\n", " {\\frac{\\lambda_{\\rm max}}{\\lambda_{\\rm min}} + 1}\n", " \\right)\n", " %\n", " \\| e^{(k)} \\|_A \n", "\\end{align}\n", "\n", "Note: this proves the theorem because $\\kappa(A) = \\frac{\\lambda_{\\rm max}}{\\lambda_{\\rm min}}$.\n", "\n", ":::\n", "\n", "\n", "We saw in class that \n", "\n", "\\begin{align}\n", " \\|e^{(k+1)}\\|_A^2 = \\|e^{(k)}\\|_A^2 - \\frac{\\|r^{(k)}\\|^4}{\\|r^{(k)}\\|_A^2}\n", "\\end{align}\n", "\n", "Since $e^{(k)} = x^{(k)} - x = A^{-1} (Ax^{(k)} - Ax) = - A^{-1} r^{(k)}$, we have \n", "\n", "\\begin{align}\n", " \\| e^{(k)} \\|_A^2 \n", " %\n", " &= (e^{(k)},e^{(k)})_A\n", " %\n", " = (A^{-1} r^{(k)}, A^{-1} r^{(k)})_A \\nonumber\\\\\n", " %\n", " &= (r^{(k)}, r^{(k)})_{A^{-1}}\n", " %\n", " = \\|r^{(k)}\\|_{A^{-1}}^2.\n", "\\end{align}\n", "\n", "As a result, we have \n", "\n", "\\begin{align}\n", " \\|e^{(k+1)}\\|_A^2 = \\|e^{(k)}\\|_A^2 \\left( 1 - \\frac{\\|r^{(k)}\\|^4}{\\|r^{(k)}\\|_{A^{-1}}^4 \\|r^{(k)}\\|_A^2} \\right)\n", "\\end{align}\n", "\n", "\n", "::: {#exr-6}\n", "\n", "Here, we prove the Kantorovich inequality: let $0 <\\lambda_1 \\leq\\cdots \\leq\\lambda_n$ be the eigenvalues of $A$.\n", "\n", "* Show that Kantorovich is equivalent to \n", "\n", "\\begin{align}\n", " (\\sum_i \\lambda_i y_i^2) (\\sum_j \\lambda_j^{-1} y_j^2)\n", " %\n", " \\leq \\frac{(\\lambda_1 + \\lambda_n)^2}{4\\lambda_1\\lambda_n}\n", "\\end{align}\n", "\n", "for all $\\|y\\|_2 = 1$.\n", "\n", "* Explain why $\\lambda^{-1} \\leq \\frac{1}{\\lambda_1} + \\frac{1}{\\lambda_n} - \\frac{\\lambda}{\\lambda_1\\lambda_n}$ for all $\\lambda \\in [\\lambda_1,\\lambda_n]$, \n", "* Thus show that \n", "\n", "\\begin{align}\n", " \\sum_j \\lambda_j^{-1} y_j^2\n", " %\n", " \\leq \n", " \\frac\n", " {\\lambda_1 + \\lambda_n - \\sum_j \\lambda_j y_j^2}\n", " {\\lambda_1\\lambda_n}.\n", "\\end{align}\n", "\n", "* Conclude by finding the maximum value of \n", "\n", "\\begin{align}\n", " X \n", " \\frac\n", " {\\lambda_1 + \\lambda_n - X}\n", " {\\lambda_1\\lambda_n}\n", " %\n", "\\end{align}\n", "\n", "for $X \\in [\\lambda_1,\\lambda_n]$.\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "5f7043fc", "metadata": {}, "source": [ "::: {#exr-8}\n", "# Bonus\n", "\n", "Adapt this code to include a preconditioner $P$:" ] }, { "cell_type": "code", "execution_count": null, "id": "e22ee989", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "CG (generic function with 1 method)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function CG(A,b,x0; P=I, Niter=50, tol=1e-3)\n", "\n", " X = [x0]\n", "\n", " g(x,y) = dot([x;y], A*[x;y])/2 - dot([x;y], b)\n", "\n", " r = b - A*x0\n", " v = r\n", "\n", " for i ∈ 1:Niter\n", "\n", " if (norm(r)