{
"cells": [
{
"cell_type": "markdown",
"id": "3d0a47e3-7c72-43fc-8563-d44cf886a24b",
"metadata": {},
"source": [
"---\n",
"title: \"A4: Sparse Arrays\"\n",
"subtitle: \"Assignment 4: Discrete Laplacian and Sparse Arrays\"\n",
"date: 2026-03-23\n",
"format:\n",
" html:\n",
" other-links:\n",
" - text: This notebook\n",
" href: A4.ipynb\n",
"---\n",
"\n",
":::{.callout-note}\n",
"\n",
"* Complete the following and submit to Canvas before April 1, 23:59\n",
"* Late work will recieve 0%,\n",
"* Each assignment is worth the same, \n",
"* Please get in contact with the TA/instructor in plenty of time if you need help,\n",
"* Before submitting your work, make sure to check everything runs as expected. Click **Kernel > Restart Kernel and Run All Cells**.\n",
"* Feel free to add more cells to experiment or test your answers,\n",
"* I encourage you to discuss the course material and assignment questions with your classmates. However, unless otherwise explicitly stated on the assignment, you must complete and write up your solutions on your own,\n",
"* The use of GenAI is prohibited as outlined in the course syllabus. If I suspect you of cheating, you may be asked to complete a written or oral exam on the content of this assignment.\n",
"\n",
":::"
]
},
{
"cell_type": "markdown",
"id": "7f848e4a",
"metadata": {},
"source": [
"::: {.box}\n",
"\n",
"***Name:***\n",
"\n",
"***Approx. time spent on this assignment:***\n",
"\n",
":::\n"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "454c9f57",
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra, Plots, LaTeXStrings"
]
},
{
"cell_type": "markdown",
"id": "ad62cb7b",
"metadata": {},
"source": [
"Suppose we have a uniform grid $\\{x_{j}\\}$ on $[a,b]$ so that $x_{j} = a + \\frac{b-a}{n} j$ for $j=0,\\dots,n$. We will use $h := \\frac{b-a}{n}$ to denote the grid spacing.\n",
"\n",
"**Exercise 1.** Write down the Taylor expansion of $u(x_j \\pm h)$ about $x_j$ and hence show that \n",
"\n",
"\\begin{align}\n",
" u''(x_j) = \\frac{u(x_j - h) - 2 u(x_j) + u(x_j+h)}{h^2} + O(h^2),\n",
" \\nonumber\n",
"\\end{align}\n",
"\n",
"as $h\\to0$.\n",
"\n",
"::: {.box}\n",
"\n",
"\n",
"***Your answer here***\n",
"\n",
"\n",
"\n",
"\n",
":::\n",
"\n",
"**Exercise 2.** Define $u_j := u(x_j)$. Explain why we can approximate the second derivative of $u$ at the interior points $\\{x_j\\}_{j=1}^{n-1}$ by the vector\n",
"\n",
"\\begin{align}\n",
" L\n",
" \\begin{pmatrix}\n",
" u_1 \\\\ u_2 \\\\ \\vdots \\\\ u_{n-2}\\\\ u_{n-1}\n",
" \\end{pmatrix}\n",
" %\n",
" \\quad \\text{where}\\quad \n",
" %\n",
" L := \\frac{1}{h^2} \n",
" \\begin{pmatrix}\n",
" -2 & 1 & & & \\\\\n",
" 1 & -2 & 1 & & \\\\\n",
" & 1 & -2 & \\ddots \\\\\n",
" & & \\ddots & \\ddots & 1\\\\\n",
" & & & 1 & -2 \n",
" \\end{pmatrix}.\\nonumber\n",
"\\end{align}\n",
"\n",
"::: {.box}\n",
"\n",
"\n",
"***Your answer here***\n",
"\n",
"\n",
"\n",
"\n",
":::\n",
"\n",
"\n",
"The Laplacian is $\\Delta u := \\sum_{i} \\frac{\\partial^2 u}{\\partial x_i^2}$ and reduces to $u''$ in 1 spatial dimension. $L$ is a discrete version of this and is known as a discrete Laplacian (with homogeneous Dirichlet boundary conditions). Suppose we want to solve $-u''(x) = f(x)$ on $[a,b]$ with the boundary conditions $u(a) = 0$ and $u(b) = 0$ (this e.g. models the deformation of a beam under external forces $f$, ....). By defining $f_j := f(x_j)$, we may approximate the solution of $-u''(x) = f(x)$ on the grid $\\{x_j\\}$ by solving the linear system of equations $- L u = f$. We did this last semester (in [Assignment 8](https://jack.thomaslabs.co.uk/teaching/Math5485/A8%20Solutions.html)). Since $L$ is tridiagonal, one can compute the LU decomposition of $L$ using only $O(n)$ operations (we had a presentation on this - see Thomas' Algorithm for a recap). Here is an implementation (for large $n$): "
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "5e1b70fe",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Solving the linear system of equations took 11.1537272 seconds\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdZ0AU18IG4DMzS+8dAcUuYlcExAr2gsYWNcaYGFNMM9ck12sSjYnpxhg/S6ImNyaxa+wliqCIChak2EEEld57WWbmfD/m3s1emgjb931+7Q6HmbMjzrunzTCUUgIAAGCsWG1XAAAAQJsQhAAAYNS4lStXarsOAKr3+eefX758efDgwardbUVFxeeff37//v0BAwaods/1PXr0aM2aNcXFxd27d1fJDrOyslavXp2Xl9ejRw+V7FC1srKyjhw5cuTIkRMnTlhaWrZt21Zjhz548KCzs7OVlZX0tqCg4OTJk76+vhqrAGgZBTBEFhYW7u7uKt9tTk4OISQwMFDle64vKiqKEPLCCy887S+mpaXt3bv3xo0bdbbHxsYSQmbMmKGiCqrS2bNnFTlECPm///s/jR06ISGBELJp0ybFlsWLF5uYmJSVlWmsDpTSwsLCdevWff/993l5eZo8LlBK0TUKhqlLly6dOnXSdi20Iyoq6tlnn92zZ0+d7ebm5l26dPHw8NBKrZq2fPnyioqK7du35+fnFxYWvvrqqxo79ObNm728vObNmye9raqq+uOPPxYvXmxtba2xOpSXlw8fPjwoKKioqGju3LkaOy5IZNquAIBaSF/zQZmvr29SUpK2a9Gw27dvt2nTRvMZUFFRsWPHjl9++UURe3v27DEzM1u+fLkmq7Fly5ba2tq+ffuOGDHimWee0eShgSAIQVlpaakgCPb29gzDJCYmJicnm5ubDxkyxM7OTrlYWVkZz/NSMcVGQRBKS0tNTEwUF5Sqqqrq6morKytTU9NHjx7FxcVxHBcUFOTo6KgocOHCheLi4m7duvXu3bvBKpWUlMTExJSWltrb2wcFBSn3nikOYW1tbWJikp6eHhcXV1tbGxoaamJikpOTw7Ksi4tLnR1WVFRcuXIlPz/fxsbGx8enffv2yj8VRfHWrVuPHj2qrKx0c3Pz9/c3NzdvyakkpLa2Nj4+PiMjgxDi5OTUr1+/+i0MxadzcHAYNGhQnU9XX0VFhVwut7W15ThOsZFSWlxczHGcra0tIaS8vLyiokI6OUVFRVIZ6V+B5/n8/Hxzc3N7e/s6e5b+uVmW7d69u4+PT52fNvMPowkFBQWXL1+uqKhwcXEJDAxUPqulpaVyubywsLBDhw5ShRmGqV9DNdm1a5e/v//06dMVWzZv3rxu3TrpZGrMoUOHgoKCZDLZw4cPNfbZ4W/a7psFHdKnTx9CSFJS0siRIxV/Iba2tvv371cu5u/vTwjJz89X3njnzh1CyMSJExVbPvzwQ0LIH3/88eabb7Lsfzrhra2t9+zZQyk9cOCAk5OT4ijPPfdcbW2t8g6rq6vfeecdU1NT5ZooD+RQSv/xj38QQvbu3fvKK68oDiENsdQfI6ypqfnnP/9pYWGh/Pfv7++vKLB69eo6wenk5PTrr78q76SZY4T79u1r166d8q5MTEyUKy8IwvLly5XzwMrKavXq1co7qT9GOHPmTELItWvXlItVVVURQnr06CG9HTNmTP3/5jt27KCNjBHGx8f37dtXufDgwYNTUlKUyzTzD6NBcrl88eLFMtnf37kdHBy2bt2qKODn51ento6Ojk/craoEBQXdvXtX8TYhIWH06NEaO7qE53krK6uNGzdq+LiggBYh1DVz5kwTE5MtW7a4u7uHhYVt2LDhhRdeCAgI8PLyasHevvnmm/z8/DVr1nTq1Ck6Onr16tUvvPACIWTevHmvv/76iBEjCgsLP/300507dwYGBr799tvSb4miOG3atBMnTgwcOPDNN9/09va+c+fOF1988cYbb8hksldeeUX5EJ9++mlhYeHy5ct9fHwyMzPNzMzqV4NS+uyzzx4+fLhr165Llizp3r17QUFBfHz87t27FWWuXbs2ZMiQcePGeXt7U0qvXLmydu3aBQsWeHh4NJgujbl3796cOXNsbGzWrl3bv39/hmHS0tLCwsKkhprkn//855o1a7y9vVesWNG5c+fbt2+vXLnygw8+kMvl0heIFvvoo4+6deu2fv36adOmzZ49W9oYEBDQYOG0tLTg4OCioqKFCxfOmjWL5/lff/117969w4YNi4+Pd3Z2Vi6s/Idx+vTpjRs3NucPY+HChb///ruPj89HH33Url272NjYzz777JVXXqGUSv+OX331VUFBwezZs93d3f/v//6PEKL87Uet4uPjR40a1a1bN8WWf//73+vXr9fM0RVu3bpVUVEhfdsA7dBuDoNOkf4rDh48WC6XKzZK0xZ++OEHxZanahHa2Ng8fvxYsXHx4sWEEIZhlNsE0dHR5H8bZ7/99hshZNSoUYIgKDY+fvzYxsbG2dm5srJS2iK1CC0sLB48eFDns9RpEe7cuZMQ0qNHj+LiYuViyvuvrq6us5NLly4xDDNy5EjFlua0CNetW0cI2bBhQ2MF7t27x7KspaXlw4cPFRtv3LjBcZypqWlWVpa0pWUtQkrpH3/8QQj5+OOP6xy3fotwzpw5hJC3335buZh0lHfffVexpcE/DCnGlP8w6rt06RIhxNnZWXkm5Pnz5wkhdnZ2paWl0hae5wkhnTt3bmJXLVBTU/Ppp5+uWLFi0qRJyv0NixYtOn78OKU0Ozu7zj+68r+IBly+fHnmzJn+/v4Mw0yZMmXmzJl1/nFBMzBrFOpavny5iYmJ4u2UKVMIISkpKS3b24IFC5RbDOPGjSOEeHh4LFiwQLExMDDQzs5O+RC//PILIeSrr75SdHgSQry8vJ577rn8/PyLFy8qH2L+/PkdOnRouhrbtm0jhHz++ed1hrWU91+/KTlo0KAOHTpcunSJPs2dCKXe10ePHjVWYN++faIovvTSS8rdpz179pwxY4ZcLj9w4EDzj9Uacrn84MGDpqamddqgn3zyCSFEua0sqfOHIc3paPoPQ5q5+tZbbyk3LocOHRoSElJSUvLXX3+1+kM05fvvv1+wYMHixYuPHz9++fJlaWNeXt5PP/1UXl5OCHFzc6vzj16nQ7tBycnJZ5qnrKys6V35+/vv3bs3KCioS5cuhw4d2rt3rwbWp0J96BqFuuqstm7Tpg0hJDMzUyV7c3V1JYR069ZNOYEIIW5ubklJSTU1NdKF6erVqyzLXrly5fr168rFcnNzCSH3798fNWqUYmOdIa4GSY0hqS3bGErp3r179+/f//Dhw4yMDKmlJc0MKisra/7sifHjx9vY2Hz77bdRUVGTJ08eMWKEv7+/8udNTEwkhAQGBtb5xaCgoD179kg/1YDk5OTq6upu3bq5u7srb+/Ro4ednV12dnZOTo6bm5vyduVizfnDaOKTRkREJCYmSq1PdSgvLy8qKvLy8jpw4ADDMIr+z8jISErp0KFDW7xnGxsbBweH5pRsZh9vXFxcc/6GQX0QhFCXjY2N8ltpgqIoii3bW52pktLe6hyizlEqKyulEHrzzTcb3KfyYBv5b7g2rbi4WCaTSdfuxsydO3fXrl02NjbDhw/v06ePg4MDy7K///57VlZWbW3tEw+h4OXldfbs2ffffz8qKkrq9XVxcVm0aNGyZcuk2TGlpaXkv0GiTNpSXFzc/GO1RklJSYPVIIR4eHiUlJSUlJQoB2EL/jAa+6TSWsan+qQZGRnHjx9vTkk/P7/+/fvL5fJFixYRQvbt2xcYGKhokkZGRnbr1q3pv4Smubu71/nq0BqU0sTExA8++EBVO4QWQBDCU5NWTdTpLayurlbV/s3NzTmOMzExkdZjNLM+TbO2ti4pKSkoKKgzAUQhMjJy165dfn5+ERERylf8ffv2Nb/mCgMGDDh79mxubu7Zs2fDw8P37Nnz2WefPXz4UOqhtbS0JIQUFBTU+a38/HxS76tDfao689JqDemgDdbkics5nqiVn1RZ83unpZKOjo6Ojo4VFRVHjx794osvFD89d+6c8txXrUtLSysqKurXr5+2K2LUEITw1KQ1Bnl5ecqhosKV2tJqtps3byYmJqpqyKRXr14XLlxITEwMCQlpsEBMTAwhZN68ecopWFpa+vDhwxYf1NXVddasWbNmzfrkk0+6deu2a9eurVu3mpiYSAv1bty48eyzzyqXl24C0MSdRaW2b15envLG+mdeWqsgCELT1evcubNMJnvw4EF5eblyJmVmZubl5dnb27em2STx8fGJioq6ceNGcHCw8nbpk9ZfsNgELy+vFtxuJiIioqKiQrFEPT8//9atW9IgaIvl5OSkp6c/sRjLsr6+vg3OYVYWHx9Pmte9D+qDyTLw1KRF6OfOnVNsoZSqdtK5dIeRL774on47oGWdtNJCgq+++qqxeJCGABUr0CXffPPNE+Okvvo19PT0dHJyksvllZWVhJDJkycTQrZs2SJN2ZDk5ORs376dYZjQ0NDG9lz/zBNCpCUHyqQAy87ObrqeVlZWI0eOrKys3Lx5s/L2NWvWEEJCQ0PrjOO2gDTTav369cp9y/fv3z98+LCZmZk0c0qtbt686eTk5O3tLb09f/68NEBYU1OzadOmlu2zrKysqBkKCgpqamqeuLf4+HhXV1cV9rVCS2hptiroImmWfJ0FBtKX96lTpyq2hIWFEUKcnJx2796dlpYWGRk5fvx46eJbf/nEzp07lfcmff+dMmVKnUNLzSDFuoiKigrpO/KECRPCwsLS09PT09OjoqI+/fRTb29vxVR4afnEwYMH63+WOssnqqurpcbl2LFjo6KiCgoKUlJS/vzzT8Xnio+PZxjG1tZ23759xcXFjx8/XrlypZmZmXQfHMVakeYsn/jggw9mzpx55MiR5OTksrKy+/fvL1myhBAyZMgQRZlp06YRQoYOHXrlypXCwsLz589Ln3fhwoWKMvWXT9y9e5fjOEtLy59++unBgwcxMTEvvPCCdOaVl0/k5OTIZDJbW9vvvvtu7969e/fulVYF1F8+cfXqVRMTEzMzs/Xr12dmZj5+/Pizzz5jWdbKyio5OVlRrJl/GPUJgjB8+HDp3zE+Pr6wsPDUqVNdunQhhCxbtkxRTE3LJyilX3/9da9evRRvp0+f3q5dO0rpiRMnrly5ovLDtcDkyZPHjRun7VoYOwQh/K351ztpOaBCz549T58+rcIgpJTm5+dL7Yk6fH19eZ6XyjQ/CKUd1m+CuLq6KgqsXLlSebjR3Nz8t99+q7NosjlB2OBtKvv165eWlqYoU1ZWNnXqVOUCDMPMnz9feVlbg0+fWLNmjXJDrW3btlKnrnIQUko3bNig3NvZxJ1lDh06pHyLH0KIh4fHuXPnlMu0OAgppQUFBXXG5FiWXbx4seIfkaozCB88eODq6hoXF5efn79q1arVq1e3adMmOzt7yZIloiiq/HAt0K5du08++UTbtTB2DH2aBVJg2GJiYsrLy0eMGKF8Q6zy8vKYmBhXV9c6twONiYmJiIioqanp1atXaGgoz/PR0dHOzs6K0Y779++npaX17NlTudtH2puLi0ud+2hER0dXVFSEhITU6Y67c+fO2bNns7KyLC0tPT09AwMDu3btqvhpUlLSo0eP+vTpU/+eolFRUTKZbNCgQXW2x8XFRUZG5uTkuLq6dunSZdSoUcr3OUtISLh48WJWVlbbtm0nTpzo6el55cqV0tLS4cOHS9N2amtrL168aGdn1/TshoyMjKioqIyMjKKiIg8Pjx49egwdOrR+T+P169fPnTuXl5fn7u4+cuTInj17Kv+0pKQkLi7O3d29zljazZs3T548WVpa2rVr12eeecbKyioiIsLKyqrOh62pqXn48KE05VX6VygvL7927Zqrq2udJ+2VlJQcP348KSmJYZgePXpMmDBBmuSi8FR/GA26dOnSxYsXpbMxduxYqVGoQCkNDw+3tLQMCgp64q6eVnp6+uHDhzmOmzZtmqur64ULF65du/bcc881Z7KxuhUWFjo5OUVGRg4bNkzbdTFqCEIAAE3Lz893dnY+derUSy+99OjRI+VvGKB5mCwDAKBRP/30k4uLS3R09NmzZ1966SWkoNbhHwAAQKPy8vImTJjQpk2bmJiYo0ePars6gK5RAADNksvlv/32W1VV1fPPP694PCdoEYIQAACMGsYIAQDAqCEIAQDAqCEIAQDAqCEIAQDAqCEIAQDAqOl3EH7xxRfSXQqbo/klQSUopS14dAO0hiAImAeuYbiwaJ7Kz7l+B+GaNWvKysqaWVh66DlojCiKcrlc27UwLnK5vGWPqYIWw4VF81R+zvU7CAEAAFoJQQgAAEYNQQgAAEYNQQgAAEZNNUG4ceNGFxcXGxub6dOnNzh75dq1a7169bKwsOjTp8/169eljZTSpUuX2tvb29nZvfnmm4oZhnfu3PH397ewsOjatWtkZKRKaggAANAgFQRhQkLCxx9/fPbs2by8PLlc/sknn9QpIIri7NmzFy1aVFlZ+eqrr86aNUua2LZ///79+/ffvXs3LS3twoULP//8s1R+/vz5kyZNqqio+Oyzz2bOnFlTU9P6SgIAADRIBUG4bdu2adOm9ezZ09zcfOnSpb/99ludlUyRkZFlZWWvv/46wzCLFi0qLS29cOGC9Iuvvfaau7u7g4PD4sWLt23bRgi5efPmrVu33n//fZZlZ8+e7ejoeOzYsdZXEgAAoEEqCMKkpKSePXtKr3v27FlYWFhQUKBcIDk52dfXl2VZQgjLsj4+PsnJyXV+sUePHtLG5OTkjh07Wlpa1tneIEppcXFx0X9h+TYAADwtFTyhvri42NraWnptY2NDCCksLHR2dlYuYGVlpXhra2tbWFhICCkqKlL+xaKiIinYGizcoPLy8n79+jEMI70dN27cli1bGiwpF8mgkyYHRpR7W+G+GxoiCIJcLsd9NzSpqqrK1NSU4zhtV8SIlJeXa7sKxiW7ivk8gWwIbO6tVMzNzU1MTJouo4IgdHJyKikpkV4XFxcTQlxcXOoUKC0tVbwtLi6WCjg7Oyv/orOzM8Mw9Qv7+vo2dmgbG5uUlBQHB4fm1PO1blXPXTK/GCqzecI5AdWQgtDCwkLbFTEiMpkMQah5UgMANONknljI19rYmKlwnyroGu3evXtCQoL0OiEhwc3NrU4y+fj43Lp1S2oZ8Dx/69YtHx+fOr+YmJgobfTx8Xnw4IFi6mlCQoK0vfVe7cwHuTLzzgki2oQAAPrpah7t76ji+wiqIAgXLFhw+PBhadboqlWrXnnlFWn70qVLDx48SAgJCgpq06bNV199VVJS8uWXX3p5eQUGBhJCFi5c+OOPP96+fTs1NfX7779fuHAhIaRr166BgYErVqwoKSlZv359bW3t+PHjW19JyfogrlhOP76GoUQAAL10LZ8OcFJxa0YFQditW7etW7e+/fbbffr08fHx+eijj6TtJSUl1dXVhBCGYf78889z58517dr1/Pnz+/fvlwpMnDhxyZIloaGhwcHBc+fOff7556Xtv//+e1JSUteuXffu3Xv06FGZTAX9txITluwfKdvzgO5KwY2JAQD0jEjJ9Xzaz0HFF3BGrx/a4ujo2PwxwrKyMqkr/0YhHXWSPzZGNtCFUXMFjRrGCDUPk2U0T3FhAQ24XUynnBauj69S7Tk3xlus9XJktgzhpp8RMiv1+EsAAICxuZpH1dGAMcYgJIRM8WZf785OOyNUY7gQAEBPXMmj/ghCFVrWl+1gw7wShSQEANAPaBGqGEPIv4dyd4vpNwmYOAMAoOtqBHK7iPZzQhCqlIWMHBrNrb8tHnuEwUIAAJ2WUEi72DGWKltG8DejDkJCiKcV8+dI7uUo/lYRshAAQHepaYCQIAgJIQGuzJoAbkqYUIDHPQEA6KoruQhCdXq+Mzu9PfNsOF+L4UIAAJ10NR9BqGZfDeQsOPJuDCaRAgDonGI5yaigvg4IQnViGbIzRHYuk/50B61CAADdcjWP9ndmOPXcDQxB+DdbE3JkDPfpdeFcFibOAADoEPXNlCEIwjo62TI7gmVzIvgHZchCAABdcSWPBiAINSbEg/m4HzfltFBWq+2qAAAAIYSQK7mivyuCUIPe9GUHuzHP4xG+AAA6IK2MsgzT1gpBqFnrg7gSPMIXAEAHXMmjAWprDhIEYWMUj/DdiUf4AgBo1WV1zpQhCMImOJuTw2O4f8QIV/PQQwoAoDWXc9Ei1J6eDszWIdw0PMIXAEBLakWSUEgHOiMItWeyN7uoOzs1DI/wBQDQgoRC2tGGsTZR4yEQhE+2rC/byZZZiEf4AgBonLr7RQmCsDkYQn4Zyt3DI3wBADQuJleNS+klCMJmkR7huwGP8AUA0KzLeTQQLUId4WnF7McjfAEANKighuRV0e72CEKdgUf4AgBoUkwuHejCsOrNQQThU3q+MzujAzPzDB7hCwCgdjG5orr7RQmCsAW+9OOsTMjiaEwiBQBQr5hcGuiq9pxCED41liE7gmXns+mPeIQvAIDaiJRcVfNdRiUIwpawNSGHR3OfXRfO4hG+AADqcauIulkwTmZqPxCCsIWkR/g+h0f4AgCoR0wuHaT+5iBBELZGiAezvB83GY/wBQBQg+hcOsgNQajz3vBlh7ozc8/iEb4AACoWnav2pfQSBGFr/d8grqyWfoRH+AIAqE5hDcmsoD0dEIT6wIQl+0bK9j6gO+5jEikAgGpczqUDXRhOEzmIIFQF6RG+Sy4LV/AIXwAAVbiUK2pmgJAgCFVFeoTvdDzCFwBAFaJzaJD6l9JLEIQqM9mbfcOXfSZMqOK1XRUAAH0maGopvQRBqEr/6sN2sWUWRgloFQIAtNiNQuplxTiqfym9BEGoSgwhPw/lkkvxCF8AgJa7lKOhFYQSBKGKSY/w3XhbPPoIWQgA0BKXcmkQglCveVgyf47iFkYJN/EIXwCAp3cphwZpaoCQIAjVxN+F+T6Am3Iaj/AFAHg6WZWkrJZ2U/NT6ZUhCNVlbmf22Y7MDDzCFwDgaVzKEQe5qvuh9P8DQahGX/hxNibMO3iELwBAs13M0egAIUEQqhXLkB3BXFQ23XQbrUIAgGbR8EwZgiBUNxsTcmQMtyoOj/AFAHiyKp7cKqIDnRGEhqWjDbMzRPZcBJ9SiiwEAGjKlTza04GxkGn0oAhCTQhuw6zoz00+LZTiEb4AAI27lEsHa7ZflCAINWZRd3aEBzP3LI9H+AIANOZCtoggNGQ/BHIVteRDPMIXAKAhIiXRuXSwm6aDCUGoOSYs2TdKtj+VbscjfAEA6rldTF3MGVcLTR9XNUGYmJg4YsQIDw+P0NDQx48f1y9QXl6+YMECDw8PPz+/kydPKrafOXNm4MCBHh4e8+fPLykpkTYuW7bM779CQkJUUkMd4WRGDo/m3rssXM5FDykAwP+4kK2FAUKikiAUBGHKlCmhoaE3b97s3Lnz3Llz65f517/+lZOTEx8fv2LFilmzZmVlZRFC8vPzp0+fvnTp0oSEhLKysvfff18qnJqaOnr06M2bN2/evHnt2rWtr6FO6eHA/DJUNiNcyKhAFgIA/O1CDh3iroUgJLTVTpw44eXlJYoipbS8vNzc3PzOnTvKBaqrq21sbK5duya9nTBhwtdff00p/eGHH0JCQqSNt27dsrS0LC8vp5TOmjVr06ZNzTm0g4NDYWFhM+tZWlrazJIa8HW84HewtrJW2/VQJ57nKysrtV0L41JZWcnzvLZrYVx06sKi79rvrr1XLD6xmMrPuQpahHfu3OnXrx/DMIQQKyurLl263L59W7lAenp6RUVF3759pbf9+vW7c+cOIeT27dv9+/eXNvr6+gqCkJaWJr1dt25d//79Z82aFRcX1/oa6qClfVgfe+ZlPMIXAIAQQsjjClrF0y52WmgRNnfV4p49e+pvHDBgQOfOnQsKCmxsbBQb7e3t8/PzlYvl5+dbWlpyHKcoIMVbQUGBp6enopidnZ30i3Pnzn3rrbdsbW3//PPPYcOGJSQkdOzYscFalZWVdejQQcpgQsj48eO3bNnS2EeoqKhQlNQFa/uR8REmn16pfd+X13Zd1EIQBLlcLgiYJas5VVVVpqamiv9roAG6dmHRX2cesgFOXEV5+RNLPtU5Nzc3l8mekHTNDcIjR47U3+jk5NS5c2cHB4ebN28qNpaUlDg6OioXc3R0rKysFEWRZVmpgJOTEyHEwcGhXOkzl5aWSr8YGhoqbendu/fFixf379//z3/+s8FaWVtbX79+3d7eXnpra2vbxCWAUmptbd28j6sJ1oQcGUsCDvP93EyneBvg9F0pCC0sND4DzIhxHIcg1DBdu7Dor2slwggvxtra7IklVX7OmxuEO3bsaOxHnTt33rx5s/S6pqbmwYMHnTt3Vi7g6elpamp69+5dX19fQsidO3d69+4t/WJUVJRUJjU1VRRFb2/vOju3tbWtrq5u7NAMw9jb2zs4ODTzU+iaNpbkwGhu4im+ow3TyxFfKgHAeEVl05e6aqdJoIKjjh8/vrS0dPfu3ZTSdevWde7cWRoOPHTo0MaNGwkhlpaWM2fO/Prrr3mej42NPXXq1PPPP08ImTt3bmRkZExMjCAIX3311dSpU21tbQVBOHr0aFVVlSAIhw4dOnny5Lhx41pfSZ3l58z8EMhNCRPyG417AAADV1RDHpbRvlpqD6ggCM3MzPbu3bt8+XJra+s9e/Zs375d2p6SkhIfHy+9/u6773Jzcx0dHSdOnLhhwwZpzM/Ly2vz5s3Tpk2zt7d/8OCBtFKCUrpq1SoHBwdbW9uPPvpo27Zt/v7+ra+kLpvTiZ3TiZmOR/gCgLG6mEMDXBmZlsaIGEpVNm+RUtr0AGZjBRrc/sS9EUIcHR1TUlKa2TVaVlamPKlHp4iUTDsjuFmQzUMMZ3QHY4Sah8kymqfLFxY9svSKYG3CLO/XrCRU+TlXZf4+MbcaK9DgdqOaiMUyZPsILjqXbsAjfAHA+ERl06FaWUpPCMG9RnWHtQk5PJr7Ml4Iz8TaQgAwIpU8uVFEA1wQhPeX/I8AACAASURBVEBIBxtmV7Bs7ln+Ph7hCwBGIyaX9nHU9MN4lSEIdcvwNsxnA7jJp4USubarAgCgEdrtFyUIQh30qg870oN57iyP268BgDE4ny0Oc9dmGCEIddHaQE4ukn9dwc3JAMDAyUVyNY8ORosQ6pCxZE+I7PAj+lsyJpECgCG7lke72jG2JtqsA4JQRzmakcOjuX9eEaLxCF8AMFyR2XSYVpuDBEGoy7rbM9uGyWaGC4/xCF8AMFDns8ThbRCE0LjxbZl/9GSfCRMqDfNJTQBg1HiRROfSIVqdKUMQhLrvvV5sLwfmxUjMIQUAQxNXQL2tGacnP3lJvRCEemDzEC69gq6Kw8QZADAoujBASBCEesGMIwdGy365J/6ZiiwEAMMRqQMDhARBqC/cLcih0dyii0J8AbpIAcAQCJRczKHD22g/hrRfA2imfk7MxsHcM2FCbpW2qwIA0GrxBdTDknEx13Y9EIT6ZWYH9sWuzLQzfA3uOQMAeu5cFh2hA/2iBEGodz7pz7lbMosuIgkBQL9FZlFdGCAkCEK9wxDy23AuroCuvYmJMwCgrwRKLuSIw7W9glCiE5WAp2IlI4dHc98lin+lY+IMAOglaYDQ1ULb9SCEIAj1VDtrZu9Ibn4kf7cYWQgA+ueszgwQEgSh/hrsxnwzkJscJhTVaLsqAABP6VymGIwghNZ7sSsb2o6ZFcHzGC4EAP3Bi7qyglCiK/WAlvnWn5OxZMllTCIFAL1xvYC2s2acdWAFoQRBqN84huwKloVl0K130SoEAP1wNpMGe+hKvyhBEBoAO1NyZDS3PFY4n42JMwCgByIyxRCdGSAkCELD0MWO2T5CNjuCTy1DFgKATpOLJCaXDtOZAUKCIDQYozyZD/tyk08LZbXargoAQOMu59Ju9oy9qbbroQRBaDje8mUHuzHPnxNENAsBQFdFZFLdWTghQRAalPVBXKmcfnwNk0gBQEedzRJDPHQrenSrNtBKJizZN1K25wHdcR+TSAFA51TyJDafDtGBp9IrQxAaGmdzcmQMt+SycCUPPaQAoFsu5NB+ToyVTNv1+F8IQgPUw4H5eSg3/YyQUYEsBAAdEpEpjtSxflGCIDRUoe3Yt3zZZ8KEKl7bVQEA+K/wDBqiS0vpJQhCg7W0D+tjz7x0XkCrEAB0QWENSSqhAa4IQtCgrUO5h+X0izhMnAEA7TuXJQ52Z0x1L3Z0r0agOuYcOThatvWeeCANWQgAWnYmg47SvQFCgiA0eO4W5OAo7vULQkIhukgBQJvCM+lIT53rFyUIQmPQ35nZOJibclrIrdJ2VQDAWD0qp8Vy2tsRQQhaMrMDO78LM+0MX4N7zgCANpzJpCM9WF2MQQSh8Vg5gGtjybx+EUkIAFpwJoOO0sl+UYIgNB4MIduGcwkFdM0NTJwBAI2ihERkiqN0bwWhBEFoRKxk5PAY7vsb4onHmDgDAJqTWEjtTZl21ghC0AFtrZj9o7iXzvO3i5GFAKAhYTrcL0oQhEZokCvzXQA35bRQUKPtqgCAcQhLF0cjCEGnzOvMTmvPzDzD12K4EADUrFog0bk0WCeX0kt0t2agVl8N5KxNmHeiMYkUANQrKpv2cWJsTbRdj8YhCI0Uy5AdwdyFbLrxNlqFAKBGYRniaE+dzhqdrhyolY0JOTKG+yJeOJOBiTMAoC6n0+kYHR4gJAhCI9fBhtkdInv+HJ9cgiwEANXLriKPK+hAFwQh6LBh7sznflzoaaFYru2qAIDBCcsQQzxYTqdzUEVBePbs2bFjxwYGBn7xxReC0MD8i8zMzBdeeMHPz+/FF1/Mzs6WNqampq5du/bFF19cunSpcuHCwsLXXnvNz89v9uzZqampKqkhNGFhN3ZcW2ZWOM9juBAAVOqUzveLEpUE4aNHj6ZMmfL8889v2rRp375933//ff0yM2bMsLKy+u2338zNzWfOnCltvH79emJiYm1tbVhYmHLhF198saysbNu2be3bt584caIo4vKsdmsCOIYh713GJFIAUBlKSFiGONZL14OQ0FZbsWLF9OnTpdfHjh1r3759nQKxsbHW1tbV1dWU0qqqKmtr6/j4eMVPd+7c2a9fP8XbtLQ0U1PT/Px8Sqkoim3atAkLC2vs0A4ODoWFhc2sZ2lpaTNLGqeiGuqzr3bLHUFVO+R5vrKyUlV7g+aorKzkeV7btTAuuLA0ITZP9NlXq/Ldqvycq6BFmJiY6O/vL70OCAhIS0srKSmpU6BPnz5mZmaEEHNz8969eyckJDS2t5s3b3bo0MHJyYkQwjDMwIEDmygMKmRvSo6M5pbHCpFZmDgDACpwKoPqQXOQEFlzCpWUlDSYRn5+fpaWlnl5efb29tIWBwcHQkhOTo6dnZ2iWG5urqKAVCY3N7exYz1V4fLy8r59+7Lsf+J87Nix3333XWOFKyoqGEYP/km0qA1Hfg5kZ4XTM6Nr21u1Ng4FQZDL5Q2OGYOaVFVVmZqachyn7YoYEVxYmnA8zeQ9X6G8XMXDW091zs3NzWWyJyRds4IwNTV15cqV9bf/+uuv3t7etra2lZWV0pby8nJCiHKSEULs7OwUBaQyyjFZR/3CXbt2baywpaXloUOHFHuzt7e3trZurDCltImfgmSSNVlRI86KYqOnyFp5JwgpCC0sLFRUNXgyjuMQhBqGC0tjSuQksbh2XAczi2blzFNQ+TlvVgX79u0bERHR2E/bt2+fnJwsvU5OTra2tnZ2dlYu4O3tff/+fUopwzCU0pSUlPbt2ze2N29v74cPH8rlclNTU0LI/fv3p02b1lhhlmXbt28vNUNBVd7wZW8V0zkR/JExMh2f9AwAOis8UwxyY1SeguqggjHCuXPn7t+/PysrixCyYcOGOXPmSH2Vv/zyy6VLlwghI0eO5Hn+0KFDhJCDBw+KohgcHNzY3vr37+/p6blt2zZCyPnz51NTU0NDQ1tfSXgq6wK5GoEsvYJeTQBoob/S6TgvPVmqrpIpNx9++KG9vb2Xl1dAQEB2dra0MTg4+IcffpBe//XXX25ubp06dXJ3d1fMAj1//ryDkhkzZkjbL1261LZt206dOjk7O+/fv7+J42LWqPoUVtOue2v/fa/lk0gxa1TzMGtU83BhaUzbnbV3ikR17Fnl55yhVDVTBCsqKsrKytzd3RsrwPN8Tk6Om5vbE8ctCSGCIGRnZ7u6upqYNDVO5ejomJKS0syu0bKyMhsbm+aUBMm9Ejr8GL9vpGyoe0t6SDFGqHmYLKN5uLA06FYRDT0tPJillo5RlZ9zlbVbraysmkhBQohMJvP09GxOChJCOI7z9PRsOgVB3brZMb+PkM2K4NPKsKACAJ7CyXQ6Th8WTkj0pAMXtGSMJ7OsDzc5TCir1XZVAEB/nHwsjm+LIARD8XYPdrAb89xZXkSzEACaoayWXMujITr8SPo69KaioEX/N4ir5Mm/rmISKQA82ZkMMdCVsdKHhRMSBCE8mQlL9o2UHUyj25JwA3QAeIITj+mEtvoULvpUV9AiRzNydAy39KpwIRs9pADQKErIyXQ6QX8GCAmCEJrPx575Y4TsWUwiBYDGJRRQSxnpYocgBAM1xpP5FyaRAkDjTjymE/WqOUgQhPC03unBBrliEikANOz4Y1G/BggJghBaYH0QV8njTqQAUFd+NblVRIe16F5UWoQghKcmTSI9/Ij+ikmkAKDkr3QxxIM107d7/CEIoSUczciR0dyyq0IUJpECwH8de0QntdOz5iBBEEKLSZNIZ0XwqZhECgCE8CIJy9C/AUKCIITWGO3JfNSXCz0tlGISKYDRi8qhnWwZdz183gyCEFrlTV92eBtmTgQvoFkIYNyOPRJD2+llpuhlpUGnrAvk5CL54DImkQIYtWOPaKgeDhASBCG0nowle0NkJx7Tn+9hEimAkbpXQqt40scJQQjGysGMHB3DfXxNOJeFHlIAY3TkIZ3UjtHLGEQQgqp0sWN2BsvmRPD3S5GFAEbn6CNxsre+Boq+1ht0UIgH8+kALvS0UCzXdlUAQIPyq8mNQhrcRk8bhAhCUKlXfdhxXsyz4TyP4UIAo3H8sTjSU/9uKKOAIAQV+y6Ak7FkcQwmkQIYiyMP6WT9nC8qQRCCinEM2R0iO59FN97BYCGA4asWSHimOFE/VxBK9LjqoLNsTciRMdzXifRMFv7AAAzcmQzaz4lxMtN2PVoB1ylQiw42zJ5g9uVLzO1itAsBDNnhh+IUvZ0vKtHv2oMuC3IlX/YTJ58W8qu1XRUAUA+RkqOPxCneejxASBCEoFZzO9JnOzDTzvByTCIFMESXcqi7BdPBBkEI0LjP/ThXC+bVKEwiBTBAhx7qfXOQIAhB3ViG/D6cu1lEv05AqxDA0Bx6SKe21/sc0fsPALrPUkYOj+Y23RYPpiELAQxHYiEVKemrnzfaVoYgBE3wtGIOjuZeuyBcz8ckUgADceghfUb/+0UJghA0ZoAzs3kI90yYkFGBLAQwBAdSxWn63y9KEISgSVPbs2/4slPChEpe21UBgNZJKaU5VTTIDS1CgKf0rz5sL0dm3jlBRLMQQJ8dfEineLOsIeQgghA0bvMQLr+afngNCyoA9NiBVHF6BwNJEAP5GKBHTFlyYLTsQBr9NQmTSAH0UkYFTSqhI/T2AYR1IAhBC5zMyNEx3LKrQmQWekgB9M+BNBrqzZoYSoAYyucAfdPNjtkRLJsdwSeXIAsB9MyfaeJ0g5gvKjGcTwJ6Z6QHs8qPm3RaKKzRdlUAoNlyqkhiIR3taSD9ogRBCNq1sBs7uR0z4wxfi+FCAD1xIE2c2JY147RdD9VBEIKWfePP2Zoyr1/AJFIA/bA/VZzRwXCagwRBCFrHMmRHMBdfSL/BXbkBdF5uFbmeT8d6GVR2GNSHAT1lJSNHRnMbb4t/piILAXTagTRxQlvW3ID6RQmCEHSEpxVzeAy36KJwNQ+TSAF0175UcWZHg+oXJQhC0B39nJifh3JTzwiPypGFALoop4rEFdBxhtUvShCEoFMme7Pv9WJDTwtltdquCgDU82eqONHg+kUJghB0zT96skFuzKwIXkCzEEDH7E0VnzW4flGCIAQdtH4QJ1KyOBoLKgB0SGYlvVFIx3gaYGoY4EcCfSdjyd6Rssgsuu4mJpEC6Ip9D+gUb4NaR6+AIARdZGtCjo3lvk0Ujz1CDymATtj9QJzV0TAjQ6aSveTn53/55ZdJSUkDBgxYunSppaVlnQKiKG7atOnUqVPu7u4ffPBB165dCSFVVVXh4eGxsbGZmZk//PCDhYWFVHjr1q3Xrl2TXltZWX3//fcqqSToF29r5uBoLvQ0/9c4WT8nAxyWANAjaWX0QRkd6WGY/xNVE++hoaEFBQXvvvvu1atXFy5cWL/At99+u3nz5jfffNPT03PYsGFlZWWEkIcPH3799ddpaWlbtmyRy+WKwuHh4eXl5QMGDBgwYECfPn1UUkPQR/4uzKYgbsppIaMC7UIAbdr9gE5vz8oMs0FICG21S5cuOTg4yOVySmlOTo6pqWl6erpyAZ7nPTw8wsPDpbeDBw/evHmz4qcZGRmEkOLiYsWWWbNmbdq0qTmHdnBwKCwsbGY9S0tLm1kSVILn+crKytbv59sEoe+B2jJ56/dk+CorK3me13YtjIuRXFj6/Fl7PkvUdi3+Q+XnXAX5Hhsb6+/vb2JiQghxdXXt1KlTXFyccoGMjIzMzMzBgwdLbwcPHnz16tWm93ngwIGXX375m2++KSoqan0NQa990Jsd6MLMxoIKAC25XUwLa8hgN8PsFyXNHyNMTk6uv9HV1dXOzi47O9vR0VGx0cnJKTs7W7lYTk6OlZWVmZmZosDt27ebOFZQUBDDMDY2NgcOHNi0aVN8fLyDg0ODJSsqKiZPnixlMCFk4MCBH3/8cWO7LS8vb+KgoHKCIMjlcp7nW7+rr3uRmedN3ois/W6ACvZmwKqqqkxNTTnOECf26SpjuLD8dkc2rS2pKK/SdkX+46nOubm5uSIjGtPcIJw8eXL9jStWrJgzZ46VlVVNzd9PVq2srLSyslIuJhWglDIMQwipqqqytrZu4ljvvPOO9OKFF14YOHDgjh073nrrrQZLmpubv/fee4q9dejQwcbGpok9N/1TUC0pCBVzoFrp4Dgy+Aj/74dmi3sa6jCFCshkMgSh5hn2hYUSsv8R/+cozsZGh1qEqj3nzQ3CO3fuNPajtm3b7tu3T3otiuKjR4/atWunXMDT01MUxfT09LZt2xJC0tLSpBdPxLJsly5dcnJyGivAcdzw4cMbay+CIZEWVAw+KrS3IVO8kYUAGhKTS8050tegZ26r4IIyadKkpKSk2NhYQsjhw4ctLS0DAwMJIVeuXPnrr78IIXZ2dmPHjv3ll18IITk5OceOHXv22Wcb25sUpdLr5OTksLAwaW8A3tbModHcqxeEa/kYLQTQkB33xbmdDfyrpwrWEdrb269du3bMmDHdu3e/d+/etm3bpJ6Zw4cPJyYmjhs3jhDy7bffTpgw4dSpU6mpqfPmzfPz8yOECILg4uIiiiIhpH379hYWFpmZmYIg9OjRo23bthYWFsnJyW+99dbEiRNbX0kwDH7OzJYh3DNhwqVQrp21IX9FBdAFtSLZlypenqyaFec6i6FUNV+uCwoK0tLSunXrphixq66uFgRBMV4ol8tv377t7Ozs5eWl+C3lSaEMw9jb20slk5OTBUFo3769ra1tEwd1dHRMSUlpZtdoWVmZYXfl6xrVjhEqW3dT/PmeeCFUZmeq8n3rN0yW0TzDvrAce0S/SRSiJulWEKr8nKvs4zk5OTk5OSlvMTc3V35ramrat2/fOr/VYIaZmpr26NFDVRUDw7O4J5tSRmeE8yfGykwMvM8GQJu23xefN/R+UYJ7jYKeWhvIWXDM6xfwhAoAdSmRk1Pp4swOhh8Thv8JwSBxDNkVwiUW0i/i8YQKALX4M00M9mAdzbRdD/VDEIK+spKRI2NkP98Td9xHFgKo3u/J4rzORjElDUEIeqyNJTk2hltyWYjMwoIKAFVKK6O3i+jEdkaREUbxIcGA9XBgdgXLZkXwd4qRhQAq88d9+mxH1tQ4IsI4PiUYtBAP5lt/buIpIUdXboUIoN8oIb8niy90MZaAMJbPCYbthS7si13Z0NN8JW7KDdBqF7OpKUv8XYxigJAgCMFgrOjH9nRg5pwV8LQmgFb6LVmc39WI0sGIPioYvM1DuCqeLo7G4kKAlqvkyYE08XnjmC8qQRCC4TBhyf5Rsqhs+t0NLKgAaKEDaWKgK+NhiSAE0E+2JuT4WG79LXHvA2QhQEv8miS+ZEz9ogRBCIbHy4o5NoZ7O1qIysZoIcDTSS2jN4voZCN75KdxfVowEr0cmZ3BspnhWFwI8HR+TRKf62QsywcVjOzjgtEY6cGsDuAmnhKysbgQoHkESrYl0QVG1i9KEIRgwOZ1Zl/uxk48xZfVarsqAPrgdDptY0l6ORrRNBkJghAM2Ud9WT9n5tlwvhZTZwCe5Jck8eVuxhgKxviZwahsDOJkLHn9AtbZAzQlp4pEZIqzOxljKBjjZwajImPJ7hDZzSK6MhYL7QEa9VuyONWbtTXRdj20AUEIhs9KRo6Oke1MoVvvoocUoAGUkJ/via/4GGkiyLRdAQBNcLUgJ8dxw47xbSyZSe2Mbi4AQNPOZlILjgS6Gul/DSPNfzBCnW2ZQ6NlL0fxl3MxXAjwP7bcFV811uYgQRCCUfF3Yf49TDb1DJ9cgiwE+I/cKnI6Q3y+s/HGgfF+cjBOE9syqwZw4/EUX4D/+jVJnOrN2plqux7agyAEo/NyN3Z+Fyy0ByCEEJGSLXfF17sbdRYY9YcHo7W8H+vnzMw4g4X2YOxOZ1BHMzLQaB5G3yAEIRipjYM5Sxmz4DwW2oNR+/GOuMjX2IPA2D8/GC2OITuDudQy+s/LWGgPRuphOb2UI87uaOxBYOyfH4yZhYwcHSM7mU6/xxPtwShtviPO68xaGv16cgQhGDUHM/LXOG7dLXFnCrIQjEuNQH5JEhcZ9zQZCU4BGDsvK+bkOG5JjHA6A8OFYET2PBD7OzFd7Ix6mowEQQhAfO2ZA6Nk887xV/OQhWAsNtwW3+rBabsWOgFBCEAIIUFuzC9DZVPC+CTcdAaMQEwuLaoh473QHCQEQQigMKkd86UfN+4vIbMSWQgGbv0t8U1flkUOEkIQhADKXuzKvubDjvtLKKrRdlUA1Cazkv6VLr7UFdf//8CJAPgfS/uwoz2Z0NN8Fa/tqgCox6bb4tzORn1z0ToQhAB1fRfAdbJlno3geSypAINTxZOt98S3jf5uMspwLgDqYgj5ZShHCMEN2MDwbL8vBriwWDWhDEEI0AAZS/aGyFLL6XsxuAEbGA5KyA83xX/0wpX/f+B0ADRMugFbeCb9Ih49pGAgTqVTU44Et0Fz8H8gCAEaZW9KTo2XbUsSf7qDLARDsOaGsATNwXqM/marAE1ytyCnxnPDjwkOZmSW0d+kH/RaQiG9W0zwrIn6EIQAT9DRhjkxlht9krczZcbhThygt75LFN/uwZogB+vBKQF4sl6OzKHRsvmR/MUcTCMFvfSonJ58LL7mg2t+A3BSAJol0JXZPkI27QwfX4AsBP3zw03xpa5YRN8wBCFAc432ZH4czE08JeDG3KBfCmvIb8niuz1xwW8YxggBnsK09mypnIw5KZyfxLWzxngh6IdNt8VnvFlPK/zFNgxBCPB0XuzKlsjJ6JPC+UkyNwtt1wbgSSp5suG2cG4irvaNQksZ4Kkt7snO7cyOOcnjIRWg+36+Jw5xZ33s0RxsFIIQoCVW9GNHezLjT/FltdquCkDj5CJZc0Nc1geX+qaoprH88OHDTZs2lZSUTJ48ecKECfUL1NTUbNq06c6dO76+vosWLTIzMyOE3Lt37/Dhww8ePHB1dZ0/f36nTp2kwjzPb926NS4urmPHjm+//baVlZVKKgmgWqsDuEUXhMmn+RNjZRbodgKd9Eey2N2eDHBGc7ApKviaUFxcHBgYWFtb6+fnt2DBgj179tQvM3/+/KNHjw4bNuzw4cMLFiyQNr7xxhvp6el+fn6lpaV9+vS5deuWtP3dd9/dtm3bkCFDLl68OHXq1NbXEEAdGEI2DeY8rZgZ4bwct2AD3SNQ8nWC+FFfTtsV0Xm01X744YeQkBDp9R9//NGvX786BVJSUszMzAoKCiilBQUFZmZmqamplNLa2lpFmUmTJq1YsYJSmp+fb25unpKSQimtqqqys7O7du1aY4d2cHAoLCxsZj1LS0ub/6Gg9Xier6ys1HYt1K5WoFPD+Oln+FpB21WhtLKykud5bdfCuOjyhWV7sjD8WO2Ty+kblZ9zFbQIL168GBISIr0OCQmJi4urrKxULhAdHd27d29HR0dCiKOjY69evaKjowkhMtnf3UnFxcUODg6EkNjY2DZt2nTs2JEQYm5uPmjQoIsXL7a+kgBqImPJ7hCuopa+dF4QsbwQdIZIyRfx4vJ+aA4+WbNGNgRBqKqqqr/d0tKSZdns7OyRI0dKW1xcXAghWVlZigE/Qkh2drazs7Piraura1ZWlvJ+tm/fnpqa+tJLL9Uv7OLiUqewsqqqqhdeeMHU9D83SxgwYMC7777bRGGOw9+E5giCIJfLKTWKcNgeRKae4xZG8uv9BS2OxlRVVfE8j79zTdLZC8u+NMbOhBtkL//fhokheKpzbmpqqtzoalCzgjA6OnrSpEn1t8fExPj4+JiZmcnlcmmL9MLc3Fy5mJmZWW3t31PrampqlAucPn36vffeO378uJ2dXf3Ccrm8zt6UmZiYTJ061draWnrr7e3dROHa2tomfgoqJwgCy7JGcs7NCTk2low/LS6L534I1NoMPUqpqampbl6XDZVuXlhESr69LX4fyJqbG+A8rqc65yz75P+PzTpHQ4YMKS4ubuynnp6e6enp0uvHjx/LZDI3N7c6BR4/fqx4m56e7unpKb0+e/bsvHnzDhw44OfnpyickZFBKWUYRtrhqFGjGq29TDZ16lSpT/WJWJZtzhkBVaGUGtU5tzUjJ8ayI0/wy67Rb/y1E0Xsf2nl6MZJN0/43hTR3oyM9dK5iqmEys+5Cvb1zDPPHDx4UOo73bVr16RJk6R26IULF1JSUgghI0eOzMzMvH79OiEkNjY2KytL6kq9dOnS7Nmzd+/ePXjwYMXeAgICZDJZWFgYIeTBgwdxcXENNkYBdJCdKTk9XvZXOl0RK2i7LmC8BEo+ixNX9kfHQHOpoNUcGhq6ZcsWf3//jh07Xr58WcowQsjSpUunTp36/vvv29rafvnllxMmTBg2bNj58+e/+uorqTPz5ZdfFgThgw8+kMpPmjRp5cqVJiYmq1evnjt3bnBw8MWLF5ctW+bu7t76SgJohqMZOTNBNuIYb8aJH/U1zO/joON2pYjO5mS0J9YONhejkrkMoihevXq1oKBg8ODB0lAfIeThw4c2NjbSZFFCSEpKirSgXpoRSgh5/Pix8nCgjY2NNNeGEPLo0aPExMTOnTv7+Pg0cVxHR8eUlJRmdo2WlZXZ2Ng87UeDFpMmy1hYGOPtOLMqyYjj/MJu7Ae9NZqFVVVVGCPUMF27sPAi6b6f3zqUG9HGYINQ5edcNUGoLQhCXWbMQUgIyaigI44Lb/qymnz2DYJQ83TtwrL1rrg3VQwbb4BzZBRUfs4N+WQBaJGnFRMxkRtxTJCx5C1f9JGCJlQL5PN4cW8Ivgk9HQQhgLq0tWLCJ3DBJwSOIYu6IwtB7X68I/ZzYgJcDbZTVE0QhABq1N6GiZjABR8XWIa85oMsBDUqrSXfJAjhE3BVf2o4ZQDq1cGGiZjIBR8XGEJeRRaC2nyXKIzzYns4oDn41BCEAGrX0YaJmMCFnBAIshDUI7uK/HhHjH0Gl/SWwFkD0IROtshCUKNPrwsvdmHb0HGEEgAAGnJJREFUWaM52BIIQgANkbJw5AlBpOR1zJ0B1blbTA+kiXdnmGi7IvoKQQigOZ1smYiJXMhxQaDkTaypABX55xXxX304BzNt10NvIQgBNKqjDXNuIhdyQhAoeacHshBa62wWvV1M94/C2sGWQxACaFr7/2ZhrUje64UshJYTKVkSI3wzkDXF31Er4OQBaEE7a+bcRG7zXfGrBFHbdQE99u8k0daETO+AK3mr4PQBaIeXFRM5UfZHsrjyOp7ZBC1RIicrYoUfBqFTtLUQhABa08aSnJsoO5hGl15BFsJTWxUnTGzL9nPCkonWQhACaJOrBYmYKAvPpO9EC3r8IBjQuLvF9Pdk8Qs/NAdVAEEIoGVOZiR8giw2n74SJYgIQ2ied6KFj/pyrkb6lDMVQxACaJ+dKTk9XpZaRp8/J9Ri9gw8yZ+pYnYVlqKqDM4jgE6wkpHjY2VltXRGuFCNEUNoXHktWXJZ3BjEyXD9VhGcSABdYc6RA6NkFhyZdIovr9V2bUBXfRYnjGjDDHXHHBmVQRAC6BATluwI5jraMqNO8kU12q4N6J4bhfS3ZHG1P+bIqBKCEEC3cAzZPIQb6sYMP85nVWq7NqBLREpeuyB8PgBzZFQMQQigcxhCVgdwszuyQ4/xD8owkRT+Y/NdkWPJQjzGS9Vwr1EAHfVhX9bJnAw/Jhwfy/V2xICQscuooJ/ECpGTZPhTUDl8swDQXa/5sN8HsqNP8lHZaBcau0UXxbd6cN3tkYOqhyAE0GkzO7A7RshmhPOHH2KBofHalSKmldN/9cEVWy1wWgF03ShP5sRY2aKLws/3kIXGKLeKLIkR/j2Mw7OW1ATnFUAPDHBmzk+SfZ0gropDFhqdNy4JL3Zl/ZzRKaouCEIA/dDZlrkYKjv8UHz9goD7cxuPXSni3WK6sj8WDqoRghBAb7hZkLMTZWnldGqYUMlruzagfhkV9N0Y4ffhnBlyUJ0QhAD6xMaEHB0jczInwcf53Cpt1wbUiRKy4Lzwli/XH52iaoYgBNAzJiz59zBunBcTdJS/V4JOUoO1/pZYWkuWYaao+uEUA+gfhpBPB3Af9WWHH8MSQ8N0s4h+Hif8MQKPmNAEnGMAffVSV3b7CNmMcH7HfUwlNShVPJkTIawO4DrbolNUExCEAHpslCcTMUG2PFb89Doebm84/hEj9HZk5nfB9VlDcKIB9FsPByZ6suyvdPH5s3iiryHY+0AMz6Q/DsE8Uc1BEALoPWlZhUjI+DNcDqaS6rP7pfTtaGFPCGdrou2qGBMEIYAhMOfIzmButAcNOkbjC9BLqpeqeDIzXPikH9ZLaBqCEMBAMIR82Ev8diAZ+xf/Zyqmz+ifNy8JPRyYN3xxWdY0PI8QwKBMb890seemhgk3iugn/Tm0LPTFT3fEa/k0ejKuyVqArx4AhqafE3Nliiwik04LE8pqtV0baIaLOXTldeHAKM4KOagNCEIAA+RqQc5MkLlbksDDfBLuPqPb0ivos+HCtuEyrBrUFgQhgGEyZcmPg7klvdihx/ijjzBkqKMqeTIlTPhHL3acF1JQa9AOBzBkL3djezkyM8OFK3l0ZX8OY4Y6RaRk3jmhlwPzfi+0SbQJZx/AwPm7MNeekV3KoRP+4vOrtV0bULLsqpBfTbcMxdp5LUMQAhg+F3NyerxsgDMz4BAfnYshQ53w0x3x8EN6YLTMFJdhbcO/AIBR4Bjy5UBuYxA3NYxfexM3JtWyIw/FVXHiiXGck5m2qwIIQgCjMqkdc3mKbM8DcWqYUFij7doYq4s59JULwuExXEcbjNnqBAQhgHHxtmaiJsk625L+B/mLOWgZalpCIZ1+ht8+QuaH+6jpDAQhgNExYcl3AdzGwdzMcH5VnCggDTUlqYRO+EvYEMSN9kQK6hAEIYCRmtiWufaM7Hy2GHycf1SOMFS7B2V09EnhCz92RgdceHWLatYR8jy/a9eue/fu+fn5PfPMMw2WuXDhwunTp93d3efNm2djY0MIqa6uDgsLi4+Pt7CwmDBhgq+vr1QyPDw8JSVFem1mZjZ//nyVVBIA6vCwZE6Nk625IQ48zP8QyM3phAu0uqSV0ZEnhGV92Be74iTrHNX8k8ybN2/Tpk12dnYffvjhhx9+WL/Arl27pk+fbmFhERYWNnTo0NraWkLIRx99tGbNmpqamsePHwcEBBw+fFgqvHXr1j/++CM2NjY2NjYhIUElNQSABrEM+aA3+9c42edx4pyzmEGjFqllNPiE8H4v9vXuSEGdRFvt7t27FhYWhYWF0msrK6uioiLlAqIodu/efffu3ZRSnue7d+++b98+SmlpaamizKpVq0JCQqTXs2bN2rRpU3MO7eDgIB23OZQPBxrA83xlZaW2a2FcKisreZ5v4e/W0nejea+dtSceiaqtlWF74oXlXrHYblftj7cFzdTHGKj8Yq6Crydnz54NCAhwcHAghHTr1s3Nze3y5cvKBbKzs+/cuTN+/HhCCMdxY8aMiYiIIIRIHaQShmHMzP5eUHP16tUNGzaEhYWJIu6RCKAJFjKyNpD7fQT35iXhlSihFI+tUIWEQhp8XFjZH21BnaaCMcLs7Gw3NzfFWzc3t8zMTOUCWVlZFhYWtra20lt3d/eYmBjlAo8fP167du3u3bult23atCkuLk5KSlq3bp2Xl9fp06dNTEwaPHR1dfWKFSsUCdq9e/e5c+c2Vs/q6urG9gPqIAiCXC5nGMyO05zq6mpRFDmu5bfsGuRALk8gH8WxPfeL6/3FsZ6YRPMETVxYLuYyc84z6/zp1HZ8NW5upzpPdTE3MTF54v+I5gahcutNYe3atQsXLmRZVrndJghCnaNyHCcIgnIBmezv4+bn50+YMOHdd98dNWqUYrfSi8rKyl69eu3cubOx+TIsy9rb21tYWEhvXV1dm/jAHMe15gIBLYBzrmHcf7VmJ/Yc2TiInM0ir19igx6S1QOJs7mqKmiAGjvhhx6RN6PJ70PJSA98F1Sxp/ojb8538eYGYW5ubv2NpqamhBAPD4/Tp08rNmZlZXl4eCgXa9OmjVwuLygocHJyIoRkZmYqChQVFY0dOzY0NPTjjz+uv39LS8uAgICkpKTGamVqarpkyRKpV/aJTExM0CLUJJZlKaU455rE83xzvv82x5h25KYH+fia0P+I+G0AN68zevYa1uCF5Yeb4pob4qnxXD8npKDqqfxi3tw/bouGSP/fxo4dGxsbK3WHXrt2raysLCgoiBCSnp6enJxMCHF1dfXz8zt48CAhpKam5uTJkxMmTCCElJSUjBs3btiwYV9++aXiQJRSnuel1+Xl5dHR0T4+Pir8wADQTJYy8n0gd2ysbO0NcdQJPOC3WXiRLLoo/DtJvBiKFNQbKhgjbNu27WuvvRYcHDx27NgDBw6sXLnS0tKSEPLjjz8mJiYePXqUELJq1aq5c+dev349Pj6+U6dOY8aMIYR8/PHHcXFxtra2o0ePJoS4ublt3769tra2Q4cOwcHB5ubmZ86c8fHxmTNnTusrCQAtM8CZuTJFtuG2OPgo/3p3dlkfzhKPMW1EfjV5Npy3MiEXQ2U26ArRHwylqvmWFxkZmZSUNGDAgP79+0tbHjx4UF5e3rt3b+ltSkpKZGSki4vL+PHjpTHC27dvK0+rMTc3HzJkCCEkMTExISFBEISuXbtKjcvGODo6pqSkNLNrtKysrMGRTlATabKMYgQXNKCqqsrU1FRN47KZlfS9y2J0Dv0uAPdG+ZviwhKbT2eEC891YlYN4Fg0BdVJ5RdzlQWhViAIdRmCUPPUGoSSyCz6TrTgZEa+D+T6ouvvvxeWH++IK68LPw3mprbHVwS1U/nFHH0cAPAUhrdhrk+Vbb0rjv//9u4/qIlrXwD42d3wG/kVBAKGwOOHP0BSioIiSimgoIWLVqVFLC32jc92fFNnWqevT/vmldfpnV7bubYXe1vgTa3aH4hWRIowipA+axF811DFIDWpPyDBiAQBY5LdPe+PvXdfBiiCSiDZ7+evsydns985bM53N3tyOEmvlpIli8hgd0GnQ4OZKD7NXLuHz+aIIr0E3RX2Cy5eAACTQxHoX+aTqg1Os91Q3BF6VxszYJ7umKbJGS1eVu8c4oHO5UIWtGOQCAEAj8LbGf1xMfW3daKe+yj6sOVP7ex9erpjsqFhGv3rOealJuaTxfSfl1Au8HNZewaJEADw6KQexH+voJrWiFr1OLLS8udLrFEA6bD+Fl54hL5nRu3rROlBzMN3ADMbJEIAwOOa70NUplN1WaIfdTii0vKndnbQQZcq7R7GLzQyr51lPltGfZlK+bo8fBcw80EiBAA8GXI/4kgG1ZAtutiHI76z7GpjdMbpjunJMdLov/7Gyo/S0d7o0vOiVXPgiaDjgEQIAHiSYn2JQ2nUz38QGcwopspSrGAu9tnxb7QQQjSLKjrZuYfp9ru4NU/0XgLlBtPtHQskQgDAk/dPs4i/JFNdG52ivYncBmb5Cfrra6zJ3p6m0Sw68Csbc4Q+9Ct7OJ2qTKfCZ8GNoAOCCxsAwFTxc0Fvy8k3F5I1N9i/XmF3/MwURJCvRJNxfjM9nQzT6Mur7Ee/sGGe6K8pVJpkpgcMHgckQgDA1BKRaG0YuTaM1AziL6+yuQ2MjzMqjCI3hhOhnjMuwVwdwJ+r2K+62Gck5Ndp1JKAGRcheOIgEQIAbCR8FvGfCdR/PI1+1OGvr7EJx5goL2JdOPmHUCLKe5rzjcGMjmjY/V3sr/fwy1FkW55INvOSNJgikAgBADZFEihVQqRKqNJk6owWH/2NTfuB9RShbCmxMoRcISE8bDgs9RrRiRvssevsjzqcGUK+uZDMlpJOMHdCYGDRbTBVYNFt27PBottTASN0sQ+fvIUbbrFtd3CcH7EskEgOJBJnE1OxkOldEzp3Gzf1sKd78PUhnBlC5smINaHko/3jJBhYbA8W3QYAOBoCoXgxES8m/k1OGmn0sx7/jw6Xqdh//hE7kUjuRyz0I+b5EFFeRIQXkrgTk8qNDxikGcSdA/iKAbXfxf97B2vv46QAYkUQ+ZdkMnE2IYL7P8GDRAgAmEHcRChNQqRJCO7HXTeGsPIuvtyPFFpcrmI1g/iuCQW5E8HuSOyCvJ0JTyc0ywnxyczEoPs0MphRvwn3GlGvEfebUNgsIsoLLfAlVkuJf3+KnO9DUPD4D1iBRAgAmLlCPYlQTyIn9P9rHjBIdx9rjajvARow42Ea3bMghv37q84U8hAhLyfk60IGuKEgNyLIHUHWA+ODRAgAsCeuFAqbRYT9/QkR5DjwBMC34wAAAAQNEiEAAABBE0oiHB4eLi0tne4ohKWjo6Ompma6oxCWmpqajo6O6Y5CWEpLS+/fvz/dUQjIgwcPPvnkkyf7nkJJhDqd7osvvpjuKITl/PnztbW10x2FsNTW1p4/f366oxCWzz//XKfTTXcUAnLnzp3PPvvsyb6nUBIhAAAAMCZIhAAAAAQNEiEAAABBs++1Rl1dXSUSCUk+PJ3TNK3T6ebMmWODqABnaGjIaDTOnj17ugMREL1e7+bm5unpOd2BCMjNmzclEolIBL/JthGGYXp6eqRS6QTbFxQUlJSUjN/Gvv94165dM5lME2xsMplcXFymNB5gjWVZhmGcnB5pJWPwSCwWC0VRE7k0BE8KDCy2N6k+l0gkD21j33eEAAAAwGOCK0cAAACCBokQAACAoEEiBAAAIGiQCAEAAAiafc8atXbx4sVvv/3W2dm5qKgoIiJidIO+vr7y8vLe3t7s7OzMzEy+vrGxsba2dvbs2cXFxQEBATYM2e4pFIrq6mqxWFxcXBwUFDTi1d7e3pqaGpVK5efnl5+fz/9Rqqqq7t69y5WDgoJyc3NtGrQ9M5vNFRUVXV1dTz31VGFh4ejZoT/88MOtW7e4speX1wsvvMCVBwYGysrKenp6nn322eeee86mQdu5wcHBsrKyW7durVixIi8vb8SrBoOhsrLSumb58uXz58+/fv16fX09X5mdnT3x6f4C19/ff+HCBbVavXjx4vj4+DHbXLly5dChQyzLbtq0KSYmhqtkWfbAgQNKpTIqKmrLli3Ozs4TP6iD3BG2trauWLHCx8fHYrEkJiZev359RAOTyZScnHz58uXw8PCioqIDBw5w9ZWVlS+++KJMJuvq6lq6dCksnjtx1dXV69atk0qlN27cSEpKunfv3ogGW7dubWxsDA4O1mq1cXFxra2tXH1JSUlTU5NarVar1d3d3TYP3I4VFBR89913UVFRe/fu3bFjx+gGe/fura2t5fr25s2bXCXDMM8880xLS0tERMT27dth9fmJwxhnZmYqFIrIyMi33nrro48+GtHAYrGo/+HSpUtbt27t6+tDCCmVypKSEv4lo9E4HeHbpQ0bNrz99tvvvfdeXV3dmA06OzuXLFlCUZSbmxs3qnP1b7zxxt69e6OioiorKzdt2jS5o2KHsHHjxt27d3Pll156aefOnSMaHDx4UC6XsyyLMa6qqpo3bx5XlsvlBw8e5NokJSVVVFTYMGr7lpSUVFZWxpVTU1NLS0tHNDAajXy5qKho27ZtXDkuLu7MmTM2idGhqFQqNzc3g8GAMdZoNK6urnq9fkSblStXfvPNNyMqjx8/HhERQdM0xri+vl4qlXJl8FANDQ1z5syxWCwYY4VCERgYaDKZfq/x/v37o6OjuYGluro6JSXFdoE6EIZhMMZr1659//33x2ywbdu21157jSvv2LFjy5YtGGO9Xu/q6qpWqzHGAwMD7u7uKpVq4gd1kDvC5uZm/tvOzMzM5ubmEQ0UCkVGRgZBEFwDlUrV29trMBiUSmVGRsY4O4IxPXjw4Pz58+P3uaurK182mUzWy51UV1d//PHHDQ0NNgjVYSgUisTERG9vb4RQWFiYTCZraWkZ3ezUqVN79uw5fvw4y7JcTXNzc3p6OkVRCKG0tDStVnvt2jVbRm6/mpub09LSuFVjli1bNjQ01NnZ+XuNKyoqiouLuUEGIaTX6/fs2VNeXt7T02OjcB3CQ5eDGHO0b2lpkUql4eHhCCEvL6/ExESFQjGJgz5qtDMITdN6vZ5fyisgIECr1Y5oo9Vq+QZeXl6urq5arVar1ZIk6e/vz9UHBgbCKTtBOp0OY8w/Uh2/6xQKRV1d3euvv85tyuVyiqK0Wu2rr77KP8QCD2V9DiOEAgICRvf53LlzPTw87ty5s3PnzpUrVzIMgxDS6XT8jk5OTn5+fqM/IGBM1l3HjRW/d56r1epz585t3ryZ2/Tw8IiPjzcYDCdPnlywYMG5c+dsFLEAWH8Q+NF+xKdjsoO5I0yWIUmSoiiaprlNmqZHPyYViUR8A+7u29nZ2cnJiVsGjLtYtlgssFTSBHELp/FdOk7X/fLLLxs2bKioqJDJZFzNV199xRXefPPNqKios2fPLlu2bOpDtnsikYhLbByLxTL6POf/YemuXbvmzp3LPcedyI5gTBPvuvLy8uzs7ODgYG4zPT09PT2dK7/zzju7du06ffr0VEcrENaDOT/aP+ZJ7gh3hCRJBgUF8dMuuru7+dORFxISwl8g9Pb2WiyW4OBgiURCEARf393dPZFV6QBCKDAwkKIo6z4fs+s6OjpWrVr16aefPv/882O+SWRkpEajmdpYHUVISIj13KIxz3Oep6dnXFwc17fWOw4NDQ0MDIyzI7Bm3XUmk6mvr2/MrqNpev/+/cXFxWO+SXJyslqtnsIoBcZ6MOc/BZP6dIzmCIkQIZSbm1tVVYUQwhhXVVXl5ORw5cbGxsHBQYRQTk5ObW0tNym0qqoqJSXF19fXw8MjPT2d29FkMtXU1MBU/gkSiUTZ2dmHDx9GCFkslurqaq7rjEZjY2OjxWJBCHV1da1aterDDz/cuHEjv6PZbOYv3Lq6ulQqFT/7GYwvKyurvb2dG1JbWlqGhoZSUlIQQhqNRqlUIoQYhuF6HiGk0+laW1u5vs3Jyamvr+em9R49enThwoX83TkYX05OzunTp/v7+xFCx48fDwsLi46ORghduXLF+mHhyZMnGYZZvXo1X2M9TfTEiROxsbE2jNoBGQwGfhZCbm4uN/IghA4fPsyN9ikpKcPDw9xX0BqNpr29PSsraxIHeJzpPTPHb7/9FhISsm7duoyMjJiYmP7+foyx2WxGCLW1tWGMGYZZs2ZNfHz85s2bxWJxc3Mzt+NPP/0kFosLCwsXLVqUmZkJs+kmrq2tTSwWFxQULFmyJDU11Ww2Y4y7uroQQrdv38YYp6WleXp6JvzD9u3bMcYXL16USqXr169fv369l5fX6Pm9YBy7d+8ODQ3lfrW5b98+vjIrKwtjrNVqAwMD8/Ly8vPzxWLxK6+8wu+Yn58fExNTVFTk7+9fV1c3PdHbp6Kionnz5r388sv+/v7ff/89X7l161a+TV5e3ogzOT8/PzU1tbCwMCEhITQ0tKOjw6ZB27MPPvggISHBx8cnJCQkISHh2LFjGOPGxkYXFxeuQW9vb0RExOrVq3NycsLDw7VaLVe/b9++wMDA4uJimUz27rvvTuqgjvPfJwYGBk6dOuXi4pKRkcHNV8QYt7a2xsbGuru7I4RYlm1qatLr9cuXL7e+a9bpdM3NzWKxOC0tjXtYCCbo9u3bTU1Nvr6+/Mw6k8mkVCqffvppkUjU2dk5NDTEN/b29o6MjMQYX758WaVSkSQpl8vHXPoAjOPChQtXr16Vy+ULFizgarq7u4eHh7k7lc7Ozo6ODpqmY2Nj58+fz++FMVYoFFqtNjk5OTQ0dHpCt08Y47Nnz968eXPp0qVhYWFcpUajIQiC31QqlTKZzMfHh99rYGCgpaWlr68vKCho6dKl1jOowfhu3Lih1+v5TZlM5u/vPzg42NnZuWjRIq5yeHj41KlTGOOMjAzr6egdHR1KpTI6OjohIWFSB3WcRAgAAAA8Agd5RggAAAA8GkiEAAAABA0SIQAAAEGDRAgAAEDQIBECAAAQNEiEAAAABA0SIQAAAEGDRAgAAEDQIBECAAAQNEiEAAAABA0SIQAAAEH7P35GkmeO10NGAAAAAElFTkSuQmCC",
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
""
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 10000\n",
"f(x) = -x^3\n",
"x = [j/n for j ∈ 1:n-1] \n",
"L_dense = diagm(-1 => ones(n-2), 0 => -2ones(n-1), 1 => ones(n-2))\n",
"t_dense = @elapsed u = ( -L_dense \\ f.(x) )/n^2\n",
"println(\"Solving the linear system of equations took $t_dense seconds\")\n",
"plot(x, u; title=L\"numerical solution of $-u\\prime\\prime = f$\", legend=false)"
]
},
{
"cell_type": "markdown",
"id": "0bf0b16e",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 3.** Print ```L_dense``` in the following cell. How many entries of this matrix are being stored in memory? "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e46a1a9a",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "7e622db4",
"metadata": {},
"source": [
"::: {.box}\n",
"\n",
"\n",
"***Your answer here***\n",
"\n",
"\n",
"\n",
"\n",
":::"
]
},
{
"cell_type": "markdown",
"id": "ecc94871",
"metadata": {},
"source": [
"We can use ```Base.summarysize()``` to find out how much memory this is using:"
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "182af5f5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Size: 762.79 MB\n"
]
}
],
"source": [
"bytes = Base.summarysize(L_dense)\n",
"println(\"Size: $(round(bytes/1024^2,digits=2)) MB\")"
]
},
{
"cell_type": "markdown",
"id": "97078db9",
"metadata": {},
"source": [
"The problem here is that we are storing lots of zeros (that we don't need to!). Instead, we should use a \"sparse format\" where we only store nonzero entries (and where they are). To do this, you need to use the package ```SparseArrays```: "
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "3cbf3a97",
"metadata": {},
"outputs": [],
"source": [
"#import Pkg; Pkg.add(\"SparseArrays\")\n",
"using SparseArrays"
]
},
{
"cell_type": "markdown",
"id": "34c14c62",
"metadata": {},
"source": [
"Here, we can instead define ```L``` to be the discrete Laplacian but as a sparse matrix:"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "47bc9d96",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"9999×9999 SparseMatrixCSC{Float64, Int64} with 29995 stored entries:\n",
"⎡⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤\n",
"⎢⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⎥\n",
"⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⎥\n",
"⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⎦"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L = spdiagm(-1 => ones(n-2), 0 => -2ones(n-1), 1 => ones(n-2))"
]
},
{
"cell_type": "markdown",
"id": "b73c6861",
"metadata": {},
"source": [
"```L``` represents the same matrix as ```L_dense``` but requires a lot less memory to represent it:"
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "195d62eb",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Size: 0.53 MB\n"
]
}
],
"source": [
"bytes = Base.summarysize(L)\n",
"println(\"Size: $(round(bytes/1024^2,digits=2)) MB\")"
]
},
{
"cell_type": "markdown",
"id": "732da1ea",
"metadata": {},
"source": [
"**Exercise 4.** As a percentage, how much memory have we saved by storing $L$ in a sparse format?\n",
"\n",
"::: {.box}\n",
"\n",
"\n",
"***Your answer here***\n",
"\n",
"\n",
"\n",
"\n",
":::"
]
},
{
"cell_type": "markdown",
"id": "077b141e",
"metadata": {},
"source": [
"It turns out that Julia is storing ```L``` in the following Compressed Sparse Column (CSC) format where three vectors are stored:\n",
"\n",
"* ```nzval``` (non-zero values) : the actual numbers that are in the matrix (length is the number of nonzero entries), \n",
"* ```rowval``` (row indices) : integer vector of row indices corresponding to each of the entries in ```nzval``` (has the same length as ```nzval```), \n",
"* ```colptr``` (column pointers) : integer vector of length $n+1$ (where the matrix is $n\\times n$) that tells you which entries in ```nzval``` and ```rowval``` correspond to which column: ```colptr[1] = 1``` and ```colptr[k+1] = colptr[k] + (the number of nonzero entries in column k)```.\n",
"\n",
"That is, if $L_{ij} = a$ then there exists $k$ such that \n",
"\n",
"* ```nzval[k] = a```, \n",
"* ```rowval[k] = i```, and \n",
"* ```colptr[j] ≤ k < colptr[j+1]```"
]
},
{
"cell_type": "markdown",
"id": "96bc23fa",
"metadata": {},
"source": [
"**Exercise 5.** Write down a valid (```nzval```, ```rowval```, ```colptr```) representation of the matrix \n",
"\n",
"\\begin{align}\n",
" \\begin{pmatrix}\n",
" 1 & 0 & 0 & 0 \\\\\n",
" 0 & 1 & 0 & 1 \\\\\n",
" 2 & 0 & 1 & 0 \\\\\n",
" 0 & 0 & 0 & 3\n",
" \\end{pmatrix}\n",
"\\end{align} \n",
"\n",
"::: {.box}\n",
"\n",
"\n",
"***Your answer here***\n",
"\n",
"\n",
"\n",
"\n",
":::"
]
},
{
"cell_type": "markdown",
"id": "641f6b98",
"metadata": {},
"source": [
"Notice that, if $A = [u_1 | u_2 | \\cdots | u_n]$ where $u_j \\in \\mathbb R^n$ is column $j$ of $A$, then \n",
"\n",
"\\begin{align}\n",
" (Ax)_i = \\sum_{j=1}^n A_{ij} x_j = \\sum_{j=1}^n (u_{j})_i x_j \n",
"\\end{align}\n",
"\n",
"This leads to an efficent algorithm for doing matrix-vector multiplication on sparse matrices:\n",
"\n",
"**Exercise 6.** Use the following code as a starting point to implement matrix-vector multiplication (you only need to add code between \"##########\" )"
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "c2e8d88a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"your implementation is incorrect: multiply_sparse(L,y) ̸= Ly\n",
"matrix vector multiplication in 29995 mulitplications and additions\n"
]
}
],
"source": [
"function multiply_sparse( A::SparseMatrixCSC, x::Vector )\n",
" evals = 0\n",
" m,n = size(A)\n",
" b = zeros(m)\n",
" nzval, rowval, colptr = A.nzval, A.rowval, A.colptr\n",
"\n",
" # sum over the columns\n",
" for j ∈ 1:n\n",
"\n",
" # if x[j] = 0, we don't do anything because (u_j)_i x_j = 0\n",
" if (x[j] != 0)\n",
"\n",
" for k ∈ colptr[j]:(colptr[j+1] - 1)\n",
"\n",
" ##### COMPLETE THIS ######### \n",
" # get row index i \n",
" # get value v = A[i,j] by only using nzval \n",
" # add v * x[j] to b[i] \n",
" ############################\n",
"\n",
" evals += 1\n",
" end\n",
" end\n",
" end\n",
" return b, evals\n",
"end\n",
"\n",
"# don't change this\n",
"y = randn(n-1) \n",
"b, evals = multiply_sparse( L, y )\n",
"b == L*y ? println(\"success!\") : println(\"your implementation is incorrect: multiply_sparse(L,y) ̸= Ly\")\n",
"println(\"matrix vector multiplication in $evals mulitplications and additions\")"
]
},
{
"cell_type": "markdown",
"id": "e46292d1",
"metadata": {},
"source": [
"**Exercsie 7.** How many multiplication and additions are required to do the matrix multiplication $L x$ if $L$ is stored in a dense format? What is the percentage saving in the number of operations required?\n",
"\n",
"::: {.box}\n",
"\n",
"***Your answer here***\n",
"\n",
":::"
]
},
{
"cell_type": "markdown",
"id": "7eced612",
"metadata": {},
"source": [
"Here, we can see that sparse matrix-vector operations are also much faster:"
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "baf6a410",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" time taken\n",
"dense matrix-vector multiplication = 209.47 ms\n",
"sparse matrix-vector multiplication (using your function) = 0.1 ms\n",
"sparse matrix-vector multiplication (built-in function) = 0.23 ms \n",
"\n",
"• sparse matrix-vector multiplication is 908× faster!\n",
"• your code above (without any fancy optimisation)\n",
" is 2173× faster than the naive method!\n"
]
}
],
"source": [
"# | echo: false \n",
"s_dense = @elapsed L_dense * y\n",
"s = @elapsed multiply_sparse( L, y )\n",
"\n",
"# the built-in method uses an optimised version of your function\n",
"s′ = @elapsed L * y;\n",
"\n",
"println(\" time taken\")\n",
"println(\"dense matrix-vector multiplication = $(round( s_dense*1e3, digits=2)) ms\")\n",
"println(\"sparse matrix-vector multiplication (using your function) = $(round( s*1e3, digits=2)) ms\")\n",
"println(\"sparse matrix-vector multiplication (built-in function) = $(round( s′*1e3, digits=2)) ms \\n\")\n",
"println(\"• sparse matrix-vector multiplication is $(round(Int, s_dense/s′))× faster!\")\n",
"println(\"• your code above (without any fancy optimisation)\\n is $(round(Int, s_dense/s))× faster than the naive method!\" )"
]
},
{
"cell_type": "markdown",
"id": "9ba363f2",
"metadata": {},
"source": [
"Here, we numerically solve $-u''(x) = f(x)$ using the exact same code as above but now use the sparse format for $L$:"
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "b0a1ed4d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Solving the linear system of equations took 0.0127439 seconds\n",
"That is 875× faster than before!\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdZ0AU18IG4DMzS+8dAcUuYlcExAr2gsYWNcaYGFNMM9ck12sSjYnpxhg/S6ImNyaxa+wliqCIChak2EEEld57WWbmfD/m3s1emgjb931+7Q6HmbMjzrunzTCUUgIAAGCsWG1XAAAAQJsQhAAAYNS4lStXarsOAKr3+eefX758efDgwardbUVFxeeff37//v0BAwaods/1PXr0aM2aNcXFxd27d1fJDrOyslavXp2Xl9ejRw+V7FC1srKyjhw5cuTIkRMnTlhaWrZt21Zjhz548KCzs7OVlZX0tqCg4OTJk76+vhqrAGgZBTBEFhYW7u7uKt9tTk4OISQwMFDle64vKiqKEPLCCy887S+mpaXt3bv3xo0bdbbHxsYSQmbMmKGiCqrS2bNnFTlECPm///s/jR06ISGBELJp0ybFlsWLF5uYmJSVlWmsDpTSwsLCdevWff/993l5eZo8LlBK0TUKhqlLly6dOnXSdi20Iyoq6tlnn92zZ0+d7ebm5l26dPHw8NBKrZq2fPnyioqK7du35+fnFxYWvvrqqxo79ObNm728vObNmye9raqq+uOPPxYvXmxtba2xOpSXlw8fPjwoKKioqGju3LkaOy5IZNquAIBaSF/zQZmvr29SUpK2a9Gw27dvt2nTRvMZUFFRsWPHjl9++UURe3v27DEzM1u+fLkmq7Fly5ba2tq+ffuOGDHimWee0eShgSAIQVlpaakgCPb29gzDJCYmJicnm5ubDxkyxM7OTrlYWVkZz/NSMcVGQRBKS0tNTEwUF5Sqqqrq6morKytTU9NHjx7FxcVxHBcUFOTo6KgocOHCheLi4m7duvXu3bvBKpWUlMTExJSWltrb2wcFBSn3nikOYW1tbWJikp6eHhcXV1tbGxoaamJikpOTw7Ksi4tLnR1WVFRcuXIlPz/fxsbGx8enffv2yj8VRfHWrVuPHj2qrKx0c3Pz9/c3NzdvyakkpLa2Nj4+PiMjgxDi5OTUr1+/+i0MxadzcHAYNGhQnU9XX0VFhVwut7W15ThOsZFSWlxczHGcra0tIaS8vLyiokI6OUVFRVIZ6V+B5/n8/Hxzc3N7e/s6e5b+uVmW7d69u4+PT52fNvMPowkFBQWXL1+uqKhwcXEJDAxUPqulpaVyubywsLBDhw5ShRmGqV9DNdm1a5e/v//06dMVWzZv3rxu3TrpZGrMoUOHgoKCZDLZw4cPNfbZ4W/a7psFHdKnTx9CSFJS0siRIxV/Iba2tvv371cu5u/vTwjJz89X3njnzh1CyMSJExVbPvzwQ0LIH3/88eabb7Lsfzrhra2t9+zZQyk9cOCAk5OT4ijPPfdcbW2t8g6rq6vfeecdU1NT5ZooD+RQSv/xj38QQvbu3fvKK68oDiENsdQfI6ypqfnnP/9pYWGh/Pfv7++vKLB69eo6wenk5PTrr78q76SZY4T79u1r166d8q5MTEyUKy8IwvLly5XzwMrKavXq1co7qT9GOHPmTELItWvXlItVVVURQnr06CG9HTNmTP3/5jt27KCNjBHGx8f37dtXufDgwYNTUlKUyzTzD6NBcrl88eLFMtnf37kdHBy2bt2qKODn51ento6Ojk/craoEBQXdvXtX8TYhIWH06NEaO7qE53krK6uNGzdq+LiggBYh1DVz5kwTE5MtW7a4u7uHhYVt2LDhhRdeCAgI8PLyasHevvnmm/z8/DVr1nTq1Ck6Onr16tUvvPACIWTevHmvv/76iBEjCgsLP/300507dwYGBr799tvSb4miOG3atBMnTgwcOPDNN9/09va+c+fOF1988cYbb8hksldeeUX5EJ9++mlhYeHy5ct9fHwyMzPNzMzqV4NS+uyzzx4+fLhr165Llizp3r17QUFBfHz87t27FWWuXbs2ZMiQcePGeXt7U0qvXLmydu3aBQsWeHh4NJgujbl3796cOXNsbGzWrl3bv39/hmHS0tLCwsKkhprkn//855o1a7y9vVesWNG5c+fbt2+vXLnygw8+kMvl0heIFvvoo4+6deu2fv36adOmzZ49W9oYEBDQYOG0tLTg4OCioqKFCxfOmjWL5/lff/117969w4YNi4+Pd3Z2Vi6s/Idx+vTpjRs3NucPY+HChb///ruPj89HH33Url272NjYzz777JVXXqGUSv+OX331VUFBwezZs93d3f/v//6PEKL87Uet4uPjR40a1a1bN8WWf//73+vXr9fM0RVu3bpVUVEhfdsA7dBuDoNOkf4rDh48WC6XKzZK0xZ++OEHxZanahHa2Ng8fvxYsXHx4sWEEIZhlNsE0dHR5H8bZ7/99hshZNSoUYIgKDY+fvzYxsbG2dm5srJS2iK1CC0sLB48eFDns9RpEe7cuZMQ0qNHj+LiYuViyvuvrq6us5NLly4xDDNy5EjFlua0CNetW0cI2bBhQ2MF7t27x7KspaXlw4cPFRtv3LjBcZypqWlWVpa0pWUtQkrpH3/8QQj5+OOP6xy3fotwzpw5hJC3335buZh0lHfffVexpcE/DCnGlP8w6rt06RIhxNnZWXkm5Pnz5wkhdnZ2paWl0hae5wkhnTt3bmJXLVBTU/Ppp5+uWLFi0qRJyv0NixYtOn78OKU0Ozu7zj+68r+IBly+fHnmzJn+/v4Mw0yZMmXmzJl1/nFBMzBrFOpavny5iYmJ4u2UKVMIISkpKS3b24IFC5RbDOPGjSOEeHh4LFiwQLExMDDQzs5O+RC//PILIeSrr75SdHgSQry8vJ577rn8/PyLFy8qH2L+/PkdOnRouhrbtm0jhHz++ed1hrWU91+/KTlo0KAOHTpcunSJPs2dCKXe10ePHjVWYN++faIovvTSS8rdpz179pwxY4ZcLj9w4EDzj9Uacrn84MGDpqamddqgn3zyCSFEua0sqfOHIc3paPoPQ5q5+tZbbyk3LocOHRoSElJSUvLXX3+1+kM05fvvv1+wYMHixYuPHz9++fJlaWNeXt5PP/1UXl5OCHFzc6vzj16nQ7tBycnJZ5qnrKys6V35+/vv3bs3KCioS5cuhw4d2rt3rwbWp0J96BqFuuqstm7Tpg0hJDMzUyV7c3V1JYR069ZNOYEIIW5ubklJSTU1NdKF6erVqyzLXrly5fr168rFcnNzCSH3798fNWqUYmOdIa4GSY0hqS3bGErp3r179+/f//Dhw4yMDKmlJc0MKisra/7sifHjx9vY2Hz77bdRUVGTJ08eMWKEv7+/8udNTEwkhAQGBtb5xaCgoD179kg/1YDk5OTq6upu3bq5u7srb+/Ro4ednV12dnZOTo6bm5vyduVizfnDaOKTRkREJCYmSq1PdSgvLy8qKvLy8jpw4ADDMIr+z8jISErp0KFDW7xnGxsbBweH5pRsZh9vXFxcc/6GQX0QhFCXjY2N8ltpgqIoii3bW52pktLe6hyizlEqKyulEHrzzTcb3KfyYBv5b7g2rbi4WCaTSdfuxsydO3fXrl02NjbDhw/v06ePg4MDy7K///57VlZWbW3tEw+h4OXldfbs2ffffz8qKkrq9XVxcVm0aNGyZcuk2TGlpaXkv0GiTNpSXFzc/GO1RklJSYPVIIR4eHiUlJSUlJQoB2EL/jAa+6TSWsan+qQZGRnHjx9vTkk/P7/+/fvL5fJFixYRQvbt2xcYGKhokkZGRnbr1q3pv4Smubu71/nq0BqU0sTExA8++EBVO4QWQBDCU5NWTdTpLayurlbV/s3NzTmOMzExkdZjNLM+TbO2ti4pKSkoKKgzAUQhMjJy165dfn5+ERERylf8ffv2Nb/mCgMGDDh79mxubu7Zs2fDw8P37Nnz2WefPXz4UOqhtbS0JIQUFBTU+a38/HxS76tDfao689JqDemgDdbkics5nqiVn1RZ83unpZKOjo6Ojo4VFRVHjx794osvFD89d+6c8txXrUtLSysqKurXr5+2K2LUEITw1KQ1Bnl5ecqhosKV2tJqtps3byYmJqpqyKRXr14XLlxITEwMCQlpsEBMTAwhZN68ecopWFpa+vDhwxYf1NXVddasWbNmzfrkk0+6deu2a9eurVu3mpiYSAv1bty48eyzzyqXl24C0MSdRaW2b15envLG+mdeWqsgCELT1evcubNMJnvw4EF5eblyJmVmZubl5dnb27em2STx8fGJioq6ceNGcHCw8nbpk9ZfsNgELy+vFtxuJiIioqKiQrFEPT8//9atW9IgaIvl5OSkp6c/sRjLsr6+vg3OYVYWHx9Pmte9D+qDyTLw1KRF6OfOnVNsoZSqdtK5dIeRL774on47oGWdtNJCgq+++qqxeJCGABUr0CXffPPNE+Okvvo19PT0dHJyksvllZWVhJDJkycTQrZs2SJN2ZDk5ORs376dYZjQ0NDG9lz/zBNCpCUHyqQAy87ObrqeVlZWI0eOrKys3Lx5s/L2NWvWEEJCQ0PrjOO2gDTTav369cp9y/fv3z98+LCZmZk0c0qtbt686eTk5O3tLb09f/68NEBYU1OzadOmlu2zrKysqBkKCgpqamqeuLf4+HhXV1cV9rVCS2hptiroImmWfJ0FBtKX96lTpyq2hIWFEUKcnJx2796dlpYWGRk5fvx46eJbf/nEzp07lfcmff+dMmVKnUNLzSDFuoiKigrpO/KECRPCwsLS09PT09OjoqI+/fRTb29vxVR4afnEwYMH63+WOssnqqurpcbl2LFjo6KiCgoKUlJS/vzzT8Xnio+PZxjG1tZ23759xcXFjx8/XrlypZmZmXQfHMVakeYsn/jggw9mzpx55MiR5OTksrKy+/fvL1myhBAyZMgQRZlp06YRQoYOHXrlypXCwsLz589Ln3fhwoWKMvWXT9y9e5fjOEtLy59++unBgwcxMTEvvPCCdOaVl0/k5OTIZDJbW9vvvvtu7969e/fulVYF1F8+cfXqVRMTEzMzs/Xr12dmZj5+/Pizzz5jWdbKyio5OVlRrJl/GPUJgjB8+HDp3zE+Pr6wsPDUqVNdunQhhCxbtkxRTE3LJyilX3/9da9evRRvp0+f3q5dO0rpiRMnrly5ovLDtcDkyZPHjRun7VoYOwQh/K351ztpOaBCz549T58+rcIgpJTm5+dL7Yk6fH19eZ6XyjQ/CKUd1m+CuLq6KgqsXLlSebjR3Nz8t99+q7NosjlB2OBtKvv165eWlqYoU1ZWNnXqVOUCDMPMnz9feVlbg0+fWLNmjXJDrW3btlKnrnIQUko3bNig3NvZxJ1lDh06pHyLH0KIh4fHuXPnlMu0OAgppQUFBXXG5FiWXbx4seIfkaozCB88eODq6hoXF5efn79q1arVq1e3adMmOzt7yZIloiiq/HAt0K5du08++UTbtTB2DH2aBVJg2GJiYsrLy0eMGKF8Q6zy8vKYmBhXV9c6twONiYmJiIioqanp1atXaGgoz/PR0dHOzs6K0Y779++npaX17NlTudtH2puLi0ud+2hER0dXVFSEhITU6Y67c+fO2bNns7KyLC0tPT09AwMDu3btqvhpUlLSo0eP+vTpU/+eolFRUTKZbNCgQXW2x8XFRUZG5uTkuLq6dunSZdSoUcr3OUtISLh48WJWVlbbtm0nTpzo6el55cqV0tLS4cOHS9N2amtrL168aGdn1/TshoyMjKioqIyMjKKiIg8Pjx49egwdOrR+T+P169fPnTuXl5fn7u4+cuTInj17Kv+0pKQkLi7O3d29zljazZs3T548WVpa2rVr12eeecbKyioiIsLKyqrOh62pqXn48KE05VX6VygvL7927Zqrq2udJ+2VlJQcP348KSmJYZgePXpMmDBBmuSi8FR/GA26dOnSxYsXpbMxduxYqVGoQCkNDw+3tLQMCgp64q6eVnp6+uHDhzmOmzZtmqur64ULF65du/bcc881Z7KxuhUWFjo5OUVGRg4bNkzbdTFqCEIAAE3Lz893dnY+derUSy+99OjRI+VvGKB5mCwDAKBRP/30k4uLS3R09NmzZ1966SWkoNbhHwAAQKPy8vImTJjQpk2bmJiYo0ePars6gK5RAADNksvlv/32W1VV1fPPP694PCdoEYIQAACMGsYIAQDAqCEIAQDAqCEIAQDAqCEIAQDAqCEIAQDAqOl3EH7xxRfSXQqbo/klQSUopS14dAO0hiAImAeuYbiwaJ7Kz7l+B+GaNWvKysqaWVh66DlojCiKcrlc27UwLnK5vGWPqYIWw4VF81R+zvU7CAEAAFoJQQgAAEYNQQgAAEYNQQgAAEZNNUG4ceNGFxcXGxub6dOnNzh75dq1a7169bKwsOjTp8/169eljZTSpUuX2tvb29nZvfnmm4oZhnfu3PH397ewsOjatWtkZKRKaggAANAgFQRhQkLCxx9/fPbs2by8PLlc/sknn9QpIIri7NmzFy1aVFlZ+eqrr86aNUua2LZ///79+/ffvXs3LS3twoULP//8s1R+/vz5kyZNqqio+Oyzz2bOnFlTU9P6SgIAADRIBUG4bdu2adOm9ezZ09zcfOnSpb/99ludlUyRkZFlZWWvv/46wzCLFi0qLS29cOGC9Iuvvfaau7u7g4PD4sWLt23bRgi5efPmrVu33n//fZZlZ8+e7ejoeOzYsdZXEgAAoEEqCMKkpKSePXtKr3v27FlYWFhQUKBcIDk52dfXl2VZQgjLsj4+PsnJyXV+sUePHtLG5OTkjh07Wlpa1tneIEppcXFx0X9h+TYAADwtFTyhvri42NraWnptY2NDCCksLHR2dlYuYGVlpXhra2tbWFhICCkqKlL+xaKiIinYGizcoPLy8n79+jEMI70dN27cli1bGiwpF8mgkyYHRpR7W+G+GxoiCIJcLsd9NzSpqqrK1NSU4zhtV8SIlJeXa7sKxiW7ivk8gWwIbO6tVMzNzU1MTJouo4IgdHJyKikpkV4XFxcTQlxcXOoUKC0tVbwtLi6WCjg7Oyv/orOzM8Mw9Qv7+vo2dmgbG5uUlBQHB4fm1PO1blXPXTK/GCqzecI5AdWQgtDCwkLbFTEiMpkMQah5UgMANONknljI19rYmKlwnyroGu3evXtCQoL0OiEhwc3NrU4y+fj43Lp1S2oZ8Dx/69YtHx+fOr+YmJgobfTx8Xnw4IFi6mlCQoK0vfVe7cwHuTLzzgki2oQAAPrpah7t76ji+wiqIAgXLFhw+PBhadboqlWrXnnlFWn70qVLDx48SAgJCgpq06bNV199VVJS8uWXX3p5eQUGBhJCFi5c+OOPP96+fTs1NfX7779fuHAhIaRr166BgYErVqwoKSlZv359bW3t+PHjW19JyfogrlhOP76GoUQAAL10LZ8OcFJxa0YFQditW7etW7e+/fbbffr08fHx+eijj6TtJSUl1dXVhBCGYf78889z58517dr1/Pnz+/fvlwpMnDhxyZIloaGhwcHBc+fOff7556Xtv//+e1JSUteuXffu3Xv06FGZTAX9txITluwfKdvzgO5KwY2JAQD0jEjJ9Xzaz0HFF3BGrx/a4ujo2PwxwrKyMqkr/0YhHXWSPzZGNtCFUXMFjRrGCDUPk2U0T3FhAQ24XUynnBauj69S7Tk3xlus9XJktgzhpp8RMiv1+EsAAICxuZpH1dGAMcYgJIRM8WZf785OOyNUY7gQAEBPXMmj/ghCFVrWl+1gw7wShSQEANAPaBGqGEPIv4dyd4vpNwmYOAMAoOtqBHK7iPZzQhCqlIWMHBrNrb8tHnuEwUIAAJ2WUEi72DGWKltG8DejDkJCiKcV8+dI7uUo/lYRshAAQHepaYCQIAgJIQGuzJoAbkqYUIDHPQEA6KoruQhCdXq+Mzu9PfNsOF+L4UIAAJ10NR9BqGZfDeQsOPJuDCaRAgDonGI5yaigvg4IQnViGbIzRHYuk/50B61CAADdcjWP9ndmOPXcDQxB+DdbE3JkDPfpdeFcFibOAADoEPXNlCEIwjo62TI7gmVzIvgHZchCAABdcSWPBiAINSbEg/m4HzfltFBWq+2qAAAAIYSQK7mivyuCUIPe9GUHuzHP4xG+AAA6IK2MsgzT1gpBqFnrg7gSPMIXAEAHXMmjAWprDhIEYWMUj/DdiUf4AgBo1WV1zpQhCMImOJuTw2O4f8QIV/PQQwoAoDWXc9Ei1J6eDszWIdw0PMIXAEBLakWSUEgHOiMItWeyN7uoOzs1DI/wBQDQgoRC2tGGsTZR4yEQhE+2rC/byZZZiEf4AgBonLr7RQmCsDkYQn4Zyt3DI3wBADQuJleNS+klCMJmkR7huwGP8AUA0KzLeTQQLUId4WnF7McjfAEANKighuRV0e72CEKdgUf4AgBoUkwuHejCsOrNQQThU3q+MzujAzPzDB7hCwCgdjG5orr7RQmCsAW+9OOsTMjiaEwiBQBQr5hcGuiq9pxCED41liE7gmXns+mPeIQvAIDaiJRcVfNdRiUIwpawNSGHR3OfXRfO4hG+AADqcauIulkwTmZqPxCCsIWkR/g+h0f4AgCoR0wuHaT+5iBBELZGiAezvB83GY/wBQBQg+hcOsgNQajz3vBlh7ozc8/iEb4AACoWnav2pfQSBGFr/d8grqyWfoRH+AIAqE5hDcmsoD0dEIT6wIQl+0bK9j6gO+5jEikAgGpczqUDXRhOEzmIIFQF6RG+Sy4LV/AIXwAAVbiUK2pmgJAgCFVFeoTvdDzCFwBAFaJzaJD6l9JLEIQqM9mbfcOXfSZMqOK1XRUAAH0maGopvQRBqEr/6sN2sWUWRgloFQIAtNiNQuplxTiqfym9BEGoSgwhPw/lkkvxCF8AgJa7lKOhFYQSBKGKSY/w3XhbPPoIWQgA0BKXcmkQglCveVgyf47iFkYJN/EIXwCAp3cphwZpaoCQIAjVxN+F+T6Am3Iaj/AFAHg6WZWkrJZ2U/NT6ZUhCNVlbmf22Y7MDDzCFwDgaVzKEQe5qvuh9P8DQahGX/hxNibMO3iELwBAs13M0egAIUEQqhXLkB3BXFQ23XQbrUIAgGbR8EwZgiBUNxsTcmQMtyoOj/AFAHiyKp7cKqIDnRGEhqWjDbMzRPZcBJ9SiiwEAGjKlTza04GxkGn0oAhCTQhuw6zoz00+LZTiEb4AAI27lEsHa7ZflCAINWZRd3aEBzP3LI9H+AIANOZCtoggNGQ/BHIVteRDPMIXAKAhIiXRuXSwm6aDCUGoOSYs2TdKtj+VbscjfAEA6rldTF3MGVcLTR9XNUGYmJg4YsQIDw+P0NDQx48f1y9QXl6+YMECDw8PPz+/kydPKrafOXNm4MCBHh4e8+fPLykpkTYuW7bM779CQkJUUkMd4WRGDo/m3rssXM5FDykAwP+4kK2FAUKikiAUBGHKlCmhoaE3b97s3Lnz3Llz65f517/+lZOTEx8fv2LFilmzZmVlZRFC8vPzp0+fvnTp0oSEhLKysvfff18qnJqaOnr06M2bN2/evHnt2rWtr6FO6eHA/DJUNiNcyKhAFgIA/O1CDh3iroUgJLTVTpw44eXlJYoipbS8vNzc3PzOnTvKBaqrq21sbK5duya9nTBhwtdff00p/eGHH0JCQqSNt27dsrS0LC8vp5TOmjVr06ZNzTm0g4NDYWFhM+tZWlrazJIa8HW84HewtrJW2/VQJ57nKysrtV0L41JZWcnzvLZrYVx06sKi79rvrr1XLD6xmMrPuQpahHfu3OnXrx/DMIQQKyurLl263L59W7lAenp6RUVF3759pbf9+vW7c+cOIeT27dv9+/eXNvr6+gqCkJaWJr1dt25d//79Z82aFRcX1/oa6qClfVgfe+ZlPMIXAIAQQsjjClrF0y52WmgRNnfV4p49e+pvHDBgQOfOnQsKCmxsbBQb7e3t8/PzlYvl5+dbWlpyHKcoIMVbQUGBp6enopidnZ30i3Pnzn3rrbdsbW3//PPPYcOGJSQkdOzYscFalZWVdejQQcpgQsj48eO3bNnS2EeoqKhQlNQFa/uR8REmn16pfd+X13Zd1EIQBLlcLgiYJas5VVVVpqamiv9roAG6dmHRX2cesgFOXEV5+RNLPtU5Nzc3l8mekHTNDcIjR47U3+jk5NS5c2cHB4ebN28qNpaUlDg6OioXc3R0rKysFEWRZVmpgJOTEyHEwcGhXOkzl5aWSr8YGhoqbendu/fFixf379//z3/+s8FaWVtbX79+3d7eXnpra2vbxCWAUmptbd28j6sJ1oQcGUsCDvP93EyneBvg9F0pCC0sND4DzIhxHIcg1DBdu7Dor2slwggvxtra7IklVX7OmxuEO3bsaOxHnTt33rx5s/S6pqbmwYMHnTt3Vi7g6elpamp69+5dX19fQsidO3d69+4t/WJUVJRUJjU1VRRFb2/vOju3tbWtrq5u7NAMw9jb2zs4ODTzU+iaNpbkwGhu4im+ow3TyxFfKgHAeEVl05e6aqdJoIKjjh8/vrS0dPfu3ZTSdevWde7cWRoOPHTo0MaNGwkhlpaWM2fO/Prrr3mej42NPXXq1PPPP08ImTt3bmRkZExMjCAIX3311dSpU21tbQVBOHr0aFVVlSAIhw4dOnny5Lhx41pfSZ3l58z8EMhNCRPyG417AAADV1RDHpbRvlpqD6ggCM3MzPbu3bt8+XJra+s9e/Zs375d2p6SkhIfHy+9/u6773Jzcx0dHSdOnLhhwwZpzM/Ly2vz5s3Tpk2zt7d/8OCBtFKCUrpq1SoHBwdbW9uPPvpo27Zt/v7+ra+kLpvTiZ3TiZmOR/gCgLG6mEMDXBmZlsaIGEpVNm+RUtr0AGZjBRrc/sS9EUIcHR1TUlKa2TVaVlamPKlHp4iUTDsjuFmQzUMMZ3QHY4Sah8kymqfLFxY9svSKYG3CLO/XrCRU+TlXZf4+MbcaK9DgdqOaiMUyZPsILjqXbsAjfAHA+ERl06FaWUpPCMG9RnWHtQk5PJr7Ml4Iz8TaQgAwIpU8uVFEA1wQhPeX/I8AACAASURBVEBIBxtmV7Bs7ln+Ph7hCwBGIyaX9nHU9MN4lSEIdcvwNsxnA7jJp4USubarAgCgEdrtFyUIQh30qg870oN57iyP268BgDE4ny0Oc9dmGCEIddHaQE4ukn9dwc3JAMDAyUVyNY8ORosQ6pCxZE+I7PAj+lsyJpECgCG7lke72jG2JtqsA4JQRzmakcOjuX9eEaLxCF8AMFyR2XSYVpuDBEGoy7rbM9uGyWaGC4/xCF8AMFDns8ThbRCE0LjxbZl/9GSfCRMqDfNJTQBg1HiRROfSIVqdKUMQhLrvvV5sLwfmxUjMIQUAQxNXQL2tGacnP3lJvRCEemDzEC69gq6Kw8QZADAoujBASBCEesGMIwdGy365J/6ZiiwEAMMRqQMDhARBqC/cLcih0dyii0J8AbpIAcAQCJRczKHD22g/hrRfA2imfk7MxsHcM2FCbpW2qwIA0GrxBdTDknEx13Y9EIT6ZWYH9sWuzLQzfA3uOQMAeu5cFh2hA/2iBEGodz7pz7lbMosuIgkBQL9FZlFdGCAkCEK9wxDy23AuroCuvYmJMwCgrwRKLuSIw7W9glCiE5WAp2IlI4dHc98lin+lY+IMAOglaYDQ1ULb9SCEIAj1VDtrZu9Ibn4kf7cYWQgA+ueszgwQEgSh/hrsxnwzkJscJhTVaLsqAABP6VymGIwghNZ7sSsb2o6ZFcHzGC4EAP3Bi7qyglCiK/WAlvnWn5OxZMllTCIFAL1xvYC2s2acdWAFoQRBqN84huwKloVl0K130SoEAP1wNpMGe+hKvyhBEBoAO1NyZDS3PFY4n42JMwCgByIyxRCdGSAkCELD0MWO2T5CNjuCTy1DFgKATpOLJCaXDtOZAUKCIDQYozyZD/tyk08LZbXargoAQOMu59Ju9oy9qbbroQRBaDje8mUHuzHPnxNENAsBQFdFZFLdWTghQRAalPVBXKmcfnwNk0gBQEedzRJDPHQrenSrNtBKJizZN1K25wHdcR+TSAFA51TyJDafDtGBp9IrQxAaGmdzcmQMt+SycCUPPaQAoFsu5NB+ToyVTNv1+F8IQgPUw4H5eSg3/YyQUYEsBAAdEpEpjtSxflGCIDRUoe3Yt3zZZ8KEKl7bVQEA+K/wDBqiS0vpJQhCg7W0D+tjz7x0XkCrEAB0QWENSSqhAa4IQtCgrUO5h+X0izhMnAEA7TuXJQ52Z0x1L3Z0r0agOuYcOThatvWeeCANWQgAWnYmg47SvQFCgiA0eO4W5OAo7vULQkIhukgBQJvCM+lIT53rFyUIQmPQ35nZOJibclrIrdJ2VQDAWD0qp8Vy2tsRQQhaMrMDO78LM+0MX4N7zgCANpzJpCM9WF2MQQSh8Vg5gGtjybx+EUkIAFpwJoOO0sl+UYIgNB4MIduGcwkFdM0NTJwBAI2ihERkiqN0bwWhBEFoRKxk5PAY7vsb4onHmDgDAJqTWEjtTZl21ghC0AFtrZj9o7iXzvO3i5GFAKAhYTrcL0oQhEZokCvzXQA35bRQUKPtqgCAcQhLF0cjCEGnzOvMTmvPzDzD12K4EADUrFog0bk0WCeX0kt0t2agVl8N5KxNmHeiMYkUANQrKpv2cWJsTbRdj8YhCI0Uy5AdwdyFbLrxNlqFAKBGYRniaE+dzhqdrhyolY0JOTKG+yJeOJOBiTMAoC6n0+kYHR4gJAhCI9fBhtkdInv+HJ9cgiwEANXLriKPK+hAFwQh6LBh7sznflzoaaFYru2qAIDBCcsQQzxYTqdzUEVBePbs2bFjxwYGBn7xxReC0MD8i8zMzBdeeMHPz+/FF1/Mzs6WNqampq5du/bFF19cunSpcuHCwsLXXnvNz89v9uzZqampKqkhNGFhN3ZcW2ZWOM9juBAAVOqUzveLEpUE4aNHj6ZMmfL8889v2rRp375933//ff0yM2bMsLKy+u2338zNzWfOnCltvH79emJiYm1tbVhYmHLhF198saysbNu2be3bt584caIo4vKsdmsCOIYh713GJFIAUBlKSFiGONZL14OQ0FZbsWLF9OnTpdfHjh1r3759nQKxsbHW1tbV1dWU0qqqKmtr6/j4eMVPd+7c2a9fP8XbtLQ0U1PT/Px8Sqkoim3atAkLC2vs0A4ODoWFhc2sZ2lpaTNLGqeiGuqzr3bLHUFVO+R5vrKyUlV7g+aorKzkeV7btTAuuLA0ITZP9NlXq/Ldqvycq6BFmJiY6O/vL70OCAhIS0srKSmpU6BPnz5mZmaEEHNz8969eyckJDS2t5s3b3bo0MHJyYkQwjDMwIEDmygMKmRvSo6M5pbHCpFZmDgDACpwKoPqQXOQEFlzCpWUlDSYRn5+fpaWlnl5efb29tIWBwcHQkhOTo6dnZ2iWG5urqKAVCY3N7exYz1V4fLy8r59+7Lsf+J87Nix3333XWOFKyoqGEYP/km0qA1Hfg5kZ4XTM6Nr21u1Ng4FQZDL5Q2OGYOaVFVVmZqachyn7YoYEVxYmnA8zeQ9X6G8XMXDW091zs3NzWWyJyRds4IwNTV15cqV9bf/+uuv3t7etra2lZWV0pby8nJCiHKSEULs7OwUBaQyyjFZR/3CXbt2baywpaXloUOHFHuzt7e3trZurDCltImfgmSSNVlRI86KYqOnyFp5JwgpCC0sLFRUNXgyjuMQhBqGC0tjSuQksbh2XAczi2blzFNQ+TlvVgX79u0bERHR2E/bt2+fnJwsvU5OTra2tnZ2dlYu4O3tff/+fUopwzCU0pSUlPbt2ze2N29v74cPH8rlclNTU0LI/fv3p02b1lhhlmXbt28vNUNBVd7wZW8V0zkR/JExMh2f9AwAOis8UwxyY1SeguqggjHCuXPn7t+/PysrixCyYcOGOXPmSH2Vv/zyy6VLlwghI0eO5Hn+0KFDhJCDBw+KohgcHNzY3vr37+/p6blt2zZCyPnz51NTU0NDQ1tfSXgq6wK5GoEsvYJeTQBoob/S6TgvPVmqrpIpNx9++KG9vb2Xl1dAQEB2dra0MTg4+IcffpBe//XXX25ubp06dXJ3d1fMAj1//ryDkhkzZkjbL1261LZt206dOjk7O+/fv7+J42LWqPoUVtOue2v/fa/lk0gxa1TzMGtU83BhaUzbnbV3ikR17Fnl55yhVDVTBCsqKsrKytzd3RsrwPN8Tk6Om5vbE8ctCSGCIGRnZ7u6upqYNDVO5ejomJKS0syu0bKyMhsbm+aUBMm9Ejr8GL9vpGyoe0t6SDFGqHmYLKN5uLA06FYRDT0tPJillo5RlZ9zlbVbraysmkhBQohMJvP09GxOChJCOI7z9PRsOgVB3brZMb+PkM2K4NPKsKACAJ7CyXQ6Th8WTkj0pAMXtGSMJ7OsDzc5TCir1XZVAEB/nHwsjm+LIARD8XYPdrAb89xZXkSzEACaoayWXMujITr8SPo69KaioEX/N4ir5Mm/rmISKQA82ZkMMdCVsdKHhRMSBCE8mQlL9o2UHUyj25JwA3QAeIITj+mEtvoULvpUV9AiRzNydAy39KpwIRs9pADQKErIyXQ6QX8GCAmCEJrPx575Y4TsWUwiBYDGJRRQSxnpYocgBAM1xpP5FyaRAkDjTjymE/WqOUgQhPC03unBBrliEikANOz4Y1G/BggJghBaYH0QV8njTqQAUFd+NblVRIe16F5UWoQghKcmTSI9/Ij+ikmkAKDkr3QxxIM107d7/CEIoSUczciR0dyyq0IUJpECwH8de0QntdOz5iBBEEKLSZNIZ0XwqZhECgCE8CIJy9C/AUKCIITWGO3JfNSXCz0tlGISKYDRi8qhnWwZdz183gyCEFrlTV92eBtmTgQvoFkIYNyOPRJD2+llpuhlpUGnrAvk5CL54DImkQIYtWOPaKgeDhASBCG0nowle0NkJx7Tn+9hEimAkbpXQqt40scJQQjGysGMHB3DfXxNOJeFHlIAY3TkIZ3UjtHLGEQQgqp0sWN2BsvmRPD3S5GFAEbn6CNxsre+Boq+1ht0UIgH8+kALvS0UCzXdlUAQIPyq8mNQhrcRk8bhAhCUKlXfdhxXsyz4TyP4UIAo3H8sTjSU/9uKKOAIAQV+y6Ak7FkcQwmkQIYiyMP6WT9nC8qQRCCinEM2R0iO59FN97BYCGA4asWSHimOFE/VxBK9LjqoLNsTciRMdzXifRMFv7AAAzcmQzaz4lxMtN2PVoB1ylQiw42zJ5g9uVLzO1itAsBDNnhh+IUvZ0vKtHv2oMuC3IlX/YTJ58W8qu1XRUAUA+RkqOPxCneejxASBCEoFZzO9JnOzDTzvByTCIFMESXcqi7BdPBBkEI0LjP/ThXC+bVKEwiBTBAhx7qfXOQIAhB3ViG/D6cu1lEv05AqxDA0Bx6SKe21/sc0fsPALrPUkYOj+Y23RYPpiELAQxHYiEVKemrnzfaVoYgBE3wtGIOjuZeuyBcz8ckUgADceghfUb/+0UJghA0ZoAzs3kI90yYkFGBLAQwBAdSxWn63y9KEISgSVPbs2/4slPChEpe21UBgNZJKaU5VTTIDS1CgKf0rz5sL0dm3jlBRLMQQJ8dfEineLOsIeQgghA0bvMQLr+afngNCyoA9NiBVHF6BwNJEAP5GKBHTFlyYLTsQBr9NQmTSAH0UkYFTSqhI/T2AYR1IAhBC5zMyNEx3LKrQmQWekgB9M+BNBrqzZoYSoAYyucAfdPNjtkRLJsdwSeXIAsB9MyfaeJ0g5gvKjGcTwJ6Z6QHs8qPm3RaKKzRdlUAoNlyqkhiIR3taSD9ogRBCNq1sBs7uR0z4wxfi+FCAD1xIE2c2JY147RdD9VBEIKWfePP2Zoyr1/AJFIA/bA/VZzRwXCagwRBCFrHMmRHMBdfSL/BXbkBdF5uFbmeT8d6GVR2GNSHAT1lJSNHRnMbb4t/piILAXTagTRxQlvW3ID6RQmCEHSEpxVzeAy36KJwNQ+TSAF0175UcWZHg+oXJQhC0B39nJifh3JTzwiPypGFALoop4rEFdBxhtUvShCEoFMme7Pv9WJDTwtltdquCgDU82eqONHg+kUJghB0zT96skFuzKwIXkCzEEDH7E0VnzW4flGCIAQdtH4QJ1KyOBoLKgB0SGYlvVFIx3gaYGoY4EcCfSdjyd6Rssgsuu4mJpEC6Ip9D+gUb4NaR6+AIARdZGtCjo3lvk0Ujz1CDymATtj9QJzV0TAjQ6aSveTn53/55ZdJSUkDBgxYunSppaVlnQKiKG7atOnUqVPu7u4ffPBB165dCSFVVVXh4eGxsbGZmZk//PCDhYWFVHjr1q3Xrl2TXltZWX3//fcqqSToF29r5uBoLvQ0/9c4WT8nAxyWANAjaWX0QRkd6WGY/xNVE++hoaEFBQXvvvvu1atXFy5cWL/At99+u3nz5jfffNPT03PYsGFlZWWEkIcPH3799ddpaWlbtmyRy+WKwuHh4eXl5QMGDBgwYECfPn1UUkPQR/4uzKYgbsppIaMC7UIAbdr9gE5vz8oMs0FICG21S5cuOTg4yOVySmlOTo6pqWl6erpyAZ7nPTw8wsPDpbeDBw/evHmz4qcZGRmEkOLiYsWWWbNmbdq0qTmHdnBwKCwsbGY9S0tLm1kSVILn+crKytbv59sEoe+B2jJ56/dk+CorK3me13YtjIuRXFj6/Fl7PkvUdi3+Q+XnXAX5Hhsb6+/vb2JiQghxdXXt1KlTXFyccoGMjIzMzMzBgwdLbwcPHnz16tWm93ngwIGXX375m2++KSoqan0NQa990Jsd6MLMxoIKAC25XUwLa8hgN8PsFyXNHyNMTk6uv9HV1dXOzi47O9vR0VGx0cnJKTs7W7lYTk6OlZWVmZmZosDt27ebOFZQUBDDMDY2NgcOHNi0aVN8fLyDg0ODJSsqKiZPnixlMCFk4MCBH3/8cWO7LS8vb+KgoHKCIMjlcp7nW7+rr3uRmedN3ois/W6ACvZmwKqqqkxNTTnOECf26SpjuLD8dkc2rS2pKK/SdkX+46nOubm5uSIjGtPcIJw8eXL9jStWrJgzZ46VlVVNzd9PVq2srLSyslIuJhWglDIMQwipqqqytrZu4ljvvPOO9OKFF14YOHDgjh073nrrrQZLmpubv/fee4q9dejQwcbGpok9N/1TUC0pCBVzoFrp4Dgy+Aj/74dmi3sa6jCFCshkMgSh5hn2hYUSsv8R/+cozsZGh1qEqj3nzQ3CO3fuNPajtm3b7tu3T3otiuKjR4/atWunXMDT01MUxfT09LZt2xJC0tLSpBdPxLJsly5dcnJyGivAcdzw4cMbay+CIZEWVAw+KrS3IVO8kYUAGhKTS8050tegZ26r4IIyadKkpKSk2NhYQsjhw4ctLS0DAwMJIVeuXPnrr78IIXZ2dmPHjv3ll18IITk5OceOHXv22Wcb25sUpdLr5OTksLAwaW8A3tbModHcqxeEa/kYLQTQkB33xbmdDfyrpwrWEdrb269du3bMmDHdu3e/d+/etm3bpJ6Zw4cPJyYmjhs3jhDy7bffTpgw4dSpU6mpqfPmzfPz8yOECILg4uIiiiIhpH379hYWFpmZmYIg9OjRo23bthYWFsnJyW+99dbEiRNbX0kwDH7OzJYh3DNhwqVQrp21IX9FBdAFtSLZlypenqyaFec6i6FUNV+uCwoK0tLSunXrphixq66uFgRBMV4ol8tv377t7Ozs5eWl+C3lSaEMw9jb20slk5OTBUFo3769ra1tEwd1dHRMSUlpZtdoWVmZYXfl6xrVjhEqW3dT/PmeeCFUZmeq8n3rN0yW0TzDvrAce0S/SRSiJulWEKr8nKvs4zk5OTk5OSlvMTc3V35ramrat2/fOr/VYIaZmpr26NFDVRUDw7O4J5tSRmeE8yfGykwMvM8GQJu23xefN/R+UYJ7jYKeWhvIWXDM6xfwhAoAdSmRk1Pp4swOhh8Thv8JwSBxDNkVwiUW0i/i8YQKALX4M00M9mAdzbRdD/VDEIK+spKRI2NkP98Td9xHFgKo3u/J4rzORjElDUEIeqyNJTk2hltyWYjMwoIKAFVKK6O3i+jEdkaREUbxIcGA9XBgdgXLZkXwd4qRhQAq88d9+mxH1tQ4IsI4PiUYtBAP5lt/buIpIUdXboUIoN8oIb8niy90MZaAMJbPCYbthS7si13Z0NN8JW7KDdBqF7OpKUv8XYxigJAgCMFgrOjH9nRg5pwV8LQmgFb6LVmc39WI0sGIPioYvM1DuCqeLo7G4kKAlqvkyYE08XnjmC8qQRCC4TBhyf5Rsqhs+t0NLKgAaKEDaWKgK+NhiSAE0E+2JuT4WG79LXHvA2QhQEv8miS+ZEz9ogRBCIbHy4o5NoZ7O1qIysZoIcDTSS2jN4voZCN75KdxfVowEr0cmZ3BspnhWFwI8HR+TRKf62QsywcVjOzjgtEY6cGsDuAmnhKysbgQoHkESrYl0QVG1i9KEIRgwOZ1Zl/uxk48xZfVarsqAPrgdDptY0l6ORrRNBkJghAM2Ud9WT9n5tlwvhZTZwCe5Jck8eVuxhgKxviZwahsDOJkLHn9AtbZAzQlp4pEZIqzOxljKBjjZwajImPJ7hDZzSK6MhYL7QEa9VuyONWbtTXRdj20AUEIhs9KRo6Oke1MoVvvoocUoAGUkJ/via/4GGkiyLRdAQBNcLUgJ8dxw47xbSyZSe2Mbi4AQNPOZlILjgS6Gul/DSPNfzBCnW2ZQ6NlL0fxl3MxXAjwP7bcFV811uYgQRCCUfF3Yf49TDb1DJ9cgiwE+I/cKnI6Q3y+s/HGgfF+cjBOE9syqwZw4/EUX4D/+jVJnOrN2plqux7agyAEo/NyN3Z+Fyy0ByCEEJGSLXfF17sbdRYY9YcHo7W8H+vnzMw4g4X2YOxOZ1BHMzLQaB5G3yAEIRipjYM5Sxmz4DwW2oNR+/GOuMjX2IPA2D8/GC2OITuDudQy+s/LWGgPRuphOb2UI87uaOxBYOyfH4yZhYwcHSM7mU6/xxPtwShtviPO68xaGv16cgQhGDUHM/LXOG7dLXFnCrIQjEuNQH5JEhcZ9zQZCU4BGDsvK+bkOG5JjHA6A8OFYET2PBD7OzFd7Ix6mowEQQhAfO2ZA6Nk887xV/OQhWAsNtwW3+rBabsWOgFBCEAIIUFuzC9DZVPC+CTcdAaMQEwuLaoh473QHCQEQQigMKkd86UfN+4vIbMSWQgGbv0t8U1flkUOEkIQhADKXuzKvubDjvtLKKrRdlUA1Cazkv6VLr7UFdf//8CJAPgfS/uwoz2Z0NN8Fa/tqgCox6bb4tzORn1z0ToQhAB1fRfAdbJlno3geSypAINTxZOt98S3jf5uMspwLgDqYgj5ZShHCMEN2MDwbL8vBriwWDWhDEEI0AAZS/aGyFLL6XsxuAEbGA5KyA83xX/0wpX/f+B0ADRMugFbeCb9Ih49pGAgTqVTU44Et0Fz8H8gCAEaZW9KTo2XbUsSf7qDLARDsOaGsATNwXqM/marAE1ytyCnxnPDjwkOZmSW0d+kH/RaQiG9W0zwrIn6EIQAT9DRhjkxlht9krczZcbhThygt75LFN/uwZogB+vBKQF4sl6OzKHRsvmR/MUcTCMFvfSonJ58LL7mg2t+A3BSAJol0JXZPkI27QwfX4AsBP3zw03xpa5YRN8wBCFAc432ZH4czE08JeDG3KBfCmvIb8niuz1xwW8YxggBnsK09mypnIw5KZyfxLWzxngh6IdNt8VnvFlPK/zFNgxBCPB0XuzKlsjJ6JPC+UkyNwtt1wbgSSp5suG2cG4irvaNQksZ4Kkt7snO7cyOOcnjIRWg+36+Jw5xZ33s0RxsFIIQoCVW9GNHezLjT/FltdquCkDj5CJZc0Nc1geX+qaoprH88OHDTZs2lZSUTJ48ecKECfUL1NTUbNq06c6dO76+vosWLTIzMyOE3Lt37/Dhww8ePHB1dZ0/f36nTp2kwjzPb926NS4urmPHjm+//baVlZVKKgmgWqsDuEUXhMmn+RNjZRbodgKd9Eey2N2eDHBGc7ApKviaUFxcHBgYWFtb6+fnt2DBgj179tQvM3/+/KNHjw4bNuzw4cMLFiyQNr7xxhvp6el+fn6lpaV9+vS5deuWtP3dd9/dtm3bkCFDLl68OHXq1NbXEEAdGEI2DeY8rZgZ4bwct2AD3SNQ8nWC+FFfTtsV0Xm01X744YeQkBDp9R9//NGvX786BVJSUszMzAoKCiilBQUFZmZmqamplNLa2lpFmUmTJq1YsYJSmp+fb25unpKSQimtqqqys7O7du1aY4d2cHAoLCxsZj1LS0ub/6Gg9Xier6ys1HYt1K5WoFPD+Oln+FpB21WhtLKykud5bdfCuOjyhWV7sjD8WO2Ty+kblZ9zFbQIL168GBISIr0OCQmJi4urrKxULhAdHd27d29HR0dCiKOjY69evaKjowkhMtnf3UnFxcUODg6EkNjY2DZt2nTs2JEQYm5uPmjQoIsXL7a+kgBqImPJ7hCuopa+dF4QsbwQdIZIyRfx4vJ+aA4+WbNGNgRBqKqqqr/d0tKSZdns7OyRI0dKW1xcXAghWVlZigE/Qkh2drazs7Piraura1ZWlvJ+tm/fnpqa+tJLL9Uv7OLiUqewsqqqqhdeeMHU9D83SxgwYMC7777bRGGOw9+E5giCIJfLKTWKcNgeRKae4xZG8uv9BS2OxlRVVfE8j79zTdLZC8u+NMbOhBtkL//fhokheKpzbmpqqtzoalCzgjA6OnrSpEn1t8fExPj4+JiZmcnlcmmL9MLc3Fy5mJmZWW3t31PrampqlAucPn36vffeO378uJ2dXf3Ccrm8zt6UmZiYTJ061draWnrr7e3dROHa2tomfgoqJwgCy7JGcs7NCTk2low/LS6L534I1NoMPUqpqampbl6XDZVuXlhESr69LX4fyJqbG+A8rqc65yz75P+PzTpHQ4YMKS4ubuynnp6e6enp0uvHjx/LZDI3N7c6BR4/fqx4m56e7unpKb0+e/bsvHnzDhw44OfnpyickZFBKWUYRtrhqFGjGq29TDZ16lSpT/WJWJZtzhkBVaGUGtU5tzUjJ8ayI0/wy67Rb/y1E0Xsf2nl6MZJN0/43hTR3oyM9dK5iqmEys+5Cvb1zDPPHDx4UOo73bVr16RJk6R26IULF1JSUgghI0eOzMzMvH79OiEkNjY2KytL6kq9dOnS7Nmzd+/ePXjwYMXeAgICZDJZWFgYIeTBgwdxcXENNkYBdJCdKTk9XvZXOl0RK2i7LmC8BEo+ixNX9kfHQHOpoNUcGhq6ZcsWf3//jh07Xr58WcowQsjSpUunTp36/vvv29rafvnllxMmTBg2bNj58+e/+uorqTPz5ZdfFgThgw8+kMpPmjRp5cqVJiYmq1evnjt3bnBw8MWLF5ctW+bu7t76SgJohqMZOTNBNuIYb8aJH/U1zO/joON2pYjO5mS0J9YONhejkrkMoihevXq1oKBg8ODB0lAfIeThw4c2NjbSZFFCSEpKirSgXpoRSgh5/Pix8nCgjY2NNNeGEPLo0aPExMTOnTv7+Pg0cVxHR8eUlJRmdo2WlZXZ2Ng87UeDFpMmy1hYGOPtOLMqyYjj/MJu7Ae9NZqFVVVVGCPUMF27sPAi6b6f3zqUG9HGYINQ5edcNUGoLQhCXWbMQUgIyaigI44Lb/qymnz2DYJQ83TtwrL1rrg3VQwbb4BzZBRUfs4N+WQBaJGnFRMxkRtxTJCx5C1f9JGCJlQL5PN4cW8Ivgk9HQQhgLq0tWLCJ3DBJwSOIYu6IwtB7X68I/ZzYgJcDbZTVE0QhABq1N6GiZjABR8XWIa85oMsBDUqrSXfJAjhE3BVf2o4ZQDq1cGGiZjIBR8XGEJeRRaC2nyXKIzzYns4oDn41BCEAGrX0YaJmMCFnBAIshDUI7uK/HhHjH0Gl/SWwFkD0IROtshCUKNPrwsvdmHb0HGEEgAAGnJJREFUWaM52BIIQgANkbJw5AlBpOR1zJ0B1blbTA+kiXdnmGi7IvoKQQigOZ1smYiJXMhxQaDkTaypABX55xXxX304BzNt10NvIQgBNKqjDXNuIhdyQhAoeacHshBa62wWvV1M94/C2sGWQxACaFr7/2ZhrUje64UshJYTKVkSI3wzkDXF31Er4OQBaEE7a+bcRG7zXfGrBFHbdQE99u8k0daETO+AK3mr4PQBaIeXFRM5UfZHsrjyOp7ZBC1RIicrYoUfBqFTtLUQhABa08aSnJsoO5hGl15BFsJTWxUnTGzL9nPCkonWQhACaJOrBYmYKAvPpO9EC3r8IBjQuLvF9Pdk8Qs/NAdVAEEIoGVOZiR8giw2n74SJYgIQ2ied6KFj/pyrkb6lDMVQxACaJ+dKTk9XpZaRp8/J9Ri9gw8yZ+pYnYVlqKqDM4jgE6wkpHjY2VltXRGuFCNEUNoXHktWXJZ3BjEyXD9VhGcSABdYc6RA6NkFhyZdIovr9V2bUBXfRYnjGjDDHXHHBmVQRAC6BATluwI5jraMqNO8kU12q4N6J4bhfS3ZHG1P+bIqBKCEEC3cAzZPIQb6sYMP85nVWq7NqBLREpeuyB8PgBzZFQMQQigcxhCVgdwszuyQ4/xD8owkRT+Y/NdkWPJQjzGS9Vwr1EAHfVhX9bJnAw/Jhwfy/V2xICQscuooJ/ECpGTZPhTUDl8swDQXa/5sN8HsqNP8lHZaBcau0UXxbd6cN3tkYOqhyAE0GkzO7A7RshmhPOHH2KBofHalSKmldN/9cEVWy1wWgF03ShP5sRY2aKLws/3kIXGKLeKLIkR/j2Mw7OW1ATnFUAPDHBmzk+SfZ0gropDFhqdNy4JL3Zl/ZzRKaouCEIA/dDZlrkYKjv8UHz9goD7cxuPXSni3WK6sj8WDqoRghBAb7hZkLMTZWnldGqYUMlruzagfhkV9N0Y4ffhnBlyUJ0QhAD6xMaEHB0jczInwcf53Cpt1wbUiRKy4Lzwli/XH52iaoYgBNAzJiz59zBunBcTdJS/V4JOUoO1/pZYWkuWYaao+uEUA+gfhpBPB3Af9WWHH8MSQ8N0s4h+Hif8MQKPmNAEnGMAffVSV3b7CNmMcH7HfUwlNShVPJkTIawO4DrbolNUExCEAHpslCcTMUG2PFb89Doebm84/hEj9HZk5nfB9VlDcKIB9FsPByZ6suyvdPH5s3iiryHY+0AMz6Q/DsE8Uc1BEALoPWlZhUjI+DNcDqaS6rP7pfTtaGFPCGdrou2qGBMEIYAhMOfIzmButAcNOkbjC9BLqpeqeDIzXPikH9ZLaBqCEMBAMIR82Ev8diAZ+xf/Zyqmz+ifNy8JPRyYN3xxWdY0PI8QwKBMb890seemhgk3iugn/Tm0LPTFT3fEa/k0ejKuyVqArx4AhqafE3Nliiwik04LE8pqtV0baIaLOXTldeHAKM4KOagNCEIAA+RqQc5MkLlbksDDfBLuPqPb0ivos+HCtuEyrBrUFgQhgGEyZcmPg7klvdihx/ijjzBkqKMqeTIlTPhHL3acF1JQa9AOBzBkL3djezkyM8OFK3l0ZX8OY4Y6RaRk3jmhlwPzfi+0SbQJZx/AwPm7MNeekV3KoRP+4vOrtV0bULLsqpBfTbcMxdp5LUMQAhg+F3NyerxsgDMz4BAfnYshQ53w0x3x8EN6YLTMFJdhbcO/AIBR4Bjy5UBuYxA3NYxfexM3JtWyIw/FVXHiiXGck5m2qwIIQgCjMqkdc3mKbM8DcWqYUFij7doYq4s59JULwuExXEcbjNnqBAQhgHHxtmaiJsk625L+B/mLOWgZalpCIZ1+ht8+QuaH+6jpDAQhgNExYcl3AdzGwdzMcH5VnCggDTUlqYRO+EvYEMSN9kQK6hAEIYCRmtiWufaM7Hy2GHycf1SOMFS7B2V09EnhCz92RgdceHWLatYR8jy/a9eue/fu+fn5PfPMMw2WuXDhwunTp93d3efNm2djY0MIqa6uDgsLi4+Pt7CwmDBhgq+vr1QyPDw8JSVFem1mZjZ//nyVVBIA6vCwZE6Nk625IQ48zP8QyM3phAu0uqSV0ZEnhGV92Be74iTrHNX8k8ybN2/Tpk12dnYffvjhhx9+WL/Arl27pk+fbmFhERYWNnTo0NraWkLIRx99tGbNmpqamsePHwcEBBw+fFgqvHXr1j/++CM2NjY2NjYhIUElNQSABrEM+aA3+9c42edx4pyzmEGjFqllNPiE8H4v9vXuSEGdRFvt7t27FhYWhYWF0msrK6uioiLlAqIodu/efffu3ZRSnue7d+++b98+SmlpaamizKpVq0JCQqTXs2bN2rRpU3MO7eDgIB23OZQPBxrA83xlZaW2a2FcKisreZ5v4e/W0nejea+dtSceiaqtlWF74oXlXrHYblftj7cFzdTHGKj8Yq6Crydnz54NCAhwcHAghHTr1s3Nze3y5cvKBbKzs+/cuTN+/HhCCMdxY8aMiYiIIIRIHaQShmHMzP5eUHP16tUNGzaEhYWJIu6RCKAJFjKyNpD7fQT35iXhlSihFI+tUIWEQhp8XFjZH21BnaaCMcLs7Gw3NzfFWzc3t8zMTOUCWVlZFhYWtra20lt3d/eYmBjlAo8fP167du3u3bult23atCkuLk5KSlq3bp2Xl9fp06dNTEwaPHR1dfWKFSsUCdq9e/e5c+c2Vs/q6urG9gPqIAiCXC5nGMyO05zq6mpRFDmu5bfsGuRALk8gH8WxPfeL6/3FsZ6YRPMETVxYLuYyc84z6/zp1HZ8NW5upzpPdTE3MTF54v+I5gahcutNYe3atQsXLmRZVrndJghCnaNyHCcIgnIBmezv4+bn50+YMOHdd98dNWqUYrfSi8rKyl69eu3cubOx+TIsy9rb21tYWEhvXV1dm/jAHMe15gIBLYBzrmHcf7VmJ/Yc2TiInM0ir19igx6S1QOJs7mqKmiAGjvhhx6RN6PJ70PJSA98F1Sxp/ojb8538eYGYW5ubv2NpqamhBAPD4/Tp08rNmZlZXl4eCgXa9OmjVwuLygocHJyIoRkZmYqChQVFY0dOzY0NPTjjz+uv39LS8uAgICkpKTGamVqarpkyRKpV/aJTExM0CLUJJZlKaU455rE83xzvv82x5h25KYH+fia0P+I+G0AN68zevYa1uCF5Yeb4pob4qnxXD8npKDqqfxi3tw/bouGSP/fxo4dGxsbK3WHXrt2raysLCgoiBCSnp6enJxMCHF1dfXz8zt48CAhpKam5uTJkxMmTCCElJSUjBs3btiwYV9++aXiQJRSnuel1+Xl5dHR0T4+Pir8wADQTJYy8n0gd2ysbO0NcdQJPOC3WXiRLLoo/DtJvBiKFNQbKhgjbNu27WuvvRYcHDx27NgDBw6sXLnS0tKSEPLjjz8mJiYePXqUELJq1aq5c+dev349Pj6+U6dOY8aMIYR8/PHHcXFxtra2o0ePJoS4ublt3769tra2Q4cOwcHB5ubmZ86c8fHxmTNnTusrCQAtM8CZuTJFtuG2OPgo/3p3dlkfzhKPMW1EfjV5Npy3MiEXQ2U26ArRHwylqvmWFxkZmZSUNGDAgP79+0tbHjx4UF5e3rt3b+ltSkpKZGSki4vL+PHjpTHC27dvK0+rMTc3HzJkCCEkMTExISFBEISuXbtKjcvGODo6pqSkNLNrtKysrMGRTlATabKMYgQXNKCqqsrU1FRN47KZlfS9y2J0Dv0uAPdG+ZviwhKbT2eEC891YlYN4Fg0BdVJ5RdzlQWhViAIdRmCUPPUGoSSyCz6TrTgZEa+D+T6ouvvvxeWH++IK68LPw3mprbHVwS1U/nFHH0cAPAUhrdhrk+Vbb0rjv//9u4/qIlrXwD42d3wG/kVBAKGwOOHP0BSioIiSimgoIWLVqVFLC32jc92fFNnWqevT/vmldfpnV7bubYXe1vgTa3aH4hWRIowipA+axF811DFIDWpPyDBiAQBY5LdPe+PvXdfBiiCSiDZ7+evsydns985bM53N3tyOEmvlpIli8hgd0GnQ4OZKD7NXLuHz+aIIr0E3RX2Cy5eAACTQxHoX+aTqg1Os91Q3BF6VxszYJ7umKbJGS1eVu8c4oHO5UIWtGOQCAEAj8LbGf1xMfW3daKe+yj6sOVP7ex9erpjsqFhGv3rOealJuaTxfSfl1Au8HNZewaJEADw6KQexH+voJrWiFr1OLLS8udLrFEA6bD+Fl54hL5nRu3rROlBzMN3ADMbJEIAwOOa70NUplN1WaIfdTii0vKndnbQQZcq7R7GLzQyr51lPltGfZlK+bo8fBcw80EiBAA8GXI/4kgG1ZAtutiHI76z7GpjdMbpjunJMdLov/7Gyo/S0d7o0vOiVXPgiaDjgEQIAHiSYn2JQ2nUz38QGcwopspSrGAu9tnxb7QQQjSLKjrZuYfp9ru4NU/0XgLlBtPtHQskQgDAk/dPs4i/JFNdG52ivYncBmb5Cfrra6zJ3p6m0Sw68Csbc4Q+9Ct7OJ2qTKfCZ8GNoAOCCxsAwFTxc0Fvy8k3F5I1N9i/XmF3/MwURJCvRJNxfjM9nQzT6Mur7Ee/sGGe6K8pVJpkpgcMHgckQgDA1BKRaG0YuTaM1AziL6+yuQ2MjzMqjCI3hhOhnjMuwVwdwJ+r2K+62Gck5Ndp1JKAGRcheOIgEQIAbCR8FvGfCdR/PI1+1OGvr7EJx5goL2JdOPmHUCLKe5rzjcGMjmjY/V3sr/fwy1FkW55INvOSNJgikAgBADZFEihVQqRKqNJk6owWH/2NTfuB9RShbCmxMoRcISE8bDgs9RrRiRvssevsjzqcGUK+uZDMlpJOMHdCYGDRbTBVYNFt27PBottTASN0sQ+fvIUbbrFtd3CcH7EskEgOJBJnE1OxkOldEzp3Gzf1sKd78PUhnBlC5smINaHko/3jJBhYbA8W3QYAOBoCoXgxES8m/k1OGmn0sx7/jw6Xqdh//hE7kUjuRyz0I+b5EFFeRIQXkrgTk8qNDxikGcSdA/iKAbXfxf97B2vv46QAYkUQ+ZdkMnE2IYL7P8GDRAgAmEHcRChNQqRJCO7HXTeGsPIuvtyPFFpcrmI1g/iuCQW5E8HuSOyCvJ0JTyc0ywnxyczEoPs0MphRvwn3GlGvEfebUNgsIsoLLfAlVkuJf3+KnO9DUPD4D1iBRAgAmLlCPYlQTyIn9P9rHjBIdx9rjajvARow42Ea3bMghv37q84U8hAhLyfk60IGuKEgNyLIHUHWA+ODRAgAsCeuFAqbRYT9/QkR5DjwBMC34wAAAAQNEiEAAABBE0oiHB4eLi0tne4ohKWjo6Ompma6oxCWmpqajo6O6Y5CWEpLS+/fvz/dUQjIgwcPPvnkkyf7nkJJhDqd7osvvpjuKITl/PnztbW10x2FsNTW1p4/f366oxCWzz//XKfTTXcUAnLnzp3PPvvsyb6nUBIhAAAAMCZIhAAAAAQNEiEAAABBs++1Rl1dXSUSCUk+PJ3TNK3T6ebMmWODqABnaGjIaDTOnj17ugMREL1e7+bm5unpOd2BCMjNmzclEolIBL/JthGGYXp6eqRS6QTbFxQUlJSUjN/Gvv94165dM5lME2xsMplcXFymNB5gjWVZhmGcnB5pJWPwSCwWC0VRE7k0BE8KDCy2N6k+l0gkD21j33eEAAAAwGOCK0cAAACCBokQAACAoEEiBAAAIGiQCAEAAAiafc8atXbx4sVvv/3W2dm5qKgoIiJidIO+vr7y8vLe3t7s7OzMzEy+vrGxsba2dvbs2cXFxQEBATYM2e4pFIrq6mqxWFxcXBwUFDTi1d7e3pqaGpVK5efnl5+fz/9Rqqqq7t69y5WDgoJyc3NtGrQ9M5vNFRUVXV1dTz31VGFh4ejZoT/88MOtW7e4speX1wsvvMCVBwYGysrKenp6nn322eeee86mQdu5wcHBsrKyW7durVixIi8vb8SrBoOhsrLSumb58uXz58+/fv16fX09X5mdnT3x6f4C19/ff+HCBbVavXjx4vj4+DHbXLly5dChQyzLbtq0KSYmhqtkWfbAgQNKpTIqKmrLli3Ozs4TP6iD3BG2trauWLHCx8fHYrEkJiZev359RAOTyZScnHz58uXw8PCioqIDBw5w9ZWVlS+++KJMJuvq6lq6dCksnjtx1dXV69atk0qlN27cSEpKunfv3ogGW7dubWxsDA4O1mq1cXFxra2tXH1JSUlTU5NarVar1d3d3TYP3I4VFBR89913UVFRe/fu3bFjx+gGe/fura2t5fr25s2bXCXDMM8880xLS0tERMT27dth9fmJwxhnZmYqFIrIyMi33nrro48+GtHAYrGo/+HSpUtbt27t6+tDCCmVypKSEv4lo9E4HeHbpQ0bNrz99tvvvfdeXV3dmA06OzuXLFlCUZSbmxs3qnP1b7zxxt69e6OioiorKzdt2jS5o2KHsHHjxt27d3Pll156aefOnSMaHDx4UC6XsyyLMa6qqpo3bx5XlsvlBw8e5NokJSVVVFTYMGr7lpSUVFZWxpVTU1NLS0tHNDAajXy5qKho27ZtXDkuLu7MmTM2idGhqFQqNzc3g8GAMdZoNK6urnq9fkSblStXfvPNNyMqjx8/HhERQdM0xri+vl4qlXJl8FANDQ1z5syxWCwYY4VCERgYaDKZfq/x/v37o6OjuYGluro6JSXFdoE6EIZhMMZr1659//33x2ywbdu21157jSvv2LFjy5YtGGO9Xu/q6qpWqzHGAwMD7u7uKpVq4gd1kDvC5uZm/tvOzMzM5ubmEQ0UCkVGRgZBEFwDlUrV29trMBiUSmVGRsY4O4IxPXjw4Pz58+P3uaurK182mUzWy51UV1d//PHHDQ0NNgjVYSgUisTERG9vb4RQWFiYTCZraWkZ3ezUqVN79uw5fvw4y7JcTXNzc3p6OkVRCKG0tDStVnvt2jVbRm6/mpub09LSuFVjli1bNjQ01NnZ+XuNKyoqiouLuUEGIaTX6/fs2VNeXt7T02OjcB3CQ5eDGHO0b2lpkUql4eHhCCEvL6/ExESFQjGJgz5qtDMITdN6vZ5fyisgIECr1Y5oo9Vq+QZeXl6urq5arVar1ZIk6e/vz9UHBgbCKTtBOp0OY8w/Uh2/6xQKRV1d3euvv85tyuVyiqK0Wu2rr77KP8QCD2V9DiOEAgICRvf53LlzPTw87ty5s3PnzpUrVzIMgxDS6XT8jk5OTn5+fqM/IGBM1l3HjRW/d56r1epz585t3ryZ2/Tw8IiPjzcYDCdPnlywYMG5c+dsFLEAWH8Q+NF+xKdjsoO5I0yWIUmSoiiaprlNmqZHPyYViUR8A+7u29nZ2cnJiVsGjLtYtlgssFTSBHELp/FdOk7X/fLLLxs2bKioqJDJZFzNV199xRXefPPNqKios2fPLlu2bOpDtnsikYhLbByLxTL6POf/YemuXbvmzp3LPcedyI5gTBPvuvLy8uzs7ODgYG4zPT09PT2dK7/zzju7du06ffr0VEcrENaDOT/aP+ZJ7gh3hCRJBgUF8dMuuru7+dORFxISwl8g9Pb2WiyW4OBgiURCEARf393dPZFV6QBCKDAwkKIo6z4fs+s6OjpWrVr16aefPv/882O+SWRkpEajmdpYHUVISIj13KIxz3Oep6dnXFwc17fWOw4NDQ0MDIyzI7Bm3XUmk6mvr2/MrqNpev/+/cXFxWO+SXJyslqtnsIoBcZ6MOc/BZP6dIzmCIkQIZSbm1tVVYUQwhhXVVXl5ORw5cbGxsHBQYRQTk5ObW0tNym0qqoqJSXF19fXw8MjPT2d29FkMtXU1MBU/gkSiUTZ2dmHDx9GCFkslurqaq7rjEZjY2OjxWJBCHV1da1aterDDz/cuHEjv6PZbOYv3Lq6ulQqFT/7GYwvKyurvb2dG1JbWlqGhoZSUlIQQhqNRqlUIoQYhuF6HiGk0+laW1u5vs3Jyamvr+em9R49enThwoX83TkYX05OzunTp/v7+xFCx48fDwsLi46ORghduXLF+mHhyZMnGYZZvXo1X2M9TfTEiROxsbE2jNoBGQwGfhZCbm4uN/IghA4fPsyN9ikpKcPDw9xX0BqNpr29PSsraxIHeJzpPTPHb7/9FhISsm7duoyMjJiYmP7+foyx2WxGCLW1tWGMGYZZs2ZNfHz85s2bxWJxc3Mzt+NPP/0kFosLCwsXLVqUmZkJs+kmrq2tTSwWFxQULFmyJDU11Ww2Y4y7uroQQrdv38YYp6WleXp6JvzD9u3bMcYXL16USqXr169fv369l5fX6Pm9YBy7d+8ODQ3lfrW5b98+vjIrKwtjrNVqAwMD8/Ly8vPzxWLxK6+8wu+Yn58fExNTVFTk7+9fV1c3PdHbp6Kionnz5r388sv+/v7ff/89X7l161a+TV5e3ogzOT8/PzU1tbCwMCEhITQ0tKOjw6ZB27MPPvggISHBx8cnJCQkISHh2LFjGOPGxkYXFxeuQW9vb0RExOrVq3NycsLDw7VaLVe/b9++wMDA4uJimUz27rvvTuqgjvPfJwYGBk6dOuXi4pKRkcHNV8QYt7a2xsbGuru7I4RYlm1qatLr9cuXL7e+a9bpdM3NzWKxOC0tjXtYCCbo9u3bTU1Nvr6+/Mw6k8mkVCqffvppkUjU2dk5NDTEN/b29o6MjMQYX758WaVSkSQpl8vHXPoAjOPChQtXr16Vy+ULFizgarq7u4eHh7k7lc7Ozo6ODpqmY2Nj58+fz++FMVYoFFqtNjk5OTQ0dHpCt08Y47Nnz968eXPp0qVhYWFcpUajIQiC31QqlTKZzMfHh99rYGCgpaWlr68vKCho6dKl1jOowfhu3Lih1+v5TZlM5u/vPzg42NnZuWjRIq5yeHj41KlTGOOMjAzr6egdHR1KpTI6OjohIWFSB3WcRAgAAAA8Agd5RggAAAA8GkiEAAAABA0SIQAAAEGDRAgAAEDQIBECAAAQNEiEAAAABA0SIQAAAEGDRAgAAEDQIBECAAAQNEiEAAAABA0SIQAAAEH7P35GkmeO10NGAAAAAElFTkSuQmCC",
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
""
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = @elapsed u = ( -L \\ f.(x) )/n^2\n",
"println(\"Solving the linear system of equations took $t seconds\\nThat is $(round(Int, t_dense/t))× faster than before!\")\n",
"plot(x, u; title=L\"numerical solution of $-u\\prime\\prime = f$\", legend=false)"
]
},
{
"cell_type": "markdown",
"id": "c3789505",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 8.** Explain why storing matrices in a sparse format is also convienent for implementating iterative methods.\n",
"\n",
"::: {.box}\n",
"\n",
"***Your answer here***\n",
"\n",
":::\n",
"\n",
"It turns out that, when solving the $-\\Delta u = f$ in higher dimensions, it becomes prohibitively expensive to solve using direct methods. If we have $n+1$ grid points per dimension (we use a uniform grid $\\{(x_{j_1}, \\cdots, x_{j_d})\\}$ where $x_j = \\frac{j}{n}$ and $0 \\leq j \\leq n$), then the size of the linear system we must solve is $(n-1)^d \\times (n-1)^d$ and the computational cost involved in a LU decomposition is $O(n^{3d-2})$. Solving the Laplace equation numerically on $[0,1]^2$ would make a nice presentation/poster/paper - see [here](https://jack.thomaslabs.co.uk/teaching/Math5486/Possion_d=2.html) for some hints!"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.11",
"language": "julia",
"name": "julia-1.11"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}