{
"cells": [
{
"cell_type": "markdown",
"id": "81b1d2f0",
"metadata": {},
"source": [
"---\n",
"title: \"Jacobi & Gauss-Seidel Methods\"\n",
"subtitle: \"Lecture 8: Introduction to Chapter 2: Iterative Methods in Matrix Algebra\"\n",
"date: 2026-02-18\n",
"abstract-title: Overview\n",
"abstract: | \n",
" (1) Reminders on matrix analysis, \n",
" (2) Jacobi & Gauss-Seidel\n",
"format:\n",
" html:\n",
" other-links:\n",
" - text: This notebook\n",
" href: L8.ipynb\n",
"---"
]
},
{
"cell_type": "markdown",
"id": "63cf0867",
"metadata": {},
"source": [
"::: {.callout-note}\n",
"\n",
"I encourage you to play around with the juptyer notebook for this lecture - you can copy the code with the ```this notebook``` button on the side of this page.\n",
"\n",
":::\n",
"\n",
"::: {.callout-warning}\n",
"\n",
"These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.\n",
"\n",
":::\n",
"\n",
"::: {.callout-tip}\n",
"\n",
"This lecture is mostly based on @Burden sections 7.1-7.3.\n",
"\n",
":::"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "db6fe346",
"metadata": {},
"outputs": [],
"source": [
"#| echo: false\n",
"\n",
"# IF YOU ARE READING THIS, THANK YOU\n",
"# FOR DOWNLOADING THE LECTURE NOTES.\n",
"# THE FIRST PERSON TO LET ME KNOW,\n",
"# WINS SOME CHOCOLATES OF YOUR CHOICE\n",
"\n",
"using Plots, LaTeXStrings"
]
},
{
"cell_type": "markdown",
"id": "5ea01014",
"metadata": {},
"source": [
"## Recap: Matrix Analysis\n",
"\n",
"@Burden 7.1 - 7.2\n",
"\n",
"We will define $\\|\\cdot\\|_{p}$ to be the induced matrix norm from the vector norm $|\\cdot|_p$: that is,\n",
"\n",
"\\begin{align}\n",
" |x|_p := \\begin{cases}\n",
" \\big(\\sum_{i=1}^n|x_i|^p\\big)^{\\frac1p}, &\\text{if } 1\\leq p < \\infty\\\\\n",
" %\n",
" \\max_{1\\leq i\\leq n} |x_i|, &\\text{if } p= \\infty.\n",
" \\end{cases}\n",
"\\end{align}\n",
"\n",
"and\n",
"\n",
"\\begin{align}\n",
" \\|A\\|_p := \\max_{|x|_p = 1} |Ax|_p\n",
" %\n",
" = \\max_{ x \\not= 0} \\frac{|Ax|_p}{|x|_p}.\n",
"\\end{align}\n",
"\n",
"Notice that $|Ax|_p \\leq \\|A\\|_p |x|_p$.\n",
"\n",
"You can use ```LinearAlgebra``` to work with matrices in Julia:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "607e643d",
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra"
]
},
{
"cell_type": "markdown",
"id": "93e0c919",
"metadata": {},
"source": [
"::: {#exr-1}\n",
"\n",
"Compute $\\| A \\|_2$ and $\\| A\\|_{\\infty}$ when \n",
"\n",
"\\begin{align}\n",
"A = \\begin{pmatrix}\n",
" 1 & 0 \\\\ -1 & 1 \n",
"\\end{pmatrix}.\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "3a9a5d8a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Matrix{Int64}:\n",
" 1 0\n",
" -1 1"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [1 0; -1 1]"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c20bf1af",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1.6180339887498951, 2.0)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"opnorm(A), opnorm(A,Inf)"
]
},
{
"cell_type": "markdown",
"id": "d4ce2d4d",
"metadata": {},
"source": [
":::\n",
"\n",
"::: {.callout-tip}\n",
"\n",
"You can show \n",
"\n",
"\\begin{align}\n",
" \\| A \\|_\\infty = \\max_{1\\leq i \\leq n} \\sum_{j=1}^n |A_{ij}|\n",
"\\end{align}\n",
"\n",
":::"
]
},
{
"cell_type": "markdown",
"id": "76cc7786",
"metadata": {},
"source": [
"::: {.callout-warning}\n",
"\n",
"In Julia, ```norm(A)``` is the *Frobenius norm*:\n",
"\n",
"\\begin{align}\n",
" \\| A \\|_{\\mathrm{F}} := \\sqrt{\\sum_{ij} |A_{ij}|^2}\n",
"\\end{align}\n",
"\n",
"which is not an induced matrix norm.\n",
"\n",
":::"
]
},
{
"cell_type": "markdown",
"id": "9d12fea3",
"metadata": {},
"source": [
"::: {#exr-2}\n",
"\n",
"Find the eigenvalues and eigenvectors of \n",
"\n",
"\\begin{align}\n",
"A = \\begin{pmatrix}\n",
" 0 & 1 \\\\ -1 & 0 \n",
"\\end{pmatrix}.\n",
"\\end{align}\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "5f3994a5",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Matrix{Int64}:\n",
" 0 1\n",
" -1 0"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [0 1 ; -1 0]"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "7c0976b0",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}\n",
"values:\n",
"2-element Vector{ComplexF64}:\n",
" 0.0 - 1.0im\n",
" 0.0 + 1.0im\n",
"vectors:\n",
"2×2 Matrix{ComplexF64}:\n",
" 0.707107-0.0im 0.707107+0.0im\n",
" 0.0-0.707107im 0.0+0.707107im"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"λ, v = eigen(A)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "c44eebba",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.0 + 0.0im, 0.0 + 0.0im)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"det(A - λ[1]I), det(A - λ[2]I)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "acfc6e32",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(ComplexF64[0.0 + 0.0im, 0.0 + 0.0im], ComplexF64[0.0 + 0.0im, 0.0 + 0.0im])"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(A - λ[1]I) * v[:,1], (A - λ[2]I) * v[:,2]"
]
},
{
"cell_type": "markdown",
"id": "a9bdbe87",
"metadata": {},
"source": [
":::"
]
},
{
"cell_type": "markdown",
"id": "1b41f72f",
"metadata": {},
"source": [
"We will write $\\sigma(A) := \\{ \\lambda : \\exists x \\not= 0 \\text{ s.t. } A x = \\lambda x \\}$ for the set of eigenvalues of $A$. The *spectral radius* $\\rho(A)$ is defined as\n",
"\n",
"\\begin{align}\n",
" \\rho(A) := \\max_{\\lambda \\in \\sigma(A)} |\\lambda|.\n",
"\\end{align}\n",
"\n",
"::: {#lem-rho}\n",
"\n",
"$\\rho(A) \\leq \\|A\\|$ for every submultiplicative matrix norm.\n",
"\n",
":::\n",
"\n",
"
\n",
"Proof. \n",
"\n",
"Here, we prove the result for the induced matrix norm - see [Appendix: Linear Algebra](App2.html) for a more general proof. \n",
"\n",
"Let $|\\cdot|$ be any vector norm and $\\|\\cdot\\|$ the corresponding induced matrix norm. For all $\\lambda \\in \\sigma(A)$, there exists a normalised eigenvector $x$ (i.e. $Ax = \\lambda x$ and $|x| = 1$). Therefore, \n",
"\n",
"\\begin{align}\n",
" |\\lambda| &= \\left|\\lambda x \\right| \n",
" %\n",
" = \\left| Ax \\right| \n",
" %\n",
" \\leq \\| A \\|.\n",
"\\end{align}\n",
"\n",
"
\n",
"\n",
"::: {#thm-rho}\n",
"\n",
"If $A$ is normal (i.e. $A^\\star A = A A^\\star$ where $A^\\star$ is the conjugate transpose of $A$), then $\\rho(A) = \\|A\\|_2$. \n",
"\n",
"For fixed $A$ and $\\epsilon > 0$, there exists a submultiplicative matrix norm $\\| \\cdot \\|_\\epsilon$ such that \n",
"\n",
"\\begin{align}\n",
" \\rho(A) \\leq \\| A \\|_{\\epsilon} \\leq \\rho(A) + \\epsilon.\n",
"\\end{align}\n",
"\n",
":::\n",
"\n",
"I've put the proof here [Appendix: Linear Algebra](App2.html) for completeness."
]
},
{
"cell_type": "markdown",
"id": "b8ee6891",
"metadata": {},
"source": [
"## Jacobi Method\n",
"\n",
"7.3 @Burden\n",
"\n",
"Implementation of Jacobi:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "d365a393",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Jacobi (generic function with 1 method)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function Jacobi(A,b,x; Niter=100, tol=1e-10) \n",
" \n",
" X = [x]\n",
"\n",
" L = -tril(A, -1)\n",
" U = -triu(A, 1)\n",
" D = Diagonal(diag(A))\n",
"\n",
" if (any(iszero, diag(D)))\n",
" @warn \"diagonal of A is not invertible\"\n",
" return X\n",
" end\n",
"\n",
" for i ∈ 1:Niter\n",
" push!( X, inv(D)*( b + (L+U)*X[end] ) ) \n",
" r = norm(X[end]-X[end-1],Inf)/norm(X[end],Inf)\n",
" if (r < tol)\n",
" return X\n",
" end \n",
" end\n",
"\n",
" @warn \"max iterations exceeded\"\n",
" return X\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "d22adb57",
"metadata": {},
"source": [
"::: {#exm-Jacobi}\n",
"# Example 1 of @Burden\n",
"\n",
"Applying the Jacobi iteration to $Ax = b$ with"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "33aad3a1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Int64}:\n",
" 10 -1 2 0\n",
" -1 11 -1 3\n",
" 2 -1 10 -1\n",
" 0 3 -1 8"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b = [6, 25, -11, 15];\n",
"A = [10 -1 2 0; \n",
" -1 11 -1 3; \n",
" 2 -1 10 -1; \n",
" 0 3 -1 8]"
]
},
{
"cell_type": "markdown",
"id": "c0fce944",
"metadata": {},
"source": [
"gives:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "049866ca",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Vector{Float64}:\n",
" 1.0000000000277243\n",
" 1.9999999999534517\n",
" -0.9999999999639979\n",
" 0.9999999999487015"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X = Jacobi(A,b,[0. ,0., 0., 0.])\n",
"X[end]"
]
},
{
"cell_type": "markdown",
"id": "d2f20078",
"metadata": {},
"source": [
":::\n",
"\n",
"::: {#exr-1}\n",
"\n",
"Verify that\n",
"\n",
"\\begin{align}\n",
" x = \\begin{pmatrix}\n",
" 1 \\\\ 2 \\\\ -1 \\\\ 1\n",
" \\end{pmatrix}\n",
"\\end{align}\n",
"\n",
"solves $Ax = b$ with $A$ and $b$ given by @exm-Jacobi.\n",
"\n",
":::\n",
"\n",
"Here is an implementation of the Gauss-Siedel iteration:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "7abf11a9",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"GaussSiedel (generic function with 1 method)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function GaussSiedel(A,b,x; Niter=100, tol=1e-10) \n",
" \n",
" X = [x]\n",
"\n",
" L = -tril(A, -1)\n",
" U = -triu(A, 1)\n",
" D = Diagonal(diag(A))\n",
"\n",
" if (any(iszero, diag(D)))\n",
" @warn \"diagonal of A is not invertible\"\n",
" return X\n",
" end\n",
"\n",
" for i ∈ 1:Niter\n",
" push!( X, inv(D-L)*( b + U*X[end] ) )\n",
" r = norm(X[end]-X[end-1],Inf)/norm(X[end],Inf)\n",
" if (r < tol)\n",
" return X\n",
" end \n",
" end\n",
"\n",
" @warn \"max iterations exceeded\"\n",
" return X\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "666c4471",
"metadata": {},
"source": [
"::: {#exm-GS}\n",
"\n",
"Returning to @exm-Jacobi, we see that Gauss-Seidel converges in fewer iterations:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "3779e8d4",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Vector{Float64}:\n",
" 0.9999999999812351\n",
" 1.9999999999852736\n",
" -0.9999999999921585\n",
" 1.0000000000065026"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Y = GaussSiedel(A,b,[0.,0., 0., 0.])\n",
"Y[end]"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "b094ab30",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"# iterations for Jacobi = 29 \n",
"# iterations for Gauss-Siedel = 12\n"
]
}
],
"source": [
"println( \"# iterations for Jacobi = $(length(X)) \\n# iterations for Gauss-Siedel = $(length(Y))\")"
]
},
{
"cell_type": "markdown",
"id": "781372bf",
"metadata": {},
"source": [
"In fact, we have the following error curves:"
]
},
{
"cell_type": "code",
"execution_count": 61,
"id": "4e6e5cc2",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdeSBU6xsH8O+MwSAhkXZtlK3FkkpdFS23pl371W2jXZtoV91ulBSt6rbp3haq26Joue0rCiW0b78WSYnszPn9oWyNpZg5M+b5/HOdcx5zvs0dHmfmPe/LYRgGhBBCiLzish2AEEIIYRM1QkIIIXKNGiEhhBC5Ro2QEEKIXKNGSAghRK5RIySEECLXqBESQgiRa9QICSGEyDVqhIQQQuQaNUJCCCFyrVo1wqdPn27cuJHtFKwRCoVsR5AKeXl5bEdgHz0JABiGoR8K0IuhAqpVI3z48OHJkyfZTsGazMxM+rEHkJ6eznYE9tGTAEAoFGZmZrKdgn30YihXtWqEhBBCyI+iRkgIIUSuUSMkhBAi16gREkIIkWvUCAkhhMg1OW2Ee/fute/4S5/uPa9evcp2FkIIIWzisR2ABSMGODw/fsyNyUkD5/eL513W+UyfMZ3tUIQQQtghd1eE6enpV08cvcTkDABGgbkqzF3t6sZ2KEIIIayRuyvCGzdutARPCbmX0eUJmnXFBYXc12yHIoQQwhrZuCJkGGbOnDk2NjZWVlY3btyozENZWlo+Qq4Q6Ijrw3GgH44lMUvfv6+qpIQQQmSMbDTCCxcuxMbGXr169Z9//pk4cWJlHqpmzZqG1jYCjlIYct8iwwsddNU616+PadNA05MRQogcko1GeO7cuf79+wNo0aIFh8N5/bpSb2aevnah0/KlM5oYjmtp+ouW8hOlgTtm3tu1C1pa2LevihITQgiREbLRCBMTE2vVqpX/de3atd9X+q3MBYsWRDyNvxh3V9XfHxoajv/0/HjjgUAAR0f06oWXLyudmBBCiIyQjUaoqan55cuX/K9TU1O1tLSq7KEdHGBggJ49lfva/+35vzdv0L492rWDhwdo2npCCJEHLDfCiIiIP/74Y+TIkVu3bi26/8GDB4MGDbKwsJg2bVpKSoqFhcW1a9cAfP78OSkpqUGDBlUZwscHp05h9Gj8+quu8udlyxAZiadP0aQJ1qypyvMQQgiRQiw3wvPnzyckJCQkJNy+fbtgZ05OTs+ePdu0abN79+63b99OmjRp4MCBz549mzhxYs+ePZctW8bjVeldH0ZGmDgRxsawtcXAgcjKatgQAQGYMgULFkBfH7duAYBQKDx79uzOnTuLRiWEECLrOAzDsJ0Brq6uycnJ27dvz988evTo3LlzHz9+DOD169dNmzZ99eqVtrZ2fHy8jo6Orq5uaY8TGBi4ePHiMWPGFOwZOHBg8+bNK5ojL09hxAgoK+cFBIDDAfDmDUaO5N26xbG0zE5N7f36tUlKSistrQs9e+rs3r3up/+9YpKens7n87lc2Xi7W3xSU1PV1dXZTsEyehIA5OXlZWVlqaqqsh2EZXL+YuDxeBwOp5wayUT5IXfv3rWyssr/un79+rq6urGxsba2tsbGxmV/Y15eXk5OzqdPnwr2ZGRk/MCi7RyOcNcuxd69uUuW5C5bBkBPD+fPZ4eEcEeP/pyWFgDUB5CUNCkkZPTVq1c7duz4E/868REKhbRCPeh5AEBPAoBvTwI9D/QklEsaG+H79+81NTULNmvVqpWQkFCRb9TS0jI0NFxTmU/2lJURHIxOnRTq18f0rxOQDhgAK6tZFy6sL6hKTu4VEXG3a9euP38iMcjLy1NWVqYrwuzsbGVlZbZTsIyeBAB5eXkA6HmgF0O5pPGXpoaGRnp6esFmampq0b4oRrm5aNcOCgoICYGnJ44cKTjSrFld4E7BppLSk6ZNq3TADiGEEJZIYyPU19d/9OhR/tdpaWlv377V19eXxIl5PFhbY+lSNGmC48cxaRKuX88/4u4+qU6dBcAlIFFB4YBQGHLjxq95eZIIRQghRKyksREOGjQoOjo6LCwMwPbt201MTAwNDSV07pUrERiI6GiYm2PPHgwejAcPADRr1uzKlf0ODvvbtRszadLd2NhTd++qtGqFyZORkiKhaIQQQsSCYZWfn59WEStWrMjfv3fvXi0trWbNmjVp0uTOnTsVfLSQkJBevXpVNtOWLYyNDSMUMgzD/P0307Qp8/bt91W5uYyrK1OjBlO7NhMQ8LWcXWlpaXl5eWynYF9KSgrbEdhHTwLDMLm5uWlpaWynYB+9GMolFbdPiJSdnZ2YmFi3bt2Kj/4IDQ319fUNCQmp1ImFQnTogJkzMWIEACxfjuPHcfEiatT4vvb4cYwdC01NaGtj0yZYWlbqzJVEt0/kk/PB4vnoSQDdPvENvRjKJb2/NJWUlOrXr8/Cr3UuF+vXY86cr296LlkCa2v074/s7O9r+/XD9evg86Gmhr594eiIDx8knZcQQkhlSG8jZFOHDujRA6tWfd309UXNmjf79RvQsaOVvv74AQNeFpmW29AQN2+idm3UqweGgakpnj1jJzUhhJCfII33EUoFLy8IBF+/VlCInj9/dqdOu3NzWwCXX74cHBl57u5dDQ2N/OPq6ggMxOrVWL8e3t6o2plQCSGEiBVdEZaiTh2EhRVs7dq4cVVurgHAAX5hmCHv34ecOlW0nMOBmxv278e8efDx+brzf/9DcrIkQxNCCPlh1AgrJPHt27pFNutmZb0XtTiwrS1u3sSRIxg5EmlpOHUKLVti+3bQ9EaEECK1qBFWSMeePY+oqRVsHq5Vq1Mp86s1bIhLl6CqCgsLdOmC2bP/mz//cZ06L//++4GkwhJCCPkB9BlhebZsQWqq8+zZI0JCht27Z5SZeTEjo7uNjbm5eWnfwefjr7+wbRvatp3P5X5IT58MKI8ZU3vnztcHDtQvffEMQgghLKArwvIIBFizhvfqVdB//y367792//yz58CBRZGRSEsr+/tGjvzC559PT98GtAOMhUKlyMhrpqbYvJneKSWEEClCjbA8DRpg5kzMnQvA1NRUIBA0GjwYNjb488+yv+/FixdKSq2AgnWwtPj8Nf/9hzNn8Py5eCMTQgipOGqEFeDqipgYFJ2wZu1abN+ePw1paZo1a8Ywd4GcbztefPlSMzcXR4+iaVMxhiWEEPJDqBFWgJISNmyAiwuysr7u0dODu3vBgoUi8fl8N7eJ2tr9gENc7k4dnUHTpnn++iuGDi28437vXmzYgNxcMecnhBBSOmqEFdOjBwwN4edXuGfGDCQkFF2z8Htz5kw+e/bPxYsfr1mTFh0dvGqV5aNHMDeHpSWcnZGYCGtrHDsGc3NcuSL2fwEhhBCRaNRoha1fjw4dMHIk6tcHAB4PmzZh1Cj07Ikid1aU0LZt27Zt2xZsqqnBzQ3jxuGPP2BsjDlzcOIEzp3Db78hfxxNw4YS+JcQQggpRFeEFdasGVxccPx44R4bG9jYFE5JWmE6OvD1xZUruH0bhoZ4+xb37sHcHG3awMND5OTehBBCxIUa4Y9YuBCTJxfbs3Yttm0re9RMaQwNERiIoCDs2wdraxgb4/p13LgBc3M8e5bx7NmzXPrwkBBCxI8aYeXo6cHNrexRM2Vr3x4XL8LXFytWYNw4LF7MNGnib2XVs0MHt4YNLffvL+szSEIIIZVHjbDSXFzKHTVTLjs7REXB2RkCQUBoaNyHD5cSEgLfvbvq4rL18ePnVRSUEEKICNQIf0qXLkhM/Po1j4eNGzFrVrlzzZSNy4WjI9q0OZ2TM/XbbfhqHz8u79JF69y5yuYlhBBSGmqEP8XSEgsWFG527vxzo2a+p6amDGQVbPJ410aNuuvsDIGA5qMhhBCxoEb4Uzw8EBKCW7cK91Ri1ExRzs5DtbSWAukAgGc5Ofs1NVvGxsLODlZW8PBAZmYlz0AIIaQYaoQ/RV0df/6JqVML58+u9KiZfAJB71WrBA0adNHVtTQ2Hn/kyJbgYJ1hwzB6NO7cwdOnMDVF8SWBCSGEVAo1wp/1229QU8OuXYV78kfN/PtvJR/Y2fn3V68iEhLCY2LO9+9veeUK2rRBu3Z4+hQBAfDzw8KF9DYpIYRUGWqEP4vDwcaNWLgQSUlf9+SPmpk5s5KjZkrg8eDhgR07MGIEXFxgZ4fISOjrV+EZCCFErlEjrARTUwwejGXLCvd07oxOnapk1EwJ+f3v0SN07oynTwv3HzuG4OAqPxshhMgRaoSVs2IFLl0qtn6Ej0+VjJr5nq4uTp7EiBGwtsaBA1931qkDV1f07YsnT6r8hIQQIheoEVZOrVqIjgavyNzlVTRqRiQOBy4uOHUKixfD0RFpabC2RnQ0unSBtTVWrFDMyBDHaQkhpDqjRigGVTRqpjQWFrhzBwwDS0tER0NJCfPmISoKjx9zTUy4YjstIYRUT9QIxUA8o2aKUlfH3r1wd4edHXx9wTBIT39kY7Nl9OjzCxcKe/dGcrKYzkwIIdUNNcIqsn8/Vqwo3MwfNePpKdZzOjoiLAwHDqBNm50dOjjNmaPs7R2WnGxta5vA4Yj1zIQQUn1QI6wi3btj40bExxfu8fGBv784Rs0U1aQJQkLSHj/emJQUmpc3MTPT/e3bv44dc9LQEOtpCSGk+qBGWEV0dTF/frExMnp6jKvrgSFDxvbtO2XYsLCwMDGd+enTB6qqloDytx1mL168KTg6fz7s7BAXJ6aTE0KIzKNGWHWmTcP790XHyLg+eHAlNnb2yZOjAwPde/c+IZ5xLA0aNOBwHhfZ8Sk7W6Xg08mVK+HoiG7d4OKClBRxnJ8QQmQbNcKqw+Nh06aCMTKpqanXjh3bKBSaAh2BIx8/es6bJ47T6urqdu7cRE1tCfASiNLQGNq27dxmzeDlhezsr6s7xcQAgJERAgLAMOJIQQghsooaYZWysUGHDlizBsDLly8NudyCMSuaQLbYBpEePLh15cr6HTq49O697vjxJWfO9Dt3DlevwtQUQUFgGGhrw9cXR49i82bY2uLePTEFIYQQ2UONsKp5e2PTJjx71rx582iGyf62+ymgqaMjpnPyeDwXF+dz5/4JDt7VpUtnACYmOHECmzbB0xPW1rh4EQAsLHDtGoYNg50dXr0SUxZCCJEx1AirWoMGWLAAgYHKyspTFy/uU7v2PxzOFgWFIVzu6vHjJZzFzg4REZg7FxMmwN4e0dFQUMCUKUhIQMOGEs5CCCFSild+CflRs2bl/3fC9OkdunU7e/Kkqrr66UaNdMaORffuMDaWZBYOBw4OGDAAu3ahVy907ozVq4stXnHpEjQ10bq1JEMRQogUoStC8TI2Np45b57T5Mk6ffpgzRr064cPHyQfQ1ERTk549Ajm5rCwgIsLkpOxY8deff32AwfOsbRMHjny46dPks9FCCHso0YoQWPGoH9/jBhRbLUKCapRA25uiI5GZib09YOnT//3xYv/Pn1am5Pz/vjxMwYGjK8vhEJWohFCCGuoEYpTr154/brYHm9v8PmYPZulQABQvz78/dGq1f6MjFVADQCAAfBo4cJLgYGwssLNmyymI4QQSaNGKE5WVnB1LbaHy8W+fTh/Hlu3spSpQNq3LggAubk1atd+ffUqpk7FwIGYPBmZmSxmI4QQyaFGKE7z5+Pmza/3LhRQV8eJE1i+HOfPs5MKADB4sH2NGtu/bWXk5BwSCjtzOBg7FnFx0NMT38oZhBAiXagRipOKCtaswbRpyMkptr9JExw4gFGj8PhxKd8pdrNnT/7113c6Op11dMbWqdNx3ryZixY1cnFBVhY0NbF0KbS12YpGCCESRY1QzAYPRr162LKl5P4uXbBoEQQCfP7MRixwudyDB7fGxR09f37Oixc3V61yuHsX79/DwqLkvDMeHpg4EYmJrMQkhBCxo0Yofps24Y8/8PZtyf1Tp6JrVwwbhrw8NmIBgLa2tomJibKyMgBNTezfD1dXdO/+dbHffK6uqF8fpqbw9WUxKSGEiAs1QvFr0QJjxmDRIhGH/PyQm4v58yWeqVSOjrh1C0FB6NXra+9WU4OHBy5dwqlTXydpI4SQ6oQaoUQsWYLwcHz5UnI/j4fAQBw9ir/+YiOWaE2a4OJFdOgAc3OcOPF1p6EhTp/G8uUYNQqOjkhIYDUiIYRUHWqEEqGujrt3UaOGiEO1auH4cSxYgEuXJB6rVDwePDxw5Ahmz4azc+EIUoEAMTHQ04OpKd68KfMhCCFERlAjlAItW2LvXgwbhqdP2Y5SjLU1bt8GAEtLREZ+3VmjBlavRmws6tVjMRohhFQZaoTSoWdPuLpi0CBpu32vZk34+2PZMvTqBQ8PCIU4e/bc5Mnz/fw8X31byenePbx7x25MQgj5edQIJSs4GAsXij40Zw6srPDbb1I43aeDA8LDcfEiGjde5OCwZ+vW7n/80djcfMCdO3cA3LgBU1P4+JS8W5IQQmQCNULJ6tQJO3ciOlr00U2b8PFjyvz5G9aunT95ctCBA0KpaYqNGiEwMDE5+dLnzwGAHcOMSEw8OGnSYgBOTrhxA//9BxMTnD7NdlBCCPlB1AglS0sLy5Zh6tTC2/SKUlR87+/f1ceH5+7efevWCCcnh+7dGZGVbHj8+JGycjuA821H87dv33/9qjlOnoS3NyZPhkCAFy/YykgIIT+MGqHETZiAnBzs2yfyoM/69QuFwsm5uXaAV2pq7ejoS1IzmrR58+ZcbiRQ0JifpKXpPHtWWCAQ4P59tG0LCwusXk133xNCZAM1QonjcrFpE+bNQ0rK9wcfREdbFXk7tH1yclxMjATDlUVXV3fIkE6amuOBSxzOodq1h40cudzSEs7OhROwqahg+XLcvImHD5GUxGpcQgipGGqEbLCwQM+eWLny+yMGZmZ3OAXvPSJCQ8PQyEiCycqxefOqv/92GDv2hKtrXHj4oY0bLeLiwOfD2BheXsjI+FrWrBn++gu6uqxmJYSQiuFIz0dQlRcaGurr6xsSEsJ2kAp4/x6mprh0CS1bFt397t27Xy0spicktMzNDeVy7zRrdvzBA06R1liG9PR0Pp/P5bLwx82DB1i8GDduYPFijBsHHq/Y0TVrkJODOXOgrCyJMKmpqerq6pI4kxSjJwFAXl5eVlaWqqoq20FYRi+GctEVIUt0dbF0KXbtKrFbT0/v3N27HxYv3j9yZPNly45++sR58oSVgD/E0BCBgTh8GAcOwNQUQUHFBgMNHYrwcJiaQib+RCGEyBu6IpRu3t44ehSXLkFBodxaFq8Iizp3Dq6u4PPh5YUuXb7uTEtLCw5OX75cp2lT+PmhSRMxBqC/f0FPAgC6IvyGXgzloitC6TZ7NpSVsXo12zl+gJ0dbt/G7NkYOxb29oiKEjo6zmjSpJuLi/OnT1b16sVbW8PdXdqm0CGEyC9qhNKNy0VAANavR1gY21F+AJcLBwfExcHBAV26bNq/v2Zi4q2EhCNv3545dswxNPTD06cwM8OHD2wHJYQQaoRSYciQsqbbrl8fmzdjzBikp0swUxVQUoKTE5o1C83NnfFtn+aXL/2fP78SGIh//4WWFpvxCCEkHzVCKWBhgdmzyyoYPBjt2sHNTVKBqpKqqjKQVWRHlrKyMgAzs8LPPd+8kbkuTwipPqgRSoHZsxEfj1OnyqrZvBnBwTh5UlKZqoyT01ANjUXfeuHDjIyT4eGdSkyhunMnjIxw5Agb+Qghco8aoRRQUsKGDXBxQVZWqTUaGggIwMSJMrc2/JgxwxcssKxXr4OuroWZ2dRTp3ZdvKjRuTOK3hWyaBH27IGHB7p3R2wse1kJIXKJGqF0sLeHsTHWrSurpnNnjB4NZ2dJZaoy8+ZNe/36TkJCRHT02Z49zc6fx9Ch6NAB27YV1vzyC+7cQb9++OUXuLggNZW9uIQQOUONUGr4+WHtWrx8WVbNihV48QI7dkgqk1hwOHBxwYUL2LIFQ4YUTknK48HFBffu4dMnGBsjMJDVlIQQuUGNUGo0aoQpU+DuXlaNsjL++Qfu7nj4UFKxxMXYGDdvonlzmJkhOLhwv54eAgKwbx+2by+cy5sQQsSHGqE0cXdHZGQ5v/6NjLBoEUaNqgbrwSsrw9MTQUGYORPOzsVusbexwdmz0NFhLxwhRG5QI5QmKiqIiyv/1/+MGdDRwapVEskkdh074s4dALCwQESE6Jpt27Bvn+jFjAkhpJKoEcogDgd//YXNm3HzJttRqkbNmvD3x4oV6NsX7u4irnXNzbF2LWxtce8eG/kIIdUaNULZVK8etmzB6NHVaXjlkCGIjkZMDGxscOHC886dB9ata1m/vrmHx5p27ZjwcIwfj5494exMS/4SQqoSNUKpdOkS5s4tp2bgQNjYwNVVIoEkpE4dnDiB33/PsbcfcfWq+7t34W/e3PTxeblmzSYuF46OuH8ffD6MjODrixJ35RNCyM+hRiiV2rXDgQO4daucso0bcf48jh+XSCYJ4XBgbR1Ts6Yp0B4AoJia6rVzZ1D+US0t+PoiOBj79qFLF6SksJiUEFJNUCOUSurq8PTE1KnlXPXUqIFdu+DsjHfvJJVMEtLT07ncomvIKWZmZhctsLTEjRuYPx98voSjEUKqIWqE0mrUKNSogZ07yynr1AnjxuH336vTkMo2bdooKV0Bvt5GwuPtTEj4xd292OehXC769IGS0tfN9HR6p5QQ8pOoEUorDgcbN2LRovJHhixbxnz6tH/s2CkODrPGjImKipJIPjFSU1P755919ev3qlNnlI5ON1vb81FRSz99QosW8PISPSHrwoXo0AHh4RLPSgiRfdQIpZiJCYYMgYdHOWU83sz69W8EBMw6dWrQ339Pt7M7d+aMJOKJU9euXV68CLt8eWlsbODZswcNDVX8/XHuHG7fRsuW2Lat5PWfjw+mTkX//nB0RGIih6XUhBCZRI1Quq1YgUOHEB1dRklSUlLk1at+DNMa+AU4nJS0ctYsiQUUHwUFBQMDg9q1axfsMTFBYCD+/ht796J162ITs3E4cHREXBxq1YK1tdrWrcjLYyEzIUQWUSOUblpa+PNPbN1aRsnz58+NinxAqAt8SU4WfzLWdOqEK1ewbh0WLIC9fbHJaDQ0sH49QkPTjxyBpSWuX2cvJSFEdvDYDkDKM3Ysxo4t47iBgcFtDif32//Lh4BOvXqSicYiOzvcuYOdOzFgADp2xB9/wMAAr169Wrfur5cv340Y0UVJaeS4cZxr16CtzXZWQoh0oytCmaeurv7brFn9tbWDgK1c7jBFRS8ZX6epgng8ODnh4UOYm8PGBsOGxbdr13/9+taHD49zcYn5++9R8fHUBQkh5aNGWB3MmD//j7NnHy5cyKxbd6ZNG9Ny78SvRlRV4eaGBw8QEfHnhw9bGWYQ0D41dVVERF50kc9Wg4Jw9SqLMQkh0oveGpUdEyZg7ly0bCnyYNu2bQ0NDfl8PrdbN3Ttil9/Rf36Eg7IIi0tKCg8BVoX7ElPb/v48ePWrb/uqVEDI0fC1harV0NPj6WUhBCpRFeEssPUFNOnl19mYgJnZ8ycKf5A0sXEpBWHU7gcR17etYYNjQs2e/dGfDyaNoWxMTw8kJ0t6iEIIXKJGqHsmDYNiYk4cqT8ysWLcf8+jh4VfyYpsn794rp1Z/F4W4FTmprj9PX1hw1reeFCYYGqKjw8cPMmwsJgZgbZv9mSEFI1qBHKDgUFbNyIWbOKLeUukrIytm7FtGmo1vdRlNCoUaOYmP9Wr86dMuVKYOCI+PgN27ZhzBg4OyM9vbCsRQucOoWVK+HkhFGjkJnJXmJCiHSgRihTbGzQsSPWrCm/sksX9O6NRYvEn0mKaGlpzZo1zdNzgb29PQB7+68L+VpYFLvdEMDgwYiNhY0NzVBKCKFGKHO8vbFxIx49qlDlsWO4dk38maSXhgb8/bFsGQQCuLsjJ6fwkKoqJk+Gqmrp30wIkQ/UCGVN/fqYO7dC6/FqaMDHBxMmiJ6mWp44OCAqCnFx6NQJ8fGiayZPxtChePVKsskIIVKAGqEMmj0bjx/j+fPyKx0cYGiI1avFHknq1amDY8cwbRq6dIGXl4h3RH18YG4Oc3N4eNBfDoTIF2qEMkhJCTEx0NevUPHmzdi4EbGx4o0kIxwdERaGkBD06FHy4k9FBW5uuHMHT5/CxASnTrEUkRAicdQIq7t69bB4MSZNqk4r91aGvj4uXIBAAHNzbNsGhmEOHAgaNWrG3Lker169atAAAQHYuhVz50IgwLNnbMclhIgfNUI5MGUKhEJs3852DmnB4cDFBRcuwN8f9epNdHK6tm/fcB8fUwuLATExMQC6d0dkJKyt0b69XN2BQoicokYoy+7cQUWWHuRysXUrFi7E69fizyQzjI3x99+PUlISUlPXAx0ZZvD79ztnzFiRf1RZ+esTpqnJbkxCiNhRI5RlrVrh2DEUnT2lNCYmmDRJDuddK9vz5495vLZFdpg+e/a8aIGiYuHXFy/iyRPJ5CKESBQ1QlmmogJvb0yfXuz+uNIsWiSH866VrWXLlsrKYUV2RGRkNC+t2z17BmtrLF2KjAyJhCOESAo1Qhk3aBAaNcLmzeVXyuW8a2Vr0qSJQGCsqTkOCFVQ2FmnjvNvv3lYW8PZGe/elSweOxbR0Xj2DAYGCAhgIy4hRDyoEco+X1+sXIm3b8uvlMt518q2Y8fagwdHzJhxbeXKlHv3zqxZ0yI+HlpaMDGBuzs+fy5WXK8eAgKwZw+8vNC7d4Wm9yGESD9qhLKvRQv8/jsWLqxQMc279p0ePex9fVe4uc3U0dEBoK0NT09ERuLTJxgYwMur5MTc3bohKgp2dujYEQsXIi+PndiEkKpCjbBaWLwYZ89yw8LKr6R51yqmYUP4++PiRdy+DQMDbNtWrOEpKmLOHERHIy+PPjIkROZRI6wW1NXh7c3bsqVCxTTvWoW1aoXAQAQF4cABmJoiKOjrtAQMw0RGRsbF/efqmlSjBtspCSGVQ42wuhg2LHvXrooW07xrP6J9e5w/Dz8/eHqiQwecPJnarp19jx7rhgw5a2TUe9eu/flls2dj3jx8+cJuWELID6NGKJe+zbv26ePHuLi47OxstgPJAPPNHIgAACAASURBVDs7RERgzhyMHLk0Onrihw8Bycme799fdnNbn5iYCGDhQmRloWVLBATQfHaEyBJqhHJKOGnS1Ph4gb7+Klvb9o0aHdm/n+1EMoDDgYMDdHRuMczAb/v4WVk9bt++DUBbG76+OHoUmzbB1vbrmsCEEOlHjbDacXHB3bvlVm319a2VlnY1NTXg/ftrCQneLi6vaC2+itHQ0AA+FtmRULt27YINCwvcuIHx49GjB1xcSt6AQQiRQtQIqx0TE0yZUu57c2cPH3ZOT8//WhVwSEm5cvmy+MNVB+7uTlpazkAiIORyj2ZmPvDyavPpU2EBlwtHR9y9i7Q0GBsjNZW9rISQCqBGWO2MH4/cXPzzT9lVqqqq6UU203g8VTU1seaqNhwcBvj7/25sPKxRI+vhw889enSoXj1eu3Y4f75YmY4O/voL4eFQV2cpKCGkYqgRVjtcLjZtgptb2e/KjZg2bYGmZhoAIA44WrOmra2tRPJVBw4OA2Nizr94EfbPPxsbNdLx9YW/P37/Hc7OSEsrVlm3buHX9+7R9HaESCNqhNWRuTl698Yff5RR0nfAAIGXV7fGjS319OZoau4aOlSTFhyqhB49vo6OMTNDae8xBwejVSvs2kVjSgmRLtQIqykvL+zdW/bIxTFOTreePw9/+/ZUdLTpP//gxQuJpauWNDTg7w8/P4wcCRcXEVP3zJ+P4GBs24aOHXHnDhsRCSGiUCOsprS1sWBBhZbtBdCoEaZNg6urmDPJhT59cOcO/vc/mJuL6Hbm5rh+HZMno08fODoiMZGNiISQ4qgRVl9TpuDjx4rezubujqgohIaKOZNc0NXF4cNYuhS//gp395KLRXI4cHREbCw0NGBqir/+YiklIeQbaoTVF4+HO3dgalqhYmVleHtjxgyajLuqODggKgqxsbCxQXw8ADx8+DAgIOD06dO5ublaWtiwAaGh+O+/kuNrCCESRo2QfNOvHwwM4OfHdo7qQ08Px49j6lR06QJ7e59OnaZOmJA8bNhZI6POHz9+BNCmDfbvB923Qgi7eGwHINLE1xdWVhg2DI0asR2l+nB0RLNmL7p2PZGT8x/A/fwZX74cnTlzWUCAb4nK5ctRuzacnaGgwEpSQuQUXRHKgQcPMGVKhSqbNcOUKXBzE3MguZOYGKmkZFfw45aX1/vmzdvflw0ZgqAgWFri+nXJ5iNEvlEjlANNmuD8eZw8WaHiBQtw61bJWVJI5dSvX5/Pf15kx3OGqVtiEA0AIyNcuIBlyzByJBwdkZAguYSEyDNqhHJASQkbNoi+te17KipYuxbTp5cc7EgqwdzcvFGj50pKO4CPQLSm5oSGDWe2aFFy4ft8AgFiY9G0KczM4OuL3Fw2EhMiT6gRygd7e5iawsenQsUDB0JfHxs2iDmTHOFyuZcv/zt9+v/MzIb36LEmNNT7/PlO+/dj3z6YmSEoqGS9qio8PHDxIk6cgIUFMjLYCE2I3OAw1Wi6p9DQUF9f35CQELaDsCM9PZ3P53O5pfxx8/IlzM0REYHGjct/rMeP0aEDoqNRr17VhpSA1NRUddmZ6PrcObi6gs/HqlUQOdtrWBgsLcHh/NjDytaTICZ5eXlZWVmqqqpsB2EZvRjKRVeEcqNRI0ydWtGBMM2bY+JEzJsn5kwEdna4fRuzZ2PCBNjbIzKyZIGVVWEXfPuW3rEmpOpRI5Qnbm4ID6/oQJhFi3DtGi5eFG8kAnC5cHBAXBwcHNC3L4YOxePHhUczMjISvg2bWbUKbdviwgV2chJSXVEjlCcqKli3DuvWVahYVRVr1mDaNLoGkQxFRTg54eFDmJujY0c4O+PFi+whQyY2bmzbuvX4Jk2sbt0K9/PDqlUYPx4CAV6+ZDsxIdUFNUI5068fDh+uaPGQIWjYEJs3izMQKUZNDW5uiI+HlhZatvzz2LGWiYm3EhKCnz8/PGjQpIyMDIEA9+/D3Bxt28LDg2bEI6QKUCOUP0pKP1Ds54eVK/H2rdjSEBFq1YKnJ/T0zuXmTvu2r2FWlk1kZCQAFRV4eCAsDBERaNMGZ8+ymJSQ6kA2GmFWVpaKioqFhYWFhcWsCi4tRKpEixYYOxbu7mznkEdKSlxAWLApFOYoKioWbDZrhuBgrF6NZcvo/gpCKkVEI8zNzT137tz//vc/yacpg76+fkRERERExLoKfsRFyrZwIcLDK1S5dCkuXcKlS2IOREr67beBNWp4fOuFd1NSIgIDW2dnF6sRCHD1KlRUWIhHSLUhohG+f//e3t7+2bNnkk9ThqSkpFmzZnl7e3/+/JntLNWCsTEmTRIxr8n3VFXh5YVp02iOEwlbsMBlwgTlOnXMdXWt27Wbd/Hi3kePlCwty1rdftMmyOtttIT8PBGNsHbt2qqqql++fJF8mvT09D3FBQYGAlBQUPDy8ho0aFB6erqtrW1eRX59k7KNHImaNSu6LOywYdDVxdatYs5EiuFyuevWLX/3LvLNm2u3b4fa2BgePYpFi76u9ytymEyrVnBxwcCBePlSNj71IEQaiPhpUVJScnZ2Xrt2bZbER6QJhcLPxaWkpADg8Xhjx47t3LnzkiVLFBQUHj16JOFg1dPGjViyBB8+VKh4wwZ4eNCoGVYoFFmWKX+930ePYGEh4r3tbt1w/z5sbWFrq+rhgcxMieYkREaJXo9QQ0MjKiqqRYsWPXr0qF27dtFDnp6eP3GajIyMLVu23L59+/3793v37tXT08vfLxQKly9ffvDgQT6fP2fOnNGjR8+YMeP7bxcKhfkzh6Wnp797905LS+snMpCSjI0xbBiWLKnQDRJGRvj9dyxciJ07xZ+MlEVPD4cPIygI/fphzBgsWwZl5cKjiopwcUGvXmkrV9YwMcGqVXBwYC8rITKBEcXAwECrFCLry5WYmDh27Nj8cS5Pnz4t2L958+ZWrVrFx8dfunRJS0vr1q1bIr89ODi4Y8eOI0aMMDQ09PT0LO0sISEhvXr1+rmE1UBaWlpeXt6Pfc+nT0zdusydOxUqTklhGjRgrl//iWySlJKSwnYECUlIYIYMYZo3Zy5dKnko/0k4c4Zp2ZLp25fJzmYhHutyc3PT0tLYTsE++fmJ+GkSnXQ7KyuLz+c/ffq0SZMm+XvatGkzd+7c0aNHA5g5c2ZaWtr27dtFfm9SUlJiYmKDBg1q1KhR2uPv379/4cKFgwcPLtgzYsQIY2PjKv1HSK9yJt0uhcLu3dzdu3MuXKjIvM6J27YtXLQoRk1NQVFxhLPz5NmzOT86G7T4ydsUwydPcmfM4PXqJVy9OldN7evOgichOxvHj3MHDRL+4OuiOqBJt/PJ209ECUpKSuX+mhL91qhkMAwTFxfXrl27/M127dptLX04hra2tra2dtkPqKCgoKioWPSNUxUVlR9tDLKL+80PfRfz+++c3bsVbt1iOnYsu1IoFI7aunV+SkqvlJRMYKan5zYud/KcOZWILBY/8STINIEANja5CxcqWFoq+fvnde6cd+DAoeDgKwYGDSdPHlOnTp2hQ1EwGiAzE3w+q3EliGEYeXsxiERPQrlKbYTZ2dnHjh2Liop69epV3bp1TUxMBg8eXLV/W6WkpGRnZ9esWTN/U1NTMzExsTIPWLNmzaZNmy5YsKAq0skeRUVFRUXFn3nFX7rEq8B0Mw8ePGjw4UMvAAAfWJeaar9jxwzpu9c+/3lgO4VE6epi+3acPIlx43gM89vnz41SUx0VFJ7u3Cm4fv1IwRswAMaOxYcP8PODgQGLeSWEy+UKhUJ5ezF8Tw5/In6U6Eb47t27nj173r17V0FBQVtb++PHj7m5uUuXLj19+nSLFi2q6tzq6uo8Hi8tLS1/MzU1tVatWlX14OQHVGzStZSUFA1h4UQnfCCbZrqUJn36ICDgdu/e3MxMLwB5eR3fvWvo6up56JB/Qc3u3di0CR07YtQo/PEH5PgNM0IKib56mDx58suXL//555/8JWAyMzOPHz+ek5OT/2FelZ2by23cuHF8fHz+Znx8fNOmTavw8UnVMjU1vcnjJX3bDOLxLG1s2AxEvvPu3SOGaVdkh3ls7IOiBTweXFwQE4NPn9CqFQICJByQEGkkohGmp6cHBwdv2LBh5MiR+RfUCgoKAoFgz549YWFhL168+LkzJScnf/r0CcDnz58/ffqUP0jH0dFx/fr1WVlZCQkJe/bscXR0rMS/hVTO//6HiRPLOM7n83327u1Rr95EHZ2B6uq7VFVX0S32UsbQ0FBdPazIjltZWa3evy9ZpqeHgADs34+1a9G1K2JiJBiREOkjohF++vQpNze3TZs2Jfbnj2p5//1PVcV07dq1U6dOTZs2HTx4sIWFRWZmJoC5c+fWrl27Tp06hoaGv/32W+/evX/uwUkVqFsX4eE4dKiMEls7u+tPn045fdrnzp2Qhg1rXr8usXSkItq2bdupk6qGhgtwjsf7S09vbs+eC4yM4OEhYmLuzp1x+zYGDcKAASgxhSkh8uX7Oyryl3pYt25dif1BQUEcDuft27dVfg9HZmZmbm5u5R+H7iP84fsIS7h6lWnUiPnypULFJ04wrVoxOTmVOqMYyPldU0Kh8MiRo+PHz/HyWp+UlMQwzPPnzG+/MQ0bMv7+TFX8nMkMuo8wn5z/RFSEgoeHR4nWqKCg8OLFC29vb6FQqKurq6io+OrVq4CAgDlz5nTv3t3JyanKmzGPx6uS0b2PHz++detW1X6QKUNycnJ4PF6lbuxr1AgREYiNRdeu5RcbGOD4caSnw9Ly588oBtnZ2cpFp1qRMxwOp1Wrll27dujWrYuKigoATU0MHAgbG6xbB29v6OqijHtrDxwAAF1dScUVJ4Zh8vLyaMCknP9EVIjI9piWljZgwIASlb/88ktiYqJk+/SPoSvCyl4RMgzz9i2jo8M8eFCh4vv3mbp1meTkyp60StHfv0zpT8LZs4yZGdO9OxMRIfob//6b0dVlZs6Utv+rP4OuCPPRT0S5RF+Hqaqq/vvvv1FRURs2bFi+fPn69etv3rx58eLFEvOOkmpITw9z52Lu3AoVGxmhb1/8+aeYM5EqY2eHyEg4OqJ/fwwdiidPCg8xDJOVlTVqFOLiIBTC2BjbtqHI/TKEVF/f98b8Wa1DQ0Ml3pUri64Iq+CKkGGYrCymZUsmOLhCxQkJjI4O8+hRFZy3itDfv0wFnoS0NMbTk6ldm3FyYl6/zhk3bpaOTmtd3U5mZt3j4uIYhrl9m+nQgbG0ZEqZAFgG0BVhPvqJKJeIK0JVVdXk5GS1glkLibxRUoKvL7y9K1Ssq4uZMyGvs/nILlVVuLkhPh5aWmje/M+9e7UTEyPfv7969+763r0dc3Jy2rXD1auYOBH9+mHGDLo0JNWZiEaorq7evXv348ePSz4NkRY9eiA0tKLFs2cjPBxXrogzEBELbW14ekJbOyQnxw3IH2ZlkpZmFhsbC4DLxcSJiI1F48bUCEl1JnqKtZkzZ44fPz45OVkgENStW7foQERzc3NJZSOsqvgwMz4ff/4JFxdERIDm9pVBvOK/BvLyuHl5eQWbtWpB+mZWJ6Qqif61NX78+ISEhO3bt/fr18/S0tKiCAnnI7Jh+HCoqmLfPrZzkJ8hENjx+Ru+bT1KTo7891+TnBzRxaNHY9IkJCWJPkqILBJ9RRgYGJhNU00QAN7esLZGuXOKcjhYuxZDhmDgQNCny7LG23vxu3fTL1605HI1NTWzfHx2+PkpWVpi+3YR94hu3IilS9GqFRYuxLRpUFBgIzEhVUpEI0xLSzt16tTw4cO/n2WNyJ3mzeHsjKgolHtXcvv2sLGBjw8WL5ZIMlJllJSUAgP9c3Jy0tLSNDU1Afz6K4KC0K8fhg7FypUouhi2piZ8fTFuHKZPx5492LgR5S1kSYi0E/HWaGpqqpeXVxatsEMADBgAfX1s3Fih4jVr4OuLV6/EnImIhaKiYn4XzOfggNhYZGaidWucPl2yuHVrXLqE2bPh4IAJE1C5hUQJYZmIRqirq1u3bt0nRW+1JfLM1xd//om3b8uvbNAATk5YulT8mYgkaGnB3x87dmD6dAwdWrLbcTgYPRpxcahZE61bIzeXpZSEVJqIRsjlcteuXbt48eLIyEjJByJSp3lzjB9f0TsFFyzA6dOIiBBzJiI5traIikLTpl/nmimhZk34+ODJk5JDTwmRIaJfvLt27fr48aO5uXmDBg10i8+/G0G/4+TQkiUwMsKVK+jcuZzKGjXg4YG5c3HxoiSCEYlQVYWnJ4YPx8SJ+PdfbNkCff1iBSoqhV+fPw9jY9SpI9mIhFSC6NsnGjdubG9vP2TIEGtr66bFSTgfkQqqqli1Ci4uKHJ7WanGj0dqKo4eFX8sIlFt2uDGDfTqBSsreHnh4cMn3boNbdDAsmXLLocPF/7vjoiAqSn8/OjNUiI72J7jrSrRXKNVM9eoSEIh07UrExJSoeLz55mmTZnMTHGFKRPNrMiI+Ul48oTp2jVFSckCuAUwwIdatXqFhp4tKHj4kOnVi2nZkjlzRnwpykdzjeajn4hy0TwgpGI4HISGokePChV37QojI2zeLOZMhB1Nm2L27Mtcbm/ACgCg/fHjWh+fPQUFLVogJASrV8PJCQIBjSMm0q5YI3RxcSk6xai/v3/RsaOnT5+2srKSXDQibZSUfmAGNW9vrFqFDx/EGYiw5vPnzwyjVWRHrY8fk0vUCASIiYGZGdq1g4+PJNMR8mOK/V47evTo/fv3878WCoWTJk0KDw8vOJqUlFR0k5CyGBpi+HAsX852DiIWnTp1rFnzX+DrPGwKCgGvX3e9dKlkmZoaVq7E9et4+ZI+MiTSi94aJT8oKQnOzmCY8iuXLcPBg/j2pxWpTvT19ZcvH6+r21Fb20VH51c7u6jNm2eMGwd7e8TGlixu0QLr19P9FUR6USMkP0hLC9HR2Lu3QpXu7nB3F38mwoJJk8Y8fHju+PFh4eFbQkP3DRjAi4tD37745Rc4O5c118zkyfDyAk1mTKQHNULyg7hcbNoEd3d8/lx+8bRpePgQZ86IPxZhgYaGRseOHRs3bpy/qaQEFxfEx4PPh7ExvLwgcqJGV1dcvQozM3pdEGlBjZD8OHNz9OmDFSvKr1RUhJcXZs2iD4jkh7Y2fH1x+fLXbhcUVLKgaVOcOIGNG+HiAoEAz5+zEJKQojhMkQ97GjduzOPxCv6+u3DhgrGxccHMMgkJCbGxsUxFPhxiSWhoqK+vb0hICNtB2JGens7n87mSWRo3KQnGxjh7Fqam5dbuatVqa2KiUEmpbuPGa3bvNjQ0FGu01NRUdXV1sZ5C+knJk3DuHGbPho4O1q5FwWI2QqEwJydHWVk5OxtbtmDlSkyZAnd38PlVfPa8vLysrCxVVdUqflxZIyUvBmlW7Jdmfgt88U3Tpk0zMjIKNjMzM2lmGfKVtjYWLsSsWeUWHjl4MOT16wtJSeFv3/5x8+Zoe/svX75IICCRBnZ2iIrCmDH49Vc4OuLly+zffpuup2feuHH3du16vnz52MUFYWG4exetW1do2iJCxKHYQK7Lly+zlYPInilTsGMHDh3CkCFlVO3buNE7NTX/b3IzoMfnz9euXevZs6dkMhLWcblwdMSgQfD2hqHhstzcRrm5kQASEiJ79frt4cNr+vrcI0fw7Bmt8UtYQ58Rkp+loAA/P3h5lV2VlpZWZFVX1MjNpStCOZQ/GXutWudyc+d829c2NbXZo0eP8jeaNCksjolBRoakExJ5Ro2QVEKXLrhypewS+4EDt6up5X+dBhyrUaMjrWgurxQUGIBTsJmXx80T9X7oX3/BxARFJrkiRLyoEZLKKW+Ew4z58x/Z2XXV0Rlbu7YNl+u+YEHdunUlE41Im169flFW3vJtKzY5Of7gQcPvb7FYvx47dmDBAtjZIS5OshGJXKJGSMSLx+PtPHo08P792efP31yyZEBYGNuJCGv8/Fb8+uttHR2rOnXsjIymnDu35/59BVNT/PdfyUpbW0RGQiBA585wcQG9m07EimY9IlVhyxYYGqJbt9KO6+jo6OjooGlTtGiBqKjCofREnvD5/CNHdmRnZ2dkZGhoaACwtcWJExg/Hl26wMcHtWsXFisqwsUFQ4bA1RUmJli/HgMGsJacVG90RUiqQosWcHJCZmY5ZWpqcHfHkiUSyUSklJKSUn4XzCcQ4P591KsHIyNs21ZyFtv69bFvH3bvxt69EAolHZXICWqEpCrY2cHMrEJr7UyahPv3cf26+DMRmaGmBk9PnDmDHTtgayvic0FbWxw+/AOLgBHyQ+iVRarI+vVYtw4vXpRTpqSExYtpJm7yvTZtcOMGRo1C585wdxc9T2m+BQtw8KAEk5HqjhohqSKNGmHaNMybV37lb78hKQmnT4s/E5ExXC6cnHD3Lp4+hYkJzp37uj8vL+/58+cZ3+4uFAiwahW6dxex5BMhP4EaIak6bm6IiCi/wykoYMUKuLnRZz5EpHr1EBgIHx9MmIChQ7F9+78NGlhYW89t3LjL5MnuDMN06IA7dzBmDLp3h4sLUlLYTkxkHDVCUnX4fPj4YPZs5OSUUzloEFRUcPiwRGIRmSQQIDYW2trPJ03yfvfuakLCocTEsH37MrdvD8C3mdvyV302MkJAQIXWiiZEJGqEpEr17w99fRw9Wn7l8uVYtIiWZyJlUFVF+/aXFRSGA/mTE3FSUqbv31+4vEytWvD1xaFD8PNDt27UC8lPovsISVU7erRCw/vs7dGgAfbswfjx4s9EZBWfr8zjZRV5iyFLQUG5RI21NcLCcPUqOBwQ8hPoipBUNUXFiq4j4OmJpUtpfmVShu7du2loHATyRyOn8/mLb98e6uWF7OxiZVwuunQp3ExMpKtD8gOoERL2WFrCwgL+/mznINJLR0fn2LHNxsZjdXUtGzTo4uPT5+bNPleuwNQUwcGlfte0abCxQVQUXSGSCqFGSMQjPR0TJ5a/1urKlfD0pGF/pAxWVpYxMecTEsJfvYqYPHmcoSGCg7FpE9zdYW//dbxMCfv3Y+xY9O3LnTVL6dMniScmsoYaIREPVVU8eYLt28spMzZGjx5Yv14imUj1YWeHyEg4OMDODs7O+PCh2FEuFxMm4N69PABGRvjrL7pVh5SFGiERmw0bsHRpyV9R31u+HH5+eP9eIplI9aGoCCcn3L8PPh9GRvD1LTkGWUsL69Zlh4Rg925YWuLGDZaCEqlHjZCIjbExhg/H4sXllOnrY/hwrF4tkUykusm/g+LSJYSGwtQUIYX3ViArKys1NbVNG1y5ghkzMHs2jaAholEjJOK0YgWOH0d4eDllixZh9268eiWRTKQaatUKISFYvRrTp0MgQExMWt++jvr6Xc3MRpiYdH306OGYMbhxg+6vIKJRIyTiVLMmVqyAi0s5f4rr6WHiRKxYIalYpHoSCBAXBzs7WFi4h4TYJibeev/+/P37vr/+OoYp/gpctw43b7IVk0gdaoREzMaOBcMgIKCcMjc3HD2K+HiJZCLVVv5yvlpaN4XCcd/2mX350ujp06dFyxo1wpAhGD8eiYmSz0ikDjVCImYcDjZswKpV5ZRpamLWLCxbJpFMpJrjFZ8yKzdXqFB8kofBg/HgARo2hIkJfH3Lv82HVG/UCIn4WVggLKz8spkzcfUqIiPFH4hUcz17/qKsvOHbVtjnz2/XrdNPTS1Wo6YGDw9cuYKQEJib4+pViackUoMaIZGImjXLr1FRgbs7Fi0SfxpSzW3c+MfAgbE6Oua6uh0tLDzCwvZmZqJVKxHv0BsYIDQU8+djxAiMG0fDSuUUNUIiTZyc8OABLl5kOweRbXw+f//+La9f34yPPxkefqpt2yb+/ti3D6tXQyDA8+cl64cNQ1wcunenYaVyihohkaCAgGL3eX1PURFLl9JFIakSXC5XWblwqYouXRAZCTs7WFnBwwNZWcWKa9TAqFGFm/SpoVyhRkgkqFkzODsjLa2smlGj8OULTp2SVCYiR/LHlN65g5gYmJnh3LlSK/v0wciRePNGguEIe6gREgnq1AldupQzgpTLhYcH3N1pdkgiJg0a4NAheHtj4kQ4Ooqe3e/ff2FgAFNTEdeOpPqhRkgky9sb27bh4cOyagYMQI0aCAqSVCYijwQC3L+Ppk1hagpf369/dz148CAoKOjWrVsqKvDwQFgYIiJgZobQULbjEnGiRkgkS08Prq6YMaPsqnQPj5VTpvQwMxvWteuF8+clE43IG1VVeHjg7FkEBsLKCsOGLercedaYMc/69NnaqZMgJyenWTMEB8PbG1OmQCDAixdsJybiQY2QSNysWXj9uoxlVRmGGb58uWZy8sF795ZdvPjnkCGnjh+XZEAiV8zMcPUq+vYNO3QoLjHxVEbGvKSkXZGRtmvXbs4vyL92NDeHnR27SYm4UCMkEsfjYf16zJyJzEyRxx8+fKjy8OFUoVALaAn8/enT+nKXsCCkEjgcaGiEMUy/gj0ZGf3PnCmcjTT/ndJHj9gIR8SPGiFhQ/fuaNsW+/aJPPjmzZvGRVaW0wWSk5MllYzIqQYN9Pj8/xXZ8UpNTa+M+v37RdyPSGQUNULCkr17i923VYSZmdllHi/n2+Z5Dse4dWuJ5SLyqXfvXnp6x7ncU0AWEKmu7n7z5sTSxpQCSEyEpSVWrCjtfQ0iS6gREpbw+Shys3NR2traExYtstPRWaek5Mbjza9de9W2bRJOR+RNjRo1rl07OmJESPPm3eztvS9f9n/0yKhePZiZwcsL2dkl62fMQFQUHj2CiUkZn3cT2cBhqtHkeqGhob6+viFlz11SfaWnp/P5fC63mvxx8/z58+vXr2s9fdr94EGlu3crPvlVamqqurq6WLNJP3oSAOTl5WVlZamqqlbmQR4+xKxZePIE69ahd28RBWfOYMYMGBhg/Xo0bVqZU4kLvRjKVU1+aRJZlZuLiRORk/P9EX19CQ+dxQAAIABJREFU/ZEjR/ZetEhJSYlu4yJsMTDAyZNYswbTp0MgQPFlDQGgRw/cvYuOHdG+Pa0tLauoERJW8Xh48wYbNpRVM3cuvLwkFYgQEQQCxMV9nafUxQUlVnRSUoK7OyIjoajIUj5SOdQICdt8fbFqFd6+LbVg6FD873+4cUOCmQgpKX+e0nv3kJkJIyMEBBSu2fTu3bvnz5/Xr8+4u7MakfwsaoSEbc2bY8IEzJ9faoGCAmbNwpo1EsxEiGh168LfH0eOYOtWtG+P06eTrKx6t2491tratXnzDjExMQWV/fph4cJyZpgnUoIaIZECixfj4kVcvlxqwfjxuHkTsbESzERIqSwtce0apk1D//4uERHT3r8PSUgIevo0oH//CQU1/v549QqGhiJWAybShhohkQKqqvD0xMyZpa4Cx+dj8mT4+Eg2FiGl4nDg6AgtrQcM0+fbPoO0tPovX77M36hbFwEB2LcPa9eiWzfcv89WUlI+aoREOgwfDk1NbN9easHUqfj3X7x6JcFMhJRDQaHYZk5OuoqKStE9Xbrg9m307w9bW7i4ICVFovFIBVEjJFLDzw8rV5Z6tFYtODqWM76UEMkSCLrz+asABgCHc+LLF8yZo5OYWKyGx4OLC2Ji8PkzLC3ZyUnKRo2QSA0TE9y7V1bBnDnYuRM07yiRGr6+y8eO/Vy3roWenkXv3kEPHwaUWOCwQJ062L0b4eEsBSVlokZIpImmZllHGzRAnz7YulVSaQgph5KS0ubNnm/e3H77NuLkyYDGjXU8PHDmDIKCYGWFsLCS9TVrFn594QK9UyotqBESmeLuDl9fZGSwnYOQUpmZ4coVzJiB/v3h7IzPn0WXnTqFVq2wdy+q0TSXsooaIZE+x47hxAnRh1q1goUFDUgnUi5/TGlcHPh8GBsXu/u+wJo1OHECmzejfXt6y5Rl1AiJ9GnQoKw/pN3csHp1qTdaECI1NDXh64vjx7F5M2xtRdxB0a4drl/HtGkQCODoiKQkNlISaoREGpmbQyDAsmWij9rYoG5dHDki2UyE/KT8bjd+POzsvs5TKhQKDxwIdHFZ7O+/IzMzw9ER9+9DVRWmpti7l+24cokaIZFKq1Zh375SB5HOm4c//6SPVois4HLh6IiYmPx5SoWGhv2dnKL8/DrNnPnJ1NQ2JSVFWxtbt+LECdy8yXZWuUSNkEilWrWweDGmThXd7QQC5OTg/HmJxyLk52lrw98fU6eefPrUIDX1T6BXZubcFy+cV6/elF9gbo5Nm9jNKKeoERJpNWkSUlMRFCTiEIdDazMRGZWVFScUWhds5uZ2CA8XPYnumDHYvr3k/YhEHKgREmmloIBNmzB3rugJ/EeNwsOHuH1b4rEIqRRT0+ZqalFFdkTyeC1EvvExezb27IG1NY0pFTtqhESKdewIW1vs2CHiUP7qcLQ2E5E1AoGgWbMbKireQASPt7NOnbXv3k3t1k3EB+KtW+PKlfw1LuDoiPfv2YgrH6gREum2fTsmTBB9yMkJFy7g8WPJBiKkUhQVFcPCQlatquHgsHPx4qT4+P/Cw7XHjkXPnnB2xocPxYrz70d89AgFM7fRfUPiQI2QSDdlZaiqij6kpgYnJ1qbicgcZWVlF5dJgYGblyxx1dTUzB9Tev/+17vvv+92amrw8MC5c/j3X9jaspO5eqNGSGTZ9Ok4cADv3rGdg5DK0tKCry9On8bhw7CywrVrJQtMTXHhgugPCkglUSMkMmLGDGRmltypq4sRI7BxIxuBCKl6bdrg8mV4eGD0aAwdWnL9TQ4HBgaFm/HxyM2VcMDqiRohkRFv3sDbW8T+uXOxbRtSUyUeiBBxEQhw/z6MjNC2LTw8kJUFALm5uXfu3Ll161bGt0nnPTxgYYGrV9mMWj1QIyQyYt06+Pri+fOS+5s0QffuZS1tT4gMUlWFhwdu3kRsLExNsXPn85YtbXr23Ny3b0CLFjZhYeEADhzAihUYPRoCQclrR/JDqBESGdGwIWbMgKuriEPu7li3DtnZEs9EiHg1b47AQGzejGnTpj95svHDh78+fNj0+vWxoUOn5hcIBIiNhbk52rSBhwf9EPwkaoREdsybh+honD5dcn/r1jAywr59bGQiROy6d2dq1kwALL7taJCZWffNmzf5G/nXjtev48YNmJvj0iW2YsowaoREdigrw9sbs2YhJ6fkITc3eHnRbFSkWuJwOAoKDFA4/UxW1qeaRVe7BwwNcfo0li3Dhg0Szyf7qBESmdKvH5o0gZ9fyf3dukFTE8HBbGQiROyGD++jpuYOZANCRcVNubmNBgyoER9fsmzQIBw6xEY+GUeNkMia9f9v787jocz/AIB/xjDuI1euROTuQHSQW6VSkV93tkO3q4NUW8l26dqI7rLp2G3b2m6i3ZC0KSIKtZEj95kMY4z5/fFsk6SSYx7G5/3qj3m+z9fzfDyN+czzfb7HQQgMbKPb+Nq1sHMnGQEh1O327Nns6yunrm6hqjp6yZKC/Pxjjo5gYQH+/m2MKuLw9YW//uJilL0WJkLU2wweDJmZwM/futzZGSor3925k5eXx8J5qBBvoVKpW7asy87+Jzf3UWjoLikpUS8vSEmB7GwwMIDbt9v+KSsrWLIEfvhBGPuUfh0mQtQLSUt/XtbU3LxUXNxhypSt1tYmqqp3IyO5HxdC3KSkBOHhcPw4rFsHjo5tjC2aOBGeP4fhw1nGxh/HI6LPYSJEPOLnn35Sy8x80Nj4R1nZX4WFGxYsqKioIDsohLqdjQ2kpoKdHZiatpHthIVh9erG5ORv3Dv2cZgIUa919y5cvszZunPlyko6nXjdD2Ainf7o0SOSIkOIq4hFyTjZ7vPWEBUVCA+H4GDw9oY5c8gIsWfDRIh6LQUFWLUKPtz2CQoKtvwq3MDPLyQkREpcCJGCyHYHDsDKleDoCHl5UFlZuXr1Fju7uUuX+hYWFjo4QFraF5c168swEaJey8AAZs2CTZuIrVnLlm2UkCAm1kgH+EtUdOTIkSRGhxApiHlKjY1h+PD3GhoTQkJ0EhNDT540NzGZXFxcLCgINjYfK1dVkRdoT4KJEPVmAQFw8yYkJgLA/CVLhvr5jVZWHiEn50OjnQ0MFBUVJTs+hEggLAz+/rB+/eV372Y0Nc0BGMBmTykq8g0Obr2Gk5MTODrC69ekhNmDYCJEvZmEBPz0E7i7E3PKeG3YkFRQcO/164gjR3SPHiU7OITIVFmZ19w8mLPJZmtnZOS2qhMdDXZ2MHIkeHnB+/fcja8nwUSIerkFC4CfH8LDPyn84QeoqGijzwBCfcbYsYYSEh+H01Mo0Q0Nhq36lBK9bJ49g6oq0NVt/WfUd2AiRL0chQIhIbBxI1RXfyykUmH7dli/HmcfRX3WpEmTRo8u79dvCUC4hMRqA4O7MjKLDQzg1q3WNYnxiOfOwd69YGsLL16QES6pMBGi3s/ICKZMgdDQTwqdnEBUFC5eJCkmhEhGoVAiIy9cuTJ3x46SX3+1T0mJPHeOdvQo+PqCoyNkZ7eub2kJyckwaRKsXElGuKSisNnsb9fqJSIjI4OCgiIiIsgOhBx0Ol1ISIiPr09+uWlsBAYDxMUBoLa2VlxcHAAgLg4WLIDMTKDRSA6P6z5ehD6MxWIxGAwRERGyAyFZqzcDkwmHD8NPP8HcubBjB4iJkRhaT9EnPzQR76HR4PPPfQsL0NKCU617yiHUlxHPBdPS2vVc8MAB+HyNC96DiRDxtD17ICAAamvJjgOhnkVREcLD4fJlCA0Fa2tIS2u7mqgoWFiAry+P9ynFRIh4i68v5cNEawAAQ4eClRWuVYpQm0xN4eFDWLgQxo8HLy+oqQEAePfu3fXr169cuVJZWblsGWRkAIMBOjoQHg489CTtE5gIEW8pLqbt3ftJyc6dcPAg4ATcCLWFjw9cXSE9HQBASwt8fJJ1dKznzn06f36mnt64v/+OkZGBoCC4ehVCQ8HK6ov3jr0aJkLEW/bsEThzBl69+liirg4uLhAYSF5MCPV00tIQFASRkRAaurqo6PL791vp9I0lJRE//LCWqDBiBCQkwMyZYGcHGzaQG2zXw0SIeIuCQuPq1eDu/knh5s1w+jTg4qQIfdWQIU3i4o0Aah8K5JhMhdLSUmKDSoWVK+H5c9DVJSm+boOJEPGaxhUroKgIrl//WKSoCMuWQUAAeUEh1Avw8/NTqY0ArA8FbDq9VFKyX8s6srLg6sr90LoXJkLEc/j54eBBWL0aGho+Fvr6wvXrkJFBXlgI9QJLlsyUkFgOUARQLizsIyFhNWaMwMOHX6xvagqenp9M69QbYSJEvMjGBoyMoGWvGUlJWLcOtmwhLyaEegF/f5+ff7Y0NnYbNmzurl2Dc3N3+/nB9Ong6gplZW3Uv3MHKBTQ0oKgoF48oSEmQsSj9u+H4OBPRj95eMCjR/DPP+TFhFBPR6FQFi2a9+TJrZSUO15ey6hU6v/+BxkZ0K8f6OtDUBCwWJ/U79cPgoIgIgJ++w3MzCApiaS4OwcTIeJRqqqQmvrJ/FFCQrB5M/j5kRcTQr2SpCQEBUF0NFy6BKambXyZNDaGhARYuhQmT4YVK6CykowoOwETIeJdSkqtSxYtgtJSiI4mIxqEerdhw+D+ffDyAicncHWF8vJP9lIosHAhZGSAgADY2ZEUYkdhIkR9CZUKAQHg49OLn2YgRB4KBVxd/2sp1dNr47mglBQEB0NyMknxdRQmQsTrHj6EX3/9uDl9OggLwx9/kBcQQr2blBQEBUFU1H8tpY8eQXp6uoXFdFVVkzFjpiZ/mgZ//bX1vWMPhIkQ8To5OfDygqKi/zYpFNi9G378EZhMUsNCqHcbPhzu3wdPT3B0fGtquuD+/a35+Y8fPgx0cFjx+vVrTrX0dNDXh8OHW/ey6VEwESJep6kJS5bA+vUfSywtQU0NTp8mLyaEeAHRUurpeaWhwR1gKAAA6JSWrj958mMbzI4dEBcH167BiBHw4AFZkX4DJkLUB2zaBHFxEBf3sWTXLvjpJ2i5TgVCqEOqqirYbLkWBfJv337SGKqtDXfuQEAAzJ0Lrq5QUsLlAL8NEyHqA0REIDAQ3N2hqem/EmNjGDMGl2dCqPMcHMz79bvE2aRQLubmjv0wQelHjo6Qng4KCjBkCBw9ytUIvwkTIeobZs4EOTk4fvxjyY4dsG9f7xvxhFAPY2dnN2OGgpycnaTkZjm5CS4uzdbW04cObWP0vZgY7NkDcXE9bko2CpuHVlqMjIwMCgqKiIggOxBy0Ol0ISEhPr6+/uWmtrZWXFy8jR0vXoC1NaSng9yHZpxly0BaGnbt4mZ43PHFi9CXsFgsBoMhIiJCdiAk486bIS8v78WLF9ra2urq6gCQlQVeXlBSAocOgbl5d5+8s3ruh2ZdXR3z0359paWlDS2nUUbou+jpwdy5cODAx5KtW+HECSgoIC8mhHiEqqrqhAkTiCwIANraEBkJAQEwfz7MmPG1NdDGj4cDB0juxM1P5sm/IDY2dvny5QUFBb/99tukSZMAgMVizZw5s6mp6c2bNwEBAVOmTOnC09HpdAaD0YUHJEt9fb2goGDvuiMUExMTEBDg3vl274aamo+bSkrg5la0ceOT//1PRkZm1KhRvevqIdTDOTqCrS3s2QOGhuDuDhs3Ao3Wuk5oKHh5wbFjEBwM48eTEWXPTITGxsapqakLFy7klNy8eVNAQOCPP/4oKSkxNzefPHlyF35gKSsrUyiUrjoaaj82my0qKlrAzRsyGu1juygAAJyRlz+2d+/4K1cKhYQ2KChci4+XkpLiXjwI8ToREfD3h3nzwNsbhgyBoCCYMOGTCpqacOsW3LgBK1aAvj6EhMDAgdwOsicmQrGWEyUDAEBCQoKtrS0A9O/fX1xcvKCgQFVVtatOV1dXV1dXx9X7EgQAANXV1ZyGFFKUlpYe2b07rrmZVlcHdXVXamq2enoGhYeTGBJCPElTE27ehBs3YNUq0NODQ4dATQ1ycnL+/PNmczPbyWmSo6OGnR0EBoKREXh4wIYNICjIvfDITIQpKSmxsbEtS4YOHWptbf15zcrKSgkJCeK1pKRkRUVFFyZC1Bdt3QqrV6ekpNg2NnKaahybmvb22BG/CPV+jo4wfjwcOQIjR4KtbUR0dEBFxXIASmDg/GPH/Jydp/j7w/z54O0NtrYQH8+9wMhMhIKCgtLS0i1LREVF26wpIyNT8+HRTlVVlaysbLcHh3hbaSls3Sq3YEFRi5aAMgDpfv1IDAohnkejgZcXuLiAtvaWuro7ANIAUF4+1dPT2tl5CgBoaMCNG9yenrRbEmFKSsqpU6dSU1MHDhx49uxZTnlOTs7SpUuTkpIGDx58+PBhY2NjXV3d9hxw7NixYWFhS5Ysefv2bX19vbKycneEjfqQHTtAT2/YwoW5ysq/VVdPb2oqAVhCo7lv3kx2ZAjxPmnpelFR4bo6zo2QZFOTTHV1NecJfcubndhYGDWqe1tKu6WPXGlpqby8vKGhYUZGRsvy+fPnGxkZFRQU/PDDD9OmTWviTPPxqYyMjBEjRkRGRnp7e7u4uADAhAkThIWF7e3tJ06cGBQUhF37UGdJS8PWrXweHpfv3UtauNBSQ2OpoaGnvPwknug/jFAPJywsTKW+B+D8uTFra8uoVMk2K584AUOGQLeOD+/GAfXnzp07ePDgkydPiM2MjAxjY+OysjKi/XPQoEHBwcGTJ09u/wFra2tFRESoVOqXKpw/f97Pz298ix64ixcvNjQ0/PphJSQksLMMKYjOMkWcdSG6SHuHDzc30ywsWF5erP/9jyjge/pUYOrUxgcP2AMGdG1I3IcD6gEH1H/QM98MoaGntm+Prq7eBECRkNg9aNDYsjLPbdua5s5tY5WKmBi+NWsE1NXZ+/cz1dS+L2fRaLRv3jtx7xlhZmampqYm5yng0KFDMzMzvysRfvP/UlhYWFJScsSIEZwSBQUFLmS4hoaG0tJSGRmZLz3jRF/R5f9BAgIC7T1mUJDA7Nl8U6YA0VHZ1JTt4UFzc2u+exd6eavDd1wE3sXHx9fc3IzXoWe+Gby9lxsYaB0/fpLNZru5uY0fbx8by/bw4L9wgT84uFlH55PK9vaQnNx85AjF0lJw+XK2r2+zkFB7T9Se0XHcS4SVlZUtx0VISkqWd/XzUBERkQEDBixfvrxrD/sV9+7dCwjYfv9+HIvVRKFQhg83Xr9+3cyZM7kWwFfU19cnJiZaWlq2Ks/MzKTRaIMGDSIlqs995Ra/wwds7zHNzMDamhoYCDt3/leyaRP8/Tc1KAjWrevaqLjsOy4CT8PrAD34Iowfbz9+vD1n08YGUlLg8GGwtKTOnQs7dkDLkXRUKqxeDS4usGYNZeRIanp6V0bCva+90tLStbW1nM2amhoZGRmunb07HDp0yNbW7v79IhZrG8CvbPbe1FT+WbNmeXh4dOaw69evX7t2befDKywsdHJy+rw8NDT03LlznT8+jwgMhJMnP/ZR4+ODs2dh/35ISSE1LIT6IgEB8PKCZ8+gqgp0deHzMb0DBsClS3DvXhefl3uJUFNTMzs7u76+nth8/vz54MGDuXb2LpeUlOTl5c1mz2exUgE2AswCWNvcnACwLSQk5Ndff/32Ib6gvr6e3p3r5AUGBvr4+HTf8XsZBQVIToaW38lUVODQIZgzB1crRIgUSkoQHg5nz8KePWBnB5mZrSt8Oj1UF+iWRFhfX5+dnV1WVsZgMLKzs4uLiwFgyJAh+vr6gYGBLBYrPDz8/fv3EydO7I6zc8fBgwepVAWAowAtG98pAFv4+Ebv2/dz50/BYDAWLVqkr6+vpaW1ePHiiooKoryqqmrlypV6enq6urpbtmwhCk+fPm1qaqqrq+vp6fn+/XvOQU6dOqWnp2dsbHz58mWi5OTJk5cuXQLEoaICrZ4iuLiAsTH4+pIUEEIIrKzg6VNwdARzc/Dyghafal2vWxJhcnKyvb19SEgInU63t7fftGkTUX727NmoqChJScl9+/ZdvnyZ9vn0qz1DdTVUVX3jX0xMQlOTA0AbT2ybm52fPn1SXMz45kG+2WN36tSpCQkJT548odFoxG0cm80mRp7ExcUlJCQQjwCvX7++ZcuW48ePP3jwIDc3193dnfjxd+/ePXz4MDExMSgoyM3NLT09HQBevXr15s2bLr1gvCgkBG7fhhs3yI4Dob7r85bSa9duDRo0qn9/k0GDRl+/frurTtQtnWXMzMxev379ebm2tnZCQkJ3nLELFRSAsfG31wSprq4G+NIEN/JsNltHp5qPr/9XjkClwqVLYGX1xQqCgoLjxo2Lj4+vrq7W1tYOCgoCgOfPn6empkZFRQkKCgIAMQVreHi4h4fH8OHDAWDv3r36+vonTpwAABaLtXv3bjExMXNz81mzZp0/f34XL6691zWePYP0dJgz579NSUk4exb+9z8wMQEFBVIjQ6hPI1pKY2LAzS0tL28Pk3kbQBqgYvFi53v3VA0MDDp/ip446Ta5VFSgpOTb1bS0+r96lfeFnW/4+QXKyqQ72WM5NzfX0tLS2tpaXV29rq6usrISAPLz81VUVAQ/nWWhsLBQTU2NeK2urt7U1FRaWgoAIiIinOnoBg4c2Gp+A/SJfv3AywvMzD5OfW9mBm5usGABRES0bjtFCHGXlRU4O1/bu3cNMSsbgEx5+erz56/t2tUFibB3D5Yi0YQJdlTqLYCyz/Y0UqnnLS2tOj9w5+LFi/b29mFhYVu2bBk3bhxRqKqqWlBQ0GqBYmVlZU5rZ3Z2toCAQP/+/QGATqeXlf0X4Zs3b3Bquq8ZMAC8vKBVf90tW6CmBg4fJikmhNBHDQ0MgJb3AIJ0etcs1Y6JsIPWrVtHo7H4+JwAilsUv6dQ5jc3/+vvv6Xzp5CWlk5NTS0rK8vOzvb39ycK9fT0DA0N3d3dS0pKysvLo6OjAWDBggWHDh1KSkoqKSnx8fGZN28ePz8/AFCp1PXr11dXV8fFxV28eHHevHmdj4qX+fpCWhpERn4s4eeH8+dh2zZISyMvLIQQAMCMGeOlpY8BEA+umNLSR2fMmPCNn2kfTIQdpKqqevXqFRGRdD4+DQAXAD8AVypVjUq9euzYUXNz8w4fub6+npihxtXV1cjIyNzc3NXVde3atVZWVgBAoVD+/PNPQUFBOzs7Kyurhw8fAsCkSZN27tzp4eFhZ2enqakZEhICAMLCwo6OjhYWFmPHjvXx8QkLC9PT0wMAHR2dnjOavmeh0WD/fvD0hJYzjg4aBLt3w5w50NA13z0RQh1jbm6+fr29nJxJ//7/k5Mz2bBhgpmZWZccuRvnGuW+yMjIoKCgiO+cnJVGo3V4rtGCgoKDBw/evh319u1bOTk5GxsLb29vIt90DJ1ONzExWbNmzeLFizt8kN6CmGu0qqqqaw/b2ZkVJ08GK6vWM8vMmgUqKrBvXydj45qeOb0kl+FcowQeezMwmcy3b98qKyt34bxx2FmmU1RUVPbt29dVH4/Xr1+fM2fO2LFj53D6LiLuCwqC0aNhzhxQUvpYePQoDB8Otrbg4EBeZAghEBAQ4PQN7CrYNNoFqqurX79+TXTU7AwbG5uSkpKIiAhhYeEuCQx1hIYGLF0KrcaZSElBeDgsWwYfpjVACPEMTISdcvv27dGmpjLS0pqamv3799fX0fnll186fDQxMTHO+hUVFRWc6ei4hsViVVRUsFgfl0GJiYlpf4K/evXqVyaHy83NTUxM7GyI3LF1K2zY0LrQwgJmzgQ3NzICQgh1I0yEHRcYGDh58mRGcvIBNvsywFEAhVevFi5cuHjx4g4/eWWz2Xv37lVSUtLU1FRWVlZTU9uxY0fXht2mxsbGRYsWSUtLm5iYSElJTZo0iSj38fFJaff00wsWLOCM1vhcVFTUnj17uiBWLhAQ+KRdlGPnTnj7Fk6d4npACKFuhM8IO+iff/7ZsGHDUjY7lMXiLHCyrLl5H4DP6dOWlpaurq4dOOymTZt+//33q1evmpqaAkBWVlZAQEDXRf1FwcHBz549y8vLk5SUbGxsjImJIcofPHhAjMRAAAACAnDmDFhagrk5aGuTHQ1CqGvgHWEHBQcHK1OpQQCtlvlaB2DBxxd04EAHjllZWXngwIGjR48SWRAAtLW1z58/T7w+ffq0oaGhoqLiiBEjbt26RRTu37//+PHjxOuMjIxp06YBAJvN9vPzU1NTU1VVHTVqVG5uLgAcO3ZMW1tbRUVl2LBh8fHxrU6dlZVlamoqKSkJADQajTN+f+nSpY8fPyaOGRwcPGzYMDU1tUWLFtXU1BAVzpw5o6urq6mpebjFqPOGhgZfX18tLa3Bgwdv3ry5qampA1ejR9i/H1qtmqmrC/7+ac7Oq2bOnG1jE7pvH/Ob0/EhhHo2/LLfljdvoLn561X+iYmZ0NQk2Nauqc3N6549a8jIEBJsc38LAwdCiwUzHz9+3NzcbG1t3WZdbW3tO3fuyMvLJyQkODk5paWlycvLFxQUENkLAOh0OjGtdnR0dERExLNnzyQkJDIyMsTFxQsLCzdu3JiSkjJgwICCgoKWTwEJEydOJLqq2tjYWFlZcSZme/78OZHzjh07dvHixYiICFlZ2bVr13p7e4eFhT158mTt2rV///23tra2p6cnZ71Jb2/vurq6pKQkFos1Y8aMkJAQb2/vb1yKnqmwEDZuhA9fNQiPDA29//13z4sXigCXHj2aeevWlS5fHg0hxEWYCD9TUAD29t9MhDUlJV9aVlgGgM1m1zg4CH1zVejTp6HFCvLv3r2TkpLirCU9ZcqUwsJCADh37pyOjo6Zmdk///wTGRnZ0NAgJSWVlJTk8IWu/DQaraqqKiYmZty4cbq6ugBQXFzMZDL//vvvadOmqaioAACbzWYwGABAoVAEBQWdnJwiIyPPnDnj5+f35s0bd3f3gwcPtjzmsWPM4dd6AAAT70lEQVTHPDw8GhoaCgoK5s+fP27cuLCwsN9++23evHlDhw4FgO3btxP3pkwm8/Tp07GxscTzwjlz5oSFhfXWROjvD7q6kJgIH+7RASDQx+dMY6MWAABsoNNnpaenp6d3ycy/CCFSYCL8jIoKvHr1zVoK2to5L1+2uSsbgCYgIPPvv/CdT9cUFBTKy8vr6+uJ4RMhISFNTU0GBgZEV8y5c+fm5+c7ODhISUnx8/N/ZRy6lZXVjz/+uHv37jlz5kyZMuXIkSMKCgoXLlw4dOjQqlWrxowZc+zYMSqVamxsDABUKpVYMNLS0pJY1CkhIcHW1tbR0ZFY2oJQUFBw4sQJzoLDJiYmDAajpKSEyIIAICcnR0wFXlJS0tTUxFkoEQDU1dW/6zr0IOLisGMHrFoFjx4B33/PEQoLCzVaVNFhMHJzczERItR7YSLsoPGTJp14/bqIxVL8tLwB4ByVamdn14E+JiYmJrKysmfPnl26dCkAqKqqAgCFQgGAmpqaS5cu1dbWCgoKstnsAx+eQYqKitbV1RGvCwoKOIdaunTp0qVLCwsLZ8+effz4cR8fn8mTJ0+ePLmmpsbLy8vf3//MmTNf6uE5ZswYNTW1/Pz8loUDBw5cvXr1jBkzWhaqqKjk5ORwzk7cYiooKNBotNDQUC0tre+9Aj2RqyucOgVnzsDChUSB/tCh99+8sWKzAYANECMouGjIEFJDRAh1CibCDlq7du3pEyem1Nf/wWJ9WLYHKgEWUCh5AL9+mCP7uwgJCYWEhLi5uVVVVdnb24uKiqakpDCZTH5+fmFhYSEhoYiIiNGjR4eGhr59+5b4EVNTUy8vL1dXVyaTuXfvXqIwPj7+3bt3RkZGjY2NDAZDTk7u1atXjx8/trS0pFKpjY2NSp+NDdi9e7e4uLiZmZmYmNjVq1fz8/MtLCxa/b5+fn5ycnIGBgb5+fnJyclubm6urq5jxoyZOnWqvr6+r68vkfv5+fm9vb0XL14cHBysrKyclZVVVlbm7OzcgQvSI1AoEBQEkyaBkxNISQFAQGjolORkl7IyZQbjD2lp+0WLiK8sCKFeChNhBykrK1+7edPFyUmrpsaOzR7IZpcCRFGpTfz8v5w6ZdrikdJ3mTFjhqqq6qFDhy5cuCAiIqKjoxMbG0s0P166dGn79u11dXWzZ8/evHkz8eHr6OiYkpLyww8/KCkp7dix448//gAAGo128uTJf//9V1xc3NnZ2dXVtaCgICIiYt++fVQq1cLCYtu2ba3OO3r06LCwsPDwcAaDoa2tfffuXWJibisrK3l5eQCYPXu2qKjo/v378/LyFBQUZs6cCQC6urrEYr8NDQ1r166l0WjEvI47d+48evSoj49PaWmphobG8uXLAUBNTa3Dl4VkhoYwZQr4+8PBgwCgrKwcl5Fx7erV8uLibTY2hoaGZMeHEOoUnHS7U5Nul5aWhoaG3rl9u7S4WKpfP0tbWw8PD1zboZ166KTbbaqoAH19ePwYBgzo4iN3Dx6bZ7ljcNJtAr4ZvgnvCDtFXl5+27Ztn99gIV4jIwOPH7c93QxCqJfDRIhQ+/SSe0GE0PfCmWUQQgj1aZgIEfoe//4LnVhgBCHUA2EiROh7iIuDry98YS4FhFBvhIkQoe/Rvz/4+YG7O9lxIIS6DHaWAQBITk7GxYa4jzNJdy/j6Qm//ALXrsHUqWSHghDqAvjpD5aWlqtWrSI7ii7AZrOJ+dh6kbFjx5Idwvfj54eDB2HhQrC3hz4/Rg0hHoCJEKKjo8kOoWvQ6XQhISE+Pmzu7n42NjByJOzfD5s3kx0KQqiz8EMToQ45eBB4Y1ZxhPo8TIQIdYiSEsycSXYQCKEugIkQIYRQn4aJECGEUJ+GiZB3PHnypKSkhOwoSFZXVxcXF0d2FOSLiYmpr68nOwqSFRUVJScnkx0F+SIjI3lplaHugImQd4SGhsbExJAdBcmePXsWEBBAdhTk27p16/Pnz8mOgmR//fXXkSNHyI6CfG5ubu/evSM7ih4NEyFPwe99CHHgnwNqJ0yECCGE+jRMhAghhPo0Ci+1Hpw5c2bdunXDhw8nOxBypKWlycnJKSgokB0ImWpqal6+fGliYkJ2ICRLTEzU1dUVFxcnOxAyFRUVVVZW6uvrkx0IyWJiYsaOHUulUskOhBxOTk4rV678eh2eSoTV1dXXrl1TVlYmOxByFBcXS0pKCgsLkx0ImZqamoqKigb0+dXk8/LylJWV++xnH4FOp9fW1vbv35/sQEiWk5Ojrq5OdhSkUVdX19DQ+HodnkqECCGE0PfCZ4QIIYT6NEyECCGE+jRMhAghhPo0TIQIIYT6NFyYlxdkZmYWFBRwNm1tbXvdUvUdxmKxXrx4UVJS0uq3zs/Pv3fvnqysrL29vYCAAIkRckdVVVVaWlr//v21tbWJEjqdnpCQwKmgpaWlqqpKUnRckpmZ+fTpU0FBwbFjx8rJyXHKi4uL7969Ky4uPmHCBEFBQRIj5IKGhob4+PjCwkJ1dXVzc3Pij4LJZMbGxnLqtKcjZZ+CvUZ5gbu7e2RkJKeHdFRUVB9JhCkpKWZmZiIiIuXl5Uwmk5//vy928fHxU6dOnTZt2osXL4SFhaOioji7eNLy5cvDwsKEhYUXLFhw8OBBovDly5cGBgaWlpbE5ooVK5ydncmLsdvt3Lnz6NGjZmZmdDo9Njb22rVrxO+emppqY2MzefLkvLy89+/fx8XF8fYQIzk5OV1dXQ0NjQcPHqipqd26dUtAQKC8vFxeXt7W1paoM3fu3AULFpAaZg/DRr3fqlWrAgICyI6CBLW1tUVFRa9evQIAJpPJKbe2tt6/fz+bzWYwGFpaWlevXiUvRm7Iy8tjMBgrVqzw8vLiFGZlZcnKypIYFZe9efOG8x748ccfra2tidcuLi6bNm1is9lNTU0mJianT58mLUSu+Pfff4kXtbW18vLyN27cYLPZZWVlfHx8pMbVo+EzQh7x9u3biIiIzMxMsgPhKjExsc9n0qmrq4uJiZk+fToA0Gg0R0fHmzdvkhEd9wwYMIBGo31ezmKx/vrrrwcPHtTV1XE/Ki4bOHAg575fUVGRwWAAAJvNvnnzJvFmoFKpTk5OPP9m4LR5iomJiYuLNzY2cnbFxsbGxsbW1taSFFrPhYmQF1Cp1LS0tCNHjpibmzs5OTGZTLIjIlNhYSEAcCYYUlZWfvv2LakRkUZKSiooKMjd3V1TU/PBgwdkh8MltbW1P//8s5ubGwBUVlY2NDSoqKgQu/rUm+HChQtMJtPe3p7YVFJS2rdvn6+vr5qa2p07d8iNrafh5QcnfceBAweIybSqq6tNTExOnTq1fPlysoMiDYvFAgDOU1IqldrU1ERqROTQ0NDIzs4mXgcEBCxZsuTFixfkhsQFDAZjxowZo0aNIp6B9dk3w4MHD7y9vf/8809ivlkZGZnc3Fw+Pj4ACA0NXbBgQWFhYR/pSdAeeEfICzhTSkpJSTk4ODx9+pTceMiloKDAZrPLy8uJzZKSEkVFRXJDIkXLiUbnzJmTmZnZ0NBAYjxcwGQyZ82aJSoqGhYWRnzKy8rK8vPzl5WVERVKSkqUlJRIjZEbEhMTnZ2dz507Z2ZmRpRQKBQiCwLAnDlziouLi4qKyAuwx8FEyFPYbPbTp095vpf810lJSQ0fPjwqKgoA2Gx2dHS0tbU12UGRLCkpSU5OTkhIiOxAuhGLxXJ1dWUwGOfPn+c8LOTj47O0tOS0BEZFRVlZWZEWIlekpKRMnTr15MmT48aNa7NCUlKSkJBQy+ElCJtGeYG1tbWZmZm4uHh0dHR+fv6yZcvIjohLGhoavLy83r17BwArVqwQFxc/cOAAAGzatGnlypXFxcWpqak1NTUzZ84kO9LudePGjZs3b8bHx/Px8S1btmzatGkODg779u3LysrS0tJ6+/bt6dOniSvDw7Zv337p0qV58+Z5enoCgISExN69ewFgw4YNLi4u9fX1ubm5z58/v3DhAtmRdiMmk2lra6uoqHjz5k2iW9D06dPHjRt38uTJ+Ph4PT29srKyU6dObdu2rS8Mrm0/TIS8wNfX98mTJ7W1tbNnzyaahsiOiEuoVKqxsTEAEPd8nDseFxcXBQWFyMhIIyOj4OBg3h43BgBKSkrGxsbEpQAAoift9OnTb9++nZ+fr6ysfP/+/WHDhpEaY7ezs7Nr2QYuIiJCvLC1tb179+61a9cGDRr0+PFjaWlpkgLkBj4+vl27drUsIVahcnBwaGpqysvLk5GRiYiIGDlyJEkB9lA4oB4hhFCfhs8IEUII9WmYCBFCCPVpmAgRQgj1aZgIEUII9WmYCBFCCPVpmAgRQgj1aZgIEeIpt2/fbrkEa5d79OhRYmIi8frp06cnTpxosxqdTr906RIx1wFCPRwmQoS6BpPJHDRo0NatW4nN3Nzc48ePV1VVdd8Zy8rKjh8/3mo5BX9//0OHDnXfGcePH8+ZuvP27dsrVqxos6aIiMj+/ft37tzZTZEg1IUwESLUNdhsdk5ODmey79TU1GXLlhFrQnWT3NzcZcuWtVpTYsWKFbNnz+6mM/7000+ampqTJk1qT+VNmzYFBQX1nWWPUO+FiRAhbmOxWMXFxV9ZH5XJZBYXF7dcV7KqqopOp7fn4AsXLiTWoW2ptra2uLi4ubn5Sz9VU1PzzeO/f//+zJkzxPJG7eHg4CAlJfWltlOEeg5MhAh1vQsXLsydOxcAxowZIy0tLS0tnZ6eDgAsFisgIEBBQUFRUVFSUnLMmDFEOWHw4ME//vjjpk2bpKWlFRUVo6Oj6+vrZ8yY0a9fP2lpaVFRUQ0NjV9++YWo/Pfff9va2gKAs7MzcYq7d+8CwMSJE5csWcI55pMnT0aNGiUhIaGoqKioqPjzzz9zdl2+fFlaWvr+/ftWVlZSUlLi4uLm5uYFBQVf+qWuXLny7t07Z2fnL1VITk7W0NCYPXs2sd4TPz//lClTwsPDO3gREeIWTIQIdT0rKytfX18AOHDgwO+///7777+rqakBgIeHx+7du319fVNTU+/duwcAtra2nNbU6urqEydOxMXFXbhwITY2Vl9fv7GxUUhI6OzZs8+fP3/48KG5ufmiRYuIHxw2bBgxvfLmzZuJUxgaGgJAeXk558FkTk6Ora1tbW1tREREUlLS1KlT16xZw8mFjY2NVVVVrq6uEydOTExMDAsLS0tLW7t27Zd+qXv37qmqqn5pPT9iuauxY8eGh4dzZj8fPXp0Tk5OTk5O5y8pQt2IjRDqCgwGAwBWrlxJbF67dg0A0tPTORWysrIoFMqBAwc4JaWlpaKionv37iU2ZWVlZWVl371796VTNDc3GxgYuLm5EZuPHz8GgKioqJZ1TExMpk+fTrz28vLi5+fPzs7m7LWyspKRkWlsbGSz2cSCRDt27ODs9fX1FRQUZLFYbZ7dyMjIzs6uZcn27dupVCqbzQ4LCxMQEFi/fn1zc3PLCgkJCQDw559/fuk3QqgnwGWYEOKS6OhoNpstIyNDtGESlJWVW7aOjhs3TlxcvOVP1dfXX7x4MSsrq7KyEgDq6upev37dzjOmpKSYmpqqq6tzSmbNmhUTE/P8+fPhw4cTJQ4ODpy9enp6DAajtLSUWMiplbKyMk1Nzc/Lt2zZsmvXrqNHjy5evLjVLhkZGQAoLS1tZ8AIkQITIUJcUlJSAgDe3t6tylv2miFWj+N4/fq1hYUFi8UaN26crKwsjUYTFBSsqalp5xnz8vKMjIxalhANmxUVFZySfv36cV7TaDQAaGxsbPNo/Pz8LBarVWFzc3NQUJCxsTHxTLQVor8PrgGLejhMhAhxiaSkJAA8ffp04MCBX6rDx/fJY/sjR440NjZmZWVxlpP9559/vtLdtBVRUVHOA0gCMQSQiOR79e/fv9XRiIDv3bs3fvz4CRMm3Lhxo9XtLJFxW2V3hHoa7CyDULcQExMDgPr6ek6JpaUlAFy6dKn9B8nJyVFXV+dkwZKSkidPnrQ6BdFFs02jRo1KTExs2TJ548YNMTExAwOD9sfAYWJikpaWxv5sKW8jI6PY2NiXL19OnDix1VQyz549o1AoI0aM6MDpEOIaTIQIdQttbW1+fv7Dhw/Hx8cnJSXV19ePGDFi+vTpW7ZsIYaZ0+n0tLS0nTt33r59+0sHGT58+NOnTy9fvsxgMNLT011cXCgUCmfvwIEDxcXFT548ef/+/aSkpM+bTD09PVks1qxZs7KysioqKnbt2nXt2jVPT09Or87vYmdnV1lZmZWV9fkuPT29e/fuvXnzxsbGpmW7a0JCwrBhw+Tl5TtwOoS4BhMhQt1CWVk5JCQkJibGxsZmxIgRL1++BIBz5865ubn5+fmpqKiIiooOHTr03LlzrZoTW1qzZo2NjY2Li4uQkJChoeGoUaMmTpzI2SssLHzq1KmsrCw7O7sRI0YQXTRbGjJkyJUrV16+fKmjoyMrK7t161YvL69t27Z17DeaOHGisrLyxYsX29yrra19//79qqoqCwuLoqIiAKirq7tx44abm1vHTocQ11A+b+hACHUMi8Xi4+NredPWJjqdnpWV1dzcrKKi0p7nZ7m5uaWlpYMHD5aSkvpmZWL6mJbPGlksVkZGBp1O19HRkZCQ+OYRvmL37t3Hjh17+fJle/q/nDp1ys/P7/Xr1508KULdDRMhQqi9GhoadHV1N27c2HLymjYxmUwdHZ01a9asWrWKO7Eh1GGYCBFC36G8vJzFYn3zRraxsbGgoEBVVZWfH7umo54OEyFCCKE+DTvLIIQQ6tMwESKEEOrTMBEihBDq0zARIoQQ6tP+D2P7/y580N/VAAAAAElFTkSuQmCC",
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
""
]
},
"execution_count": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ξ = [1.,2.,-1.,1.]\n",
"\n",
"errJ = zeros(length(X))\n",
"for (i,x) ∈ enumerate(X)\n",
" errJ[i] = norm( x - ξ, Inf ) \n",
"end\n",
"\n",
"errGS = zeros(length(Y))\n",
"for (i,y) ∈ enumerate(Y)\n",
" errGS[i] = norm( y - ξ, Inf )\n",
"end\n",
"\n",
"plot( errJ, label=\"Jacobi\", m=3,\n",
" xlabel=\"Iteration (k)\", yaxis=(:log10, \"Error\"), \n",
" legend=:bottomleft, color=:blue ) \n",
"\n",
"plot!( errGS, label=\"Gauss-Siedel\", m=3, color=:red ) \n",
"\n",
"L = -tril(A, -1)\n",
"U = -triu(A, 1)\n",
"D = Diagonal(diag(A))\n",
"\n",
"ρJ = maximum( abs.( eigen(inv(D)*(L+U)).values ) )\n",
"ρGS = maximum( abs.( eigen(inv(D-L)*U).values ) )\n",
"\n",
"plot!( 1:length(X), (errJ[1]/ρJ)ρJ.^(1:length(X)), \n",
" primary=false, linestyle=:dash, color=:blue )\n",
"plot!( 1:length(Y), (errGS[1]/ρGS)ρGS.^(1:length(Y)), \n",
" primary=false, linestyle=:dash, color=:red )"
]
},
{
"cell_type": "markdown",
"id": "6bdf9c2a",
"metadata": {},
"source": [
":::\n",
"\n",
"::: {#exr-}\n",
"\n",
"What is the rate of approximation in the error curves of @exm-GS?\n",
"\n",
"::: \n",
"\n",
"## Convergence\n",
"\n",
"Both Jacobi and Gauss-Siedel can be written as \n",
"\n",
"\\begin{align}\n",
" x^{(k+1)} = T x^{(k)} + c \n",
"\\end{align}\n",
"\n",
"::: {#exr-T}\n",
"\n",
"Write down what $T$ is for *(i)* Jacobi, *(ii)* Gauss-Siedel.\n",
"\n",
":::"
]
},
{
"cell_type": "markdown",
"id": "b94ee96d",
"metadata": {},
"source": [
":::{.callout-tip}\n",
"# TO-DO\n",
"\n",
"* Assignment 3: Due next Wednesday 11:59 PM,\n",
"* Make sure you are up to date with reading from chapter 1: @Burden sections 5.1 - 5.6, @FNC [zero stability](https://fncbook.com/zerostability/) \n",
"* This lecture: sections 7.1 - 7.3 of @Burden \n",
"* (Sign-up to office hours if you would like)\n",
"\n",
"Next \n",
"\n",
"* Read: sections 7.4 - 7.5 \n",
"\n",
":::"
]
}
],
"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
}