{ "cells": [ { "cell_type": "markdown", "id": "e7b74a1b", "metadata": {}, "source": [ "---\n", "subtitle: \"My solutions (of course, you may have got the correct answers in different/better ways)\"\n", "date: 2026-03-18\n", "abstract-title: Overview\n", "format:\n", " html:\n", " other-links:\n", " - text: This notebook\n", " href: Exam-midterm-answers.ipynb\n", "---" ] }, { "cell_type": "code", "execution_count": 1, "id": "bf72ca2e", "metadata": {}, "outputs": [], "source": [ "#| include: false\n", "\n", "using Plots, LaTeXStrings, LinearAlgebra, Polynomials\n", "\n", "function splitting(A,b,x; P=0., Niter=100, tol=1e-10) \n", " \n", " X = [x]\n", "\n", " if (P==0.)\n", " P = Diagonal( diag(A) )\n", " end\n", "\n", " if ( isapprox(det(P), 0; atol=1e-12) )\n", " @warn \"P is singular\"\n", " return X\n", " end\n", "\n", " iP = inv(P)\n", " r = b - A*x\n", " for i ∈ 1:Niter\n", " push!( X, X[end] + iP*r ) \n", " r = b - A*X[end]\n", " if (norm(r) < tol)\n", " return X\n", " end \n", " end\n", "\n", " @warn \"max iterations exceeded\"\n", " return X\n", "end;\n", "\n", "function SD(A,b,x0; α=0., Niter=50, tol=1e-3, xlims=[-25,25], ylims=[-25,25])\n", "\n", " g(x,y) = dot([x;y], A*[x;y])/2 - dot([x;y], b)\n", "\n", " x_range = range(xlims[1], xlims[2], length=100)\n", " y_range = range(ylims[1], ylims[2], length=100)\n", "\n", " P = contour(x_range, y_range, g, \n", " xlims=(xlims[1],xlims[2]),\n", " ylims=(ylims[1],ylims[2]),\n", " fill=true, \n", " levels=50, \n", " xlabel=\"x\", ylabel=\"y\", legend=false,\n", " colorbar=true)\n", "\n", " if (α==0.)\n", " adaptive_step = true\n", " title!( P, \"Gradient descent (with line search)\" )\n", " else\n", " adaptive_step = false\n", " title!( P, \"Gradient descent (with fixed α=$α)\" )\n", " end \n", " \n", " r = b - A*x0\n", " anim = @animate for i ∈ 0:Niter\n", "\n", " if (i==0)\n", " scatter!(P, [x0[1]], [x0[2]], markersize=5, color=:white)\n", " else \n", " \n", " if (norm(r) **Print name:**\n", "::: \n", "\n", "---\n", "\n", "Time Allowed: **1 hr 30 mins**\n", "\n", "Answer **SECTION A** (40%) and **TWO** of the three (30% + 30%) optional sections (B, C, D). If you have answered more than than the required two optional sections, you will be given credit for your best two optional sections. As a guide, I suggest you spend around **30 minutes** on each section.\n", "\n", "Please indicate your answers to multiple choice questions clearly with an ╳ to the left of your answers.\n", "\n", "**Please label any additional paper you use with your name and page number.**\n", "\n", "Electronic devices are not needed and not permitted in this examination. You may use your own notes, limited to **2 sides of letter** sized paper. Additional paper may be requested. \n", "\n", "---" ] }, { "cell_type": "markdown", "id": "3df6a28d", "metadata": {}, "source": [ "## A: Compulsory Section **[40 points]** {.unnumbered}" ] }, { "cell_type": "markdown", "id": "c27edcf9", "metadata": {}, "source": [ "::: {.callout-note}\n", "# Useful definitions\n", "\n", "Consider the initial value problem: seek $u:[0,T] \\to \\mathbb R$ such that\n", "\n", "\\begin{align}\n", "\\tag{IVP}\n", "\\begin{split}\n", " u(0) &= u_0 \\\\\n", " u'(t) &= f\\big(t,u(t)\\big). \n", "\\end{split}\n", "\\end{align}\n", "\n", "Recall, $u_0$ is known as the *initial condition* and $(u_0, f)$ is known as the *data*.\n", "\n", "For a fixed mesh $t_j = j h = \\frac{jT}{n}$ (for $0 \\leq j \\leq n$), we use the notation $u_j \\approx u(t_j)$ for approximations coming from the numerical scheme in question.\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "eb07bb94", "metadata": {}, "source": [ "**1.**
[10 points]
\n", "\n", "Please select all the **true** statements. \n", "\n", "
\n", "\n", "\n", "There is $f$ for which there exists a solution $u : [0,T) \\to \\mathbb R$ which blows up at time $T$. \n", "\n", "\n", "::: {.box}\n", "\n", "True! We saw that the solution to $u'(t) = u(t)^3$ with $u(0) = 1$ has finite-time blow-up at $t=\\frac12$.\n", "\n", ":::\n", "\n", "\n", "If $f$ is continuous, then there exists a unique solution to $\\text{(IVP)}$.\n", "\n", "\n", "::: {.box}\n", "\n", "False! The theorem we had in class needed $f$ to also be Lipschitz in the second component with constant independent of $t$. Indeed, $f(t,u) := \\sqrt{|u|}$ is continuous on $[0,T] \\times \\mathbb R$ but there is not a unique solution to $u'(t) = \\sqrt{|u(t)|}$ with $u(0) = 0$.\n", "\n", ":::\n", "\n", "\n", "When $f(t,u) = \\cos(u)$, the problem $\\text{(IVP)}$ has a unique solution for all time.\n", "\n", "\n", "::: {.box}\n", "\n", "True. The solution is $u(t) = u_0 + \\sin t$\n", "\n", ":::\n", "\n", "\n", "$\\text{(IVP)}$ with $f(t,u) = \\sqrt{u}$ is well-posed.\n", "\n", "\n", "::: {.box}\n", "\n", "False. There are mulitple solutions when $u(0) = 0$.\n", "\n", ":::\n", "\n", "\n", "Euler's method is consistent\n", "\n", "\n", "::: {.box}\n", "\n", "True. It has order of accuracy $1$ since\n", "\n", "\\begin{align}\n", " \\tau_{j+1}(h) \n", " %\n", " &= \\frac{u(t_j+h)-u(t_j)}{h} - u'(t_j) \\nonumber\\\\\n", " %\n", " &= \\frac{u(t_j) + hu'(t_j) + \\frac{h^2}{2} u''(t_j) + O(h^3) -u(t_j)}{h} - u'(t_j) \\nonumber\\\\\n", " %\n", " &= O(h)\n", "\\end{align}\n", "\n", "as $h\\to0$.\n", "\n", ":::\n", "\n", "\n", "The order of accuracy of Runge-Kutta methods can be seen numerically by plotting the error against the number of grid points $n$ on a log-log scale \n", "\n", "\n", "::: {.box}\n", "\n", "True! Order is $O(h^p) = O(n^{-p})$ and so if we plot $y = C n^{-p}$ on a log-log, we should see $\\log y = -p \\log n + \\log C$ and so the order is the negative gradient of the error curve.\n", "\n", ":::\n", "\n", "\n", "If one uses an adaptive time stepping Runge-Kutta method, it is common to notice that the step size is larger in regions where the derivative of the exact solution is large \n", "\n", "\n", "::: {.box}\n", "\n", "False. We expect to see the opposite - more mesh points where the derivatives are large.\n", "\n", ":::\n", "\n", "\n", "If a multistep method is consistent then it is automatically convergent\n", "\n", "\n", "::: {.box}\n", "\n", "False. It also has to be zero-stable.\n", "\n", ":::\n", "\n", "\n", "A linear $m$-step method can have order of accuracy $p = m+1$ \n", "\n", "\n", "::: {.box}\n", "\n", "True. Implicit methods aim to achieve this!\n", "\n", ":::\n", "\n", "\n", "If a method is zero stable, then it produces a sequence of approximations $(u_j)$ which remains bounded as $j\\to \\infty$ \n", "\n", "\n", "::: {.box}\n", "\n", "False. Zero stability is about the limit as $h \\to 0$.\n", "\n", ":::\n", "\n", "\n", "$u$ and $v$ are $A$-orthogonal if $(Au, Av) = 0$\n", "\n", "\n", "::: {.box}\n", "\n", "False. We defined $A$-orthogonality to be $(u,v)_A := (u, Av) = 0$.\n", "\n", ":::\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "5be7c49d", "metadata": {}, "source": [ "---\n", "\n", "**2.**
[5 points]
\n", "\n", "Suppose that $f$ is continuous on $[0,T]\\times \\mathbb R$ and $f$ is Lipschitz in its second component with constant $L$ (independent of $t$) and that $u$ solves $\\text{(IVP)}$. Suppose that $v$ solves the perturbed problem\n", "\n", "\\begin{align}\n", "\\tag{IVP$\\delta$}\n", "\\begin{split}\n", " v(0) &= u_0 + \\delta \\\\\n", " v'(t) &= f\\big(t,v(t)\\big). \n", "\\end{split}\n", "\\end{align}\n", "\n", "Explain why one can not improve the general bound $\\max_{t\\in[0,T]} | u(t) - v(t) | \\leq |\\delta| e^{LT}$ for all $|\\delta|$ sufficiently small. \n", "\n", "Hint: It is useful to consider the function $f(t,v) = v$." ] }, { "cell_type": "markdown", "id": "21c0dc64", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**Answer.** \n", "\n", "If $f(t, u) = u$ then the solution to $\\text{(IVP)}$ is $u(t) = u_0 e^{t}$ and the solution to $(\\text{IVP}_\\delta)$ is $v(t) = (u_0 + \\delta)e^t$ and so \n", "\n", "\\begin{align}\n", " |u(t) - v(t)| = |\\delta| e^{t}.\n", "\\end{align}\n", "\n", "Since $f$ is Lipschitz with constant $L=1$ (since $|f(t,u)-f(t,v)|=|u-v|$), we have $\\|u-v\\|_{L^\\infty([0,T])} = |\\delta| e^{LT}$. In particular, this bound is sharp (it cannot be improved).\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "83221034", "metadata": {}, "source": [ "---\n", "\n", "**3.**
[10 points]
\n", "\n", "Write down the Runge-Kutta method corrseponding to the following Butcher tableau:\n", "\n", "\\begin{array}\n", "{c|cccc}\n", "0\\\\\n", "c & a \\\\\n", "\\hline\n", "& b_1 & b_2 .\n", "\\end{array}\n", "\n", "Show that you can choose the constants $a,b_1,b_2,c$ to obtain a order $2$ method." ] }, { "cell_type": "markdown", "id": "9abc0a80", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**Answer.** \n", "\n", "The Butcher tableau corresponds to the RK method given by\n", "\n", "\\begin{align}\n", " u_{j+1} = u_j + b_1 h f(t_j,u_j) + b_2 h f\\big(t_j + ch, u_j + a h f(t_j,u_j) \\big).\n", "\\end{align}\n", "\n", "The local truncation error is \n", "\n", "\\begin{align}\n", " \\tau_{j+1}(h) &= \\frac{u(t_{j+1}) - u(t_j)}{h} - b_1 f(t_j,u(t_j)) \n", " %\n", " \\nonumber\\\\ &\\qquad - b_2 f\\big(t_j + ch, u(t_j) + a h f(t_j,u(t_j)) \\big) \\nonumber\\\\\n", " %\n", " &= u'(t_j) + \\frac{h}2 u''(t_j) + O(h^2) \\nonumber\\\\\n", " %\n", " &\\qquad - b_1 u'(t_j) - b_2 \\Big[ f(t_j,u(t_j)) + ch \\partial_t f(t_j, u(t_j)) \\nonumber\\\\\n", " %\n", " &\\qquad\\qquad + a h f(t_j,u(t_j)) \\partial_u f(t_j,u(t_j)) + O(h^2) \\Big] \\nonumber\\\\\n", " %\n", " &= (1-b_1-b_2) u'(t_j) + \\frac{h}2 u''(t_j) \\nonumber\\\\\n", " %\n", " &\\qquad - h b_2 \\Big[ c \\partial_t f(t_j, u(t_j)) \\nonumber\\\\\n", " %\n", " &\\qquad\\qquad + a f(t_j,u(t_j)) \\partial_u f(t_j,u(t_j)) \\Big] + O(h^2)\n", "\\end{align}\n", "\n", "Therefore, if $b_1 + b_2 = 1$ and $a = c = \\frac{1}{2b_2}$, then $b_2 \\big[ c \\partial_t f(t_j, u(t_j)) + a f(t_j,u(t_j)) \\partial_u f(t_j,u(t_j)) \\big] = \\frac12 \\big[ \\partial_t f(t_j, u(t_j)) + u'(t_j) \\partial_u f(t_j,u(t_j)) \\big] = \\frac12 u''(t_j)$ and thus $\\tau_{j+1}(h) = O(h^2)$ and $h\\to0$.\n", "\n", "Therefore, the RK method has (at least) order $2$ convergence for this choice of constants.\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "39a376b1", "metadata": {}, "source": [ "---\n", "\n", "**4.**
[15 points]
" ] }, { "cell_type": "markdown", "id": "41fd8113", "metadata": {}, "source": [ "::: {.callout-note}\n", "# Hint\n", "\n", "It will be useful in the following to use the Taylor expansion: \n", "\n", "\\begin{align}\n", " f\\big(t\\pm h, u(t\\pm h)\\big) &= f\\big(t, u(t)\\big) \\pm h \\left[ \\frac{\\partial f}{\\partial t}\\big(t,u(t)\\big) + \\frac{\\partial f}{\\partial u}\\big(t,u(t)\\big) u'(t) \\right] + O(h^2)\n", " \\nonumber \n", "\\end{align}\n", "\n", "as $h \\to 0$. \n", "\n", ":::\n", "\n", "Here is a multistep method for some fixed $\\alpha \\in \\mathbb R$:\n", "\n", "\\begin{align}\n", " u_{j+1} &= (1-\\alpha) u_j + \\alpha u_{j-1} \n", " %\n", " + \\frac{h}{2} \\big( (\\alpha + 3) f( t_j, u_j ) + (\\alpha-1) f( t_{j-1}, u_{j-1} ) \\big)\n", " %\n", " \\tag{$\\alpha2$}\n", "\\end{align}\n", "\n", "1. Is this an implicit or explicit method? \n", " \n", "Recall that the local truncation error of $(\\alpha2)$ is given by \n", "\n", "\\begin{align}\n", " \\tau_{j+1}(h) := \\frac\n", " {u(t_{j+1}) - (1-\\alpha) u(t_j) - \\alpha u(t_{j-1})}\n", " {h}\n", " %\n", " - \\frac{1}2 \\Big( \n", " (\\alpha+3) f\\big(t_j,u(t_j)\\big) + (\\alpha-1) f\\big(t_{j-1}, u(t_{j-1})\\big) \n", " \\Big).\n", " \\nonumber\n", "\\end{align}\n", "\n", "2. Show that $(\\alpha2)$ has order of accuracy at least $p = 2$, \n", "3. For what $\\alpha$ is $(\\alpha2)$ zero stable?" ] }, { "cell_type": "markdown", "id": "45a9fa67", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**Answer.** \n", "\n", "**1.** This is an explicit method (there is no $u_{j+1}$ on the right hand side), \n", "\n", "**2.** Using the Taylor expansion of $u$ about $t_j$, we have \n", "\n", "\\begin{align}\n", " &\\frac{1}{h} \\big[ u(t_{j+1}) -(1-\\alpha) u(t_j) - \\alpha u(t_{j-1}) \\big]\n", " \\nonumber\\\\\n", " %\n", " &= \\frac{1}{h} \\Big[ u(t_{j}) + h u'(t_j) + \\tfrac{h^2}2 u''(t_j) -(1-\\alpha) u(t_j) \\nonumber\\\\\n", " %\n", " &\\qquad - \\alpha \\big[u(t_{j}) - hu'(t_{j}) + \\tfrac{h^2}2u''(u_j) \\big] \\Big] + O(h^2) \\nonumber\\\\\n", " %\n", " &= (1+\\alpha) u'(t_j) + \\frac{1-\\alpha}{2} h u''(t_j) + O(h^2) \\nonumber\n", "\\end{align}\n", "\n", "Moreover, \n", "\n", "\\begin{align}\n", " &\\frac{1}2 \\Big( \n", " (\\alpha+3) f\\big(t_j,u(t_j)\\big) + (\\alpha-1) f\\big(t_{j-1}, u(t_{j-1})\\big) \n", " \\Big) \\nonumber\\\\\n", " %\n", " &= \\frac{1}2 \\Big( \n", " (\\alpha+3) u'(t_j) \\nonumber\\\\ &\\qquad + (\\alpha-1) \\big[ \n", " f\\big(t_{j}, u(t_{j})\\big) \n", " %\n", " - h \\partial_t f\\big(t_{j}, u(t_{j})\\big) \\nonumber\\\\\n", " %\n", " &\\qquad\\qquad - h \\partial_u f\\big(t_{j}, u(t_{j})\\big) u'(t_j)\n", " %\n", " + O(h^2) \n", " \\big] \n", " \\Big) \\nonumber\\\\\n", " %\n", " &= (\\alpha+1) u'(t_j) - \\frac{\\alpha-1}{2}h u''(t_j) \n", " + O(h^2).\n", "\\end{align}\n", "\n", "Taking the difference of these terms, we get $\\tau_{j+1}(h) = O(h^2)$. \n", "\n", "**Extra:** for what $\\alpha$ is this method 3rd order accurate?\n", "\n", "**3.** This method is zero stable if it satisfies the root condition: the roots $\\zeta$ of the characteristic polynomial $\\rho(z) = z^2 + (\\alpha-1)z - \\alpha$ satisfy $|\\zeta|\\leq1$ and if $|\\zeta|=1$ then $\\zeta$ is simple. We have $\\rho(z) = (z-1)(z+\\alpha)$ so the roots are $\\zeta = 1$ and $\\zeta= -\\alpha$. Therefore, this method is zero stable if $|\\alpha|\\leq 1$ and $\\alpha \\not= 1$. \n", "\n", "As a result, the method $(\\alpha2)$ is convergent for $|\\alpha|\\leq 1$ and $\\alpha \\not= 1$ (we saw this in Example 7.4 in the lecture notes)\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "61b5657d", "metadata": {}, "source": [ "Choose ***TWO*** sections from B, C, D (each worth ***30 points***)\n", "\n", "## B: Simple Mixing **[30 points]** {.unnumbered} \n", "\n", "::: {.callout-note}\n", "# Useful facts and definitions\n", "\n", "Recall that splitting methods for solving $Ax=b$ take the form \n", "\n", "\\begin{align}\n", " P x^{(k+1)} = N x^{(k)} + b\n", " \\tag{splitting}\n", "\\end{align}\n", "\n", "where $A = P - N$. This can also be written as $x^{(k+1)} = P^{-1} N x^{(k)} + P^{-1} b$.\n", "\n", "In this question, we consider the following method: for fixed $\\alpha > 0$, define\n", "\n", "\\begin{align}\n", " x^{(k+1)} = (I-\\alpha A)x^{(k)} + \\alpha b.\n", " \\tag{$\\star$}\n", "\\end{align}\n", "\n", "It is useful to remember that, when $A$ is symmetric, $A$ is SPD if and only if all eigenvalues of $A$ are strictly positive, \n", "\n", ":::\n", "\n", "1. Suppose that $x^{(k)}$ defined by $(\\text{splitting})$ converges to some $x$ as $k\\to\\infty$. Show that $x$ is a solution to $Ax=b$, \n", "2. For $(\\text{splitting})$ to be practical, what is the key property that the preconditioner $P$ must satisfy?\n", "3. Which condition leads to convergence $x^{(k)} \\to x$? \n", "4. For which choice of preconditioner $P$ gives $(\\star)$? \n", "5. Show that, if $A$ is symmetric and positive definite (SPD), then $(\\star)$ converges for all $0 < \\alpha < \\frac{2}{\\rho(A)}$. \n", "\n", "For fixed $\\beta\\in \\mathbb R$, consider the matrix \n", "\n", "\\begin{align}\n", "A_{\\beta} := \\begin{pmatrix}\n", "1 & \\beta & \\beta \\\\\n", "\\beta & 1 & \\beta \\\\\n", "\\beta & \\beta & 1\n", "\\end{pmatrix}\\nonumber.\n", "\\end{align} \n", "\n", "which has $\\sigma(A_{\\beta}) = \\{1-\\beta, 1 + 2\\beta\\}$.\n", "\n", "\n", "6. Hence, show that for $-\\frac12 < \\beta < 1$, there exists $\\alpha$ sufficiently small such that $(\\star)$ converges when applied to $A_\\beta$." ] }, { "cell_type": "markdown", "id": "697c0191", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**Answer.** \n", "\n", "**1.** If $x^{(k)} \\to x$ as $k \\to \\infty$, then we have $Px = Nx + b$ and so $(P-N)x = Ax = b$. \n", "\n", "**2.** $P$ must be easy to invert,\n", "\n", "**3.** We require $\\rho( P^{-1} N ) < 1$, \n", "\n", "**4.** $P = \\alpha^{-1}I$ gives $N = \\alpha^{-1}I - A$ and thus $P^{-1} N = \\alpha(\\alpha^{-1}I - A) = I - \\alpha A$ as required.\n", "\n", "**5.** The eigenvalues of $I - \\alpha A$ are $1 - \\alpha \\lambda$ for $\\lambda \\in \\sigma(A)$. If $A$ is SPD then all the eigenvalues are strictly positive and so \n", "\n", "\\begin{align}\n", " \\rho(I-\\alpha A) < 1 \\qquad \\iff \\qquad &-1 < 1 - \\alpha \\lambda < 1 &&\\text{for all } \\lambda \\in \\sigma(A) \\nonumber\\\\\n", " %\n", " \\iff\\qquad \n", " &0 < \\alpha \\lambda < 2 &&\\text{for all } \\lambda \\in \\sigma(A) \\nonumber\\\\\n", " %\n", " \\iff\\qquad\n", " &0 < \\alpha < \\frac{2}{\\lambda} &&\\text{for all } \\lambda \\in \\sigma(A) \\nonumber\\\\\n", " %\n", " \\iff\\qquad \n", " &0 < \\alpha < \\frac{2}{\\rho(A)}.\n", "\\end{align}\n", "\n", "By question 3, under this condition, $(\\star)$ is convergent.\n", "\n", "**6.** $A_\\beta$ is symmetric. When $1-\\beta > 0$ and $1+2\\beta > 0$ we have that $A_\\beta$ is SPD and thus $(\\star)$ is convergent when applied to $A_\\beta$ for all $\\alpha$ sufficiently small. We have $1-\\beta > 0$ and $1+2\\beta > 0$ if and only if $-\\frac12<\\beta < 1$.\n", "\n", "**Extra:** Show that $(\\star)$ does not converge for any $\\alpha$ if $\\beta$ is outside this range.\n", "\n", ":::" ] }, { "cell_type": "code", "execution_count": 2, "id": "bcd26d53", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd1xT1/sH8JMJIQkz7CFDNogyVVQUERW3IIhaqVpFtM5at3XvVq0iWHDP1tGquHBPEAQZylL2kBU2Cdn390f8Ub6oFSEQQp73yz/C4eTkySXy4dx77r04DMMQAAAAIK/w0i4AAAAAkCYIQgAAAHINghAAAIBcgyAEAAAg1yAIAQAAyDUIQgAAAHKN2JZOQqEwOjr6wYMHTCbT1tY2ODhYRUXl024JCQnHjh0TCAQzZ84cPHiwpEsFAAAAJK9NM8KysrK1a9dqa2uPGDHi0aNHQ4YM4fF4rfqkp6d7enpaWFg4OzuPGzcuNja2E6oFAAAAJAzXlhPqRSIRHv8xMjkcjrq6+pMnT1xcXFr2mT9/PpFIDA0NRQht3rz57du3ly5d6oyKAQAAAAlq04ywOQURQvX19TweT1NTs1Wfly9fenh4iB8PHToUZoQAAABkQpuOETbDMCw4OHjGjBnGxsatvlVWVsZgMMSPGQxGeXl5y3lkS/v27bty5Yqurq74SxwOt2fPHgMDg2+uXTYJhUICgSDtKmQPbLd2EIlEOBwOh8NJuxAZAx+2dsAwDMOwz/7Ol5SSwgJKQyUh5Qml/6jsWpaNo2tbnkUkEr/6X+AbghDDsEWLFlVWVp49e/bT7yooKDQfOORyuQoKCl/aIlVVVVpaWv7+/uIvyWSyjo6O/Hzs2Gw2nU6XdhWyB7ZbO/B4PDKZLD//uSQFPmztIBAIBAIBhULpvJe4c+tWwV9hC52NE+9GHS4nX737oC3Passfgt8QhD///HN8fPy9e/eoVOqn3zU0NCwsLBQ/Liws/I8ZHo1Gs7a2bg5CeYPH4zv1j6aeCrZbO+D/n7QLkTGw0dqhsz9smIA/nlixiEsYG5VOU1E9/udlCb5WWwdat27d/fv379y50/LEiaqqqosXL4ofT548+fz58yKRCCF09uxZX19fSZUIAABAnmECftXxrQpKSlfiUpNzCp+/TrWwsJDg+G0KwqysrB07dpSVlbm4uJiZmZmZmUVHR4vbAwICxH3mzp3LYrFcXV0HDRqUkZGxdOlSCVYJAABAPmF8XlXkRrwCReO7VQjfKfv527Rr1MTEJCcnp2WLtrY2QsjJySk3N1fcQqfTY2Nj4+PjhUKhq6srmUyWeK0AAADkCsbnMSM3EmjK6jNWdlIKojYGIZlMNjU1/bRdQUHBxMSk+UsCgTBgwACJlQYAAECOdU0Kom89fQIAAADoAhiPyzy6iUBTVp+xCnXy2iVYGQUAAKB7wXhc5tGNBJpKF6QggiAEAADQTcQ8fz7U1XGAg+3ZHwOJalrq33VFCiLYNQoAAKA7EIlEwUEzTgwzppKIM2+98N4ShrrqokgyE4QYhqWkpAiFQmkX0lFsNltJSUnaVcgeKW43Y2NjDQ0Nqbw0APKjtrZWk0LUoSkihOz1tfILCrR1dLrmpWUmCLOysgYOHGhjYyPtQjoKwzC49mM7SGu7VVZWTpw48ffff+/6lwZArqjRqHwO5+SbYhqJmFLD7dOnT5e9tMwEoVAoNDU1TUhIkHYhQL4cPnw4IyND2lUA0MNhPC4z8pfz6378p4bE4XBu/zC3Uy9b2orMBCEAAIAeCeNxmJEbiRo6hgFLF0tjxw+sGgUAACA1GI/DjPiFqKGjFrC0y1bHtAIzQgAAANLxMQUZemoBS6SVggiCEAAAgFR8TEFNfTX/xVJMQQS7RgEAAHQ9jMdhRmzoDimIYEYIAACgy5SWliYlJfW1syFdP0zUNOgOKYhgRth5MAyrqanh8XjSLgQAALqFmBcvxgzuf3/3qtH9nZJrBd0kBREEYSfx8/NTU1NTV1f/+++/pV0LAAB0C38c+G2nm8Hyfvq7h1meSynoJimIIAg7ybx581JTU/v27SvtQgAAoLug0+nlLA5CqIzFU+lOly2EIOyQ2bNnHzx4UPwYwzAXF5dHjx4hhLy9vY2MjKRaGgAAdCMiDmteL3J4WqXv7axjxYI1m7ZKu6J/yXYQvijHsusxhNCbaiypCkMIFTRiT8swDKFqLrpfgjUJEF+EHnzAKpok07+VadOmhYWFYRiGEHr27FlFRcWQIUO66t0DAIBsEHFYzPB1hjZ9Yt8VPHuT9SwhWaerLqjdFjIchAIR2p0iupyHIYSOZolC00QIoRuF2K4UIYuPEpnY7lRhTgNWzMJ2pwifl4sk0r+V4cOHI4SePXuGEIqIiAgODiYQCF24DQAAoLsTpyDZ2FrN70eEw5HJZGlX1JoMnz5BxKPr3h9T5/cBHx8stMEvtMEjhEbo40bof3x3d0d/fNDx/q3gcLjg4OCIiAhbW9tr167t2bNHYm8PAABkn4jDYoatJZvYqE4KlnYtXyTDQdhNzJo1a+vWrUZGRt7e3np6etIuBwAAugtRUyMzfB3Z1FZ14jxp1/JfIAg7SlVVdeLEibt3775161Zz461bt4qLi6uqqh48eFBfXz927FjISACAXJGVFEQyfYyw+xg9erSxsfGIESOaW8rKynJzc6dNm6ahoZGbm9vU1CTF8gAAoIuJmhqZ4WsVZCEFEcwIO04kEkVGRi5cuBCP//evitmzZ0uxJAAAkIr79+7u2rieoqi43NHAcdBQFVlIQQQzwg56/vy5trY2n88PCQmRdi0AACBNTCZz1YJ5u22pyw2wpRcfKk+YK+2K2kq+ZoQsFis1NTUrK0tdXd3e3t7ExKSDAw4aNKiyslIitQEAgEzLzc3tq0nTpipoUxU0lOk1NTXq6urSLqpN5CUIRSLRr7/+umPHjrq6uuZGT0/PP/74o3fv3h0ZWSgUFhQU9OrVC84gBADIM2tT41f5pfc0Fep5Qr4iTVZSEMnPrtEff/xx1apVdXV1vdWpY8113A3VKUTCw4cPBw4cmJ2d3ZGRa2trzczMYF4IAJBnInaj6M+957asfN97SIPruKvR96Vd0TeQixnhw4cPw8PDSQT87yPtx5l/vK4Pk82bfys5rqRywYIFd+/elW6FAAAgu0TsxpqI9SQzO/vJ8+2lXUw7yMWM8Pjx4wihuf16NacgQoihRA73cVAg4u/fv19YWCi96gAAQIaJ2I3M8DXk3n2URgdJu5Z2ku0ZIScjAeN+/RS95PhYhJCXiVardk0lhb7aKnElNa/+Oa/p4f7VccgmNgSVbnTrEAAAkC5RY11l2GpFa2fq6CA+ny/tctpJloNQJGQnPsQEX9/07Lo6hJAS6TOLWagkIkKoNiuZrSL8+isSSRQIQgAAQAiJU/DwKkUbV5VxswUCgbTLaT9ZDkI8QX3GyrZ0ND37KO/Bg/TKeltNest2IYalMxsQQvbfL9dwde2UIgEAoCdqmYLSrqWj5OIY4fjx4xFCB1/l1nD+Z/oY8Tq/rJGjr6/v7OwspdIAAECWlJWVxcTEsCrLe0wKIjkJwnnz5llaWubXskedj4lMKogtrr6dXR5yK2XH83cIob1797a8OhoAAIDPuvbPP+MG9z+2Yu4AB5s6Q7uekYJItneNtpmiouKdO3cmTJiQmpq65WlmczuJRNqzZ09gYGBHBq+vr8fhcKqqqh0uEwAAurXftm0+6WWmokC6qkL8O7d6tbTrkRS5CEKEkLGx8atXr/7888+oqKj8/Hwqlerm5jZ37twOXlamvLw8LCysT58+ioqKkioVAAC6JzKJ2MgTqiiQGgQYhaIk7XIkRl6CECFEJpNnzpw5c+ZMCY45Y8YMPB5/6tQpCY4JAADdkKix9icnox+ux2ir0vkUlRvzuu8d57+VHAUhQqimpiYxMTEjI0NTU9POzs7Ozq6DA967d08ihQEAQHcmbKhhHl7tPt4v6feLNTU1mpqa0q5IkuQlCAUCwdatW3/99Vc2m93c6Obmdvz4cRsbmw4OzuVyExMTa2pqDAwMrK2tyWRy87eSk5NLSkrU1dX79u1LoVA6+EIAAND1xClI6eehPHIaQqiHpSCSnyCcO3fuyZMnEUI0ExWqoTK/gVefVR0XFzdo0KAXL15YW1u3e+Rr164FBwfr6uoaGxsXFBQwmczIyMiRI0eyWKzRo0dXVFQ4ODiUlZW9fv26rq4OlqcCAGSLsKGGeXgVxXGYsneH1hV2Z3IRhNHR0SdPnsSTCDaLndT7aYsbBWx+xsHEmreVCxYsePToUftGTk5ODggIOHfunK+vr7glNTX1w4cPCKGTJ09yOJy0tDTx7ZkqKipwOJwk3g0AAHQReUhBJCfnEYrnggZjTJtTECFEVCJZLXTEkwmPHz/Oy8tr38iHDh0aN25ccwoihPr06TNq1CiEUHV1NZVKbZ4CamlpQRACAGTIxxR06uEpiGR9RhhTEs8Tfv1ao3FJ8QghdQftVu0kOpluplqXUXXh4cWBIwZ9dRxL9d66tP8ZRDwjFD/m8XgsFgshRCKRaDTatGnTQkND7ezsRowYMWzYsFGjRikoKLTxfQEAgHQJ62uYYasoTsOUR/TwFEQyHYQiTPSkMIYj4H61Zx2rHiGEJ31m+ksgExBCSSVvuAVfGQSPwysRKa2CkMfjKSl9PJnmypUrCxcu5HA4I0aMuHbtmpmZWVZW1j///PP06dPg4GBNTc1nz57BefcAgO7s2tV/jh363UBfL8SQYDx0DH3EVGlX1BVkOAjxOPyaAUvb0jPW+sG94nuN+XU0Y5WW7ZgQayyoRwitHrPMycmpHTWYmJi8f/9e/DgwMDAwMHD37t0xMTHiFlVV1VmzZs2aNauxsdHe3v7MmTOLFi1qx6sAAEAXePPmzW9rlu8ZaJRemfnzW+6d7WekXVEXkYtjhBMnTkQIFV57z6v7n+lj4fX3vFqOkZFRv3792jfy1KlTT58+XV5e/um3hMJ/7+tEo9F0dXV5PF77XgUAALpAUlKStx7NgE7xNtWqaGB//Qk9hQzPCNvuhx9+CA8Pf/v27eu1T/S8TWi9VAQsXsXLD9VJ5Tgcbt++fe0+qyEwMPDBgwd9+/YNCQkxNzevrq7+888/Bw4ciBDasmVLenq6p6enqqrqo0ePMjMzJ0+eLNG3BQAAkuRiY3ko64OjFj29mmVibiHtcrqOXAQhmUy+c+fOpEmTXr16lX/p34tuKyoqHjhwoOWaz2+Fw+GOHTv2+PHjGzdu3L9/38jI6MCBAx4eHgihBQsW3LhxIzU1lc1mm5qavnnzRl9fXwJvBgAAOoGwtlLtzh/bflp0KSHDwNn0+Moec0ntr5OLIEQI6evrx8bG/vPPP9evX8/Pz6fT6W5ubrNmzTI0NOz44EOHDh06dGirRm1t7Tlz5nR8cAAA6GzC2srK0JXU/qNGegWMlHYxXU9eghAhRCAQ/Pz8/Pz8pF0IAAB0I8LaysrQVbTBE2geE6Vdi3TIURAihMrLy1++fJmVlaWhoWFvb+/i4gInuQMA5Jl4LijPKYjkJwh5PN66desOHjzYculmnz59Tp482e4lo50hKysrNTXVysrK3t7+S30SEhLy8vIMDAxcXFyIxI8/wbq6uqdPn4pEIicnJwMDg+bOqampWVlZ2trabm5uzWf0s1isJ0+e8Hi8Pn36mJqatnp1DQ2NAQMGNF8inMPhPH36tLGx0dbW1tLSsrlzbm5uUlISnU4fOHAgjUYTN/L5/GfPntXU1FhaWra8uUdxcXF8fDyFQhkwYEDzyZRCofDFixcVFRVmZmYODg7NS5bKy8tjY2NJJJKbmxuDwRA3YhgWFxdXVFTUq1cvJycn8YXrEELV1dXPnz9HCLm4uOjq6ja/4uvXr3NycvT09FxdXUkkkrixvr7+6dOnQqGwX79+RkZGzZ3fvn2bkZGhpaXl5uYGt5YE8kNYU1F5eBVtyETakAnSrkWqsC63Y8eONWvWfOuz3r59a2tr2+4XnTp1KkIIh3AO+o6THQI9LUaqKakjhJSVlVNSUto9rGSFhYVpa2vPmDFDX19/+/btn3ZoamoaO3asqanplClT3Nzcrl69Km5/8OABg8Hw8vKaNGnSoEGDxI0CgeC7774zNDT09fUdNGjQsWPHxO2vXr3S1dX18PDw9fW1t7dvHnzRokU6OjqTJ08eNmzY3r17xY2ZmZnGxsYDBgzw8/MzMTHh8/ni9k2bNjEYjIkTJ44YMaL5p1lUVGRtbe3o6Ojv729paVlaWipuP3jwoIaGxrhx40aPHh0cHCxuZDKZLi4utra2AQEBdnZ2qamp4vYzZ85oaGiMGTNm7NixAQEB4saGhgZPT08LCwt/f38nJ6cHDx6I26OiojQ0NLy9vSdOnOjt7S1u5HK5kyZNMjY29vPzGzBgwMWLF8XtT58+1dTU9PT0nDx5spubm7hRKBTOmjVLX1/f19d3yJAh4eHhrbZ5aGjowoUL2/YD7HbYbLZAIJB2FbKnvr5e2iV0BUF1eenW7xueXJXIaHw+n81mS2SoricXQRgVFYUQUiRRDvgefbbsjfjf3R/jBpp6IITc3d3bN2yzwsLCiooKDMOampra/VFgsVhqamoxMTEYhmVmZiopKTGZzFZ9NmzY4OHh0dTUJP5SJBJhGFZfX6+pqXnlypWWjRiGHTx4sG/fvs3/pcXtPB7PxMTkyJEjrTpfuHDB1NS0+RWb2x0dHbdt29aq8d69e5qamsXFxa3avb29lyxZ0lytuD0xMZFOp2dlZbXqPG3atBkzZgiFwpbtOTk5VCo1ISGhVeclS5b4+PjweLyW7RUVFSoqKnfv3m3Vedu2bQMHDmz+QYjbWSyWjo7O+fPnW3X+448/bG1ta2trW7U3gyCUQz07CCsqKjIyMvhVZRJMQUzGg1AuTqg/c+YMQmia8ywnI7fmRgpJaf3IHYokyosXL7Kzs9s3cm1traen5/jx452cnEJDQ+fPn3/r1q32DfX06VNlZeUBAwYghCwtLa2srKKjoz99IytXriwuLk5MTORyueIDnHfu3NHS0ho9erR4z2HzUc8zZ84sW7aspqYmPj6+qalJ3P7ixQsOhzNz5sz4+Pi8vLyWnRcsWMDlcl++fNnY2ChuT0tLy8zMXLJkSUJCwvv371t2nj17NolEio2NraurE7eXl5ffu3dv7dq1ycnJ6enpGIaJ28+fPz9lyhQGgxETE1NdXS1u5HK5ly5dWrt2bXp6empqqlAoFLdfunRp+PDh5ubmMTExLe/XcebMmTVr1uTk5CQlJfH5fHH79evXbWxsBg4c+PLlyw8fPrTs/PPPP5eWliYkJDRvpfv379Pp9IkTJ8bFxRUWFrbsvGTJkvr6+ri4ODabDceMQc92POKIz0Dn1YETvFz7kgeOk/c9ov9Pto8RVibWCnmir3ZLeZWKEHLt5d6qna6obKNj/7oo/vnfcbTB6l8ZBYdT7U1VZJBbtu3atUtZWfnhw4cfPnxwdnYmEAjh4eGtnjdz5syrV6+2ahw8ePDNmzdbthQXF7c8l8PQ0LC4uLhlBx6PV1RUdPDgwcbGRgzDKisro6OjTUxMsrOziUSii4uLqalpfHz8rFmzdu7ciRDKycm5cOFCaGgolUrNycm5efOmvb19dnY2nU7v379/r169kpKSRo8eHRERIe4sFArPnDmjpaWVmpp65coVd3f37OxsNTU1T09PLS2t9PR0R0fHv/76i0AgZGdnf/jwYdiwYcbGxnFxcWfOnBk9enROTo6ysrKfnx+NRsvNzTUwMLhx44aiomJ2dnZ9fb27u7uZmVlsbOzhw4enTp1aUFAgEomWLl2KYVh5ebmiomJ0dLSqqmp2dnZVVZWbm5u5ufmLFy+2bdsWEhJSU1NTXV29fft2DofDYrGampqio6P19PSys7MFAoGLi0vv3r1fvny5dOnStWvXCoXCvLy8P/74o66ujkAgfPjw4c6dO+bm5tnZ2QoKCq6uriYmJgkJCYGBgb/99pv4jf/9999HjhxRVVXNysqKiorqVseMAZCsg3v3/DPKkkTA743Pf95IHC/teroJGQ5CTIQ1FDaJBNhXe3I5XIQQmUD69FskAhkh1FjObizmfHUcRXVSqyB88eJFSEgIQkhPT8/IyMjMzOzT29AfP348MjKyVeOn17Lh8/nNC0AQQkQikc//nxtrcDgcoVBoYmIizto5c+Zs2LDh7NmzbDY7PT397du3FhYWRUVFVlZWU6dOdXBwYLFYNBrt1q1bOBxu9erVK1asiI6OZrPZ7969e/nypZubW3V1tfiom5eXF4vF4vF4r1+/xuPxv/766+LFixMTE9lsdklJyZEjR8aOHctms21tbS9dujR16lQ2m81isZKTk0kk0qlTp4KDgwsLC1ksVl1d3Zw5c4KCgvh8vqur69GjR3/88Uc2m52Xl5eWlqakpHTz5s3p06dPmjSJxWIJhcIRI0asWLFCJBKNGDFi//79mzdvZrPZmZmZ7969U1dXj4+P9/Dw8Pf3b2pqQgj16dNn9+7dGIb5+/tv27YtLCyMzWa/ffs2KyvL0NDw3bt3ffr0CQgIMDAw4PP5enp6t2/fRgiFhISsXbv20qVLbDY7IyMjJSXF1ta2tLTUwsIiICDA1dWVxWKRyeSEhAQcDvfLL78sX7683TenBEAGYELx1EGIJ8D+j2YyHIQ4PM50ku7X+yFkecwivzQvreyNudb/3IleIBK8q0hHCLnPcO7dT68dNTQ2NtLpdPFjPT294cOHf9rn5MmTSUlJrRotLCyWLFnSskVXV7eysrL5y4qKipZrIBFCysrKNBptxIgR4i+9vb23bt0qfmKvXr0sLCwQQoaGhtbW1ikpKQ4ODnp6eiNGjBB/1r29vcU3ZdTV1VVRUXFzc0MIqauru7i4pKSkeHl5iYsXx7O3t/eaNWtEIpG4AG9vb4SQkpLS4MGDU1JSpk6dqqenZ25uLl6K6e3t/f3331dVVenp6SGExOWRSKRhw4YlJyeLX1FTU1N8jw5vb+/6+vr8/PyWnfF4vJeX18uXL8WdnZ2d1dXVEUKurq6KioqZmZlubm54PF7cGYfDtXwv1tbW4mm0hYWFkZHRmzdvzMzM1NXVW26lVatWiTvr6enZ2tqKH/fp0yclJcXV1VVXV9fLy6t5Kx06dKhNP3gAZJCwumKuQy+/66n66qoiVW3xnVMBkpOLbosv8nk6LqKioaxl+7GY0Bp2tampad++fds3spmZ2du3bxFCdXV1Dx48yMnJ+bSPsbGxwyfMzMxadevfv39OTk5JSQlCqL6+PiEhYdCgQQghoVDYfMqHh4dHfn6++HFeXp44Tjw8PCoqKthsNkKIz+cXFxeLA+yznd3d3TkcDpPJRAhhGJafn//Zzjo6Ong83tHRkU6nt2z/bGclJSVVVVULCwtdXd1POw8dOrS5MT8/H4fD6ejoaGlpWVtbf3bkgoICDMMQQkwms6GhQVdXl0gkuru7f7ZzSUmJeN7MZrPLy8s/La85dIcMGVJdXd3Q0IAQEggEhYWF/7GVAOh5hNUVlYdXBS1b+Sgl89Dft289etJ8WhGQi1WjPB7P0dERIaRCUf3Ode7mMb+u9NroaOiGEMLj8VFRUe0bFsOwR48eaWpqhoeH+/v7T5o0SVtb+8yZM+0ebe7cuW5ubpGRkcOGDfPz8xM3Hjp0yNHRUfz42bNnmpqaBw4c2Ldvn4aGRvOCycmTJ/v4+Jw6dcrX17d///7ihYJv3rxhMBi7du0KDQ3V0dH5888/xZ2Dg4OHDBly6tSpoKAgKysr8UKvgoICLS2tjRs3RkRE9OrV6/Dhw+LO69evd3R0PHHixMKFCw0MDKqrqzEMYzKZ+vr6P//88/Hjx62srDZt2iTu/Pvvv1tZWR0/fnzlypUMBqOoqAjDMDab3bt374ULF544ccLJyal5Bea5c+eMjIwiIiI2bdqkpqYmXl8jFAqdnJyCgoJOnTo1ZMgQf39/cec7d+7o6OiEhobu2rWreW0thmFeXl6TJ08+deqUj4/P8OHDxWs+4+LiNDU19+3bd+DAAQaDcePGDXHnwMBAb2/vU6dO+fv7Ozo6itegZmRkaGpqbt++PSwsTE9P7/Tp061+KLBqVA71vFWjgqry0i1BDc/a/7vuq2R61SgOw75+jE2ydu7c2dDQsGPHjm96VlpaWkBAgHj61Q7l5eVTpkx59uxZy0YqlRoWFjZz5sz2jSkWFxcXFRVFpVKXLVuWlJRUVFTk7+/fvqGEQuHJkyeTk5Otra1/+OEHMpmMEEpJSUlPTw8M/HiT6NevX1+4cEFBQcHPz695Isvj8Y4fP56enm5tbT1r1qzmU8LT09PFK2bHjx8vXo8qfpWzZ88mJiaamprOmTOneddubm7uyZMnORzOqFGjPD09xY0Yhl26dOn58+eGhoazZ8/W0NAQt3/48OHYsWO1tbXDhg0bO3Zs81u4cePG3bt3dXV1g4KCmmdXVVVVkZGRZWVl7u7ufn5+zUcmHj58GBUVpa6uPn369Obz+hsaGiIjI/Pz811cXKZNm9Z83DQ2Nvby5cs0Gi0gIMDGxkbcyOFwjh49+u7dO3t7+6CgIPEWQwglJydfuHCBSCROnjy5+U6TfD7/5MmTb968sbCwmD17dvMdlbOysk6fPi0UCseOHSuehbd0+PDhjIyM0NDQb/1pdgdNTU1kMrnlsWfQFg0NDc3/L3oAQXU58/Aq2jA/2qCxX+/d7lcRCPh8/qeLJGSCvAQhQgjDsNu3b0dFReXn51OpVDc3t6CgIC0trXYPCOQBBKEc6klB2DUpiGQ8CGV4scy3wuFwPj4+Pj4+0i4EAAC6gqCypPLwavqIAJp756agrJOLxTIAACBvBJUllYdXQQq2hRzNCAEAoMf78/y565f+tLO0nEYs0Rw9g+o+RtoVyQCYEQIAQA8RfefO6Z2/zKXX1D+L2p/VACnYRjAjBACAHuLZg7vTequZq9NMVKlTH+VKuxyZATNCAADoIVwc+lxM/1BU33Ty7Qc399ZXVwZfIkszwpqaGvEVogHoMi9evBBf8g2Abk5QWeKa/aDcP3Dvqy1Bx+8AACAASURBVFR7N5+dGzZKuyKZITNBaGRkNHny5MTERGkX0lF8Ph+ubNQO0tpudDp90qRJXf+6AHwTQUVxZdgaFZ+Z81xHzJN2MTJHZoKQTqf3jAsi96RzdbsSbDcAvkRQUVwZtlplzPdKLl7SrkUmyUwQAgAA+NT/p+AsJZfP3P0GtAUslgEAAFklKC+CFOw4mBECAIBMEpQXVYavURk7W8nZU9q1yLa2BmF1dXVsbGxSUpKmpmZwcPCnHUpKSloew/P19XVxcZFMjQAAAP5fXV2dQCBQEbAhBSWlrUF4+vTpS5cu4XA4IpH42SAsKyuLjIxcuXKl+EsFBQWJ1QgAAAAhhNDurZsvnT5OxiE7NcqBw2GQghLR1iBcunTp0qVLjxw58ueff36pj4qKyqpVqyRUGAAAgP/BYrH+Onn82jhrHEJz72UVKxtYSLuknkGSxwjr6+tXrlxJo9F8fHycnZ0lODIAAAChUEjCYeIbW1MUFfh8vpQL6ikkFoRUKnXKlCkGBgYFBQWenp6hoaFfuvN7bm7us2fPMjMzxV/icLjt27cbGBhIqpJujsViNd+iHbQdbLd2gBvztk+3/bBhFUWODKVZ0RlUCoWkb2pkZNTY2Cjtoj4S35hXKBRKu5DWlJSU8PivnB8hsSC0srIKDw8XP7a3t9+4ceOXglBLS8va2nrq1KniL8lkcq9eveTnmKJQKFRSUpJ2FbIHtls74HA4CMJ26J4fNn5JbsOpbb8eDi+h63O5XFtbW2lX9D+67R3qv5qCqJNOn3BwcPjw4QOGYZ/9q4pGo1lbW/v7+3fGS3d/eDy+LT8Y0Apst3bA/z9pFyJjuuFG45fkVkdsUJkUrOQ4tHseF5TpD1tHi7506VJlZSVCiMlkilswDDt79qyjo2P33LcAAACyhV+SwzyyTpyC0q6lZ2rrjPDu3bshISH19fWNjY1mZmY+Pj7iswanT59+7949Dw+PnTt3Xr161dzcvKCggM/nX7lypTPLBgAAucAvyWH+sUF1yiJKn4HSrqXHamsQuru737t3r/lLKpUqfvDu3TsdHR2E0N69e3/44YeSkhJNTU0bGxu4wQIAAHTQxxT0+xFSsFO1NQipVKqpqemn7cbGxuIHeDze2tra2tpaUpUBAIA845fkMP9YrzplEcUeUrBzwbVGAQCgu8Aw7NyZ07GPHw506utZHg8p2DVkcoUPAAD0SEePhN/4fbs3K+tK6K+3yCaQgl0DZoQAANBd3LtxfbG9jrGqkrIC6ejrt3OkXY+cgBkhAAB0Fw42ln9nldZx+f/kMJ3cYDrYRSAIAQCgW+AXZ39HLlOwdVucVK8+dOLin1ZIuyJ5AbtGAQBA+vjF2cyIDZqBy3ba9Zd2LXIHghAAAKSMV/S+KvIXVf8lFEhBaYAgBAAAaRKnoFrAEkVbSEHpgGOEAAAgNZCC3QEEIQAASAekYDcBu0YBAKBLcblcAoEgKn7PPLpJPXC5oq2btCuSdxCEAADQdVYtXXz3+lW+gD+7j8GC3aGQgt0BBCEAAHSRzMzM1Ee3o8bbCETYuKupiywcpV0RQAiOEQIAQJdpbGxUJ+MRQkQ8jkpR5HK50q4IIAQzQgAA6DK2qgpl1fVrX+TVC5HzkGE0Gk3aFQGEIAgBAKBr8PLSak/tuHHzRiKziUaj9e3bV9oVgY8gCAEAoNPx8tKqjm9T/26lgkW/QdIuBrQCxwgBAKBzcXPTmMe2qn+3SsGin7RrAZ8BQQgAAJ2Im5tWdXyrxszVChawL7Sbgl2jAADQWT6mYNBqBXNIwe4LghAAADoFN/dt9YntkILdHwQhAABIjEgkOnf2TNab1FEDnExfX9cIWqNg7iDtosBXQBACAIDEbFy9svTJjSG61J8unNx/KFQPUlAWQBACAIDERN+8ccmrFwGHE4iw+6/fDpwk7YJAG8CqUQAAkBhTA72HeZUCEfbwA8uuL1xKVDZAEAIAgGRws9+stVG+x1cLeFRoPTZgsq+vtCsCbQK7RgEAQAI4Wa+rT++yXLjlfO8+0q4FfBsIQgAA6ChOZmL1md0aszYo9LaXdi3gm8GuUQAA6BBIQVkHM0IAAGg/TmZizblfGXM3k42tpV0LaCcIQgAAaCdxCmrM2Ug2tpJ2LaD9YNcoAAC0B6RgjwEzQgAAaCsMw+bPmhnz5DEBhzYPsRy9KRRSsAeAIAQAgLaKjo7mZSTcnmDLZPOCn+VPgBTsEWDXKAAAtFUVk2mohEcIqVNIPKFI2uUAyYAZIQAAtJVnb93dGR/YQpRZx5s8dZq0ywGSAUEIAABtwsl4hW4effjw/vN3RQH6+i4uLtKuCEgGBCEAAHwdJ+NVzfl9GnM3kY0sJ9o4SbscIEkQhAAA8BWc9Fc1Fz6moLRrAZIHQQgAAP/l/1NwM9nIQtq1gE4BQQgAAF/ESY+vubAfUrBng9MnAADg8yAF5QTMCAEA4F8CgeDin3/W19eNd7RSenhOY95msiGkYA8HQQgAAP/6YUYgvSjNlEacuKfwzp3bkILyAIIQAAD+lZqYcH2MJUKoiovFZBdP6SftgkDng2OEAADwLzIRn1fLFoiwxKomU1NTaZcDugLMCAEA4KOmlGfbB5uvf1Xc2MQNmhvs5AQnzssFCEIAAEAIoaaUZ7VXwodsCr9BZ9DpdGmXA7oO7BoFAICPKcgI2UHSM5F2LaCrwYwQACDvmpKf1f4dzgjZQdI1lnYtQApgRggAkGtNyc9q/zkCKSjPYEYIAJBfTUlPa6/+wZi/HVJQnsGMEAAgpyAFgRjMCAEAcoTH480K9M94k6KipLRnqHnfVQcgBQHMCAEAcuToH0eMqrKv+1gutqDvzWiAFAQIghAAIFdKiwqtVRUQQjaa9NJKprTLAd0C7BoFAMiR8c42c04fLWrgPi5jf79kpbTLAd0CBCEAQF6wXz82SHtw5frN52nvdtjZubi4SLsi0C1AEAIA5AL79eO6a5GMkJ26Okbmru7SLgd0IxCEAICejxV3t/7mScaCXSQdI2nXArodWCwDAOjhWHF362+dghQEXwIzQgBAT8aKi66/dZqxYCdJG1IQfB7MCAEAPRakIGgLmBECAHoUPp9/9epVHpfrbagseHxZc9FeIkNP2kWBbg2CEADQo4wb4WklrKbg0MH8iifxiZCC4KsgCAEAPUdFRYWwunyFpylCqIAjyiqvcYAFMuBrIAgBAD2HiopKaXUdiy8k4XG5tU3a2trSrgjIAAhCAEDPIUh6tGKwTeC9HBFCS1eu1dHRkXZFQAZAEAIAeghW7O36exdmHDr/PUNX2rUAWQJBCADoCcQpqLlwN1EDUhB8GziPEAAg8yAFQUfAjBAAINtYMbca7v8FKQjaDYIQACDDxCnIWLibqAHrYkA7wa5RAICsghQEEgEzQgCALKmvrw+eOSMzI83WSG9rfyPDpb9BCoIOauuMMC8v78KFC+vXr793796X+sTHxwcHB8+ZM+fJkycSKg8AAP7H1vVr3bhFUaPNrYTVpzlakIKg49oahOvWrTt9+vTFixdfvHjx2Q5paWleXl42Njb9+/efMGFCTEyM5IoEAICP8rPfu+jQEUKueqr5hUXSLgf0BG0NwvPnz9++fdvZ2flLHQ4ePDhz5swlS5bMnTt3+fLl+/btk0h9ly9fHj58+KFDhyQyGgBA1o13tV//OPOfzNItr0r8g2ZJuxzQpZKSkhYsWHDhwgXJDiuxxTIvX74cOnSo+LGHh8fLly87PubmTZsWz11qKex3cNvhyRMnd3xAAIBMa3z89wh85Zbw43yvGQfOXR7tM0baFYGuc/fu3VHDfFCa4o6Ve4KCvpfgyIRNmza1vfeVK1c0NDSaA6+lzZs3f/fddyYmJgghFov122+/rV+/HofDfdrz1KlT0dHRDx8+/Ouvv/7666+LFy86ODgoKiryPvHD93PXeG4fYTl6kOnQsFsHcN4LD6WJvBicGwXC9QkiezqnsJ7/YwxGx/NU8LxZz7BKFt+ayl0eJ3pUwh+oxj34VtgN+w9Rqb9bSug+9chK/x9e4KqahN2nHpnoH5ou8tLk3CrCukk9HezPfHCl4untfL+ttnY2EcIBJBXNTqpHScgKjiNJ/f3KXP/DGZi3Fq/zxv9lyY8L7OeMtfP1svTZe3nb/IXBn6bGp4hE4meTqFOCMDQ01MfHp3fv3gih4uLis2fPrlu37rODxMfHYxgWGBhoZ2dnZ2fXr18/Z2dnCoVC+kRkRASDrG2tbZde9qawJmf5T3NNVYm2DJKKIkGDgnfWJqsrkZRI+P46JAaVpEQiOGsRdOhkOhlvp04wUSUrKxBMVAjdrX9vukhLhdJ96pGV/nghf4C+YvepRyb6GyiJ7LUU1JSI3aSejvS3y4givr6b6bvTqbdOZ9dDJ4pUqRTZ2j5S708j4YyoqI82RfLj03D0Iq5yQrUwr7K84UNffec85vvnxY9+XrXi09T4FB7/9R2fOAzD2h6E06ZNs7Cw+Gx2Dho0KCgoaO7cuQihq1evrl69OjMz87OD7Ny5s6GhYceOHV99udzcXDfH/gRExAjYoUkbzMyGWU43phlR2l5wN9TQ0ECn06VdheyB7dYOTU1NZDKZQCBIu5COanj8N+vFTc0f9xBUNLri5eDD9u0EAgGfz6dQOvT7+cnjJ6uXr+HxeKvWrfTz86t9x2Km1FW/bVDSVSjRKT7JPh33y0NuPRcR0OkLJ0eMGCGp4jt0HiGTybx///7UqVMRQpMnTz537tycOXPwePyZM2d8fX07XpypqWllbQVC6Lf4sIqKsn4vr6VFjtfur2E0UgtP/MpUFwDQMzQ8usKKudVlKQikRSgUhsxesMs7jEKirFg1XylWx7RPL0ZfFfXR1N3Jv/MEvMNjdjOCOuUz0NbFMtu2bTMzM4uKijp48KCZmdnZs2cRQu/evQsMDBR3mDt3blNTk7Oz84ABA7KyspYtWybBKoPsp0Y35Sh4mxkQQtkldcn7shsKmyQ4PgCge4IUlBNCjijnaSEdr6qupEEhKVnoWqtPIdn80KvEsCj48fLeqib7vbYxlDrrM9DWGWFISMi0adOav2QwGAghJyengoICcQudTo+JiUlMTBQKhS4uLkSiJK9Zw6CojzL1/Aer+cF7LOH+dpzn5oxjBYy+KsZjtJ/GPK2trR05cmQHp+QAgG6Cw+FER0fT6XRnrIr98g6kYE+Sm5v76tUrV1dX8cpKEU9U+57FTKmrTm9QNlYi00kXU05TiNTsxgwH5z4nUi9EZd9ZM2Cpi26/Tq3q244RSkTbjxG2VMetnxEVEjFqn/Kb1/X3LqjO2lnwULQhbDVPyNWkaqfWJjx5+VgmshAOP7QPbLd2kMVjhDwez3OAqwtVwGRxeCJ0/nFc16cgfNjaoS3HCB89fLRs3k+DjYY/K3ywbcMOC5KdOP8YfVXU7ehERUJpaemCH39sYrPXbl17qeYWCU9cP3C5OkWts4uXmWuNqigoTzAffebtxZUDFmFCYd2Jtcbzd2Ruehs+6RxCKPLV7y9evPDy8pJ2mQCADklMTLRQEPzkaIAQmnI7k0NUpEq7JCApEQePLh+4wZRh7mIwMDL86NGwo2a+egSFf4/Q+Uwey7FDeD3cuIAJey7sm+sShP/amQ8SITNBiBCaaj1p+vX5RfUlhoPGIhyqPbJWiOOzeI1KZGpBeR6FR5N2gQCAjtLQ0Mhn1oowfY5A1MATKioqSrsi0FGCJmH12wZmSh2ulFyCFZkyzIvqCkycDDUdVVt2q6qqqhc1GAy1QggJ8jgOIquuSUEkW0FII1OnWE84kXr+l0E/09zH4gjE6XHvvz87CUNoaP9h5AS1bNYHkwk6BDLcWwoAWaVblDLQRMfnWhqBSNy6d59s7dcFVy5d2bV1t4KCwu79u9yc+lenNTCT6+pyWSqmVEZflX1XdgUGBJ5PO6apw7h47K9Wz2URm6oqqxiVbIIikZPfaGpq2mVly1IQIoT8LMdOux6cXZPbW82U2n/Ulcy5VyZZaFLJcx68VN6LFyYLk3/LMZ+qr2yiJO1KAQDfrOHBRVbc3c2X7m5TVpd2LeCblZeXb1+3c6d3KJvHnu3/w4mgy+pWdC1nVcuZhs3zk0cvHn72udG5D8Nen1i9e+3Vw5fZXG7onlAtLa0uq1zGglCRqBho43s89fwOj/UCgYBAIBqpUBBCTppKeSU5o74bVfWmPvNUoWY/lV6jtfEwNQRAdtRHn2tKeqq5aC+B3umLI4DECdjC1OgMU7oFlUyjkmnqygzLpToa2l9f6MTmN+1/dSSrOnuf11YzVePVU37qgmpbkbEgRAhNMB91MeNqOjPLhmGpb2oWllysRyHcf1eyztwYIaRhr6zSm5p/o/z13mzzAH2MwX/y5ImxsbGDg4O0CwcAfFHDg4vs1481F+6GFOzmMjIysrKy3N3dNTU1kfj43//v/9TT65VV8/ZOxnUWv0FBndSWFMyqzt7y/FcbhmXEqH2KRIXOL//zvu1aoxLx/PlzHo83fPjw9j2dgCdQSIp/Z90Yaeo50c+/QECqYxhvnjeDGH2CYueGV6LjSXh1WzpFk/zyWIpfyCTuO+z02dONnHq3AW6SfSPtw+PxFBSk9vOWXbDd2kEgEBAIhLZca1G66u+cZSc+1Fy4m9A99ojCh+1LTp04vX7xxoY03pbfNroaD6x7zs27Vibkihh9VUwma2s6qUydOTWjJlXHjrF3/x4ymfwfQ2EIu5J547e4sPmOs4LsA4h4ac7KZOY8wpaEmPD7Gz8ud13QT9u+uZEVF11/6zRjwU6StpG4JfRAWOHNivH2fgKRYNndOUnpiR0tXRLgFKX2ge3WDjJxHmH9nbPspKeaC3cTlLvLXBA+bF/Sv+/ArYMPUEhKj9/frSKXrdu4Xt2ahifhUdvOI2xqapq/ZEFM7Au3/m46U0waRexfBq3Qo+l0Vflf1N3/VPwsAo4w0z7gaMqZlo1Ut5HKPjOZYWv45YXiFk0djQpOKUKoilVJFJKEPJEUagUAfFn97TPdLQXBp/iNgrLY6rdH8okshfL6UoRQGeuDxTATRh9lcQq20ZadW+O5SXqrrV6yk7Kuvj08cnd3SEEko0GIEBrey4PFb3r5IaFlI9VtpLJPEDNsDb+sECHk6+fL1WoIuTZ9y7OVa2dtTNqTXfu+UUr1AgA+KisrCxg/xsXWcsOMiU1vYrQW/wop2H2EHQpztncZM2Jsdna2gCWsSKhNP1qQsONddVqDlrNq+KWD+xO3Lrg+o0ghe/bc2d86eHJaMr2PBkJIvZ+WqIJPwHWXfRWyt1hGDI/Dze4zLSLptKuuU8uTLqlu3gghZthqxoJdJB2jC1fOYxgmvitjdXrD+wsl6nZ04zE6La9lAADoSot+mDWJzBw43Hjdk6QXY6ZMpipLuyLwUVxc3KWIf/YMO5JfnTtt9MyD/sfUbOja/dWsvjcS3/BHC6kmvk1o/qX6TQrrSzjmqOivdwwPvcYXVUsWb+mEd9BOMpwHgw37kwnkp0UxAoGgZTvVzVt57CzmkbWCyhKEUPMPTN2G3u/n3iIulrQ3u/Y9SwoVAwAQysvNHWSoTsTjhhqqZ757L+1ywEe8BkH8zSQHTRcyUcFCy7pJyHLdbGU53UDDTrnVbe/akYLRuQ9/vLtq4fchx7dEjCJ5/LEpbPrUaV9/WleR1RkhQgiHcL2KtSbMH6tMoY8fNf7w/tDmb1FdR+CIpMrQlQrfrX/5Ls/ExMTa2hohRKQQzAP1azIb3/9ZrGZFNxmv8y43Ky8vb/DgwXBsHICu4WlrtvFpxhBD9cis6ohtE6VdjtxJSUkpLS318PAQL2zhs4Q1GQ3M5LqGAra9er8jOYe1aNq5de8d3fpJ5LavLD77t/iw3Jr8A17bTVV7ISs0etTojg8rWTK5arSZkZWx0UpbggKhJOLdxf3nHR0dW3634P7VsTNmeZrrpVY1+YcsCVm0pPlbQo4oL6rs6OnIp+/vW2nbvSp98eD5/a65kAEsSGsf2G7t0A1XjdbfPMVKi3vI6Jf17v2kgMB+/Tr39jrt04M/bFt+2fro6tNeqqapVQmXw69x3gkbCthq1nSGg4qqJQ1PxGVlZZ0/fcHI2PC7oO/++/yHVj67ajS5/M2O2AMehgPn9Z1JIpAk/W4kRoZnhEKhEOERQYGAECKpk2tqalp1uPWuZIadwXQ7fb5QNCX8cMsgJCjie0/Ru7326u9jThDxRM007cuXLi9YuKBL3wAAcqYu6hgn67X2wl3T4biglFw8dyls/Fk8Dn86PuLm3zenz5tmPcsIR/h35mdpabl5+6Z2j49hWGpqqpKSkomZ6Yk35+/kPFjVf7GrnuPXnylVMhyEBAJhlOfIR0efCugYVsBzd3dv1UFNTT2LK0IIVTXxFEmf+WOEqqxUw67WpGmV15dZcHp1RdEAyKu660c575I1F+zCK/XMyVZ3xqvjM1PrmSn1oiaMxW2kKypXCSqtx5qp20jyZyESiUaOH5XPLuazeCqGqsNDRkeM3q/R+XcT7DgZDkKEUEToH/Hx8dkfcs5zojg4riL6nzu2+E2Z8s+f5ybefMsXCHYMt+OX5pN0jVt22H9439ygeXgRwdTE1K6pf+apot5T9IhK3Wg/EgA9Q93Nk5ysJM0FOyEFOxWGYZs3bIm6GmVnb7f/8D4agV6d3lCd1iDe/2kwjHHA7Ldli+YREcnF3bndl/f6krS0tAJ2se5sM4RQ1vbEzW4raRTZuDuebAchQsjV1dUVuVYnsk6kXljmMr/lt4hE4l/XbvD5fBKJ1JT0lHlkHWP+9pZZ6NbfLTUrRdxByBXlXS9L+i279xR9NSvZ+OEBIBPqbpzgZCZoLtiJhz2inezC+Qtp99795nX0afaDOaPnr/berG6nrDtIo3n/p7etd1r2W/EvPYm/Oh8vaGhs0EEIYRhJRCASZCZfZPj0iZaC7Kc+KXyRV1f46bfEP29KvyGqE4OZR9bxS/M/24GggO89Rc98qn7O5Q/v/yoRcISdXjQAcuBjCoZACnY6bi0/7laim/5gIp442MyzsC7PdZOVeYC+mhWt5VFA9P+/9CQrqfzN/sJII6Ne+Xve5u5IXRT8owzdVFlmEvu/0cm0abZ+EUmndg7d8KU+lH5DEEKfzgtbUjWnOa7sXXi38vWu92Z+esqWlOfPn6urq8PNKwBoOzabHRsbq6+vr5sdy8lMgOOCkhUfH8/lct3d3cWXU+fW8qtS65kpdewK7kDzwXvP7sTjcLHFT8dPGdcq/zqJQCQ88/ZiVPadJf2CPSYPrKmpIZFINJos7VfrIUGIEJpkMebau9sJpcnOun2/1IfSbwjC4f47C/FkvPFYbXVbetrZ3MWnf+itaVXJLus71P7XA3s7q3QAepDa2lov9/6u6qSMsmovG5M1F25DCkpQUOD3Vdl1FKLSVtL241vOVL9pYFdw1a3pBp6aala0/gRr/eGMG1dvTpzkEzgtsAvqKagr2vriNx2aduTI/Up4CkJITU0GVse00nOCkIQnzu37XXjSiUidA/gvX/iA0ncwJhIyw9cyQnZ8KQsRQsomSizXcvMbVgv6r0AILbg+Q/CrgEjsOZsLgE5yIyrKR5s0r4++QKQ3+c77dZCCklNbW5uRlLXPJxIhtD16XXJC6qDxbiq9qTj8v7/xBg0eNGjwoC4oBkPYjey7kclnZtoF+FmNE59H2AWv2xl61G/2oUbulzOj7uY9GmXq+R/dlByHIhzuq1mooq5cK6hBCHEFXG4TF+P1sK0FQKeg0+mVLA5CqIEnIJHhrn6SwanmMVPqSxLK62rqhCIhHo9vwNXa+pqqmnTpHsiamppFKxa/SXszYcJE/gBCNac2bOQeA7peV9bQGXrar/aFTrM3PN011GigIvG/jtMq9fPA4fFpe5f9+p6TV1w8fdacBYuXturj4uJi6mwUcm26QCRYPOOnpD3ZZr56GvZwwB+AL8OwQYKS0zw0+WZGkwj3e8RRaRckY7hc7uqf1sQ8jxk8dPCOPdsxFq7qTT0zpa6pgqtmTTf3MVqmuGT+vkA8juAX6GtiYtLF5c1bFJyumqc2V+/omRPjFcYfXnqAiO8J55v1tCC01rDoo2n9V8a1IPuA/+5JcRj8S9xyfy3MdaDOimOh1vYOw4YNa9Un9I9DXC6XRCLh8fi6HNb7v0qqUutNJ+nCuYYAfAaG1f5zRFiQeelFEhdHUFRUbMfVmeXc7u17eOm4nUMOX0g8ucp341THIIa9cq/R2sqmSuL9n/Ntg38ImSMSib7p+meSkvo2VXu5OY6A03DRUa5S6hkpiHrM6RMtzesXdCUrqqqp9RXXPlXMrBnai6FEInjpUVOSXn+2j4KCgnhplooZ1XFFbyKV8HrP+6o39RIuGgBZh2G1f4fzCrMYC3fhlWgUCgVS8FtxqnjxDxIHGw8nExU8TL2K+bmuGy3N/PRaHQUkEolSScF0ZhbBnFJ4+V1dRlX9g8oJY8Z3fQ2dpAcGoQ5Vy8fM63jqua/2HDBo8P7EomeFVWdSizz6WH21P56MN52oa/29Uf7N8sxTRXyWsKys7PHjx59e5hQA+YJhtf8c4RW9Z8zfgVekSruabk0kEsXGxqakpDS3cKp5H55WpR7KTTmQM8h2yImEsOTihFPJRybNmNAy/6RIIBKeSL2w/umO8L2HV41d6lJre+K3Y+4DW1/VUnb1tF2jYjPtAr6LCsmtLTBV/a8riO4P/+NYxB+v32Xu85un/fgMz9KcbGjx1cHpxkr9fjIrjK6MXHAm/On+vgbOr0sXXblxycrq61EKQA8kngsWvWfM345XVJJ2Nd2aUCj08RpDbVJp5DaY9jH+Zc5WZkpdUyVPzYpm4KmpZk1zw1ubX+719MGzH7fMnzBhgrTrRQih/Lqi7TH71BRVI0cf0KCoecwbKO2KJK9nBqESiTLddsqh+lEF2AAAIABJREFUxMj9w7f9RzcymRzy4yLxY066Q1XERo25m8lGX89CPAlvPFY7as/lDZ679VQMEgpfHjn0x4HD+yVTPQAyRJyCxdmQgm2RlJREYdGXD9qAEFp4cWb1qFrjMTrKJlTUYuLn6+fr6+crtRJbEJ8gcSzl3ByH6eN6j5R2OZ2oB+4aFZtgMaq6qSa+9PNH/j6laOOqFrisKnIjr/BdG5+iwqDXNdUghGqbqgkcKeyyB0DKMKz2ShivOJsRvA1S8L81VXCL7lUWXayuqqlCCAlFQgGZZ+FrqGz6PykodaWlpYfDw6KiosoaK5bdX38758Fh7909OwVRT50RIoQIOMK8vkGhicd+H2yiTFNuy71JFW1c1QKXV0W2dV64ZdfmKRP8KTiqiCDc7Xwo/WhBb399snKP3aQAiOXl5a1bvqS89EOQu8MoIxXN+dtxCpSvP00+VFdXr16+JiM9wzdg8tKflnKqedVvG5gpddwavoa98vCQQS+4TiHXpgtFgmUrl0plzct/KCsr6+85kOKuyr3Gxh/Brd+5Mcg+AI/rsfOlZrJ9h/r/xuVyLYZYCzh8Apdw6siJYUNbnx3xWZz0VzUX9mnM3UQ2smxL//r6emVlZZEAK7xdXpFQazpJl9FX5T/69+CbX3cq2G7t0El3qPdwcVxqQjZVU1oYnRZx7Y6NQ3e8y3xHdOTDNiNgpjmrT3/jQb8/3entNMbdxIPhoKzhoKLcS6l55tfU1EQkEjvjytcddPbs2W33ftX2NkII5WxLLs74zG0MvuSzd6iXFT056i9evKhgTO29sp/BEqtFPy9u47MUbVzUpi2vitzEK8xqS39lZWWEEJ6IMx6nYz27V2F0RebJQn6joP11A9C9NdRUu+ipalDIQ0200t5lS7ucboRdxk19lTrM3JtKpg0x8aqgF7pssDSZoKtsrNRy/yeFQumGKYgQqlVsrM2uQhjWVM5Spf3XH/Q9TE8OwkY2C08lIoQIFCKXx2v7ExWtP2ZhU156XFzcu3dtPWpI70Xpt6I3zYjyek92WWw1Qig5OTk5ObkdxQPQTWGYnorS6TdFzwqrbhc3DBgwQNoFdTUejxcTE5Ofn9/cwi7jFkZXvN79Pi0iv7+d+7H4wyklr//JujAucEy3Ov73Hxp5rB2xB2LJyeMG+BTuSONdqjl79Iy0i+o6PfmAlr/flP2H95cz88uzS1cuWPFNz1W0dqH5Lx7lOcyQoVrWJHQfN3nzzj1teSKOgDPw1FSzor+7UBy8bB6fz0M4pGpMP3/56+c1AtDdYVjt5cOHpg6/2KT+qrz85OaFBgYG0q6pSzU2NnoO8jJW6l1QkzN1xlRf50Bmcp2QJ9KwU+7tr69srNRXcOBoxLG0lFe7w3fIyl8JCaXJe+IOueo6HvU5QBmviH6VdkFdricHoYaGRmp8Snx8fLowO13wHkMY7lv+PEtgciy01Tb2N8YQGn/prw1bd7T97hNUPUUNX4WqSOYun1CE0IYHSwsKCnr1+q+TGgHo7jCs5q8DAuYH4yV7V5Fl5p6rknUj6oaT6oDpjnMEIkFI6IxJEQEWgQY0o38PjJFIpJCF86VY4TfhCnknUy88LHi2qv9iJx35vetqT941ihBSVFQcMmTIvKHfs/lNDwuef9NzKRRKgxCHEBIIRQKB4FtXHCjRKByMjSEMIVRXV08Q9OS/OUDPh2E1fx4QMEsZ87bi5DIFWR84BbfKS+/UNrDrEUJcAUdRlWw8TqdlCsqEwsLCV69eCQSCNGbmnFtLPjSWHfU5IM8piHr2jLAZHodf4bZw3ZPt/fWcqKS2nu3k5uZGs+znfye2ns0JdjTmF2SRjb/h2jEMBmPi1PELz32HEBrpMbrkVANxdLXOAPX2vAEApAvDai6HCipLGME9PwWLiopWLVtdVVkVsnT+xEkT2WVcZkpdZVKdiC/SsFMO2jDtyeLoFXeC6zl1e3/fLe1iv9mvB34LPRWmqK3EKWc5rBq4wv3HwYb9pV2U9PXk0yda2fXyoDKZtsBx9jc9q66ujkKhiHLf1Jz7VWPOxm/KQoQQi8VCCFGp1Mbipvd/liiokXRGq6jrqX7TIADB6RPtIpnTJzCs5sI+AbOUEbxVHs4XHDpg2GSDGYZqxtsfrF09fpMxw4zRV4XhoEwz+Pe919TU0Gi07rny878Z25garbbFEXDF/2RvDVg/bco0SY0Mp0/Ihvn9gu7mPcqrLfimZ6moqJDJZEUrJ7XpK6qObeblZ37T06lUKpVKRQjRDCh9l5nR9CnvwksrEmq/aRAApAbDqi/sE1SVyUkKNhY3VZZUORq6adK0+xsNaTRnOq+zMB6j3TIFEUJqamqymIJ8kaBJyBEJRAghIp+gSoW/yD+SoyBUVVAJsp/6W3y4+LjdtxJnITPyF27Om/YVgCPgjEZpmQZpfXhalRaRz63ls1isly9fVldXt29AADoJm82Oi4tjVlZWX9gnrCpjzNvSk1KwqakpLi6usrKyuUV8/kPirveZJ4v0tfSvv72UVPzqWfF9z7EeUqxTsnJq80PurBgwfXDerpSi/RkmOIORI3v4hdPaTo6CECE0wdxHIBLcz3/SvqcrWjmpf7eq6vhWbnY7sxAhRNEmOyw1VelNvbnhsbON6/6fDnu4Dnv65Gm7BwRAsoqKigb2tT+ybM4IR9tHMS972FywuLjY1cFt3/JDXgO8r1+4URhdkbjzffrRAgFbaB6g77ze4vLDvxScsCzV1yf+OmZkZCTteiVAIBKeT7/y04MNEyxGX/3lr9w3OXG3Y+5cvS3xSw7JLjk6RiiWWfV+3ZPtp8eFtX3VTCuczMTqs3sZc34hm9i04+nNx7qWzftJr7L3AJPBpXUlJ/MP3Xxwo331yAk4RtgO7TtGuHHtaoM3d0ebaZWzuKvT2Pdi4jqpPKnYvGELPpk2zNy7ml21/f7ai4f+ZjioUPX+ZwVQT/qw5dTm74r9XV1R7We3hQwljc57IThGKEusNMzd9JxOpl5o9wiKVk7qM35mHtvCy0vrSCVUdSWOgI0QauKzEVfufhCg21IgK7D5QoQQmy9UUFSQdjkSgqGGwqb8qLL/Y++uA5rc3geAn3WwIMZGd0urKChlFxZ2N3aL3djY3X29Xr2K106wUMRAFEFapDdGbYzV+/7+mF+uP/Qq4hLO5y/wHs559t53e3biPafshUAkrQEAiCQ1+hZ06+6cellQ18lksvT0dKFQKEe/dAR7O3bbFLpCpVlQ1zWLxyfqmewzZvT16V3tOjgY2DauBrJLS6MxSwsOrjrM13vzMbNzj15zIhdiML+2mdLMuTN6dOp1P+8mT1CyKnxzyqFch4HmJAPdm4GHmhQEGWyKHZpWer1QVCKSHv/zoqYD+mUpKSlrlq2ViCWLVy7ya+P35fmH1xUoAgzd6PM3zR40Ljy+KLashnv8j2OaDlbJeDxeYOcgxAArKKz0Ge/fopX7ke47YAr8qeaYCBkk+hiPITtfHtrVef0v7TXzNZKDxyE+FfPh2Xo3sx0Xj58wYo2dMOGXamCz2YlvE4qLi9lsNhaDLYjjJW3Psu7ONmlrqCv7E0JNDYLw/4hmopKEjE/FvDJjY+OG76akJRAEGdJv6OzWy4gM4pgB4/aPO6VHpBm60R2HWDBsv8yGPH/9rKioSBdf3U/tP3wA05Zq0t5MUl6bdSbl2oK/NB2RbmimI3JhDt3EcvHdnLjfqeR1avrgFubGVNIAe8PEp41Z7YLBYExNTXE4HAaLsehg7DHFpvh5+ftDueJy6e8EBkGNgSD8M5sRYZXR+JUYAtHU1FT38gQKsl7mGuCMnTlutkYOzpwW+PbiVkud7PqZ1mVBBZ18dQ3AFfIQgCh+phH1NBuMDmmmiRCLwczzm3rgzQmBRNjoSjp07bb3bWEKt/rwm09BPh6/HxXVlOw1y07fUS9pe1bhozKAgs+fP79//179C5qgZgdB+H9EIzXVRuNWYAjadVrst6qrq1+/fq3YrUJB8fzDy/XpVQ9kFeKy+JyHr/ISsqo/+rRrLjuHSeSS48nnMuwLah7xS45nF+xN37xG9za+0ZQm+J2ogZwNHdqatpq6eZalzGTYwKEtWrT41RrmL152SN/wr6ePB03tHVT6SpzxluT4u+86RdfQ0I2R8Wf+1m1bn6TFsejGMpr4+t1rTfILLKQVEIR/NhoRVir6gpqO5ifevHkzcuBoZ2O3j9yUMyfOGgpNSl9WYPEYljfTbYI1lUO6MeLa1o3bJBLJ39cuUKmNXByuW1J4aZue7zanmZweclB/JCM9Pd3KyopGo2k6Lp3R7B6f+Nr46RPuZMcxnAyq7pTeu3zH2blBR9J/lzj7fdmxKKPRi0iO3j8u2cCV2XKZvIWt56EB57AY7N6ELZOixoSGhjY6vCagKa1oV5sGPT6hyII1VTrRFwQADAsf3onS14nt+rbg1aO8u9FrtrO8GBS2Mle36tDNpjg+4k5O7MxWk4KtAjQYCXx8QlfFxcXZDXdjtTalBbOu3bz+O1WR7NyNxi0rO7lRnKGcY3ixOCyWgEFQBAAgqhEhNXB0FFIBBOGf3aIbWRAFVdnC7MtFwhyxRCYGAIjlYmMPA8vOxsrNgtrv48ePycnJAIDXxcljrk3n15af6LVHs1lQ1zXr0TYzM/OqDD7dzkCYWuHe95eHRush2bkbjVtedmxtQ/qFP4XBYOZGzp4ZPYZOZrDZbOozTj6FaxFqDBeUQr8vNTV1w4qlUrE4orW9O5upbVkw/mn81g3bGAzG8qhldrZ2Vbk1vKRK3ttKPAXH8mau3btqxJjh7I8m3NqSf25d0XSw6jZp+qQHbx5iSTgihug8zWeO3+T2Fm00HZTOa9ZDo3l5eeOmjs/5lEPyZjw4cMtEj/37dYqzU8qOrTUatYjk9P1c+EujLtXV1TU1NRwOp7ZMkvFnAYqgjkPMKcbN6/uvgg6NVmmP7w6NyuXylq6Oa305RBx2QdzHp+/T6QZadDoYl8vt6N85st3qClH54cSdh0efU+Q/Yx9mXc9PJpMVFRWZmZmpaJMwrb3ZqqqqPAO9rSPdAQDZB99f2vVXa+9Wmg7qCzg0qqusrKzuXbub9S5z+eJlWxL2NG4z7npIdi2Mxi4tO6WcMVI6nc7hcAAAZCOi+xQbIw9G8q7soidlyogUaqaKi4vNqAQfE2YLY7qnpUnWpzxNR/QvFEFf3nvjzvKxMbL3tmhFwlNsxrB8FzpadWV/Pf6Jx+MtLS2b4VaZQnlNdY0AoCgAgCgj6OsxNR1RE9Gsh0brDHbtG5f39GbW/R72nX6/NpK9B2v88rJjUYYjIwsJ+gAAW9tGbmHzNQwWYx7CMvJgZJwv4L6pdBxiTjTEp6ammpqaGhnBnSOghjJhGxeVlT/KYxBxmPf8GicnJ/XHUFFR8fnzZxcXly+HGaGgbvyTghq+LXz5vjCpUlyBpaFsG5b6w9NOcXlPd7481K5n4Ot1iTgCvmeH7o6OjpoOqomAiRAAALAY7MK2M+fcW+Zn6qOU7YiIti2Mxi2bMqhfnlAKcHgnv/Z7jx7//WoBAGQjoscU2+Ln/ITo9/OvTDFlWOSV5yyLWjJw8ECl1A81bahcVnF647Gpg49lCWRy+fmrS9X/gMH1a9cXzVpia2RfIPp86VSMPBdbN//nMc2WYkz6o/+ZnVt2M4zol/b+rebYtFOxsHRrwr4yEX9D8DKX/o4CgUAulzOZsDuoNM16jrCeY8l/ZPCzNoQsV0ptRUVFQzv4n+nqAgAYfS/r5K1YCwsLoLzph4M7D7+/nD7Yd3StVLTwwZRXKS9/v05tprXTNtqs3hwhKpfxT6xH5TKjccsxeI3tatvWJ2CF/xYGmXn13UVAQ2bMnMnyYpANtWi1jvbcbAiKXs+6cyTpTF+nHiPcBxKw2tt10ek5Qu29rOo30n3QxJtzYj89CbVu//u1YbFY+f+mYGVSCRar5OlYEpOAYBEAgByVy8UoQAFcUAr9gMazIIqglZlC3tsqCV+meC4IwSBWIcYWoXDw8/9JTEx8Gv80KDDIwM5oc8IeHAa3u8sGK4aFpuNqymAi/BcBi1/YdsbiuChvjocB+XeHHTgcTkDXnkNuX8fI5Z5GFMOKfGBmppQ4FQYNGnTswPH0uA9FFQWTu856uzvbaYh5c3ugCmogNWfBa9eunTx4ysLGYvnqZQb6BtWfRLykSm5SJYGKY3kzo7atiYycasa0rAYVt8feUnUwuuXS5UtzouZT/PQ3Ht1i29tl3ojZ/Z17YX/xZBvoV8Gh0fr2vT7GF1UsazdXKbWVlJQAAAxEfN6RVYZD55JbtFHiqAuCIHl5eSwWi6ZHK37O/3Sz1DzEyDyEhcE2wbeN9oxW6RDF0CgWoPwT6wAGazh6MQan8u++ycnJEwdNntl2cQY39Vlp3Joe2xT5j91Sn8z6Mv4pFApLS0utra2VPlKiFBq82XqE9+T5iyimNGF+NecV7fpfVzUSRiPo9NCoNt6FmjXOc/iHso/PCpQz5cbhcDgcDtHGlTVhFf/cttqU50qpVgGLxdrY2NBoNIABJv6GXrPsytMEyXtyRKViJbYC6TRULuOfWAewOPVkQVSOPvznSaBlZ2tD207OPQq4Bd5z7BXPP9RlQQCAnp6era2tdmZBDaqWCARUUVVmBQCgJrPKzcFV0xE1F/BGrI+MJ0W2nRH9Yu/vHEzxLaKNKysiqvz8TtlHVa1qUSwo5bTWT96Tk/+AiyJoVVVVeno6giAqahHSWiKRKC0tTSyqqTi1AWBxhqMWKTcL8vn8zMzMul9RBK3Kqcm+XJS45qNplW1s1q1sXsbN1CtO7o4kfXjWdIPEF7wYd31mlzE9rIvZeetTHCosVi5Zoemgmgs4R/gd3mz3APPWB5NOzvObqsRqiZaORhPX8A4tp1Ao5BZtlVjzvzDAxN9Q35mW+VfhpXMxh+7vtjS0qUDL7sTdguOKzUdycvLI8D5O+pS0wtITM0f7zt+k3Cx4eP/hvdv2s+kmWH30j0PnK94JuUmVFBaR5cX0nmvvx3Shd8SeOHjS2s362PIjSmy3qSoSlGxL3F9Ww18duNCN5Qw6ztd0RM0OnCP8PqG0Zuz1GQvbzmxpouTzzCpSX9ec3Ww4fB7ZtbVya/5/UODr3HpD5916JNrF5LNuA20nRUxSYXOqB+cIG25keN9hpFIPNuNxXlm8ofvuI8p5hrVOC3uPvWGn8Vj83sdbg1oH9xkSxvJiEJlNp+entptNjsovf7xx+v1f2v90xE/p9ByhDl93ldIjUOf5Td0cv2uCwVATFsfDQwnn7irgLBxZE1fxDq8yHDaP7KayXIgBePKXcW9UDsQV8Mj7ZuT/DYYrb8EhiqCK9Z9SgUzxBRpHxNiEccyC4MZGv+DajevLo1bQ9PTmL4+8XvOAiqfs6bLRkmGu6biaNZgI/5O3kXv8mvsvTJ/gqkDPgO47t+xQVs1EaxfjyVG8Q8v/wdnEvXnXso3/rPmRX/aaUp7la5fNnxdhoW9dJirtW74391qJdXc2BtcEF5RCX0PlssleFtMOP3M2M86oFF+8saQRlZw+cfpazHWvll7zIucScMSqbCE/pbpu/HPO4tlz9o5n001wBpjOnTsr/SU0YcXFxVMjp1lMcxELJaNGj/r7weVOtsGaDgqCifC/PX36VM+azhlkCwC4vC5m64ZoJZ4RT7BwiOX43dm7dUZr25jrZ9dVV69av1FZlSuE9QkLDg0uLi52cHCQCZDMCwVJ27Mch1rQzMnKbQjSHqhcVnYsys3G/EVGbm7eZ1NT00YcU34l5srpHX+O9p788O7dmQnzxnpNo5qQWF5Mn3kORAYeADAlKGLwqIF8Pt/BwUEFL6Ipy8rKItnoERhEAoNoZGTYykDJMy9Q48BE+J9oNJq8Wg4AQKQIKkeUvtX9s6R3Iz0sHAz0JnoQZ8XeV27lCgwGg8FgAACIDKzbOOvihPKUg7lmQUYWHZrms4bNHCqVlB1bg6XoGY6IBFici4uLSCT61UoQGXr37we9HPtbG9oOZoxedHfannOOBHr9DwpDQ0NDQy06vEknFAqKL1Te4KUWk17qIUK5AUlfX19f00FBAMDHJ37Az8+vvVPb7HVvU1YnBI7sgFH25g7+waHnMvh5laJjSXltWjgrt/LvwACTtgY+8x2qP9UkbcsS5IsAAPn5+eXl5SpvGlI9VCYtO7EOQ6YqsuBPy3O53KKiorpfERnKT6lO/yP/xco0B5Lr9bTLhZX5f787G9w18NssCP0qGSK/mHZ1yq0FXhbur+4ndqMED7Prf//6PU3HBX0Bb/EfOX7gWG1trRwjn3xnvrL2IK0zbMTI6qqqbVev+LTvOZJYJHoXT/EIUGL930Vk4N3GW/PeVr4/nLv+/pLKqspKUcX4aeNmzJ6u6qYh1UFl0rJjazEkstHIhQ3JgssXr7h1+TYeR/AJ8I6avYH3tpKfUq0Y/7QJM2lLd6UfJv55+bB3oHfkkgVqiL9pe1uasu3FPjOayaHuWzl6bABA9Potmg4K+n9gIvwJMpkMAFjWbt7C2NVuLGeOnrESK4+YOi1i6jQAgLQgi3dwGQBADbkQAMDyYn4oTUZjsOs67UJQZMquYVNnTGmGx5w2DYosiCVRDEc2qC9YWVl5/e8bu3qdAAAsvDIt0TrJq2ML294mBNq/nwbjJ44bP3Gc6mJu2nJycrKysvz8/AAZc+jNqeeFr6a3HB9i1U7TcUH/CQ6NNoizoUO4c9j6+O2Iah67JJjbsyKiKi7sESXHq6L+b+GoWBL1y5ZXiBRFZHD3GZ2ESiW8QysamAURKcJLrvr4x2d5zZfbmEQj2Q8xNfE3/DoLQr/j9NnTIf06ztw337mV6+BT44g44umwfTALajl49zfU8BYDXhS9uZB2ZbBrX1XUTzC3Z0Ws5R1cDgCgeKq8X+jv77+JuXlNXGS5sCzMv9+7nZ+chprTLHXyYdhmKCkpac+WjUwmc6wl0czcrN68oFwu379n/7PHzwM7tJ8waQIiQys+CurGPy28TDuEhSy8M5WAI9i6W8OVn8q1acdmyxkuODIeNcW5fraZOWGipiOCfg4mwobCYrDLAuZE3JrX0sTTwcBOFU18yYUHlqEyCdU3RBVN1MFisf/cuvLx40c6nW5mZsZ7W5ly5BOntYFVNzYWDxeUajUejzd2QN9VLU3K8iWT7pTGp3+q1xfcEb0j8XJyT+cB5w+dqnwrCmZ1U8z/2YaZKFa+bA/a+unTJ7lcbmenkju52aqWCIRApFctwZHxqEBu62al6YigBoFDo7+Ao8ee3Xry2qdbxXKJipogmNuzJkdVxhwUJT9VURNfc3Z2NjMzAwAonhKrKal9u+PLglJIa7179y7AlN7SVL+LHRtLIAlFtV//V7kEuXvlwSCPkbZGDgM9R7zKedFqqZPnDDuzIKOv139aW1vDLKhEKEBvZz8YdW1a96m9iw9mfdr8nl3InDBugqbjghrkF3qEjx49OnnyJBaLHTdunL+/f73/WlBQsHv37rpfw8PDW7dW5V6aGhJi1e7x5+cH3hyf1SpCRU0QzO1ZEet4B5aiMqnMwbempsbExERFbX2tbkHph8Of2K0NrLqx5aisuLjYzMwMrqPRKq6ODvHZRS9MqdVimYxA1tPTAwAgUqQiXagY/3Q2cvsn5WIf94HXP14KGRGEp8L/faqVWZ69PfGAVC5bF7TUjeUExm6rqamhUqmajgtqqIb2CBMSEsLCwvz8/Hx8fLp16/b27dt6BYqLiw8fPmzwPyRSkz0qfa7flGcFL+MLElXXBMHcjjV53cGVkcEezkM7th/ar7fajlJieTF95juIeOK/I294O/uO7T2hlWfrwsJC9bQO/RQqEWNj9nLYlsueClc/r+geNvh/z/99zH/ApVlQWi52jL603rGb1ZGs7a36ek2cDDslKiSQCHe9PDzv/spQq8AD3ba6sZwU/w6zoG5p6OkTQ4YMcXBwiIqKAgDMmzevoqLi6NGjXxd49erVwIEDs7Ozf1qVTpw+8WPvuKkrH2880n2HIcXgV/+2gRvboyjqaWtxtbc7HotZ8ixvYvTBwMDARgXbSH07hPcxG+LIdn2SHcuzyYveoeEnn+DpEwAAVFLLO7TifZU8+lLm4qB1CIpM/GPIhTVXTXyNDN3peEr9np/ihHrYof9VP77Z5iyce+mfS3QaferKGXclT1uZeE9vOZ5JYqgzQi3ULE6fSEhIGD9+vOLn4ODgxYsXf1umqqpq/vz5NBqtZ8+eTXJctI6HsWsP+85bEvasD1mGASpZWoKiKAaLU+yDRgRyqVTdx0dgyCgRTwIAEHEkEb/2p+UhVUMltaUH1opwPqUiKlacCwDAYDAEPYLzGIsmPACjbR48eHD11U2bpZ61vJpVC1c+fvzY2RAuu9V5DU2ExcXFRkZfDlsxNjb+enMmBT09vcGDB9vY2OTm5nbs2HHXrl1jxoz5blVZWVlPnjxJS0tT/IrBYNatW2dhYdGY8DVngG3YgvxV22J2M3nUgIAAK6uGLg8TCoUN3K1t8Kgxw/48zaESKqoE7nihQCD4jXh/2awFM6eNn+7CcU8teh8Vtu3jX3mmnfSxRI2trmr4dWsCamtrD+w5kJ6W0W9Q385dOiMSpDpdwLuRKBQNoNkz2nSmnIw/FPVwcaWovM+AMKlU+l/fk2CPsHF+cLO9T0/BmhMABpCNqRSUbE40UfMbU2speoRyuVzTgdRHpVKx2J98cDU0EZLJZInky1LJ2trab0fAXVxc9u7dq/jZw8Nj1apV/5UITUxMXF1dhwwZoviVQCBYW1vr4lfaVny3xRuXGrcyrdm29u8TF1q1atWQv5LL5Q2cP1i8YtWIseMrKipcjJn8Q8swFCqlZejvhfwLgoKDnr5+kpGR4erqSkCI2ZeK0g8U2w80ZdrrqS2GrzX8ujUB82ctwOWpwkltAAAgAElEQVRQg017bloULX2HtZI4kjB5+iYi17Et8BQ8AODa3avv37+n0+nW1tY/qAeDwcBE2AjfvdkQFL2XG3eb+IT/pBiLxcoKxL269mg+9+RPae3Q6E+zIGh4IrSwsMjLy2vbti0AIC8v78cdOE9Pz8LCQhRFv/utSk9Pz9XVddCgQQ1sWmudOXPGaYo3yYBcaU0/cPzgMT+/hvwVFottyP8YBWtra8UnnfGUDdz9izEYDLVVh8ZH/Iv09fXrhrhdx1jzkiozzhSwW+lbdeeo/1nDX7puOk1WK38aG7+v7xkMwPRyDk8uehnslkhkGxkMnv31Kbuenp4/rQr7P6qMtwn69qIll6bsfHmITqTt6b+Z2pN89epVS0vLLl26NJ9Rip/S6ZutoUGHh4efPn0aRVEEQc6cORMeHq749/Pnz5eWlgIAuFyu4l9QFD116lTLli2b/C1ixjEVFQgAAKICgaWZaod28RxL4ykbKq8dq0lUyYFNDcHyZvoucpTVIq83ZVRmCQEAVVVVlZWVmopH15WUlNSNsgAAZCJ56cuKD0c+Ja75aM92upX6D09QGpt7y5eQQWSz6mVBSG24Nbw1T6Oj4rePdB+0o9M6W31rDoczYcKErl27NvmPuOajoatG+Xx+aGgoiURCEASDwdy/f//LQXdE4t27d4ODg+fPn3/p0iV7e/vPnz8jCHLx4sX/+sbaBFaNKhQUFPQb1r+4tBg1wO44tCvcPawhf/U7qx9lpfncfYuYPcdSW3dsXA1Kwf9QnXWx8Oy7o3Fv7mGx2O59u63bHKXqRpvSqlGJRNK7Wx8hV1Qm5G3csqGtdSAvqbIyW8i002N5M43cGXfu34wYMxGHYikE6d2N88xHLWhcFoRzhI1QUFAwK3J2Ka909tQ5iDPuTMpf3e06jfIYTMHDE61/RGuHRhuioUOjhoaGr1+/TkxMxGKxLVu2rHtrZWVlsdlsAMCWLVsmTpyYn5/PZrNdXV2VeJi71jI3N3/xMAEA8LmqYMbdRe5mrqpeP4ZnW7Amr+ftXwywGGpL9Y2R1mPoRpeMNYhrf39P+CkAwILrkwtnFyp2qIEaIuZyjIXMbkznKUKxYO6MiL/XX2W30nceaYkjfRmh2bBy6T/9XFhU4pG3eVeFepNhz0ON+g3rX9sWS/GjTVwyacCiofv7RZvSOJoOClKtX0hXOBxOMUf4NUtLS8UPGAzG2dnZ2Vn1B8xqH0uG+ezWk1c93nyw+1YGUbW9FoKJFWvqBt6+xQAF6pwv/CYORI/2ZdUMESGJBHBXtgaR1cjL3lVl3ykgATIAgIgn4ulYlzH1lxxLJRIyHgsA0CPgxWL47IpaFZYW2Xt6AwBYvqahBH+YBZuDpt9vU48Qq3YpvI+rn2zZEroKi1HtdDGBY2U8bRN33+InL99ceZFs5+QyddZsNY9ImJiYuLZyXh27AIvBGbGMyi/IKwcLmQ6aWVCqhSorK3dE7+SWcCdMHe/t7S0Tyfkp1XXjn4NGDxwyI7z0eWFOWdas+TPr/S1SWzOppd2Ia4keFux3ldLbo8Zo4hU0R1WS6hPJf8roaMmjfIoptfZNZdu19b/6Q01SQ+cIlajJzBHWI0flc+8t9zXxHO0x5AfFlDXX9ebhvcnDBy8KcHzFFRSZuR89++fv1/mrPnz4gKJoixYtylOrMy8UGrjSbXub1I3vKZduzRF269DDl+RvQjM/nLhz1/RDZD79y/yfB0NxfaRSaVJSkrm5eb0hZaS2hndgCcHMThw4MO/zZy8vr995sgjOETaQDJHfzL53PPmPAHO/gTZhO7buKC0rnRUx069hS8Eh0EzmCKGfwmFwqwIjJ92c62Lk2Maspaqbe/I2ZUgL89Zm+q3N9HvfVOHGpz/g5uam+MHAle6zwCH3Wsmb6EzHwebNvGsoFcoKsguW9e8HAGhXEpqLpo9aObze9wMCgfDt7ktIrZB3YCnBzM5g4AyAwZiYmqov6GbsVfHb3a+OGJL1ozussdO3BgCsWb5ah751Qb9PJ5/50FoGZP1VgZEbnu0sEpSoui3flq1uFQi5NeJbWSVWxoaqbu6n8BScw0Az+/6m6X/kZ14olIsRAEBZWZnatgtXv3qvTlbz5fmHl+vSyXhKUv7L0uriV6XP2/dt25Becr0sqMrAm7ujx491799j2erlGaXZi+LWbE88MMZjyLaOaxVZEGqGYCJUshYsl2Fu/Zc92qC6MwsV2rVrN2zu0kWptYk0+6gAG+HT6yptroEMXOk+kQ4AgLjVL1u7t+kbEu7t6vP+/XtNx6Vk5eXl/r4BfUPCvVy8kxLffsl/Uem8pEqWN9NvpcvF++efg/uHMrau2LzU0dHxpxUitULefpgF1eHKP1eijm8oayf+M+1y94m9PNktjvfcHWLVTtNxQZoE5wiVDwXoqseb6UTa/DbTvv2vqpjrknELuHsXMjoP1WvXU7k1N9qiKUuJnxhdXHpm8dJjeGcu37j0mxVq1Rzh2lVRkhe4ri5hOWWZxxL2Hlh2nOXFMHChYQmN+Wb5JQta2BsMmK7cLAjnCL81fd7MR8REA3djRIrwduekvkr5toxW3Wy6QqfnCGGPUPkwALPIf9Y77ocbWffU0yLe2Nx42qaqu+cET6+pp8WfoyI0Ig0AQCPRq8uFmo5GaaTVsuJ4fn5CiR5B8eoYOGOM8wgLIw/G72RBoo2r0rMg9K34ghfvqOncuHxRsbD0zufQwBBNRwRpBbhYRiUoePLaoCUj909cdWcZGU/evGZj2zaqXYetyIXcvYsAALR2vVTaVkNMmjaxd9c+SbzEdwVJU0NmZ14oVN2CUlUoLi7evX2PXCafNnuqpaWlTCjnp/77/MPEKRNGzx2RXPYqpThpy77Njaj/xYsXZ48esrG1DccUMpw89ftFKP0lQF9LK8vY+/pYtbh685R1qY7vzsf81c2n75IF3zlODmqG4NCoqtTW1tp62llOdgUAlBzOSnudqhg0UOmoi4xfwtu7kBY6gNZe87mwurr67du3Tk5Ohgyj3KslFR8FDoPM9J1oja5NbaNVCIK08vTrbzcMh8GffXf0TORFUZ5U8fwDy5OhOItKIBAkJSU5OTkpdlb6JWlpaaPCui3wYieVVuVi6GcfqmrRLxwaBQCU1vAOJ51+U/JutMfgnvadG/KYLxwabQSdHhqFPUJVyc/PZ1gYUDh6AACKmV5eXp4att3BG3JYU9Zz9y7C4HB6/t1V3dyP0en09u3bK352GGhWnibIOF9g4Eq3DdP2rmFOWq4BxijEoQsA4OmnOKEp33+CX70DN2g0Wt2r+1WPHj4cYEP3tzD0tzDsfSNdCRFDX3mR+GLu4nliiWTF0uX5xtyb2fcHuvSe5zeVjNe9s94g9dDqzyOdZmNjg/BkZS+Ly14Wl+eX29nZqaddPMuMPSO6+t55wZNrCIJoz6mhBi40nwUOqBx9syWzIl0AAJDJZCKRJvdmq66u/vpXSZWs8HHZu705RWeFBfzPn8s/FVUWfBJkenVsodxjp9xdnO7m8CrF0od5ZWb/26QQUgqZTDZw5CBJDyJuEHPE1FE8Pu9krz0j3QfBLAj9AOwRqgoej4+7+SB611YZIuMtE17OvD7Ita96msYZso2nbb4eOWbx0AlMPaqBudWlG7e14ehjPBnnONi84qMg43zB06IHB2/uJuHJLdv6Hjl1WM0n2mRlZfXvNYCCpSIE+eWYywQehZdUWf2pxsCVbh7CcnexOdP31MrFq+Ry+YET+2m0Rg7nfhciEji8uhzevfO0hBRbO/vDZ3YqsXKolFsK9HGKkRgTJ7Oehh0NyPqaDgrSdnCOUB14Iv7U2wtmtJwYaNlWbdMP7b1b7PHjsKjEPUn5bmPnjxo9Rg2NNpC8FnGzb7Ev/AwJT4p+snr+tlk/HWZU7nUbOXh0AKaTh5nPw8z76WUfls1YwfJi6DvRMDjV5mOkppq7bzHZyZvZe4JKG1JobnOE77lpB14fPz/7pHGoOZaIrb1f+e7F20bMWsE5wkaAc4TQT7AohlFBSxY8WMWiGloQ1bRvllgqo5PwAABDAqaqsko9jTYQhggIZDwRTwQA0HHMSr76whNXSMuSq0o+8pg+BgAAfYo+1RnnOMRcDU0jwiru/sVkZ19m2Hg1NNes5FZ+Pp78R2pZxkj3gSsezd+zf49YIp55a4aOfi5DagYToZo4Gdov8p+1/NHGze1WqOfL5rS580duXu9uRE3ILfpnpoEaWmw4LBbbf3D/ZVdns2kmeRU5zMRZZXZVRu4MJTaRlZV1eN8RAyP9ydMmM5lMRf7jva2sKRUbutJnzpmxZO0ib7NWrwtfnL9yTont/hekRsA7sIzs5AOzoFIIhUI9PT0AALeGd/Ld+Sf5zwe59l0WMJeAIwAA1qxYrekAIV0Ch0bV6q/UmGsZd/Z3j9YjUNXQXE5OTl5enq+jjfDYKlpQX1pQHzU02nAfPnwoLy9v06ZNTZ4443yBnhnFYaAZnvr9cbxfGq0qLy8P8gsZ6T6prIb7rPThnnFHFfmP5cU0cPky/llSUpKamurp6WloqPKdWpEaAW//EpKjp3pGROs0yaHRT58+devbvRaIKVjy2E0RjyoSetl3GeE+UInvKTg02ghwaBRqqEGufXP4eWueRG8IWabqYwsBALa2tra2tgAA6rRN3L0LAQBalQvrDq9g2OG95zrkXit+E53pMNDMwPV3P4MSHr7wMfYLsA0GANxJv8YOpnPcbTHY/zf/x+FwOBx1nLmK1Ah4+xeTHL3UnAWbqiWrl5J6GLDdjMpeFV85dvnK0UtGFO0a8IB0Dnx8Qt0iPEYhANn18rA6G8UZsI2nbRI8ihE8uqLOdhsOR8Lah5s5DbPIulSUfi5fJpIDAH76cEW9ArV8SUEs7+2OLEkc6W3B6wpReTYvA1AQE0/jellQDRAEqa2thVlQuWSIPKM4C88gAgAITJIlyRRmQej3wUSobjgMbnX7he+4Hy6mXVVruwZs4+mbBY9iBA9j1NnuL2E66PlGOhDphMQNab1D+7bzDnKzb/Ew7uG3JROeJ7g7erTzDgrr2ruqWFD4qCx5d/bb7VnColrLzuywrcFRu1dtfr3sQvHxMxdOq/+FxFz629PBpl0Lp7GhrYkwCyoDgqJxeU9HX5tm29W58HgG93Je2bm8yJkLNB0X1BTAOUJ1U0w/8GrKpt6J7GfYTZRZ7eXl1bKlyg/yVZBXcLl7FtICe9OC1fRQY+Oc3H3m/tlHkwJmVojKVz6Yd/3ErdraWjKZXFeg74ReC9tFGdPYpxOPcIxMRo0YxfJiMu2p6u/5fZeXo+1fXRz0CLjIuPQpu08GBgZqJIwmM0f4qvjt/tfHCTjCRK+RviaeJSUlycnJPj4+LBZLFc3BOcJGgHOE0C9jUY0G0npOGDeRE2gh3lcdNWf1yOEj1dAuTt/YePom7p6FqFRM7zRYDS02jowqZtFYAAAaiV4rqhXk10qlUhnh3wKimloGmQkAMKQaGQVRHAaYaSrUb6EoCuQyKgEHAGAzqFVV2vXsivYrKio6cuIok8EcP2Zcdk3e4aRTVeLqMZ5Dg60CMAADAOBwOJ07d9Z0mFDTAROhxtyJuWU7xIXpaiT1E+87ul89iRB8yYWbuXsiPxWXXvhQwDQwmDRlmnI3T/l9/fr12xXdkSfhZpZ9nDo/wmGgWb0v6bOrZi6JnuFq7P6qOOHe8DsaDPVbqLCqiy17wq0UGyPmexFuQ4cOmo5IlwgEgvadAknBTLlAFn1iW+DyruM8h4VaB2LhGVWQysBEqDHWltZPUhOBq5Egt8rO1ESdTeP0WfgRS4b5+01vZVMqkg26cf1G7CN1BvBThoaGz149jY+Pt7KycnJy+rbA2PFjAoPb5+Tk7PHfplVZHBFUcPcuWr5gzidzLy6Xu7t9eyKRqOmgdMmbN28IDlTjduYAgMzkpD3BGw2YcI80SLVgItSYeTPnPhv+LGX9OzKTSpzsUimuYpKU+UT5j71Ozw51MOthzwYAXL35USKRaNvnNZVK7dSp0w8KODg4ODg4qC2ehpBXl/P2LqL4BDG6DvfQdDC6KLfy82XerZKPhaxaC5lQSpTg9RlMTQcFNX1w1ajGUCiU65eu5abkpMWndPHqOPf+8mqJ+k6KcHJyelEqrBRLM/gCmVSqbVlQFyGCCt6+xYosqOlYdE+RoCQ6Ye/se0s97N23LY8u2/tJ+lfFuePn1LwbO9Q8wR6hVpjgNUIql0bGrt7WcS0FT/75H/w2Gxub+VGbIjasY9D0tvXxr7p7jtF5qBrabaoUI6IU70CYBRtCLBbzeDwzMzMMBlMsLD3z/sKT/Od9HHuc7X1Aj0AFnmDs8DGajhFqRmAi1BaTfcdsTdi35GHUppAVRJw6+mfhAweFDxwEvgzoLQQAwFzYOP8bEQ1mdB2m6Vh0QGxc7OiIMSQWhSIjDd446hnvVZhD1zNhB2hEPU2HBjVTcGhUW2AAZq7fFAMSc+XjzTJErs6mcXQD1rRNolexVXfVsf10EwOz4K+as3ie+SwXs6lOIlc05d67P/scjvAZDbMgpEEwEWoRLAa7JGAuBgPWxW9DUESdTf+bC+/AXNggT588WbEo8sLpE9w9Cym+ITALNlCRoIQr4GIJWAAAnoT3MXKHKRDSOJgItQsei1sduEgordn0fBei3k1/vuTC13EwF/7Uw7i4heNH2Kfev7hpxYlsAaMLHFL+OcVymIk35/Qe2y9vx4fSs7nI85qxo8ZoOi4IgnOE2oeAxa8JXBQZu2r3q0MdSAF8Pj8gIIBAIPz8L38bjm7AmrqBu3chBovV5n1nNO7q33/N9GAHWBi2MTeISExbpOl4tFyhoPiPlL/j8p6GOXQ91+cQnUhbM3J5Xl6eu7s7iUTSdHQQBBOhViLjSRtClrcbFbI5ez2Dw6QsIjyLi1fPRwaOYciesYW7ZyEql8Phvv/i7OAY++x2W3OD+5/4ri3cNR2O1pkdOTvm6hUGg7Fp++bX4ENi0es+jj3+7HO4bhTU2NjY2NhYs0FCUB04NKqlKDgy702xwwwv9mAbkYn00SP17fyCpekbT9+UGXcjwNnW18G6c3v/iooKtbWu/eRV5T2FH6i2Lv3v5iTS7Tds36npiLTL3bt3r725bb3UgzyYNWLiSDt963N9Do/1HArnAiGtBXuEWgqLxWIBFpHKsUScsFxIo6t1FzEsTX9XhmCOp7G/ueHfH4u3b964ev1GdQagteSVZdw9kXT/7luXDNB0LFrqeWoi1pwIMIDC1qMD2lC3/pqOCIJ+AiZC7bVp9YbI5QsBAcOwM4jHvmmDtlXnvsPllZWWlhQAgCWNlFFSrLZ2tZm8gsfdu1AvoAc9NFzTsWijd9wPf6T8/YGeWvm0FIvFyPLFvbuHaTooCPo5mAi119DBQwcNGCSVSgEeLH20Pupp9JKAuXisms6Wi5g9b9acaSFmjJsf8/cv762eRrWZvILH3ROp17YrzIL1ICj6vDDx1Lu/auXioa79ooKXlHYrvXbtmpWVVdeuXTUdHQT9HDyYV90ad+anVC5d8zRaikhXBy4iqWXfGQBAdnZ2cnJya3dX/N/byW5+zLBx6mn3uzR7VqriQGO9tl11azGtig7mTXjxYsL0CQKhcPTwUf5DAs+kXCTjScNbDKg7L1DXwYN5GwEezAupHAFHWBUYuSF+Z2Ts6g3By6gEddxtdnZ2dnZ2AABk2ibu3oUAAM3mQk2RV3C5eyL12nbTrSyoOqMmjTIYZ2lgQN6/+1CmQf78/jO9OfCwDUiHwVWjOgOHwS0JmGPNsJhzf1mVuFqdTWNpTONpm2o/vKi8ekyd7WoDRV+QFtgHZkGF8tqK8poKkhEFg8XoO7D6sbvBLAjpOpgIdQkWg5nrNyXYKmD63UXcGp5am1bkwtTEyqtH1dmuppw8etjdzsrdxmLryN60wN604L6ajkjzCgXFu14eHnl1qrmzZcGfmaUPPkuSBKGhoZqOC4J+Fxwa1T3D3MJJOOKMu4uXeM56F//Wzs7O399fDe1iaUzjqRu5+xZhcHhGj9FqaFFTBALB9nVrrvZ0xWAwg6++G+MaoNaHV7QDn8/X19fHYrEAgI/8zItpV18WJ/V26Ha29wF6OO2ff/4pKCoctHGgoaGhpiOFoN8FE6FOCncOq+JWBXYI5LS3kuWKRnUbtmrpKjW0q+gXcvcuBChg9GyyubCqqoqtRyLgsAAACyMmn89ns9maDkp9KioqQrqGViNCRChfvXPNCzS5TFQe7hw2z28qGf9le6O+fWEXGWo6YCLUVTXvKk072XBCLAGKnt50Vj2JEACA1WMYT9/C27+kWlB9hYsiCDpy9Ggmk6me1tWDQ8FjpeI1z7IJBGIVxcDJyUnTEanVzj07pT5Yi0CXmkLBsuXLz144285Crc+wQpCawUSoq0xNTOTFUgBALU+EIap1rhdLpRlNXhfQwrGXnREOh+t69NDT12+VvkZfU+Tlpdy9C88f3P1MRJbL5Tu7dlUMDzYTpTW8Z3mJGBoOAIDXI5iQjAMt1THwDkEaBBOhrurbp++la5cfRT2mkMmu47yjE/bO8ZuMw6gpG33mlrEN9Md6WgIAkio/p6enu7q6qqdplZKXl3L3LKSF9KMF9u6m6WDUTDER+LzgZcuebT4sPIJkiQWZlQe37td0XBCkcjAR6iosFnvm6GnFzzVS0dqn0Qtj16xqH6merY05HM6nypqKWikWg0kv4pqZmamhUVWT80u5exfSQvvT2jfxjcGEQuHps2dkUunIESP16LQn+c8vpF3hiyrCHLvOajWJRtRb3HFOUlKSo6Mjh8PRdLAQpHIwETYFVAJlXfDSw0mnI27NWx+81JppqeoWKRTK5r0HJixaABBkQdc2IPY86DtJ1Y2qVPPJggCAoM7BQnsZwGM2749uvSrElMEZ6NInyNIfi/kyCEyj0dq3b6/ZICFIbWAibCKwGGyEz2gLhtmse0uXBsxpbeqj6ha7de/RrXsPAAAiEvD2L6mMOcTU2Vz4vywYTmvfS9OxqFxhYWE5WmnZzRUAkPuperzp0B7t4I6gULPWjFYBNAc97TuvDVq04dmOi2lX1dYolkJjTVkvzk6piDmktkaVSMYv4e6NbA5ZUCqX3smJW/0mmlfEkwokcpEMLZG1cWml6bggSMNgj7Cp8TB2O9ht69KH6z8Wp1feK01LTxs/cnyfMNUeH4Gl0FhT1vH2L62IOaSvI/3Cx48fJyY8b9/S2+rJWVrogKaXBaurq/H4L29wnoh/NeP2lYybdvrWQz3C+x7qMm/xfASRR6/ZYmRkpNk4IUjj4OkT6qaeje1rZbWtwtpKOAjD1Yh74dNfe8+1adNG1Y0itULeviVEWzf9fhFKr1y51+3sqZMnN60Os6SfTSlYNn9Wr1nLlVWzNsjJyenWt7sEJ6MA0vbjO+OrXr0uSQ62DBjg0ttG9fPHTQA8faIRdPr0CTg02jSR8WRRvsC0s42eBZ3a1iD2UZwaGsWS9VhT10tyUysuHwRq/4L1S86fOh7V1qqfi+nKQKd/Xn7QdDhKtmjlYkpvI6t5btgQvXlrF/qZ+V7oe2x+m2kwC0LQd8FE2GR5tvDkPioQlQi5Twsr2EIZIldDo1iyHmvKOkluavnFPdqcC20tLRILKwAAz4ur7ZyawhOQCihA35S8e1vwnkAnAgAIdKInw7Wnfee6rdEgCPoWTIRN1omDxzvR2zPjcJvnr8dYE6fdicyvLlRDu4pcWJ3zsX9bTx9Hm87t2hYUFKih3YaTcQumcCSxQlLYzYzPJi1mL4jUdERKUC0RXM28Pfb6zO2JBwaNGVp0PIv7dx7v/OfImfM1HRoEaTu4WKbJYjKZu6J31v16O/vB1NuRo9wHD3BR+XNyWLLeX1K2J0W+q6fL03z+kjkzT/71t6obbSAZt4C7d6Fl/4kXN3bWdCyNd+DIwa27tpJJ5N1bd9Od9P/JuJVY9Ka1qc+MlhNamngBAKZ1mZCSkuLs7GxiYqLpYCFI28FE2Fx0tevgZOSw9kl0Ci9tnt9UVW9AU1xU3NqIBgBwMaIdeZ+v0rYaTlaaz923iNljNNVPh7Ngdnb2xn2bLee6ymqkfUb3G7x3TC/HLpFtZ+gRqHVlOBwOh8MRiUQajBOCdAUcGm1GbJlWB7pFG1L0x92Y9bY0BQBQUVEhFotV0dawseM3JRWfeJc/807KYF8nbZgvlHELuPsW63oWFMslVxKvYczxWAKWyCQx6IwdIVG9Hbt9nQUhCPolsEfYvBBxxBktJwaY+615Gp13Mo2bWyoTSVcsXD5p3ETlNuTbsuXFu3GPHj3a0cLNIuFC+V+7DAbNBJo7ykdWms/dt5jZY5TuZsGP/Mzb2bH3cx+ZUIwlGTUlD/JQAepk4UClwhQIQb8FPkeoblryiNKDp7Gjloyzn+SBypCc9cmfUnNVd9gQKqnlHVqON7b4nVz4O9fty4hozzHU1p0aV4OapaenDxkzlF/BDwoI3LR9y71PD29lP8Bhcd3sOnSxCWFRjcrKyv449weTyRwyeAiRSPyvekQiEZFIbDInZKmNlrxJdYtOP0cIe4TNFAHB0+g0AAAGh61FJDUSEY2sqllDDJHMmrSWd2iFRvqFOpcFAQBjp4zDhtGtrcwfnU/ouTx85PARSwPmurH+PR/YyMhoxvQZGowQgpoSOEfYTPn7+xsKGSWns/N2fXANdJ9we/bT/Beqaw5DJLMmrZFxC8rP71TPfOGu6M1eDjatXBwuLxjD7DlWV7KgFJE9zX+RWZBFNacBAKhW9H4m3ea0nvx1FoQgSLlgj7CZwuPxT+4/fvnypYGBgaOj49vSlO0v9l/LvD2z1SRTmkqOoMMQyayJq3kHl5Vf2G0wcIZK+4Xp6ekxxw9e6eFcLZENv/lhwIkOqmurcSQSSVJSkpmZmYWFBQAAQdG90AcAABVvSURBVNEUXmrsp6cPPj22oJt26t0l/lA80ZFa86xi2PWhmg4Wgpo4mAibLywW6+fnp/jZi93iSI+dMek3ptye38exx/AW4UTcf848NRqGRGFFRKkhFxYWFjrrk/FYjAGZQCGTxGIxmUxWUVuNUF1d3TYkAOFgagoFk8ZPNAu2uZ3zgIQjhli139d1sxnNBHQBjx8/Ts9I77Kmi6Ul3BcNglQLJkLoCzwWN8AlLMjKf/fLw+NuzJrqNfb2qRtPnj8J69pr1vRZGCUlLQyJwpq8nnd4Rdmf259Q7MrKyvr262doaKiUyuv42pjNySk6SAA8CWrn5qFVWRAAcPnyZZkb1qSbNSJFotdt3Rq8a3PoSiuGxddlAgMDAwMDNRUhBDUrMBFC/w+bylobtDih8NWEBRHC2hrjTma7rx2k0xjjx45TVhMYIok1cc3wQF89udiSSe64eUPci1dMJlNZ9ctKPguPr7517tStHF5Lff2wMG05cR4F6MeyzId58affnxELxAAARCpn67HGeQ7TdGgQ1KzBRAh9RxuzlrRiEi3ciMgkMQONb8bdUmIiBACgeMLHMsGVXq4AgCpZwePHj3v1Us5xgNKiXN7+Jcx+EVSf4NHBSqmyMdLS0iJmTynj8ebPnDdq5OgUXurDvPhHn5+RcKRgq4BD8/ZMGRmRvytNUiHet32vxqKEIAgAABMh9F+6hHY+d/ciM4Bdei8feMn3vT421K2/AVlfKZVjsVgUR+DWiA0pxHfF5SPMzZVSrbQwh3dgqX6/yRSfIKVU2GgDRgwkhxvpsy0Xblt6piLG2sE6yDJgc+iquoOQHt6JKy0t1dfX/8FTgBAEqQd8fAL6vuWLlk0IHMV5obdy1OI7G/6RIrJRV6ftfX2svLZCKfXvPnJ8yrOS3tc/Brg5WKfF/v4zFdLCHN6BJcw+EzWbBavE1bezY4urSvSsGDgynulqOITV+2iPnaM9Btc7DpDNZsMsCEHaAO4so266u2lFeW3lX6kx/2TeCrVqP8ZzKB2rd+nyJQRB+vfr/zvbSaASMe/ISry+scHQuT9YR/rj6yYtyOYdWMrsF0H1DWl0JL/kyLGjl65fauPbZvGCRUQisVBQHJ+fGF/wIq0sw8XIMXbDTZGpDMciCm5zX8Q+53BU8kTKT8GdZRpHd9+kGqTTO8vARKhuuv4e+5IOM26lbHpBctDDYDHUXPzLp4m/s0ObIhfiaEzDEQvBf9Tzg+um/ix48dLfkTsXs8Isyl+U2tAsTfrZSOTS1qY+ARatW5v4EHCEmpqaA4cPlnJLJoyZ4ODgoJ6ovgUTYePo+ptUI3Q6EcI5QujXGJCZET6j2+r59iSGWfe2AQAUHEnPyMhwdnZudJ0YIok1YTXvyCr+mU2GIyIBtkEf3I8fP87Kygr1cCLG7FZbFkQBmlWee+jKUXqgMcWURuxKztybfmj7QVum1dfFqFTq3Flz1BAPBEG/D84RQo3hYG6HqUDlYjkikfPyueve7ohJvyGU1jS6QgyRxJq4Wi6o4p/ZDBD5T8tvWRe1efrY/FNbevboUREwQOlZ8M69ux6tPd1aul+5egUAUCLk3si6t+ZpdL+/R618vMnUzbz8UbGoRFh6N79Pp7B6WRCCIN2CW7VqlZqbfPLkiUQi6dixo5rb1RISiYREImk6it9FJBLN2KZXNl6ofla2Ydn6sJCej/Of73h58FNVPpNE5+ixAQDl5eXPnz8nkUgNHGXC4HBU7yDh89vijCSKu3+9+cJ6123WxLEnOjr4mjCpRFy6noV/u3ZKfHVisbhTWCfjqXaU1sxjKw+9Yqdezr4JANrS1DvCe9Qoj8H92vemysm5d9O7eXVas3w1Hq+lIysymQyHw6nuXJGmqmm8SdUMQRAEQQgEgqYDaQwtfQND2m/IoCFDBg2p+7WVqXeluOp2TuzWF/sQFPHFu2+bs5nuZiD4WHFo24FuXbs1pE4MgciauLrsyCr+H9GGw+Z/O18oK82vSXosevuYIavJLhc6GdFSKyQhlo3pkCUnJ79+/TogIMDJ6d/9rAsFxe+5qY/fx0uZCIFOBADQrZhjrAd38++CAf8vMUdMiIiYENGIdiEI0jZwsYy6xcTE9OjRo2mvm0/hpU2eO7WCU2PozRbza9HL1QkPnjX8z1GZtOzY2k1XH977mM8yZu05dsqUSki79qeNpExeXkLxbEfxDsoRYyeNHFpVUdE2MGjf0RO/2um5dPnS7DXzqN5M4YvyjVs2k+yo77gf3nE/YADGg+3mznJZNnwRxosEcBjwsvZdYrLW9vl+7NmzZxYWFnC30l914cKFgQMHajoKHfPp06eioqK2bdtqOpDGaOjbG0XREydO3Lt3j8PhzJkz57tvrYSEhCNHjshkslGjRoWGhio1zqZjwYIFbm5uX/dCmp4WLJcQ+3Z/F98AAEirxEW1BRNvzvEz821j6tvC2AWH+bIW5u3bt1wuNygoqN7XAgyekGjm9/nTySvdXTL4wok9Qk8Nav8mh+e9YQ/RxlUxZOoCwKOXST+IQSQS3bx5U19fPzQ09OuNUhEUya8ujNq73nS0PdmYWm3LWL1nzeRV09tZtJniO9ZEj60o1iU25OiJY3K5bPzm8TqaBQEABw4cCA4OHjdOmbsCNXkIggwZMgQmwl91586d+Pj4Jp4It27deuTIkbVr18bHxwcGBqamptZbJvv+/fvOnTuvW7eOSqX269fv+vXr7ZQ6bQPplrkz517vcSPv+QecDHP3r5sIC5NQ+HrPq6OFgmJfE08/U5+rey/fS4wlsEikJdjnD5/Vu50+ff7sy2FiMRgnI5oAED53n35q8eIIW7cGtl5bW9u6vZ/cAYtUyb2Ou63duT6Dn51Rnp3Bz86uyDUg66M0bM3narIxtbZAOKBVn7l+U+rVwGAw5sycrZxroVHqH/KBmiedvtMalAjlcvmOHTtOnDjRqVOngQMHPn78+MKFC6NGjfq6zO7du0ePHj1jxgwAQGFh4bZt22AibM4MDAxeP3slEAhoNJriXzyM3SZ4jSivrUwsepOQ/+rPK+e91wQCDCj8O+tYzMlh/YYakP/dd7t7jx6d1q2+9KlEKsOE9Oj37YP2OTk5IyaMLCopDuvea8fm7XV9vkpxVUF10c37t0SmMovejgCAO6vvUx5td+E4OhrYh1q3dzSw0yNQS9qW9B8W/unWO0dbxyVbl6jlkkAQpKUalAg/f/5cWFgYFPRl56rg4ODnz5/XS4TPnz9fvny54uegoKCDBw8qN1BIF9VlwToGZGYX25AutiFHSHvkEjmOhKvl18aVxt++9gRFUWumpQ3T0oppgRZLhQyCXrhVTVZ1fiVf8bdCaQ2CImK5RCKXDJ0wTBJCsLR1u3L6hmQ7MPEzL6guKhAUYQDGnG5KlGFrS4QAReW1cjpG70SfPfWOkeJwOE/vP1HTVYAgSLs1aLHM8+fPu3btWllZqfh1w4YNiYmJly5d+rqMiYnJuXPnFFODKSkp3t7eYrH4u0sYevXq9eDBg7qhMAwG4+zsrG0nxqnO8+fPvby8dHT/BSXi8nlFgmICnYivxjpaOQAAABmDGmBRJgZlYrn5pRhDHMvbBACQtv+142APsVhMpVLr/jztXJLLFF+AgvL3XMkHoamFKahCQBWCEX8pkFv4qRoRyCVySwNzIwMjDbxCLZCamqqvr29qaqrpQHQJiqJxcXFwlcOvKiwsrKysdHV11XQg9fXv33/atGk/LtOgHiGFQhGLxXW/1tbWfv2RpEAmk+vKKA4E/6+FfEuWLPHx8al7tgyDwXh4eOjueoRflZuba21traxzbpuFdUAmkxUVFf2/JVpR//vBAoAGPZrRHJWUlNDp9G/frdCPDR061NbWVtNR6BihUCgUCtlstqYDqa8h/ysblH7Mzc0lEklJSYli7+C8vDwLC4t6ZSwsLPLy8hQ/f7dAnYCAgICAgIa0C0EQBEGq1qCnr1gsVnBw8KlTpwAAZWVl169fDw8PBwDweLw//vhDUSY8PPzs2bNyuRwAcOrUKUUBCIIgCNJyDX2gPjExMSwszNPT8+PHj927dz9w4AAAID4+vl27dooaBAJB586dRSIRmUwWCoWxsbEsFku1sUMQBEHQb/uFnWUEAsGrV69MTEzqzhmQSCSlpaV1o6AIgrx+/Vomk7Vq1ar5zPlBEARBOk0DW6xBEARBkPbQwOkTzVZBQUFMTEx+fr6tre23Z6UKhcL4+PiHDx/m5OSw2ezmvNIvJSXl+vXrAoHA2tr6uwV4PF5MTExGRoa1tbWO7navdAiC3L9//9GjR3Q63dDQ8NsCmZmZcXFxCQkJYrH4B2vZmpvCwsKYmJjPnz/b2Nj81wnGqampN27cyM7OZrFYenp6ao5QOyUkJNy9exdF0f96OCc1NfXmzZtlZWU2NjY6sEgehdQiPj7e0NBw7NixAQEBQUFBEomkXoGNGzeGhISMGzeuZ8+e+vr6z54900icGnfs2DE2mx0REeHk5DR16tRvC6Snp7PZ7GHDhnXt2tXNza2yslL9QWqh8PBwLy+viRMnGhkZ/fPPP98WsLW1HTRo0NixY62trQcOHIggiPqD1DYJCQmKd2X79u0DAgLEYvG3ZVasWMHhcIYNGxYeHh4ZGan+ILXQihUrrK2tIyIiLCwstmzZ8m2BrVu3mpubz5gxo3Xr1v3791d/hL8KJkI16dy586ZNm1AUlUgkbm5uFy9e/EHhmTNnDhkyRF2haRGJRGJqanrnzh0URUtKSvT09DIzM+uVGT9+/PTp01EURRCkQ4cO27dv10CgWiY+Pp7NZldVVaEoeubMGQ8Pjx8U5nK5eDw+NTVVXdFprx49eqxbtw5FUalU6uHh8eeff9YrEBsba2RkVFBQoInotBSXy6VQKOnp6SiKJiUl0el0xY1Xh8/nE4nE5ORkFEXFYrGtrW1sbKxGQm04eGKnOtTW1t67d0/xSAmBQOjdu/e1a9d+UL6mpqZ5rrl99epVbW2t4tBmNpvdrl2769ev1ytz9erVAQMGAAAwGEx4ePiPr2Qzce3ata5duyo2qejXr9+HDx9yc3P/q3BNTQ0Oh2Mymf9VoJmQSqW3b99WvCvxeHyfPn2+vZfOnTs3atQoBEHu3btXXFysiTC1zt27d52dnR0dHQEAXl5eHA4nLi7u6wI5OTlYLNbDwwMAQCQSfXx8rly5opFQGw4mQnUoKipCUdTc3Fzxq7m5eUFBwbfF3rx507lzZy8vr8+fP0dFRX1boMkrKCgwMzOr25PI3Ny8sLDw6wJSqZTH4/30SjY3BQUFddN+VCpVX1//u5dl7ty5HTp0CAgIOHPmDNx3rbi4WC6X1123795LWVlZr169CgsLO3ToUIsWLc6fP6/2MLXO1zcb+N51s7KykkgkGRkZAAAEQVJSUrT/TQofclCapUuXftt98fLyOnnypGKfgbrPdxwOJ5PJvq3B2tp64cKFeXl5UVFR58+fnzRpkqpj1jZyufzrefVvL5RcLkcQpK7Mf13J5qbedcPj8d+9LMOGDevYsePVq1dXrFjRpUsXBoOhxhi1juJd+eN7qba2trKy8s2bN3g8/sqVK+PGjRs4cOCvngLdxPz0ZmOxWHPmzOnRo8fw4cPr1jqoPcxfAxOh0kyZMmXo0KH1/lGx+NPExAQAwOVyFV2ZkpISMzOzb2swNDTs1KkTAEBfX3/+/PnNMBGampqWlpbW/VpSUuLi4vJ1ATKZbGBgwOVy7e3twX9fyebm6+smlUr5fP53L0urVq0AAD179vTy8oqJial3gExzU/euVCxO/u69ZGZm5u7urngqOiQkhM/nFxYWNvM1t9++Sb+9btHR0b1793737t3q1asPHz6s/RM9zfqrjXJZWFi4f8POzg4AQKPRWrVqdefOHUXJO3fuhISEAAAQBCkrK/v26xKXy22eUzi+vr5isTgpKQkAIBKJHj9+rDgEQCKRVFRUKMqEhoZ+eyWbuZCQkPv37yu+mD948MDU1FRx4wkEgpqamnqFJRJJZWVlM+8OAgDIZLK/v/+P35UdOnTIzMxUFEhPTycSiVq4qbSaBQUFJSUl8Xg8AEB+fn5GRobi6FmRSFRVVfV1sWnTpjk5OV2+fDksLExj4TYM7BGqydKlSydMmMDlclNSUkpLSxV9x/T0dFdX1/Lycn19/VGjRhkbG5uamubm5p4+ffrIkSOaDlkD9PT0Zs+ePWTIkMmTJ1+7di0gIMDX1xcAcPHixaVLl+bk5AAAIiMju3TpgsFgysrKYmNjt2/frumoNa9bt25GRkYDBgwICgrasWPHkiVLFI/ETZgwgcPh7Ny588mTJ1u2bGndujUWi7169aqxsXG3bvDMDrBkyZLRo0fz+fyPHz/m5+cPHz4cAJCdne3o6FhaWmpsbDxixIitW7dOnjzZ3d19z549ixcvJhKJmo5aw2xsbAYPHtyrV6+hQ4eePHlScY8BALZt2/bgwYP79+8DANasWSOXy7FY7JkzZ0aMGNG+fXtNR/0TcGcZ9fm/9u6YJbkoDOD4ma4OgohIi0g4BNniFgmJDjno5GJc0BTa+gJ+Db+BODgoKA5aCgriIChGheikayhBmegiGA0Hgnd6F+FePP/fpjg8uPzxHjlPv99/fHx0OByZTMbpdAohvr+/y+VyOp3WNG08Hnc6ncVicXJyEovFzs7OjJ7XMLVabTAYeL3eu7s7i8UihJjP56PR6Pb2Vn5gPB5XKhWLxZJMJv/ZzaSwzWaTz+cXi0UoFLq5uZFvysWfV1dX2+220WhMJhMhxMXFRTwe5yICaTAY1Ot1u92eTqddLpcQYr1el0qlVColl6R+fn4WCoWvr6/r6+u/L1Zx+/2+WCxOp1O/359IJOSh6evr6/v7ezQaFUIMh8Nms7nb7YLBYCQSMXre/yOEAAClcUYIAFAaIQQAKI0QAgCURggBAEojhAAApRFCAIDSCCEAQGmEEACgNK5YA8yu1+t1u93T01OPxzOdTh8eHoyeCDgq3CwDmFq1Ws3lcu12W9O0bDbbarXkpeQADoVfhIB5fXx83N/fPz09ybuef35+2LYBHBxnhIB5FQoFu91+eXkpX3a7XbmXCsABEULAvJ6fnwOBgFwIvlqt3t7egsHgbrczei7gqBBCwLw0TftbyNVut30+n8PhKBaLxk4FHBlCCJiXrutyQ/rLy0un03G73cvl0mazGT0XcFT41yhgaq1WazabnZ+fh8PhfD5vtVp1XZcPSwEcBCEEACiNR6MAAKURQgCA0gghAEBphBAAoDRCCABQGiEEACiNEAIAlEYIAQBKI4QAAKURQgCA0n4Bh4Rb8FcPXuIAAAAASUVORK5CYII=", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#| include: false \n", "\n", "α = 2/3\n", "A0 = [0 1 1 ; 1 0 1 ; 1 1 0]\n", "L = -tril(A0, -1); \n", "U = -triu(A0, 1); \n", "\n", "aa = range(-1/2,1,50)\n", "ρα = zeros(50)\n", "ρJ = zeros(50)\n", "ρGS = zeros(50)\n", "for (i,a) ∈ enumerate(aa)\n", " A = I + a*A0\n", " D = Diagonal(diag(A))\n", " L = -tril(A, -1); \n", " U = -triu(A, 1); \n", " ρα[i] = maximum( abs.( eigen( I-α*A ).values ) )\n", " ρJ[i] = maximum( abs.( eigen( inv(D)*(L+U) ).values ) )\n", " ρGS[i] = maximum( abs.( eigen( inv(D-L)*U ).values ) )\n", "end \n", "\n", "plot(aa, ones(length(aa)); linestyle=:dot, ylims=[0, 2])\n", "plot!(aa, ρJ; xlabel=L\"a\", label=\"J\", m=2)\n", "plot!(aa, ρGS; label=\"GS\", m=2)\n", "plot!(aa, ρα; label=\"α = $α\", m=2)" ] }, { "cell_type": "markdown", "id": "5af52bef", "metadata": {}, "source": [ "Choose ***TWO*** sections from B, C, D (each worth ***30 points***)\n", "\n", "## C: Strictly Diagonally Dominant Matrices **[30 points]** {.unnumbered}\n", "\n", "::: {.callout-note}\n", "# Useful facts and definitions\n", "\n", "Recall that $|x|_\\infty = \\max_{i}|x_i|$ and $\\|A\\|_{\\infty} := \\max_{|x|_\\infty = 1} |Ax|_{\\infty}$ define a vector norm and the induced matrix norm, respectively. In this question, it will be useful to use the fact that \n", "\n", "\\begin{align}\n", " \\rho(A) \\leq \\|A\\|_{\\infty} = \\max_{1\\leq i\\leq n} \\sum_{j=1}^n |A_{ij}|\n", " \\tag{hint}\n", "\\end{align}\n", "\n", "(which you can use without proof).\n", "\n", "We say a matrix $A \\in \\mathbb R^{n\\times n}$ is *strictly diagonally dominant (SDD)* if \n", "\n", "\\begin{align}\n", " |A_{ii}| > \\sum_{j\\not=i} |A_{ij}| \\qquad \\text{for all } i = 1,\\dots,n.\n", " \\tag{SDD}\n", "\\end{align}\n", "\n", ":::\n", "\n", "Recall that Jacobi's method for solving $Ax = b$ is given by $x^{(k+1)} = D^{-1} \\left( (L+U)x^{(k)} + b \\right)$ where $A = D - L - U$, $D$ is diagonal, $L$ is strictly lower triangular, and $U$ is strictly upper triangular. Define $T_{\\rm J} := D^{-1} (L+U)$.\n", "\n", "1. Show that, if $A$ is strictly diagonally dominant, then $\\| T_{\\rm J} \\|_\\infty < 1$,\n", "2. Hence show that Jacobi's method converges for strictly diagonally dominant matrices $A$.\n", "\n", "Recall that the Gauss-Seidel method for solving $Ax = b$ is given by $x^{(k+1)} = (D-L)^{-1} \\left( Ux^{(k)} + b \\right)$ and let $\\lambda$ be an eigenvalue of $T_{\\rm GS} := (D-L)^{-1} U$ with corresponding eigenvector $v$.\n", "\n", "3. Show that $\\lambda Dv = (\\lambda L+U) v$.\n", "4. Hence show that, if $|\\lambda| \\geq 1$, then \n", "\n", "\\begin{align}\n", " |A_{ii}| |v_i| \\leq \\sum_{j\\not=i} |A_{ij}| |v_j| \\qquad \\text{for all } i=1,\\dots,n.\\nonumber\n", "\\end{align}\n", "\n", "5. By choosing $i$ such that $|v_i| \\geq |v_j|$ for all $j$ show that, if $|\\lambda|\\geq1$, then $|A_{ii}| \\leq \\sum_{j\\not=i} |A_{ij}|$.\n", "6. Hence conclude that Gauss-Seidel converges for all strictly diagonally dominant matrices." ] }, { "cell_type": "markdown", "id": "b959e72a", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**Answer.** \n", "\n", "**1.** We will use the hint, the fact that $(L+U)_{ij} = A_{ij}$ if $i\\not=j$ and $(L+U)_{ii} = 0$, and $D^{-1}$ is diagonal with $[D^{-1}]_{ii} = (D_{ii})^{-1}$: \n", "\n", "\\begin{align}\n", " \\| T_{\\rm J} \\|_\\infty &= \\max_i \\sum_j | [T_{\\rm J}]_{ij} | \\nonumber\\\\\n", " %\n", " &= \\max_i \\sum_j | [D^{-1}]_{ii} (L+U)_{ij} | \\nonumber\\\\\n", " %\n", " &= \\max_i \\sum_{j\\not=i} \\frac{|A_{ij}|}{|A_{ii}|} \\nonumber\\\\\n", " %\n", " & < 1.\n", "\\end{align}\n", "\n", "In the final line, we are using the fact that $A$ is strictly diagonally dominant.\n", "\n", "**2.** We have $\\rho( T_{\\rm J} ) \\leq \\| T_{\\rm J} \\|_\\infty < 1$ and so the iteration is convergent.\n", "\n", "**3.** We have $T_{\\rm GS} v = (D-L)^{-1} U v = \\lambda v$ and thus $Uv = (D-L) \\lambda v$ and thus $\\lambda Dv = (\\lambda L+ U)v$.\n", "\n", "**4.** If $|\\lambda| \\geq 1$, we have $\\lambda D_{ii} v_i = \\sum_j (\\lambda L+ U)_{ij} v_j$ and $D_{ii} = A_{ii}$ and so\n", "\n", "\\begin{align}\n", " |\\lambda| |A_{ii}| |v_i| &\\leq \\sum_j \\left| (\\lambda L+ U)_{ij} v_j \\right| \\nonumber\\\\\n", " %\n", " &\\leq |\\lambda| \\sum_j \\left| (L+ U)_{ij}\\right| \\left| v_j \\right| \\nonumber\\\\\n", " %\n", " &= |\\lambda| \\sum_j \\left|A_{ij}\\right| \\left| v_j \\right|, \n", "\\end{align}\n", "\n", "as required.\n", "\n", "**5.** If $|v_i| \\geq |v_j|$ for all $j$ then $|v_j|/|v_i| \\leq 1$ and so \n", "\n", "\\begin{align}\n", " |A_{ii}| |v_i| &\\leq \\sum_j \\left|A_{ij}\\right| \\frac{\\left| v_j \\right|}{|v_i|} \\nonumber\\\\\n", " %\n", " &\\leq \\sum_j \\left|A_{ij}\\right|.\n", "\\end{align}\n", "\n", "**6.** We have seen that, if $|\\lambda|\\geq1$, then $A$ is **not** strictly diagonally dominant. Therefore, if $A$ is strictly diagonally dominant, then $|\\lambda|<1$ for all $\\lambda \\in \\sigma(T_{\\rm GS})$ (that is, $\\rho(T_{\\rm GS}) < 1$). Therefore GS converges for all strictly diagonally dominant matrices $A$. \n", "\n", "**Extra:** (possible presentation) Show that Jacobi and GS converge for irreducibly diagonally dominant matrices.\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "fa18160c", "metadata": {}, "source": [ "Choose ***TWO*** sections from B, C, D (each worth ***30 points***)\n", "\n", "## D: Power and Inverse Methods **[30 points]** {.unnumbered}\n", " \n", "::: {.callout-note}\n", "# Useful facts and definitions\n", "\n", "Suppose that $A \\in \\mathbb R^{n\\times n}$ has orthonormal eigenvectors $v_i$ with corresponding eigenvalues $\\lambda_i$ (that is, $Av_i = \\lambda_i v_i$ and $(v_i,v_j) = v_i^Tv_j = \\delta_{ij}$). We have seen in lectures that, in this case, we can write $A = \\sum_i \\lambda_i v_i v_i^T$. \n", "\n", "We first suppose that $\\lambda_1$ is a *dominant eigenvalue*: $|\\lambda_1| > |\\lambda_i|$ for all $i\\not=1$. We call the corresponding eigenvector $v_1$ a *dominant eigenvector*.\n", "\n", ":::\n", "\n", "1. Explain why $\\alpha v_1$ is also a dominant eigenvector for all $\\alpha\\not=0$. \n", "2. Suppose that $x^{(0)} = \\sum_i c_i v_i$ for some constants $c_i$. Show that \n", "\n", "\\begin{align}\n", "A^{k} x^{(0)} &= \\sum_{i=1}^n \\lambda_i^{k} c_i v_i \n", "%\n", "= \\lambda_1^{k} \\Big[ c_1 v_1 + \\sum_{i\\not=1} c_i (\\tfrac{\\lambda_i}{\\lambda_1})^{k} v_i \\Big].\\nonumber\n", "\\end{align}\n", "\n", "In particular, the sequence $(A^k x^{(0)})_{k=1}^{\\infty}$ approaches a constant multiple of the dominant eigenvector $v_1$. Therefore, on normalising $A^k x^{(0)}$, we obtain a sequence that converges to $v_1$ or $-v_1$ (which are both dominant eigenvectors) as $k \\to \\infty$. Using this sequence one is also able to approximate the dominant eigenvalue $\\lambda_1$. We implement this so-called *power method* here:" ] }, { "cell_type": "code", "execution_count": 3, "id": "a9defcfd", "metadata": {}, "outputs": [], "source": [ "function power( A; Niter=100, tol=1e-8 )\n", " n, = size(A)\n", " x = normalize( randn(n) )\n", " λ = []\n", " for k ∈ 1:Niter\n", " x′ = A*x\n", " push!( λ, dot(x,x′) )\n", " x = normalize(x′)\n", " \n", " # stopping criteria\n", " r = norm(A*x-λ[end]*x)/norm(A)\n", " if ( r < tol )\n", " println(\"converged in $k iterations\\nλ ≈ $(λ[end]), residual = $r\")\n", " break;\n", " elseif (k == Niter)\n", " println(\"did not converge in $Niter iterations, residual = $r\")\n", " end\n", " end\n", " return λ, x\n", "end;" ] }, { "cell_type": "markdown", "id": "006ba7e6", "metadata": {}, "source": [ "3. In the stopping criteria, why do use the norm of the residual $\\|Ax - \\lambda x\\|_2$ instead of computing either of the errors $\\|x - v_1\\|_2$ or $|\\lambda - \\lambda_1|$?\n", "\n", "Here, we apply the power method to a random $1000 \\times 1000$ symmetric matrix and plot the errors in the eigenvalue against the number of iterations:" ] }, { "cell_type": "code", "execution_count": 4, "id": "64287c69", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "converged in 3461 iterations\n", "λ ≈ -89.03496134424624, residual = 9.988658606610933e-9\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdZ0AU194G8DOF3otUAbEERQQRULGBLREVEVtMjCXRaLrJTXITTW5iqi0mGm/UcKNGTdQb0aAoKsUGUkREVATFQhEQ6Z2F2Z33wyT77l0sgAuzyz6/T3Ccnf3vjOyzZ/bMORTP8wQAAEBb0WIXAAAAICYEIXSpzz77bPHixS0tLZ33FFlZWcuWLduzZ0/nPcVD3bp1a9myZb/++mtXPukvv/yyePHiO3fudOWTdoH4+PjFixdHRUWJXYjK7N27d/HixVevXhW7EHgIBCF0qUOHDu3YsYPjuKffVU5OTlhYWGpqqlL7vXv3wsLC4uPjn/4p2uX+/fthYWFnzpzpyic9c+bMjh07SkpKuvJJW/vtt9927Nihwh3euHFjx44dGRkZKtynuJKSknbs2FFQUCB2IfAQCELoUkOHDg0MDGQY5ul3lZKSsmzZsiNHjjz9rjSXu7t7YGCgmZmZuGV8+OGHb731lgp36ODgEBgY6OLiosJ9AjwKK3YBoF1U22+AlStXrly5UuwqVG/y5MmTJ08WuwrQFgjC7qOlpSU+Pv769escx/Xr12/8+PH6+vqKG1y+fJmiKC8vr8bGxuPHj+fl5bm5uU2ePDk7O7u+vt7T05OiqLi4uBs3bpibmy9YsEB4VGVlZWxsbEFBgYGBgZ+fn6+vr+I+eZ6/dOmSvr7+wIEDa2trjx8/fu/ePW9v77Fjxz60yOvXrzc2Ng4ZMoSiKEKIRCK5du2aqalpv379KioqoqKi7t+/7+LiEhQUZGxs/JgXm5OTk5ubSwgpLi5OS0sTGl1dXS0tLRU3KywsPHnyZGVlZb9+/SZNmqSrq9t6VyUlJXFxccXFxcbGxqNGjRo4cOBjD/P/1BAXFyeRSAYNGhQYGPiozerq6mJjY+/evaujozN48OARI0bQ9P9cicnIyJDJZN7e3g0NDVFRUQUFBY6OjsHBwQYGBsIGycnJFy9epGl6woQJzzzzjOJjc3Nzy8vL+/fvb2RkJLTIz3JTU9OxY8dyc3N79Ojx7LPP2tnZKRVWX1+fmJiYn59fUVEh9MAcHR0VN2hpably5YqRkVH//v2rqqqioqKKioqcnJwmTZok74PW1NTk5OS0tLTIZDL5iTA0NBwwYMDjj15+fv6ZM2dKSkosLS0DAwP79Omj+K/l5eW5ubmOjo5KZd++fTsuLq6hoaF///4TJkzgOC4zM9PMzKxv376Km/E8f/HixYsXLzY0NDg7O0+cONHc3Fxxgxs3btTV1Xl6etI0HRsbm5WVZWBgMH78eMX9XLt2TTi5rf/bCP/k6empo6Mjb7l27VpxcbGurq6Xl1frs9xaWlqanp6eh4dH6xfu4OBgb2+v2N7c3Hzu3LmsrCypVOrm5jZu3Dg9Pb3H7x/agYdu4fTp071791Y8s05OTvHx8YrbmJqaWlpaxsfHy99cpk2bxvP88OHDCSHx8fHu7u5Ce9++fYWH/PzzzyYmJoq7DQwMLC4ulu+zubmZEPLMM89ERUVZWFgI27zyyiuPqlOImYaGBuHXW7duEUImTJjwxx9/KCafg4PDlStXHvN6n3322db/mX///Xee56Ojowkhr7766qZNm+TvU4SQAQMG3Lt3T3EnHMd98MEHSm9zc+bMqaure/zRlslkH3/8seI73dChQw8ePEgIWbhwoeKWBw8etLa2Vty/t7d3Tk6O4jY2NjZGRkYpKSkODg7yzZydnXNycqqqqiZNmiRvZFl227Ztio+dN28eISQxMVHeYmJiYmVllZqa2rNnT/kDjYyMDh8+rPjA1atXK31O0tHR+fjjj2UymXyb4uJiQsjw4cMjIyMVr7726NEjNTVV2ObkyZOtT8SQIUMec/QaGhoWL16sePRomn7ttdeam5vl2/znP/8hhHz77beKD/z0008VL6p7eHgIo2mCg4MVN8vOzlb6uGZmZvbbb78pbjNixAhCSEpKire3t3wzhmG+++47+TbCZ8H//ve/SvUXFhYyDOPi4iKVSnmez8vLa30J19PTU+ksC5eOjx07Jvwqk8mIwh+a3K5duwghq1atUmyMjo5WegoXF5fk5OTHHGRoFwRhd5CSkqKnp2dkZPT111+npqZmZGSsW7fOwMDAxMTk9u3b8s1MTU0NDAx69Ogxf/78Q4cOnT9/Pioqiv87CIUPzvv27UtMTDxw4ADP8/v27SOEWFhYhIWF3bx5MykpKTg4mBAi9DaEfQpBaGlpaWJi8tprrx05ciQhISE2NvZRpT40CB0dHY2MjD799NP4+PizZ8/Onj37iW+mly5d+vjjjwkh8+bNi/mbkNBCEPbq1cvExGTdunUpKSkxMTEBAQGEkJkzZyru5NVXXxWe6NChQ1lZWadPnw4KCiKEzJkz5/EHfNOmTcJTHDlyJD8//+zZs35+fkKMKQbhmTNnGIYxMDD47rvvsrOz09LSFi5cSAhxcnKqqKiQb2ZjY6Orq2tvb7948eJTp06dP39+7ty5hJCxY8fOmDHDy8srPDw8LS3thx9+0NPT09PTy8vLkz/2oUFoaGhoa2u7aNGimJiYlJSUDz74gKIoc3Pz6upq+Wbvvffe0qVLDx8+fPXq1WvXru3atUvokykGrRCEdnZ2xsbG//znP8+ePRsfH//SSy8RQvr37y9EZllZWUxMjKWlpZ6envxEPOY9WiaTTZ06lRASEBBw7Nix7OzskydP+vv7E0KWL18u36x1EG7ZskUIgIiIiPz8/JSUlEmTJgnHXDEI7927Z2Njw7LsO++8k5CQcP369R07dvTo0UPo+ck3E4LQ1dV14sSJkZGRFy9e3LBhg76+Pk3T8k9gcXFxhJDJkycrvYQ1a9YQQj7//HPh18zMTH9//61btyYkJNy6devs2bPCWXZ3d+c4Tv6oDgfhuXPndHR0TE1N165de/HixcuXL3/77bd6enrm5ub5+fmPOs7QLgjC7kC40nj8+HHFRmEcv2LnzNTUVEgOpYcLQejn56f4d9vS0iJ0KSIjI+WNHMf5+fkpvl0KQUgIeeedd9pS6kODkBCyfft2xWcR3pQVU7w14QaJTz/9VKldCEJCSFxcnLyxqqrK1NRUR0dH3u1ITEwUehWNjY3yzaRS6bBhwwghaWlpj3rehoYGCwsLmqavXbsmb6ysrBSuyioGodAvUerDCR8m/vWvf8lbbGxsCCGvvfaavEV+8J2cnOrr6+Xt77zzDiHk3//+t7zloUFICHnrrbcUnzQ0NJQQcvDgwUe9KJ7n7969q6+vP2DAAHmLEISEkE2bNskbZTLZoEGDCCEZGRnyRjs7OwMDg8fsXC48PFxIQaE7JWhsbOzTpw/LsvJ3dqUglEgkNjY2FEWlp6fLH9Xc3NyvXz+lIBRyesOGDYpPmpKSQtP00KFD5S1CEI4bN06xB/yvf/2LEPLFF1/IX2nv3r1Zli0qKlLcm4eHB0VRt27deszLfPHFFwkhin+SHQtCmUzm7u5O0/Tp06cVN9u6dSsh5I033nhMDdB2GDWq8TIzMy9duuTr66t4DY0QsmDBAmNj4+PHjytt/9FHHz10P++//77idafU1NR79+55eHgIn98FDMN8+OGHhJBDhw61cbdt0aNHj0WLFik+y/jx4wkhd+/e7fA+fX19x40bJ//VzMxs2LBhLS0t8vHrQo5+8MEHilcIhWt0hJDH3MF29uzZysrKSZMmKX6baG5uvnjxYsXNCgoKLl68aGtr+/LLLyu2Cx3Z1gfwH//4h/xnlmWFd+qlS5caGhrK24V+bVsOi3Ca5CZOnPjEB/bq1cvLyys7O7u2tlax3cjI6I033pD/SlFUW/b2KMJh/+STTxQvjerr6y9evJjjOPmHGCUpKSkPHjwIDAwcPHiwvFFHR0dpqGpTU1N4eLipqeny5csV24cOHerr65uamlpWVqbY/v777wvfVQuUXhdFUS+99BLHcb///rt8mwsXLly7dm3MmDFKX2oqCQkJETZ+zDZtkZaWdv369VGjRil9Cb1kyRI9Pb3Wf93QMRgso/EuXbpECJFIJMI7rCJ9ff3i4mKJRCL/Xp1hmP79+z90P0qDRK5fv04IGTJkiNJmPj4+hJDMzEzFRgsLC8Xvt9qrX79+SiMLbG1tCSFPc3ucm5ubUot8n8KXqenp6YSQ06dPZ2VlKW6Wn59PCBFG4jyUcGS8vLyU2hXfo8nfh8jDw0PpO0ih+56dnS2TyeSvWkdHR+kr3h49ehBClIbGCH3HJx4WQ0NDJycnxRbhtd+/f1/ewnFcWFjY/v3779y5c//+falUKv+n8vJyxS+Ghb5a67117OwIh/3w4cPChUc54ag+6rDfvHmTECL/DltO6T9tVlZWU1OTpaXlJ598orRlTU0Nz/N5eXmKX9kq/SdpfZQWLVr01Vdf/frrrx988IHQIvTYhIufcpcvX16/fn16enpBQUFdXZ28XSl3O0A4XHV1da3/ug0MDPLz8xX/F0GHIQg1XlVVFSHk5s2bD71X18LCorGxUR6EZmZmiuNHFCmN6RD+noV3XkXCm4VSp0Hpse2l2OkRCH/bwuWjTtqncNwiIiJav48IVz4ftWfhyAhBpUjpWD3qAOrr65uZmVVVVdXX18vzRldXV+neSqEA+cBRgdB94Z80P7CBgYFiR0e+N8UHzp079+DBgz179pwyZYq9vb3wRNu3b8/JyVGa7qD1kRR23rGzIxz2vXv3tv4n+WCr1hobGx+6gVKLsPPS0tKwsLCH7l9+JV+g9NJaHyVXV9fRo0efO3cuLS3Nx8enubn5jz/+MDIymjVrlnwb4atlnufHjRs3bdo0oaSbN2/u2LFD8eNFxwivKCsrq3X/m6IoU1PT5uZmpUFP0AEIQo0nvJm+8MILqr1FT9jtgwcPlNqFz8vC141ySm+7GkEYpBobG6s0wvCJHnVklHpIj9qssbGxurqaYRj5DQ9dLzEx8eDBg35+fmfPnlXMWmF4VKcyNTWtqam5du1auy4hCOlSVFSk1F5YWKj4q3DMvb29U1JSnrrSvyxcuPDcuXO7du3y8fGJjIwsKytbuHChYo/5k08+kUgkJ06ceO655+SNu3fvfuLfI0VRrT9M1NfXK/4qPNGiRYuEsULQSdCn1njC+G9h6IcKCbc3CWNGFNuFKc2Ubn4ShdC17fBsbR0+bsLlOOGalSLhGrWccIiuXr0qkUgU2y9evMjzvDACor1PrSrC1GUhISGKKVhbWytcgewAHR2dNvZ+OnbYhUv08fHxSs+iNKGdu7u7np5eRkaG4vXJpzRnzhxjY+O9e/dKJJKHXhfNyMiwtLRUTEFCiPyWykehKMrOzq715WWlC/XC4UpKSupw/dAWCEKN5+np6ePjc+PGjZ07d7b+1w6/I/j4+PTq1SszMzMyMlLeyHHc+vXrCSHCHQ7iEu7+vnfvXsceLoxh2bBhQ2lpqdI/cRwnXIt7qICAACsrq+jo6CtXrsgbKyoqlHoAjo6OI0aMePDggWI7z/OrV68mYh9A4bqu8G2o3OrVqx/zqh/P0dGxubm5Ld8aCqOivvrqK6WuDyFEIpE8ajb2AQMG+Pn55eTkbN++Xd54+/ZtpUughoaGc+bMkUgkq1atar2Tjv0tGBsbz5gxo7y8/Ndffz1x4oSLi4swZEnO2tq6tra2srJS3lJQUPDLL788cc+urq719fVnz56VtxQWFipN2u7n5zdw4MDLly/v37+/9R5UmPdaDkGo8SiK+uWXX4yMjJYsWfL6669HRkZevXo1Li5u27ZtY8eOffPNNzu2W4ZhNmzYQFHU/PnzN2/enJmZGRcXN2XKFGGE6vz581X7Kjpg4MCBBgYG4eHh77///pYtW8LCwuQ3Y7SFv7//8uXL8/PzfX19N2/eHB8ff/ny5YiIiJUrV7q4uCiGnBI9Pb2vv/6a5/kpU6bs378/JycnKipq/PjxSpeLCSHr169nWfa999779ttvr1y5kpCQMHfu3OPHj7u6uioNa+xi/v7+enp6O3bsWLt27dWrVy9evPjuu+9+//33SkNs2k64vDx79uwNGzaEhYW1HhMrFxoaOmfOnCtXrvj5+YWFhSUnJ6elpYWHh7/33nuOjo6KA1WUbNmyxdDQ8PXXX583b97GjRs//PBDPz8/YWytonXr1jk5OW3YsGHGjBnh4eEZGRnx8fG//vrrjBkzJkyY0LFXJ3QB//GPf7S0tCxcuFCpKz927NiWlpaQkJBTp07dvHlz7969AQEBrefxaU24W/Sll17as2fP+fPnt23bNnToUKUJZWia3r59u76+/ksvvfT2228fPXr06tWrsbGxW7ZsGT16tHwIDzwt0W7cAJVKT08fOnSo0sm1sbHZvHmzfBthZpnWjxXuI1ScL0Zuz549SpOWTZ48ubS0VL6BfGaZNtb5qJlllDb7/PPPCSF79ux5/N4OHDigOOOG0swyShsLE4UkJCTIW2Qy2fr165Um36Ioys/P786dO49/6i+//FJxLOXo0aMjIiJIq5lljh07pvRl2IgRI3JzcxW3EWaWUdr/22+/TQg5cuSIYmNCQgL53ztBHzWzjNLe/vzzT0LI+++/L2/5/fffFa+LmpqaHjhwQLi+J58SRT6zjNLe1q1bRwjZunWrvKWsrCw4OFg+DuvxkyG0tLSsXLlSaaAKwzCjR4+WzzPw0JllkpKShNtYCSHm5uaffvrpuXPnCCEvvPCC4mb37t1TvOdHYGJi8t5778m3ERJU6QbBnJwcQshzzz2nVLBwQ6Hwf6P17YOlpaVKf3qhoaHCRwHF+/yU7iMUjsMLL7yg+MBly5YJ1w+UZpa5cOGC4gw4Ajs7u7CwsMccZ2g7iscK9d2IcE9hfX19jx49XFxcvL29FcciCmPTe/XqpfSowsJCiUTi4uLy0EUh6uvrz507l5eXZ2ho6OfnpzSHJM/zwiyabexM5OfnSySSvn37CuNrhBv7DAwMlD4IV1ZWVlZW2tjYPH7GUUFDQ0NJSQnP88L2jY2NxcXFJiYmSgM7Hzx4UFdX5+DgoDTKrqGhITExMTc3l6ZpOzs7Ly8vpSk3H6WgoODMmTMSicTd3d3f318ikRQVFbV+3qampvPnz9+6dUtPT8/Ly2vw4MFKY4vy8vJkMpmrq6tiY1lZWU1NjZ2dnWJgNDU1FRUVGRsbywejlpSU1NTUODk5yV9Ubm4uRVFKM3LV19eXlJSYmZlZWVnJG0tLS5OTk4uKiuzt7ceOHWtiYlJcXNzY2Ojk5CREmlQqzcvL09PTUzogVVVVFRUV1tbWSp1gmUx2//79pqam1g9praam5vz58/n5+fr6+vb29oMHD1YcYVtTU1NSUmJtbd16pGhFRUVjY6OtrS3Lsrt37164cOHHH38sXHBWOqrJycmVlZVmZmbOzs6+vr6Kk3MWFRU1NTUp/Z9/1P9GQkhJSUl9fb0ws1rr1yKTyVJSUrKzsxmG8fHxGThwoPCf0NTUVD6gurS0tKqqytHRUekTQHp6+sWLF3V1dYcPH+7m5lZXV/fgwQMLCwulF87z/LVr1y5fvlxfX29jY9OrV6/BgwfjxglVQRACgEaSyWSBgYHx8fHR0dHCvfAAHYMPFACgASoqKoKCgv773/9mZmbeuXMnKioqKCgoPj7e399fmIcIoMPQIwQADVBVVWVnZ6d0L8qzzz67Z8+e1rMWALQLghAANENDQ0NSUlJBQUFVVZW5ufmwYcOeuOohQFsgCAEAQKvhO0IAANBqCEIAANBqCEIAANBqCEIAANBqCEIAANBqCEIAANBqah2EEolkzZo1bd/+aRY0B7HgrGkcYZ5isauA9sEpewy1DsKKiooff/yx7du3XuEM1BzP8x1eAw/E0tLSIqw6AhqksbERHzofRa2DEAAAoLMhCAEAQKshCAEAQKshCAEAQKshCAEAQKt1nyCUyWT379+XSqViFwIAAJqkmwThvXv3+g4ePuqld109h968eVPscgAAQGN0kyBctW5jbuAnD16NKJi28cMv1opdDgAAaIxuEoRNkmZe35gQQvSMTuU3fZ4mvVKBaRQAAODJWLELUI1P3n39VOiLEld/nbvJG7dsy6bI7Dhps4xMc6Zmu9Ij7ShK7AoBAEA9Ueo8AV1xcbGPj09RUVFbNq6trb1w4YKPj4+5ubnQklnJH7gr++MOX9dCQntRwc50oD3FdpM+cDfB83xDQ4ORkZHYhUA7NDc38zyvp6cndiHQDg0NDfr6+jSNd8CH6D5BSAipra01MTFp3Z5ZyR/N5yPzZTeq+aCe9Oze1KSetA7+P6gBBKEmQhBqIgThY3STS6OPN9CCGmhBfeRF59byh/P4tRmyl89KJzvRs3tTz/WkdfEfAwBAi2lXCPQyoZZ70AnB7KVQ1sea+jFT5vB7y5w46e4cWV2L2MUBAIAYtCsI5ZyNqeUedEwQmz1bZ6ozdeCuzHFvS3A0tztHVtNC3vrwE/v+3v39xlzOyBC7UgAA6FxacWn0Maz1yYJ+9IJ+dFkTiciT7b8je/0/0S0phS3vp9wvy33h1VezLpwVu0YAAOhEWtojbM1anyxxo6OeY9e7lcicBhOKIta9SmvqxK4LAAA6F4JQ2Ywpz9mk7WLObDH8fZnEbcKCM9IKidg1AQBAp0EQKrOzs7t05vjWQONDK+beO7DWQo94HuIO3pWJXRcAAHQKrbiP8Ckl3OeXJkj7mJJtIxlHI8xRo0q4j1AT4T5CTYT7CB8DB+XJRtlRwu0WQyK4sGyZ+n5wAACA9kMQtok+Q1YNYWKD2O03ZAFHuRvVSEMAgG4CQdgOgyyppGnsS33p0ZHcqkvSFnxvCACg+RCE7UNTZGl/OiWETSrh/SK41FJ0DQEANBuCsCNcTaiTQewn3vS0aG55krSeE7sgAADoKARhx812pTNm6FRKiOdBLrYQXUMAAI2EIHwqNgZkdyCzeQSzJF664Iy0HLfeAwBoGgShCkx2oq7PYh2MiOdBbncOhtAAAGgSBKFqGLJkjR8TPoFZd0UWHM0V1ONKKQCAZkAQqpK/DZUeyo6ypX0juE3XZLj3HgBA/SEIVUyHJh950QlT2Yg82Zij3P7Y5JFBoWOmzLyYliZ2aQAA8BDavh5hJ+lnRp2awv50uW5eyOuypb8TKRc8b2HelRRdXV2xSwMAgP+BHmFnoQgZr1dg0duD2PYjDgMk1v0KCgrELgoAAJQhCDtRnz59jEqzqIyjVPrh2uK7q/KcHjSKXRMAAPwvBGEn0tXVjY86tEwn5Q2DS1mnIvqYM4MOtWy6JpNiEA0AgNrAeoRd6mY1/2aitEJCtoxghtlgaUOsR6iRsB6hJsJ6hI+Bg9KlnjGjooPY5QPp6TEcZqIBAFAHCMKuRhGyoB+dNVvHQo+4h7dgpV8AAHEhCMVhrks2+TMnJrE7bsoCj3LXKpGGAADiQBCKyduKSgxmF7vRE6O45UnS2haxCwIA0D4IQpHRFFnQj86cpdMkJQPCMWc3AEBXQxCqBUs98vMoZt9YZv0V2YQo7kY1rpQCAHQR0YKQ5/nMzMysrCyxClBDo+2o9FA22Jkec5RbdUnaJBW7IAAALSDaXKNLlixpaGiQSCTOzs4bN24Uqwx1w9JkuQc9uzf18QXZoIPc5hHMpJ643RAAoBOJE4SZmZmZmZnJyck8z3t5eeXl5bm4uIhSiXpyMKR2BzKnivg3E6U/mZKfRjDOxohDAIBOIc6l0UuXLg0fPpwQQlHU0KFD09PTRSlDzY1zoC5NZ32sKe8/uVWXpM0ycvr06b1791ZWVopdGgBA9yFOj7C6ulo+q5aRkVFVVZUoZag/A5asGsLM70u/nSR1mPFhc2WJxN7d5uvxGQmxlpaWYlcHANAddGmPsLa2VvjBxsamrKxM+LmsrMzW1rYry9A4fUypqOdY6fVTtfO3N094v9Tz+ejoaLGLAgDoJlQchOfOnZs+fbqLi8v48eMV2y9duuTm5ta7d++ePXvGxMQEBgaeOXOmvr6+uro6JSVl5MiRqi2jW7I01CXVxYTnmwsy0zkb3GABAKASKg5ClmWnT5/+yiuvyDt/goULF77xxhulpaU//fTTiy++aGpqumLFinHjxk2cOPHbb781NTVVbRnd0u4t37vunmO7YdgsL/t40zFjMDEbAIAqdMoyTHv27Nm8efOFCxeEX9PS0saOHVtWVqarq0sI6d+//7fffjtjxown7icvL8/NzW3EiBHylrlz57744ouP2r6urs7Y2Pipy9cAPCH77tKfZbAznGWfe0qNWE1NRJ7nGxsbDQ0NxS4E2gHLMGkirV2GSV9fn2WfMBqmKwbL3L59u0+fPkIKEkL69+9/+/bttjxQV1fX2Nh4xYoV8pa+ffs+5k1TKpVqz1vqkoEkuDf55wWZ73H2h2HUjF4aeX+F8DlMe85a98CyLIJQE2lnELblJXdFEFZXVyu+05mYmLT9BgBdXd2JEye2cWOaprXqNNsbkT1j6bPF/Ovnpbtu8Zp4uyHP89p21roBmqaFEyd2IdAO9N/ELkQddcVB6dGjR3V1tfzXqqoqGxubLnheLRFgT2XMYEfZ0t5/cmszZBxm7QYAaI+uCEJ3d/dbt24Jw2d4nr906dLAgQO74Hm1hw5NPvKiL4Swp4tlvhFc8gNN/coQAKDrqTgIy8vLY2NjMzMza2pqYmNjL1++TAh55plnRowY8dFHH927d++bb74xNDRUurkCVKKPKXViEvuVLz0nTrrgjLRcInZBAACaQMVBmJubu3bt2rS0NCcnp7Vr1x48eFBo/+2338rLywMDA5OSko4ePYrr1J0n2Jm+MpO10COeB7ndOTL0DQEAHq9Tbp9QleLiYh8fn6KiojZuX1tba2Ji0qklaZBLZfxr56VGLNkykhlgrqaDaHieb2hokM+3BxoBt09oIq29faItcFC6rSHWVGIwO92FHh3JfZyK1Q0BAB4OQdidCasbXpnJFtWTQSfbsvUAACAASURBVAe56EL17f0DAIgFQdj9CasbbvRnliVI58RJSxrFLggAQJ0gCLXFFCfq+kzW3YJ4HmpZn1o7b+nb/YaMfOejT6VSXDMFAK2GINQiwuqGcZPZ71d/ua/R7dbS6O1ZLT9uDRO7LgAAMSEItY6HBeVYfZ33nk5Y3YZBIfFpV8WuCABATAhCbfT8tMmmkZ+Q7NN6Ud+c7TFpdw6mZQMA7dUVk26DuvngnTec7G1OJcXNWLPcyGvi6+elB+7K/j2CcdG0ObsBAJ4ebqgH0iIjW67LvrksfcOdXjmY0e3CywS4oV4T4YZ6TYQb6h8DBwWIDk2We9DJIWzKA94vgkvCnN0AoE0QhPCX3ibU8Uns17703FOYsxsAtAiCEP5HsDN9ZQbm7AYALYIgBGVmumSTPxP5LPPv67Kxx7jsKqQhAHRnCEJ4OGHO7lAXelQk93GqVIL5ZwCgm0IQwiMpztntcZCLwZzdANAdIQjhCYQ5u38YzixNkM6Jkz7AnN0A0L0gCKFNpjpTmTNZdwsy6FDLpmsyjKIBgG4DQQhtZciSVUOYmCD2j7uyMUe5a5UIQwDoDhCE0D6ellRCMLu0Pz0xilueJK1rIVKp9M6dO/X19WKXBgDQEQhCaDeKkAX96MxZOk1S4ra7zNl79IhFH/XxGR0bd0rs0gAA2g1BCB1kqUd+HsVMvfd7sddLJQv2liyLfPfTr8QuCgCg3RCE8FT6mBCaYQghhKLLGmVSfG8IAJoGQQhPZcmi+b0u7+yxb7Hl1imOsz7yjeBSMGc3AGgUrEcIT8XS0jL7YkJ2draTk5OZmdmBu7LpMdw0F3r9MMZUR+ziAADaAD1CeFosy3p4eJiZmRFCZrvSWbN19BniHs5h4XsA0AgIQlAxc12yyZ/57zhm/RXZ+CjuZjWulAKAWkMQQqcYaUulh7LTnOlRkdyqS9JmdA4BQF0hCKGzCHN2p4ey1yrJoIPcqSJ0DQFAHSEIoXM5GlHh45nvhtGvnJMuOCMtaxK7IACA/4UghK4Q7ExnzmIdjIh7eEtYNqbsBgA1giCELmLEkjV+THQQu+OmLPAodx0L3wOAekAQQpcabEUlBrPz+tJjIrmPU6VNWPgeAMSGIISuRlNkaX/66kydonoy6JA0thj/CQFATHgPAnHYG5Ldgcym4fS7F9ngaK6wHldKAUAcCEIQ02QnKjWo2ceaGhLBbbomKygsWvXN2h+3bGtoaBC7NADQFghCEJkBS1YNYeKnsn/eqO47euqXufb/jK8dHzJH7LoAQFsgCEEtPGNGfW6VwbqP5Ye9IBn/7p0HNegUAkDXQBCCuujt6mqUf4FI6kh5XnlNXVyZvtgVAYBWQBCCunBxcfnu47f6/DLFM/K1n7b+/PEFWXA0l1+HQTQA0LkonlffN5ri4mIfH5+ioqI2bl9bW2tiYtKpJYFq8Tzf0NBgZGTU+p9aZOT7q7J1V6T/9GTeH0Sz+MymNpqbm3me19PTE7sQaIeGhgZ9fX2axh/SQ+CggJrSoclHXvSFEPZ0scw3gkvGwvcA0DkQhKDW+phSJyaxX/nSc+KkC85IKyRiFwQA3Q6CEDRAsDN9ZSZroUc8D2HhewBQMQQhaAZh4fvDE5nNmbJxx7gbWPgeAFQEQQiaxMeaSprGhrjQoyO5VZekEszZDQBPDUEIGka+8H1mJRl0iIstRNcQAJ4KghA0kqMRdWA8s2EYvSReuuCMtBQL3wNARyEIQYMFO9PXZ7EORmRgeMumazKsfA8AHYAgBM1myJI1fkzMZHb/HVngMS6zEmEIAO2DIITuwMuSSpzGLnGjx0dxy5Ok9ZzYBQGA5kAQQjdBEbKgH315hk6lhHge5DYeSR44PNDVc9gvu34TuzQAUGus2AUAqJKdAdkdyJy4xwePep17809iZLHih5DxY0a6urqKXRoAqCn0CKEbCrBqsjQ1IhaORNewoafvnTt3xa4IANQXghC6IQMDg772lvrR6+ik3bKs0+8/8E4txSAaAHg4BCF0T3GH/9j4rN3X/avvJkb9w9csJIZbniStbRG7LABQP/iOELonfX39ZUsWCz8vsCXBzvSqS1L3cO4bX3pBP3z+A4D/h3cE0AoWemSTP7N3LLPuiiw4msvDwvcA8DcEIWiR0XZUeig7wYH2i+BWXZI2Y0EnAEAQgrbRoclyDzo5hE15wPtFcElY+B5A6yEIQRv1NqGOT2K/9qXnnpIuOCMtx8L3AFoMQQjaK9iZvjKDtdAjngex8D2A9kIQglYz0yWb/JnIZ5l/X5eNPcZlV+FKKYDWQRACkCHWVGIwO92FDjiGhe8BtA6CEIAQhYXv79QQj4NY+B5AiyAIAf6fgyG1O5D5fji9JF46J076oFHsggCg8yEIAZQJC9+7W5BBh/5a+L68vDwjI0MiwehSgG4IU6wBPIQhS1YNYWb0ol9LkIaFH7u//0vGycO4JPPCqShra2uxqwMAVUKPEOCRPC2phGC2POK7itePlj4flue3LGznbrGLAgAVQxACPA5NkR6GDOEkhBC+uelOPSN2RQCgYrg0CvAEm1evmrt0GrFy0ZdUxD97ODia2zKScTKixK4LAFQDQQjwBIEBY3IzkktLS3v27CmRUWsypL4R3Eov5q2BNIM0BNB8uDQK8GT6+vpOTk4URekzZNUQJmEqe7RA5hfBXcDC9wCaD0EI0G79zKiYIHbFYHpaNLcsQVqDhe8BNBmCEKCDZrvSWbN09BniHo45uwE0GIIQoOOEhe/3j2XWXZFNPcnl1uJKKYDmQRACPK1RdlR6KDvRkR56GAvfA2geBCGACggL36eEsBdKed8ILrEEXUMAjYEgBFAZVxMq6jn2G1/6hdPSBWekZU1iFwQAbYAgBFAx+cL37uEtYdky9A0B1ByCEED1hIXvTwax22/Ixh7jsrDwPYAaQxACdBZvKyppGvtiH3p0JPdxqrQJC98DqCUEIUAnoimytD99daZOUT0ZdJCLLuSjjp949d1//r7/vzyPbiKAWsBcowCdzt6Q7A5kjhXw8388UhUd1jz+3T/+s+9e8YOP3ntb7NIAAD1CgK4yxYkKLI9tDlpB+gfWTP36jyNRYlcEAIQgCAG6kr/3QMMrEaSliUoLL7LwuFaJq6MA4kMQAnSdt197dZmnUb9fguab3Ply1WcTo7jlSdI6zNkNICpKnb+xLy4u9vHxKSoqauP2tbW1JiYmnVoSqBbP8w0NDUZGRmIXIo4KCVmRKo0q4DcOp2e6asyn0ubmZp7n9fT0xC4E2qGhoUFfX5+mNea/WVfCQQEQjaUe+XkU81sg8680WXA0l1+nvp9KAboxBCGAyALsqYwZ7Chb2vtPbm2GTIo0BOhaCEIA8enQ5CMvOiWEPVUk843gUh4gDAG6DoIQQF30NaVOBrErB9PTY7DwPUDXQRACqJfZrnTWbCx8D9B1EIQAasdcl2zyZ/47jvnuqmx8FHezGldKAToRghBATY20pS5NZ6c506MiuVWXpBLM2Q3QOUQLwpqamsrKysrKyurqarFqAFBzLE2We9Dpoey1SuJ5iDtVhK4hgOqJNun2sGHDHBwcaJq2sbH5/fffxSoDQP05GlHh45nIfNkr56Rj7KjvhzPW+mLXBNCNiLn6RFRUFCanAGijYGd6nAP9VbrUPbzla1/m1f40JXZJAN2DmN8RLly4cOnSpZcuXRKxBgANYsSSNX5MdBC746Ys8Cj36opvHN19BgwNuJSeLnZpABqs03uE33//fWNjo2LL8uXLjY2Nt2zZ0qtXr+zs7OnTp58/f97JyamzKwHoHgZbUeeD2Q/3JmyMu8y/l1hUnj93ycs30xLErgtAU3V6ELq5uTU3Nyu26OjoEELGjh1LCHF1dQ0KCoqPj3/xxRc7uxKAboOhyDC2UNfVR0LRxLpXeV3jkx8DAI/Q8SD8+uuvExISbt26tXXr1okTJyq2b9q0qaWlZe7cuZs3b54yZUrrxzY3N1MUpaOj09zcnJaWtmDBgg6XAaCdxo0bZ/X5syX6pnrF1zjX4cHR3LaRjKMRvjcEaLeOf0cok8kWLVrEcVx9fb288cSJE1u3bk1NTc3NzU1LS/vpp58e+tiSkhIvL68xY8YMGjQoKCho5MiRHS4DQDv16NEj7czxjcPp/W8HFR3Z7GNNDYngNl3DnN0A7fa06xEOGDBg9erV06dPF359/vnn3dzcvvzyS0LIvn37Vq9efeXKlUc99onLB+bl5Q0YMCAoKEjeMmPGjNDQ0EdtX1dXZ2xs3O7XAOLR8vUIVSunhixPZapbqM1DpUMsOzEPsR6hJtLa9Qh1dXVZ9gnXPlX8HeHNmzdnzpwp/Dxo0KCbN2/yPE9RD79c88RFdHV0dPT19efMmSNv8fDweMyfX3NzM/44NQvP81KpFGdNJTx6kNjJ5Ldb/Oyz1GxX6sshlIlOpzwRRVEIQo0j/KFpYRC25SWrOAgrKyvl8WZiYiKRSJ7m8z5FUfr6+s8//3wbt2cYhmGYjj0XiILneZw11VrkRqb1IitSpR5/8huH07M6YeF7hmGEE6fyPUPnEf7QtDAI20LFB8XKyqqmpkb4ubq62sDAAFe9ALqYsPD93rHM51j4HqANVByEbm5uGRkZws9Xrlxxc3NT7f4BoI3G2FGXFRa+57CgE8AjdDwIb968mZaW1tTUdPv27bS0NGHs6JIlS3bs2JGZmVlcXLx+/folS5aorlQAaB9h4fsLIezpYplvBJeMhe8BHqbjo0aVZkfbuXPnoEGDCCE//vjjxo0bJRLJvHnzVq9e/TRfJBQXF/v4+BQVFbVx+ycOQwV1g1GjXSYyX/bmeVmgPbXRn7F8umEuGDWqibR21GhbPO3tE50KQdjtIQi7UlUz+TxNejCX/9aXXtCv42+ICEJNhCB8DBwUAG0hLHx/eCKzOVM27hh3AwvfAxBCEIQA2sbHmkqaxoa40KOx8D0AIQRBCKCFFBe+H3SIiyviCSHp6elpaWlilwYgAjEX5gUAEckXvl98Tirb+Wp9YxPNMD6WshMH94pdGkCXQo8QQKsFO9OxI8vK7hdWLNhdNm9n+v3GO3fuiF0UQJdCEAJoOytjPWNpLeF5QkhlZWVhi77YFQF0KQQhgLazsLB4fd5M6++GW2/wHztu3IzUHsuTpPWc2GUBdBXcRwhiwn2E6qOxsZHneUNDw/uN5J8p0sQH/L9HMJN6PmTpGNxHqIlwH+FjYLAMABBCiIGBgfCDnQHZHcicLuZfT5D2MyNbRzI9sfA9dGv4dAAADzHWnkoPZX2sKR8sfA/dHYIQAB7OgCWrhjAJU9nIfNnQw1xqKcIQuicEIQA8Tj8zKmYyu3wgHRLDLU+S1raIXRCAqiEIAeAJKEIW9KMzZ+o0ScmAcO5gntgFAagUghAA2sRCj/w8itk3lvnqMjXjNJWHhe+hu0AQAkA7jLajUoP5cXbEL4JbdUnajIXvQfMhCAGgfXRo8tYAPjmETXnA+0VwSVj4HjQcghAAOqK3CXV8Evu1L/18nHTBGWm5ROyCADoKQQgAHRfsTF+dyVroEc+D3O4cXCcFjYQgBICnYqZLNvkzkc8y/74uG3uMSy+q+9fXa15c+nbC+fNilwbQJghCAFCBIdZUYjA73YUePnvp6myDfZYzQ5d9mJ2dLXZdAE+GIAQA1RAWvreovCkNfIP0HVnh89LZc/FiFwXwZAhCAFCl3j0dqCvHSGUhuXQ4XDq4tEnsggCeBEEIAKr0557/zKo94XvizV9WLBk5bMiggy2brslkuMMC1BiWYQIAVbK1tf1j5zb5rzN60a8lSMNzZdtGMgMtsJwTqCP0CAGgE3laUuensa+60ROiuOVJ0jrM2Q3qB0EIAJ1LmLM7fYZOpYR4HeKOF+A6KagXBCEAdAVh4fvtY5h/pEiDo7mCesQhqAsEIQB0ncC/F773xcL3oDYQhADQpfSZvxa+P1og84vgLmDhexAbghAARNDPjIoJYlcMpqdFc8sSpDUYRAPiQRACgGhmu9JZs3T0GeIejjm7QTQIQgAQk4Ue2eTP7B/LrLsim3qSy63FlVLoaghCABDfKDsqPZSd6EgPPYyF76GrIQgBQC3o0GS5B50Swl4o5X0juMQSdA2hi2CKNQBQI64mVNRzbGS+7IXT0mH6pbc2Li4qKenXy+no/l1mZmZiVwfdE3qEAKB2gp3pKzPYK3u+TfdYUvJ+crLjtM/XbBC7KOi2EIQAoI7MdIldywPi4E4I4ewHXs8tFrsi6LYQhACgpt5busBq/1Kd0/82/uOdC33nfZwqbZKKXRN0RwhCAFBTIVOnnNu/7ZfnrNKP/Z61IrCongw6yEUXYhANqBgGywCA+nJ3d3d3dxd+3h3IHM3nlyVI/aypf49gbAzELQ26D/QIAUBjTHWmMmey7hZk0CEsfA8qgyAEAE1iyJJVQ5jYIPaPu7IxR7lrlQhDeFoIQgDQPIMsqYRgdml/eiIWvoenhiAEAI0kLHx/eYZOpYQMCOcO5WJaNuggBCEAaDBbA7I7kPktkPn0ogwL30PHIAgBQOMF2FMZM9hRtvTgQ9zaDCx8D+2DIASA7kCHJh950SkhbFyRzDeCS3mAMIS2QhACQPfR15SKDmJXDqanx2Dhe2grBCEAdDezXems2Vj4HtoKQQgA3ZC5Ltnkz/x3HLP+imzKSS63lq+pqYmKisrOzha7NFA7CEIA6LZG2lLpoeyzjrTvrgKnIQFzfzk/6sW3Nv/8i9h1gXpBEAJAd8bSZLkHvawhonbUG7WTV5Uvjfh+63axiwL1giAEgO6vr625fu09QgipKa3m9SskYhcE6gRBCADd37wX5vq3XLfZOLrnrtnB/1jT/0BLWDam7Ia/YBkmAOj+dHV1444caG5u1tXVJYRcLueXJUh/vyXbOopxN6fErg5Ehh4hAGgLIQUJIYOtqKRp7Ly+9JhIDgvfA4IQALQRTZGl/emrM3WK6onHQe7kPVwo1V4IQgDQXvaGZHcg86M/89p56Zw4aUmj2AWBGBCEAKDtJjtR12ey7hbEEwvfayUEIQAAMWDJqiFM/FT2SL5s6GHuYhnCUIsgCAEA/vKMGRU7mX1nIB18klueJK3FnN3aAUEIAPD/hIXvM2fpNEnJgHDu4F3M2d39IQgBAJRZ6pGfRzF7xzKfpcmCo7n8Olwp7c4QhAAADzfGjro8gx1lS3v/ya3NkHHoHHZTCEIAgEcSFr6/EMKeLpb5RnDJWPi+O0IQAgA8QR9T6sQk9hNvOjSGW3BGevbC5ZB5i19a+nZ+fr7YpYEKIAgBANpEWPjesKVq/JxXjvR6Za/RpHEhz4tdFKgAghAAoK3MdclC45vGA0YQ16H8wGerGZOKigqxi4KnhSAEAGiH/v376+cmkfzLJPt0ZWXVj7lmEszZreEQhAAA7WBhYXHst7BJ2ZvnVPwZf2T/tUoy6BAXV4RBNBoM6xECALSPj4/P8T92Cz+HDySR+bLF56Rj7KgNw5ke+uKWBh2BHiEAwFMJdqavz2IdjMjAcCx8r5EQhAAAT8uQJWv8mJjJ7M6bsoCjXGYl0lCTIAgBAFTDy5I6H8wucaPHR3HLk6T1nNgFQdsgCAEAVIamyIJ+9OUZOpUS4nWIO4GF7zUBBssAAKiYnQHZHcicLuZfT5D2MyNbRzI9jSixi4JHQo8QAKBTjLWn0kNZH2vKJ4LbdE0mRedQXSEIAQA6i7DwfcJUNjJfNvQwl1qKMFRHCEIAgM7Vz4yKmcwuH0iHxGDhe3WEIAQA6HR/LXw/86+F7w9g4Xt1giAEAOgiFnrk51HMvrHMqjRZcDSXh4Xv1QOCEACgS422oy7PYCc40H4R3PMbj9g949Wjn+cbH6wUuy7thdsnAAC6mg5NlnvQU5wpD89/ST44TfRN/tg1b3Famo+Pj9ilaSP0CAEAxOFqJDMz0CH6JoSQWlOXvJJysSvSUghCAABxMAwzedxoy99eNon81Cj3/FsP/HfnYBCNCHBpFABANDt/+iExMbGysnL8+C+u1+ktS5DuvCnbOpLpb46ZaLoOghAAQEwjRowQfhiiT5KmsT9dlwUc414fQK/wYvQYcUvTFrg0CgCgLliaLPeg00PZzEoy6BAXW4j7K7oCghAAQL04GFIHxjMbhtGvJkgXnJGWNoldUHeHIAQAUEfBznTmTLa3KfE42LLpmgwr33ceBCEAgJoyZMmqIUx0ELv/jizwGBa+7ywIQgAAteZlSSVOw8L3nQhBCACg7oQ5u4WF7z0PYuF7FcPtEwAAmgEL33cSMXuEtbW19fX1IhYAAKBxsPC9ynVFEBYWFg4ZMsTKymr16tXyxi+//HLKlCkTJkz4/vvvu6AGAIBuQ77w/dECLHyvAl1xadTa2vr06dPbtm2Tyf6aRu/u3bvh4eHp6elSqXTQoEHz5s2ztbXtgkoAALqNfmZUdBC7J0cWEsPNdqW/9mXyblzLysoaPXq0nZ2d2NVpkq7oEerp6ZmZmSm2JCUlBQQEMAyjq6vr7+9/4cKFLigDAKCbkS98TwhxWf6r/7x3F4bneAVMunHjhtilaRJV9ggLCwsPHz6s2GJvbx8aGtp6y4qKClNTU+Fnc3Pz8nIsPgIA0EEWemSTPxPx1s78JYeIgWmjjXvY7r0bvvlC7Lo0hip7hDRN6/8vXV3dh25paWlZXV0t/FxVVWVtba3CMgAAtJCLnTUpvU0IoUpu5kitOSzo1GZt6hEWFRVt2bLl4sWLlZWVKSkp8vba2tq33nrr5MmTNjY233zzTXBw8CuvvNKWHY4YMWL16tUcx8lksqSkpHXr1nWwfAAAIIQQsn3jmqkvvlJZW9+3T5+mETt9Irhtoxh/G9xf8WRtCsLKykqJRBIQEPDZZ58ptq9cubK8vDwzM/PChQvPP/98dna2g4ND64dzHDd8+PD79+8TQg4ePJiYmNirV6/Zs2cHBgZyHLds2TIbGxuVvBgAAK3Vr1+/G6nx8l8j82XPx0kD7akf/BkrPRHr0gAUz7d13O3Vq1eHDBnS0tIi/CqRSGxsbE6dOuXj40MICQkJGT58+IoVK9r+3HV1dTRNGxoaPmqD/Px8d3f3uXPnylumTp06adKkR21fW1trYmLS9gJAdDzPNzQ0GBkZiV0ItENzczPP83p6eHNVd1XN5Osr9KE86itvWahdvb6+Pk1r3WxiOjo6DPOEdR07PlimsLCwtrbWy8tL+HXw4MHZ2dnt2oOxsfHjN2AYRkdHRwhagYODw2NeEsMwT3zBoFZ4nsdZ0zgMwwgnTuxC4AmsDMgPw8iCvuSNJPrXHMOf/KkBFloXhBT15IvDHQ/CiooKAwMDlv1rD2ZmZqWlpR3e20PRNG1gYPD666+3cXsdHR0dHR3V1gCdiud5nDWNw/O8cOLELgTaZKgdSQ4hP1yWToim33DHwvcP0fFPBxYWFo2NjRz310To1dXVGPwJAKCGWJq82V+WFkJj4fuH6ngQOjo66uvrZ2VlCb9ev369b9++KqoKAABUzNHor4Xvl8Rj4fv/0aYglEqld+7cuXfvHs/zd+7cyc/PJ4To6+vPnTt39erVzc3NFy9ePHny5Pz58zu5WgAAeCrBzvT1WayDERkYjoXv/9KmUaPl5eVDhw6V/2pra5uYmEgIKSsrW7Ro0dmzZ01MTNauXavyICwuLvbx8SkqKmrj9hg1qnEwalQTYdSoJmpoaFAaNZpRwb+WINWhydaRzEALrb7dsB23T3Q9BGG3hyDURAhCTdQ6CAkhMp78dkv2zwvS53vT3/oxRtq6QK3WDaUFAAABTWHhe0KwQj0AgJbDwvfoEQIAgPLC9xu3hI2aPPPDf33Z1NT9R5ciCAEAgJC/F76Pn8r+Z+8fH/wefz5g3eYb7D8+6f7LOSEIAQDg/z1jRg2tvSgduYhYOUsCXj+dkCR2RZ0OQQgAAP9jyrhRZud/JkVZzMkN93qOPni3m69tiMEyAADwP2aGTq+ta9h/5LtRowaPeOGtN1Nkv+bItoxknLrpIBrcRwhiwn2Emgj3EWqih95H2EYtMvL9Vdm6K9J/ejIfeNJMt0tDXBoFAIDH0aHJR1508jQ2tkjmF8Gllqpv96ljEIQAAPBk/cyo6CD2XQ96WjS3PEla2yJ2QaqDIAQAgDahCFnQj86cpdMkJQPCuW4ziAZBCAAA7WCpR34exewdy/wrTRYczeXXafyVUgQhAAC02xg7KmMGO8qW9v6TW5shk2pyGiIIAQCgI4RBNCkhGj+IBkEIAAAd19dU4wfRIAgBAOCpCINorv89iCZc0wbRIAgBAEAFLP4eRPOZpg2iQRACAIDKaOIgGgQhAACokuIgGt8IbtHH39r39+7rPSLhfKLYpT0cghAAAFRPGEQTLEvbE5N6//3k2y/tn/fau2IX9XAIQgAA6BQUIT7sfYNenoSiiZl9jUSqnss8YBkmAADoLAEBAZYrv2zWM9Etuy3t5fv8KdnWUYyVmq1cgh4hAAB0FnNz88vxMWHP9fjj3eDiY9t6mxLPg9wBNbu/AusRgpiwHqEmwnqEmuhp1iNUraQH/JJz0t6mZNtIxlE9VvoV/6AAAID28LehLs9gR9nSQyK4sGyZOnTFEIQAANClhPsrYoPY7TdkgUe5nGqR0xBBCAAAIhhkSZ0PZmf0okdGinzrPYIQAADEwdJkuQedHMLGFMpGR3LXq8QJQwQhAACIqbcJFTOZfW0APe4Y93GqtLnLh5QiCAEAQGTC+hXpoTo3q4nPn9yFrl3aEEEIAABqwd6QHJrAfDaEnhbNLUuQ1nNd9LwIQgAAUCOzXemsWTqEEM+DXFxRV3QNEYQAAKBehKUN/z2CWXxOuuCMtELSuU+HIAQAAHUU5ERdmcla6BHPQ9yh3E4cQoNJtwEAQE2Z6pBN/syc3vySeOn+2/yqfuU7wrZKJM0fvLXUxcVFVc+CHiEAAKi1kbZUeig7cKbBcAAAD8JJREFUwJz3HB/yfWmfnxp9Rk+e2djYqKr9IwgBAEDd6TPkdccyMwsr3u95fnBwY0/vrKwsVe0cl0YBAEAD9OjRQ7+mkNy/QXQN6fxLvXv3VtWeEYQAAKABGIaJ3Lfj7RWfSiSS77b9YG5urqo9IwgBAEAzDPH2Pn8iQuW77T7fEUql0piYGLGrgPaprKxMTk4Wuwpon9u3b9+4cUPsKqB9UlNTS0tLxa5CTXWfICwuLn7nnXfErgLaJzEx8bvvvhO7Cmif8PDwffv2iV0FtM+PP/547tw5satQU90nCAGga/C8OiwqDu2GE/coCEIAANBqCEIAANBqlDp3lgsKCp555plRo0a1ZWOJRJKamtrGjUFNlJeXFxQUDB48WOxCoB3u3r0rk8n69OkjdiHQDleuXLG1tbW1tRW7kK4WGhr6xhtvPH4btQ5CQsjOnTudnJzasiXP87m5ua6urp1dEqiQRCIpKytzdHQUuxBoh+rqaqlUamlpKXYh0A7FxcUWFhb6+vpiF9LVXF1dn/ihTd2DEAAAoFPhO0IAANBqCEIAANBqCEIAANBqCEIAANBqzKpVq8SuQTXOnTt36tQpXV1dGxsbsWsBQggpLi6+ePGirq6uqampvLG5uTkqKurChQs2NjYmJiby9jt37hw+fLi0tNTV1ZWiKKFRKpVGR0efP3/ewsJChTPNw6NUVlbGxcUlJyfTNK04zr6lpeXEiRPJycnW1taKZzM3N/fw4cMlJSWurq40/denaplMFhMTEx8fb2pqamFh0dWvQcvIZLLMzMxTp06lp6ezLKv47ldTU3P48OHMzEwnJyc9PT15e3p6+okTJyQSSc+ePeWNDQ0NkZGR6enpDg4OhoaGXfoa1AHfLbz22mvPPPPMsmXLbGxsdu7cKXY5wI8ePdrIyMjAwGDbtm3yxqampuHDh48aNWrRokWWlpYXL14U2o8dO2ZlZbVkyRJvb++QkBChUSaTTZo0ydfXd/HixVZWVrGxsSK8DG2SnZ1tYmISFBS0aNEia2vrDz74QGhvaWkZM2aMv7//K6+8YmlpmZiYKLTHxMRYWlouXrzY19c3KChIJpMJ7dOmTfP29l6yZImVlVVUVJQ4L0Zr5OTkDBgw4KWXXpo/f761tfXKlSuF9qKiImdn59DQ0NDQUGdn5+LiYqH9+++/d3BwWLZsmaur6yeffCI0VlVVDRgw4LnnnnvhhRdsbW1zcnLEeTHi6Q5BmJOTY2hoWFJSwvN8dHS0g4NDS0uL2EVpu7t373IcFxAQoBiEv/32m6enp3B2vvrqq+DgYKF98ODBwseX2tpae3v7c+fO8Tx/8uRJFxeXhoYGnue3bdvm7+/f9a9Cq1RXVz948ED4OTMzkxAivHuGh4f3799fIpHwPL9+/fqJEycK2wwfPnzr1q08zzc0NLi4uMTExPA8f/bsWXt7+9raWp7nf/31V29vb1Fei3ZKTk7W1dVtamrieX7FihXPP/+80D537twVK1bwPF9XV2dmZpaamsrz/J07dwwMDIS3zQ0bNowbN074KPPmm28uWbJEtNcgku7wHeGxY8dGjRolXBMYP358Q0NDWlqa2EVpu169ejEMo9R49OjR6dOnsyxLCJk1a9bx48elUmlBQUFGRsbMmTMJIcbGxpMmTTp69Kiw8eTJkw0MDISNk5KSysrKuvx1aBFTU9MePXoIP9vZ2VEU1dzcTAg5evTotGnTdHV1CSGzZs2KjY1tamoqLS1NTk6eNWsWIcTAwGDy5MnyszZp0iRjY2NCyMyZMy9fvnzv3j3RXpKWaWhoMDMzE/6+IiMjhb8pQsjMmTOFsxMfH29hYeHr60sIcXV19fDwOHnyJCHk6NGjM2fOFL6SmDVrlrCxVukOQVhYWCi/2E3TtL29fWFhobglwUMVFhbKJ5FxdHTkOK6kpKSoqMjMzEz+faGjo6Nw+hQ3trKy0tfXx2ntMl999dWECROcnZ1Jq7NGCCkqKioqKtLT07O2tpa3y8+a/I/R2NjY1NQUZ60LzJkzJyAg4JVXXomIiBA+gCqdtdZnhzzib83R0fHBgwctLS1d/RpE1R1WqJdKpfLhFYQQlmU5jhOxHngUqVQqH1Ih/LlyHKd0+hiGEU6f4saK7dDZwsLC/vzzz4SEBOFXxRMhnCnhrD307OCPURRvvfVWVVXV9u3bP/vss+joaJqmFU/EE8+O0h8mz/NSqVRHR6fLX4doukMQ2tvb5+TkyH8tKSlxcHAQsR54FHt7+wcPHgg/l5SU0DRtZ2fHcVx1dbVEIhEGtpWUlNjb2yttXFdXV19fj9PaBXbv3v3VV1+dPn1a3nVQPBGlpaU8zzs4ONTU1DQ2NtbW1gpd+Yeetebm5srKSpy1LjBmzBhCSFBQkJWV1YULF4YPH25vby9fj17+lqh4dggh9+/fDwoKIq3+MC0tLbVtStLucGk0MDAwPj6+sbGREJKeni6RSLy9vcUuCh4iMDBQ+E6CEBIdHT1ixAhdXd1evXq5uLjExsYSQqRSaVxc3NixY4WNY2NjZTKZsHH//v3t7OxELF4bHDhwYMWKFSdPnuzbt6+8MTAwMDo6mud5Qkh0dLSvr6+xsbG9vb2bm1t0dDQhRCaTxcbGys9aXFycVColhMTExLi4uLi4uIj0arROTU2NRCIR7m8ZO3ascHYIIdHR0YGBgYQQf3//vLy83NxcQkhlZWVaWlpAQAD5+xQrbaxdRB6soyKTJk0aP378Dz/84Obm9sUXX4hdDvA///zz0qVL7e3tR48evXTp0vT0dJ7nq6urXVxcXn755TVr1lhaWsrH1oeFhTk6Om7YsCE0NNTb25vjOJ7nJRLJwIED58yZ891339nZ2e3Zs0fM16MFMjMzWZYNCAhY+rfs7Gye5+vq6vr06TN//vx169ZZW1v/+eefwva7d++2s7P77rvv5syZ4+Hh0dzczPM8x3He3t6hoaEbNmxwdHT85ZdfxHxJWmDXrl3z5s1bs2bN559/7ubmNmvWLGHwZ3Z2tpmZ2YoVK1auXGlubi6cSp7n33nnHS8vr40bN44YMWLevHlC471796ytrZcvX/7FF1+YmpqmpKSI9npE0k1Wn5BIJLt27bp79+7w4cNDQkLELgdIbGzsnTt35L8+++yzvXr1IoQ8ePBg165dNTU1wcHBQ4cOVdz+9OnTtra2L7/8snzgTHV19c6dO0tLSydOnKiNn1K71v37948cOaLYMnXqVOGSWnl5+a+//lpZWTl58uQRI0bINzhz5kxMTEyPHj1efvllMzMzobG2tnbnzp0lJSXjxo0bP358V74ELVReXn706NFbt27p6ur6+fk999xz8m8Bb926tW/fPkLICy+8IO/i8zz/xx9/pKenDxgwYN68ecIQU0JIQUHBnj17mpubZ82a5eHhIcprEVE3CUIAAICO6Q7fEQIAAHQYghAAALQaghAAALQaghAAALQagvD/2rvbkKa+OA7gZ6TL5px61fagUhZqGkUhpC/E9kIEjaQgNMIXWmKjCFYaRFgWJnsTQRJIRcLsnU9EUREkTX2hpoiJyHaVrO6ki1erzfmAc+7/4sDlon+W/4e53L6fVzu/c3bPuUP4eXd3zw8AAEIaEiEAAIQ0JEIAf3E4HGazObDlF5xOp9ls/vbtWwDXAPCHQyIE8Be73V5eXj4yMkKb79+/7+zs9OuMfX19z58/l0Z4ni8vL0dhMgAfkAgB/EWtVtfX1x84cIA2nz59eu/ePb/O2Nraev36dWkkPj6+vr4+MzPTr/MCbGvBUH0C4M8UHx9fW1u7mZFLS0sOh4NhGFr/dqOFhQWXy6VWq2nT4/HMzc3FxcVtrH68DsMwG9ewtrYmCMKuXbvoBs0bra2t0eNLay0BBCv8lQP4C8uyWq2WFtY4derUixcvRkdHGYZhGCY3N5eO4Tju9OnT0dHRWq02JibGYDAsLy/Trq6uLro1eUFBQVRUFL2y7O3tzc3NjYiIUKvVkZGRer3earXS8QaDoampSRAEOgXdXvLz589arfbt27fiqh4+fKjT6TQaTXR0dHZ29sePH8Wus2fPFhUVtbW1JScn7969W6VSVVdXYxdGCHq4IgTwF7fbzfM8TWy3b9/+9evX9PR0U1MTIUSpVBJCfvz4kZeXp1Ao2traUlNTBwYGrl27trCwQO/zud3unz9/XrhwoaysrLa2dn5+nhDC83xBQYHJZEpISGBZ9saNG4WFhTabTS6XX7lyRRCEDx8+tLa2EkJoYdXV1VVxDYSQxsZGo9FYWVlpMBhmZ2erq6vz8/NHRkb27dtHCHG5XIODgxzHPXr0KCkpyWw2P3jwQK/Xnzx5MjCfIMDWCGzxC4AgNjY2Rgh59eoVbZaUlBw9elQ6oK6uLiIi4uvXr2LkyZMnMpnMbrd7vV56GXfx4kUfUwwPDxNCLBYLbRqNRrVaLR1gs9kIIZ2dnV6v1+PxJCQk5OXlib1fvnwJDw+/dOkSbZ44cSI8PHxycpI2V1dXk5OTKyoq/tXZA2wbuCIECJh3796lpKSwLMuyLI2EhYV5vd7x8fHExEQaKS4uXveuqampjo4OjuOWl5dpCdzJyUlaYdW36elpQRDu3r0rRvbs2ZOTk9PT0yNGUlNT9+/fT1/v2LEjPT2d47j/cIoA2wASIUDAzMzM2O32kpISaTA2NnZ2dlZsajQaae/jx48vX7585MiRrKys2NhYmggdDsdmpqNPE2q1WmlQp9NNTEyITYZhpL07d+6kX8kCBDEkQoCAUalUx44d6+3t9TFGrLNK3blz59y5cy0tLbTJcdz9+/c3OV1kZCQhRBAEaVAQhJiYmH+waICgg1+NAmwRpVK5tLQkjRw/fnxoaGhqamqTR1hcXOR5PisrS4y8efPG9xRSmZmZKpXq9evXYkQQhP7+/pycnE0uACAoIRECbJGDBw9OTEy0tLQMDg6Oj48TQmpqapRKZXFxcVdX1/z8/MzMjMViOX/+/MrKyt8eQaFQpKWlPXv2jGXZxcXF9vb2hoaGdVM4nc7GxsaBgYFPnz6te7tcLjcajS9fvjSZTHNzczabrbS01O12X7161U+nDLAt4KtRgC1SVVU1NDRUU1MjCMKhQ4dGR0eTkpK6u7urqqry8/PpGLlcrtfrfTzG3tzcXFpamp6eTgjRaDTNzc1FRUVi75kzZywWi8lk4nk+Li5Oeq+RunXrlsvlqquru3nzJiFEp9N1dHQcPnz4/z9bgO1D5sXTsgB+4/F4frv5CyHk+/fvHMcplcq9e/cqFArfg1dWVqxWq0wmy8jICAv7/f+yG9fgdDqtVqtCocjIyNjM8gCCGxIhAACENNwjBACAkIZECAAAIQ2JEAAAQhoSIQAAhDQkQgAACGlIhAAAENL+Av3bNoZldQf6AAAAAElFTkSuQmCC", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# | echo: false \n", "\n", "B = randn(1000,1000)#[5.001 1 0; 1 3.2 0 ; 0 0 1]\n", "A = B + B' #create a random symmetric matrix\n", "λmax = argmax( abs, eigen(A).values )\n", "\n", "λ1, v = power(A; Niter=10000)\n", "#λ2, v = power(A; Niter=10000)\n", "#λ3, v = power(A; Niter=10000)\n", "\n", "ii1 = 1:Int(round(length(λ1)/20)):length(λ1)\n", "#ii2 = 1:Int(round(length(λ)/20)):length(λ2)\n", "#ii3 = 1:Int(round(length(λ)/20)):length(λ3)\n", "\n", "plot( ii1, @. abs( λ1[ii1] - λmax ) + 1e-20; \n", " m=2, title=\"error in the dominant eigenvalue\", yaxis=:log, xlabel=\"iteration\", legend=false )\n", "#plot!( ii2, @. abs( λ2[ii2] - λmax ) + 1e-20; m=2 )\n", "#plot!( ii3, @. abs( λ3[ii3] - λmax ) + 1e-20; m=2 )" ] }, { "cell_type": "markdown", "id": "dc8be9c0", "metadata": {}, "source": [ "4. What is the rate of convergence in this plot?\n", "\n", "One can show that, if $z \\not\\in \\sigma(A)$, then $zI - A$ is invertible and $\\sigma\\big( (zI-A)^{-1} \\big) = \\left\\{ (z- \\lambda)^{-1} : \\lambda\\in\\sigma(A) \\right\\}$. Suppose that, after possibly reordering $\\lambda_1,\\dots,\\lambda_n \\in \\sigma(A)$, we have\n", "\n", "\\begin{align}\n", " |z - \\lambda_1| < |z - \\lambda_2| \\leq \\cdots \\leq |z- \\lambda_n|\\nonumber\n", "\\end{align}\n", "\n", "5. Explain how you can apply the power method to $B := (zI-A)^{-1}$ to approximate $\\lambda_1$ (the closest eigenvalue to $z$).\n", "\n", "Here, we apply this idea to the same $1000 \\times 1000$ matrix as before to find the closest eigenvalue to zero." ] }, { "cell_type": "code", "execution_count": 5, "id": "672ee26a", "metadata": {}, "outputs": [], "source": [ "#| include: false\n", "function invpower( A, z; Niter=100, tol=1e-8 )\n", " n, = size(A)\n", " x = normalize( randn(n) )\n", " λ = []\n", " for k ∈ 1:Niter\n", " x′ = A*x\n", " push!( λ, dot(x,x′) )\n", " x = normalize(x′)\n", " \n", " # stopping criteria\n", " r = norm(A*x-λ[end]*x)/norm(A)\n", " if ( r < tol )\n", " σ = z - 1/λ[end]\n", " println(\"converged in $k iterations\\neigenvalue closest to z=$z is λ ≈ $σ, \\nresidual = $r\")\n", " break;\n", " elseif (k == Niter)\n", " println(\"did not converge in $Niter iterations, residual = $r\")\n", " end\n", " end\n", " return λ, x\n", "end;" ] }, { "cell_type": "code", "execution_count": 6, "id": "5f9f0185", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "converged in 34 iterations\n", "eigenvalue closest to z=0 is λ ≈ 0.06388010931076409, \n", "residual = 7.3399398262554825e-9\n" ] } ], "source": [ "#| echo: false\n", "\n", "z = 0;\n", "σ, v = invpower( inv( z*I - A ), z );" ] }, { "cell_type": "markdown", "id": "9977e583", "metadata": {}, "source": [ "::: {.box}\n", "\n", "**Answer.** \n", "\n", "**1.** If $v_1$ is an eigenvector of $A$ with eigenvalue $\\lambda_1$ then for all $\\alpha \\not=0$ we have $A(\\alpha v_1) = \\alpha Av_1 = \\alpha \\lambda_1 v_1 = \\lambda_1 (\\alpha v_1)$ and so $\\alpha v_1$ is an eigenvector of $A$ with eigenvalue $\\lambda$.\n", "\n", "**2.** If $x^{(0)} = \\sum_i c_i v_i$, then $Ax^{(0)} = \\sum_i c_i Av_i = \\sum_i \\lambda_i c_i v_i$, therefore, by induction, we have \n", "\n", "\\begin{align}\n", " A^k x^{(0)} &= A A^{k-1} x^{(0)} = \\sum_i \\lambda_i^{k-1} c_i Ax^{(0)} \\nonumber\\\\\n", " %\n", " &= \\sum_i \\lambda_i^{k} c_i x^{(0)} \\nonumber\\\\\n", " %\n", " &= \\lambda_1^{k} c_1 x^{(0)} + \\sum_{i\\not=1} \\lambda_i^{k} c_i x^{(0)} \\nonumber\\\\\n", " %\n", " &= \\lambda_1^{k} \\Big[ c_1 x^{(0)} + \\sum_{i\\not=1} (\\frac{\\lambda_i}{\\lambda_1})^{k} c_i x^{(0)} \\Big]. \n", "\\end{align}\n", "\n", "**3.** $(\\lambda_1, v_1)$ are unknown! We are trying to approximate these things.\n", "\n", "**4.** This is a straight line on a semi-log plot and so we have exponential decay: $O(\\rho^{-n})$ for some $\\rho > 1$ (if $y = C \\rho^{-n}$ then $\\log y = -n \\log \\rho + \\log C$). In fact, we have $\\rho = \\left|\\frac{\\lambda_1}{\\lambda_2}\\right|$ where $\\lambda_2$ is the second largest eigenvalue in absolute value.\n", "\n", "**5.** The eigenvalues of $B$ are $(z-\\lambda_i)^{-1}$ with \n", "\n", "\\begin{align}\n", " \\left|(z-\\lambda_1)^{-1}\\right| > \\left|(z-\\lambda_2)^{-1}\\right| > \\cdots.\n", "\\end{align}\n", "\n", "As a result $(z-\\lambda_1)^{-1}$ is the dominant eigenvalue of $B$ and applying the power method to $B$ converges to $\\sigma_1 \\approx (z-\\lambda_1)^{-1}$ and we can recover the eigenvalue of $A$ closest to $z$ by rearranging this equation:\n", "\n", "\\begin{align}\n", " \\lambda_1 \\approx z-\\frac{1}{\\sigma_1}.\n", "\\end{align}\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 }